Program f08chfe ! F08CHF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: dgerqf, dormrq, dtrtrs, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Real (Kind=nag_wp), Parameter :: zero = 0.0E0_nag_wp Integer, Parameter :: nb = 64, nin = 5, nout = 6 ! .. Local Scalars .. Integer :: i, info, lda, lwork, m, n ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: a(:,:), b(:), tau(:), work(:), x(:) ! .. Executable Statements .. Write (nout,*) 'F08CHF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) Read (nin,*) m, n lda = m lwork = nb*m Allocate (a(lda,n),b(m),tau(m),work(lwork),x(n)) ! Read the matrix A and the vector b from data file Read (nin,*)(a(i,1:n),i=1,m) Read (nin,*) b(1:m) ! Compute the RQ factorization of A ! The NAG name equivalent of dgerqf is f08chf Call dgerqf(m,n,a,lda,tau,work,lwork,info) ! Copy the m element vector b into elements x(n-m+1), ..., x(n) of x x(n-m+1:n) = b(1:m) ! Solve R*y2 = b, storing the result in x2 ! The NAG name equivalent of dtrtrs is f07tef Call dtrtrs('Upper','No transpose','Non-Unit',m,1,a(1,n-m+1),lda, & x(n-m+1),m,info) If (info>0) Then Write (nout,*) 'The upper triangular factor, R, of A is singular, ' Write (nout,*) 'the least squares solution could not be computed' Else x(1:n-m) = zero ! Compute the minimum-norm solution x = (Q**T)*y ! The NAG name equivalent of dormrq is f08ckf Call dormrq('Left','Transpose',n,1,m,a,lda,tau,x,n,work,lwork,info) ! Print minimum-norm solution Write (nout,*) 'Minimum-norm solution' Write (nout,99999) x(1:n) End If 99999 Format (1X,8F9.4) End Program f08chfe