* F08CSF 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=NB*NMAX) * .. Local Scalars .. INTEGER I, IFAIL, INFO, J, M, N, NRHS * .. Local Arrays .. COMPLEX *16 A(LDA,NMAX), B(LDB,NRHSMX), TAU(NMAX), + WORK(LWORK) DOUBLE PRECISION RNORM(NRHSMX) CHARACTER CLABS(1), RLABS(1) * .. External Functions .. DOUBLE PRECISION DZNRM2 EXTERNAL DZNRM2 * .. External Subroutines .. EXTERNAL X04DBF, ZGEQLF, ZTRTRS, ZUNMQL * .. Executable Statements .. WRITE (NOUT,*) 'F08CSF 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) * * Compute the QL factorization of A * CALL ZGEQLF(M,N,A,LDA,TAU,WORK,LWORK,INFO) * * Compute C = (C1) = (Q**H)*B, storing the result in B * (C2) * CALL ZUNMQL('Left','Conjugate Transpose',M,NRHS,N,A,LDA,TAU,B, + LDB,WORK,LWORK,INFO) * * Compute least-squares solutions by backsubstitution in * L*X = C2 * CALL ZTRTRS('Lower','No transpose','Non-Unit',N,NRHS,A(M-N+1,1) + ,LDA,B(M-N+1,1),LDB,INFO) * IF (INFO.GT.0) THEN WRITE (NOUT,*) + 'The lower triangular factor, L, of A is singular, ' WRITE (NOUT,*) + 'the least squares solution could not be computed' GO TO 40 END IF * * Print least-squares solution(s) * IFAIL = 0 CALL X04DBF('General',' ',N,NRHS,B(M-N+1,1),LDB,'Bracketed', + 'F7.4','Least-squares solution(s)','Integer',RLABS, + 'Integer',CLABS,80,0,IFAIL) * * Compute and print estimates of the square roots of the residual * sums of squares * DO 20 J = 1, NRHS RNORM(J) = DZNRM2(M-N,B(1,J),1) 20 CONTINUE * 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 or NRHSMX is too small, ', + 'and/or M.LT.N' END IF 40 CONTINUE * 99999 FORMAT (3X,1P,7E11.2) END