Program f08kgfe ! F08KGF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: dgebrd, dgelqf, dgeqrf, dorglq, dorgqr, dormbr, & f06qff, f06qhf, nag_wp, x04caf ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Real (Kind=nag_wp), Parameter :: zero = 0.0E0_nag_wp Integer, Parameter :: nin = 5, nout = 6 ! .. Local Scalars .. Integer :: i, ic, ifail, info, lda, ldpt, ldu, & lwork, m, n ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: a(:,:), d(:), e(:), pt(:,:), tau(:), & taup(:), tauq(:), u(:,:), work(:) ! .. Executable Statements .. Write (nout,*) 'F08KGF Example Program Results' ! Skip heading in data file Read (nin,*) Do ic = 1, 2 Read (nin,*) m, n lda = m ldpt = n ldu = m lwork = 64*(m+n) Allocate (a(lda,n),d(n),e(n-1),pt(ldpt,n),tau(n),taup(n),tauq(n), & u(ldu,n),work(lwork)) ! Read A from data file Read (nin,*)(a(i,1:n),i=1,m) If (m>=n) Then ! Compute the QR factorization of A ! The NAG name equivalent of dgeqrf is f08aef Call dgeqrf(m,n,a,lda,tau,work,lwork,info) ! Copy A to U Call f06qff('Lower',m,n,a,lda,u,ldu) ! Form Q explicitly, storing the result in U ! The NAG name equivalent of dorgqr is f08aff Call dorgqr(m,n,n,u,ldu,tau,work,lwork,info) ! Copy R to PT (used as workspace) Call f06qff('Upper',n,n,a,lda,pt,ldpt) ! Set the strictly lower triangular part of R to zero Call f06qhf('Lower',n-1,n-1,zero,zero,pt(2,1),ldpt) ! Bidiagonalize R ! The NAG name equivalent of dgebrd is f08kef Call dgebrd(n,n,pt,ldpt,d,e,tauq,taup,work,lwork,info) ! Update Q, storing the result in U ! The NAG name equivalent of dormbr is f08kgf Call dormbr('Q','Right','No transpose',m,n,n,pt,ldpt,tauq,u,ldu, & work,lwork,info) ! Print bidiagonal form and matrix Q Write (nout,*) Write (nout,*) 'Example 1: bidiagonal matrix B' Write (nout,*) 'Diagonal' Write (nout,99999) d(1:n) Write (nout,*) 'Super-diagonal' Write (nout,99999) e(1:n-1) Write (nout,*) Flush (nout) ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call x04caf('General',' ',m,n,u,ldu,'Example 1: matrix Q',ifail) Else ! Compute the LQ factorization of A ! The NAG name equivalent of dgelqf is f08ahf Call dgelqf(m,n,a,lda,tau,work,lwork,info) ! Copy A to PT Call f06qff('Upper',m,n,a,lda,pt,ldpt) ! Form Q explicitly, storing the result in PT ! The NAG name equivalent of dorglq is f08ajf Call dorglq(n,n,m,pt,ldpt,tau,work,lwork,info) ! Copy L to U (used as workspace) Call f06qff('Lower',m,m,a,lda,u,ldu) ! Set the strictly upper triangular part of L to zero Call f06qhf('Upper',m-1,m-1,zero,zero,u(1,2),ldu) ! Bidiagonalize L ! The NAG name equivalent of dgebrd is f08kef Call dgebrd(m,m,u,ldu,d,e,tauq,taup,work,lwork,info) ! Update P**T, storing the result in PT ! The NAG name equivalent of dormbr is f08kgf Call dormbr('P','Left','Transpose',m,n,m,u,ldu,taup,pt,ldpt,work, & lwork,info) ! Print bidiagonal form and matrix P**T Write (nout,*) Write (nout,*) 'Example 2: bidiagonal matrix B' Write (nout,*) 'Diagonal' Write (nout,99999) d(1:m) Write (nout,*) 'Super-diagonal' Write (nout,99999) e(1:m-1) Write (nout,*) Flush (nout) ifail = 0 Call x04caf('General',' ',m,n,pt,ldpt,'Example 2: matrix P**T', & ifail) End If Deallocate (a,d,e,pt,tau,taup,tauq,u,work) End Do 99999 Format (3X,(8F8.4)) End Program f08kgfe