* F08BHF Example Program Text * Mark 21 Release. NAG Copyright 2004. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER MMAX, NB, NMAX, NRHSMX PARAMETER (MMAX=8,NB=64,NMAX=8,NRHSMX=2) INTEGER LDA, LDB, LWORK PARAMETER (LDA=MMAX,LDB=MMAX,LWORK=2*NMAX+(NMAX+1)*NB) DOUBLE PRECISION ONE, ZERO PARAMETER (ONE=1.0D0,ZERO=0.0D0) * .. Local Scalars .. DOUBLE PRECISION TOL INTEGER I, IFAIL, INFO, J, K, M, N, NRHS * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), B(LDB,NRHSMX), RNORM(NMAX), + TAU(NMAX), WORK(LWORK) INTEGER JPVT(NMAX) * .. External Functions .. DOUBLE PRECISION DNRM2 EXTERNAL DNRM2 * .. External Subroutines .. EXTERNAL DCOPY, DGEQP3, DORMQR, DORMRZ, DTRSM, DTZRZF, + F06DBF, F06QHF, X04CAF * .. Intrinsic Functions .. INTRINSIC ABS * .. Executable Statements .. WRITE (NOUT,*) 'F08BHF Example Program Results' WRITE (NOUT,*) * Skip heading in data file READ (NIN,*) READ (NIN,*) M, N, NRHS IF (M.LE.MMAX .AND. N.LE.NMAX .AND. M.GE.N .AND. NRHS.LE.NRHSMX) + THEN * * Read A and B from data file * READ (NIN,*) ((A(I,J),J=1,N),I=1,M) READ (NIN,*) ((B(I,J),J=1,NRHS),I=1,M) * * Initialize JPVT to be zero so that all columns are free * CALL F06DBF(N,0,JPVT,1) * * Compute the QR factorization of A with column pivoting as * A = Q*(R11 R12)*(P**T) * ( 0 R22) * CALL DGEQP3(M,N,A,LDA,JPVT,TAU,WORK,LWORK,INFO) * * Compute C = (C1) = (Q**T)*B, storing the result in B * (C2) * CALL DORMQR('Left','Transpose',M,NRHS,N,A,LDA,TAU,B,LDB,WORK, + LWORK,INFO) * * Choose TOL to reflect the relative accuracy of the input data * TOL = 0.01D0 * * Determine and print the rank, K, of R relative to TOL * DO 20 K = 1, N IF (ABS(A(K,K)).LE.TOL*ABS(A(1,1))) GO TO 40 20 CONTINUE 40 K = K - 1 * WRITE (NOUT,*) 'Tolerance used to estimate the rank of A' WRITE (NOUT,99999) TOL WRITE (NOUT,*) 'Estimated rank of A' WRITE (NOUT,99998) K WRITE (NOUT,*) * * Compute the RZ factorization of the K by K part of R as * (R11 R12) = (T 0)*Z * CALL DTZRZF(K,N,A,LDA,TAU,WORK,LWORK,INFO) * * Compute least-squares solutions of triangular problems by * back substitution in T*Y1 = C1, storing the result in B * CALL DTRSM('Left','Upper','No transpose','Non-Unit',K,NRHS,ONE, + A,LDA,B,LDB) * * Compute estimates of the square roots of the residual sums of * squares (2-norm of each of the columns of C2) * DO 60 J = 1, NRHS RNORM(J) = DNRM2(M-K,B(K+1,J),1) 60 CONTINUE * * Set the remaining elements of the solutions to zero (to give * the minimum-norm solutions), Y2 = 0 * CALL F06QHF('General',N-K,NRHS,ZERO,ZERO,B(K+1,1),LDB) * * Form W = (Z**T)*Y * CALL DORMRZ('Left','Transpose',N,NRHS,K,N-K,A,LDA,TAU,B,LDB, + WORK,LWORK,INFO) * * Permute the least-squares solutions stored in B to give X = P*W * DO 100 J = 1, NRHS DO 80 I = 1, N WORK(JPVT(I)) = B(I,J) 80 CONTINUE CALL DCOPY(N,WORK,1,B(1,J),1) 100 CONTINUE * * Print least-squares solutions * IFAIL = 0 CALL X04CAF('General',' ',N,NRHS,B,LDB, + 'Least-squares solution(s)',IFAIL) * * Print the square roots of the residual sums of squares * WRITE (NOUT,*) WRITE (NOUT,*) + 'Square root(s) of the residual sum(s) of squares' WRITE (NOUT,99999) (RNORM(J),J=1,NRHS) ELSE WRITE (NOUT,*) + 'One or more of MMAX, NMAX and NRHSMX is too small, ', + 'and/or M.LT.N' END IF * 99999 FORMAT (5X,1P,6E11.2) 99998 FORMAT (1X,I8) END