* F08CHF 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) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) * .. Local Scalars .. INTEGER I, INFO, J, M, N * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), B(MMAX), TAU(MMAX), WORK(LWORK), + X(NMAX) * .. External Subroutines .. EXTERNAL DCOPY, DGERQF, DORMRQ, DTRTRS, F06FBF * .. Executable Statements .. WRITE (NOUT,*) 'F08CHF 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 DGERQF(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 DCOPY(M,B,1,X(N-M+1),1) * * Solve R*y2 = b, storing the result in x2 * CALL DTRTRS('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 F06FBF(N-M,ZERO,X,1) * * Compute the minimum-norm solution x = (Q**T)*y * 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(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 (1X,8F9.4) END