* F08CVF Example Program Text * Mark 21 Release. NAG Copyright 2004. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER MMAX, NB, NMAX PARAMETER (MMAX=8,NB=64,NMAX=8) INTEGER LDA, LWORK PARAMETER (LDA=MMAX,LWORK=NB*MMAX) COMPLEX *16 ZERO PARAMETER (ZERO=(0.0D0,0.0D0)) * .. Local Scalars .. INTEGER I, INFO, J, M, N * .. Local Arrays .. COMPLEX *16 A(LDA,NMAX), B(MMAX), TAU(MMAX), WORK(LWORK), + X(NMAX) * .. External Subroutines .. EXTERNAL F06HBF, ZCOPY, ZGERQF, ZTRTRS, ZUNMRQ * .. Executable Statements .. WRITE (NOUT,*) 'F08CVF Example Program Results' WRITE (NOUT,*) * Skip heading in data file READ (NIN,*) READ (NIN,*) M, N IF (M.LE.MMAX .AND. N.LE.NMAX .AND. M.LE.N) THEN * * Read the matrix A and the vector b from data file * READ (NIN,*) ((A(I,J),J=1,N),I=1,M) READ (NIN,*) (B(I),I=1,M) * * Compute the RQ factorization of A * CALL ZGERQF(M,N,A,LDA,TAU,WORK,LWORK,INFO) * * Copy the m-element vector b into x2, where x2 is the vector * containing the elements x(n-m+1), ..., x(n) of x * CALL ZCOPY(M,B,1,X(N-M+1),1) * * Solve R*y2 = b, storing the result in x2 * CALL ZTRTRS('Upper','No transpose','Non-Unit',M,1,A(1,N-M+1), + LDA,X(N-M+1),M,INFO) * IF (INFO.GT.0) THEN WRITE (NOUT,*) + 'The upper triangular factor, R, of A is singular, ' WRITE (NOUT,*) + 'the least squares solution could not be computed' GO TO 20 END IF * * Set y1 to zero (stored in rows 1 to (N-M) of x) * CALL F06HBF(N-M,ZERO,X,1) * * Compute minimum-norm solution x = (Q**H)*y * CALL ZUNMRQ('Left','Conjugate 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(I),I=1,N) ELSE WRITE (NOUT,*) + 'One or both of MMAX and NMAX is too small, and/or M.GT.N' END IF 20 CONTINUE * 99999 FORMAT (4(' (',F8.4,',',F8.4,')',:)) END