* F08KPF Example Program Text * Mark 21 Release. NAG Copyright 2004. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER MMAX, NB, NMAX PARAMETER (MMAX=10,NB=64,NMAX=8) INTEGER LDA, LDVT, LWORK PARAMETER (LDA=MMAX,LDVT=NMAX, + LWORK=MMAX+3*NMAX+NB*(MMAX+NMAX)) * .. Local Scalars .. DOUBLE PRECISION EPS, SERRBD INTEGER I, IFAIL, INFO, J, LWKOPT, M, N * .. Local Arrays .. COMPLEX *16 A(LDA,NMAX), DUMMY(1,1), VT(LDVT,NMAX), + WORK(LWORK) DOUBLE PRECISION RCONDU(NMAX), RCONDV(NMAX), RWORK(5*NMAX), + S(NMAX), UERRBD(NMAX), VERRBD(NMAX) * .. External Functions .. DOUBLE PRECISION X02AJF EXTERNAL X02AJF * .. External Subroutines .. EXTERNAL DDISNA, X04DAF, ZGESVD * .. Executable Statements .. WRITE (NOUT,*) 'F08KPF Example Program Results' WRITE (NOUT,*) * Skip heading in data file READ (NIN,*) READ (NIN,*) M, N IF (M.LE.MMAX .AND. N.LE.NMAX) THEN * * Read the m by n matrix A from data file * READ (NIN,*) ((A(I,J),J=1,N),I=1,M) * * Compute the singular values and left and right singular vectors * of A (A = U*S*(V**H), m.ge.n) * CALL ZGESVD('Overwrite A by U','Singular vectors (V)',M,N,A, + LDA,S,DUMMY,1,VT,LDVT,WORK,LWORK,RWORK,INFO) LWKOPT = WORK(1) * IF (INFO.EQ.0) THEN * * Print solution * WRITE (NOUT,*) 'Singular values' WRITE (NOUT,99999) (S(J),J=1,N) * IFAIL = 0 CALL X04DAF('General',' ',M,N,A,LDA, + 'Left singular vectors (first n columns of U)', + IFAIL) WRITE (NOUT,*) CALL X04DAF('General',' ',N,N,VT,LDVT,'Conjugates of '// + 'right singular vectors by row (V**H)',IFAIL) * * Get the machine precision, EPS and compute the approximate * error bound for the computed singular values. Note that for * the 2-norm, S(1) = norm(A) * EPS = X02AJF() SERRBD = EPS*S(1) * * Call DDISNA (F08FLF) to estimate reciprocal condition * numbers for the singular vectors * CALL DDISNA('Left',M,N,S,RCONDU,INFO) CALL DDISNA('Right',M,N,S,RCONDV,INFO) * * Compute the error estimates for the singular vectors * DO 20 I = 1, N UERRBD(I) = SERRBD/RCONDU(I) VERRBD(I) = SERRBD/RCONDV(I) 20 CONTINUE * * Print the approximate error bounds for the singular values * and vectors * WRITE (NOUT,*) WRITE (NOUT,*) 'Error estimate for the singular values' WRITE (NOUT,99998) SERRBD WRITE (NOUT,*) WRITE (NOUT,*) + 'Error estimates for the left singular vectors' WRITE (NOUT,99998) (UERRBD(I),I=1,N) WRITE (NOUT,*) WRITE (NOUT,*) + 'Error estimates for the right singular vectors' WRITE (NOUT,99998) (VERRBD(I),I=1,N) ELSE WRITE (NOUT,99997) 'Failure in ZGESVD. INFO =', INFO END IF * * Print workspace information * IF (LWORK.LT.LWKOPT) THEN WRITE (NOUT,*) WRITE (NOUT,99996) 'Optimum workspace required = ', LWKOPT, + 'Workspace provided = ', LWORK END IF ELSE WRITE (NOUT,*) 'MMAX and/or NMAX too small' END IF STOP * 99999 FORMAT (3X,(8F8.4)) 99998 FORMAT (4X,1P,6E11.1) 99997 FORMAT (1X,A,I4) 99996 FORMAT (1X,A,I5,/1X,A,I5) END