* F08KDF 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=10) INTEGER LDA, LDU, LWORK PARAMETER (LDA=MMAX,LDU=MMAX,LWORK=(5*MMAX+9) + *MMAX+NMAX+NB*(MMAX+NMAX)) * .. Local Scalars .. DOUBLE PRECISION EPS, SERRBD INTEGER I, IFAIL, INFO, J, LWKOPT, M, N * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), DUMMY(1,1), RCONDU(MMAX), + RCONDV(MMAX), S(MMAX), U(LDU,MMAX), UERRBD(MMAX), + VERRBD(MMAX), WORK(LWORK) INTEGER IWORK(8*MMAX) * .. External Functions .. DOUBLE PRECISION X02AJF EXTERNAL X02AJF * .. External Subroutines .. EXTERNAL DDISNA, DGESDD, X04CAF * .. Executable Statements .. WRITE (NOUT,*) 'F08KDF 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**T), m.le.n) * CALL DGESDD('Overwrite A by tranpose(V)',M,N,A,LDA,S,U,LDU, + DUMMY,1,WORK,LWORK,IWORK,INFO) LWKOPT = WORK(1) * IF (INFO.EQ.0) THEN * * Print solution * WRITE (NOUT,*) 'Singular values' WRITE (NOUT,99999) (S(J),J=1,M) * IFAIL = 0 CALL X04CAF('General',' ',M,M,U,LDU,'Left singular vectors', + IFAIL) WRITE (NOUT,*) CALL X04CAF('General',' ',M,N,A,LDA, + 'Right singular vectors by row '// + '(first m rows of V**T)',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, M 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,M) WRITE (NOUT,*) WRITE (NOUT,*) + 'Error estimates for the right singular vectors' WRITE (NOUT,99998) (VERRBD(I),I=1,M) ELSE WRITE (NOUT,99997) 'Failure in DGESDD. 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