Program f08zffe ! F08ZFF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: dgemv, dggrqf, dnrm2, dormqr, dormrq, dtrmv, & dtrtrs, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Real (Kind=nag_wp), Parameter :: one = 1.0E0_nag_wp Integer, Parameter :: nb = 64, nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: rnorm Integer :: i, info, lda, ldb, lwork, m, n, p ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: a(:,:), b(:,:), c(:), d(:), taua(:), & taub(:), work(:), x(:) ! .. Intrinsic Procedures .. Intrinsic :: min ! .. Executable Statements .. Write (nout,*) 'F08ZFF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) Read (nin,*) p, n, m lda = m ldb = p lwork = nb*(p+n) Allocate (a(lda,n),b(ldb,n),c(p),d(m),taua(min(m,n)),taub(min(p, & n)),work(lwork),x(n)) ! Read B, A, C and D from data file Read (nin,*)(b(i,1:n),i=1,p) Read (nin,*)(a(i,1:n),i=1,m) Read (nin,*) c(1:p) Read (nin,*) d(1:m) ! Compute the generalized RQ factorization of (B,A) as ! A = (0 R12)*Q, B = Z*(T11 T12 T13)*Q, where R12, T11 and T22 ! ( 0 T22 T23) ! are upper triangular ! The NAG name equivalent of dggrqf is f08zff Call dggrqf(m,p,n,a,lda,taua,b,ldb,taub,work,lwork,info) ! Compute (f1) = (Z**T)*c, storing the result in C ! (f2) ! The NAG name equivalent of dormqr is f08agf Call dormqr('Left','Transpose',p,1,min(p,n),b,ldb,taub,c,p,work,lwork, & info) ! Putting Q*x = (y1), solve R12*w = d for w, storing result in D ! (w ) ! The NAG name equivalent of dtrtrs is f07tef Call dtrtrs('Upper','No transpose','Non-unit',m,1,a(1,n-m+1),lda,d,m, & info) If (info>0) Then Write (nout,*) 'The upper triangular factor, R12, of A is singular, ' Write (nout,*) 'the least squares solution could not be computed' Else ! Form f1 - T1*w, T1 = (T12 T13), in C Call dgemv('No transpose',n-m,m,-one,b(1,n-m+1),ldb,d,1,one,c,1) ! Solve T11*y1 = f1 - T1*w for y1, storing result in C ! The NAG name equivalent of dtrtrs is f07tef Call dtrtrs('Upper','No transpose','Non-unit',n-m,1,b,ldb,c,n-m,info) If (info>0) Then Write (nout,*) & 'The upper triangular factor, T11, of B is singular, ' Write (nout,*) 'the least squares solution could not be computed' Else ! Copy y into X (first y1, then w) x(1:n-m) = c(1:n-m) x(n-m+1:n) = d(1:m) ! Compute x = (Q**T)*y ! The NAG name equivalent of dormrq is f08ckf Call dormrq('Left','Transpose',n,1,m,a,lda,taua,x,n,work,lwork,info) ! Putting w = (y2), form f2 - T22*y2 - T23*y3 ! (y3) ! d = T22*y2 ! The NAG name equivalent of dtrmv is f06pff Call dtrmv('Upper','No transpose','Non-unit',min(p,n)-n+m, & b(n-m+1,n-m+1),ldb,d,1) ! c = f2 - T22*y2 Do i = 1, min(p,n) - n + m c(n-m+i) = c(n-m+i) - d(i) End Do If (p