* F08BSF Example Program Text * Mark 16 Release. NAG Copyright 1992. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER MMAX, NMAX, LDA, LDB, LDX, NRHMAX, LWORK PARAMETER (MMAX=8,NMAX=8,LDA=MMAX,LDB=MMAX,LDX=MMAX, + NRHMAX=NMAX,LWORK=64*NMAX) COMPLEX*16 ZERO PARAMETER (ZERO=(0.0D0,0.0D0)) * .. Local Scalars .. DOUBLE PRECISION TOL INTEGER I, IFAIL, INFO, J, K, M, N, NRHS * .. Local Arrays .. COMPLEX*16 A(LDA,NMAX), B(LDB,NRHMAX), TAU(NMAX), + WORK(LWORK), X(LDX,NRHMAX) DOUBLE PRECISION RWORK(2*NMAX) INTEGER JPVT(NMAX) CHARACTER CLABS(1), RLABS(1) * .. External Subroutines .. EXTERNAL F06DBF, F06HBF, X04DBF, ZGEQPF, ZTRSV, ZUNMQR * .. Intrinsic Functions .. INTRINSIC ABS * .. Executable Statements .. WRITE (NOUT,*) 'F08BSF Example Program Results' * 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.NRHMAX) + 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 * CALL ZGEQPF(M,N,A,LDA,JPVT,TAU,WORK,RWORK,INFO) * * Choose TOL to reflect the relative accuracy of the input data * TOL = 0.01D0 * * Determine which columns of R to use * DO 20 K = 1, N IF (ABS(A(K,K)).LE.TOL*ABS(A(1,1))) GO TO 40 20 CONTINUE * * Compute C = (Q**H)*B, storing the result in B * 40 K = K - 1 * CALL ZUNMQR('Left','Conjugate Transpose',M,NRHS,K,A,LDA,TAU,B, + LDB,WORK,LWORK,INFO) * * Compute least-squares solution by backsubstitution in R*B = C * DO 60 I = 1, NRHS * CALL ZTRSV('Upper','No transpose','Non-Unit',K,A,LDA,B(1,I), + 1) * * Set the unused elements of the I-th solution vector to zero * CALL F06HBF(N-K,ZERO,B(K+1,I),1) * 60 CONTINUE * * Unscramble the least-squares solution stored in B * DO 100 I = 1, N DO 80 J = 1, NRHS X(JPVT(I),J) = B(I,J) 80 CONTINUE 100 CONTINUE * * Print least-squares solution * WRITE (NOUT,*) IFAIL = 0 * CALL X04DBF('General',' ',N,NRHS,X,LDX,'Bracketed','F7.4', + 'Least-squares solution','Integer',RLABS,'Integer', + CLABS,80,0,IFAIL) * END IF STOP END