ddrvbd.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 663 行 · 第 1/2 页
F
663 行
IF( JTYPE.EQ.1 ) THEN
*
* Zero matrix
*
CALL DLASET( 'Full', M, N, ZERO, ZERO, A, LDA )
*
ELSE IF( JTYPE.EQ.2 ) THEN
*
* Identity matrix
*
CALL DLASET( 'Full', M, N, ZERO, ONE, A, LDA )
*
ELSE
*
* (Scaled) random matrix
*
IF( JTYPE.EQ.3 )
$ ANORM = ONE
IF( JTYPE.EQ.4 )
$ ANORM = UNFL / ULP
IF( JTYPE.EQ.5 )
$ ANORM = OVFL*ULP
CALL DLATMS( M, N, 'U', ISEED, 'N', S, 4, DBLE( MNMIN ),
$ ANORM, M-1, N-1, 'N', A, LDA, WORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUT, FMT = 9996 )'Generator', IINFO, M, N,
$ JTYPE, IOLDSD
INFO = ABS( IINFO )
RETURN
END IF
END IF
*
30 CONTINUE
CALL DLACPY( 'F', M, N, A, LDA, ASAV, LDA )
*
* Do for minimal and adequate (for blocking) workspace
*
DO 130 IWS = 1, 4
*
DO 40 J = 1, 14
RESULT( J ) = -ONE
40 CONTINUE
*
* Test DGESVD: Factorize A
*
IWTMP = MAX( 3*MIN( M, N )+MAX( M, N ), 5*MIN( M, N ) )
LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
IF( IWS.EQ.4 )
$ LSWORK = LWORK
*
IF( IWS.GT.1 )
$ CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
SRNAMT = 'DGESVD'
CALL DGESVD( 'A', 'A', M, N, A, LDA, SSAV, USAV, LDU,
$ VTSAV, LDVT, WORK, LSWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUT, FMT = 9995 )'GESVD', IINFO, M, N, JTYPE,
$ LSWORK, IOLDSD
INFO = ABS( IINFO )
RETURN
END IF
*
* Do tests 1--4
*
CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
$ VTSAV, LDVT, WORK, RESULT( 1 ) )
IF( M.NE.0 .AND. N.NE.0 ) THEN
CALL DORT01( 'Columns', M, M, USAV, LDU, WORK, LWORK,
$ RESULT( 2 ) )
CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK, LWORK,
$ RESULT( 3 ) )
END IF
RESULT( 4 ) = ZERO
DO 50 I = 1, MNMIN - 1
IF( SSAV( I ).LT.SSAV( I+1 ) )
$ RESULT( 4 ) = ULPINV
IF( SSAV( I ).LT.ZERO )
$ RESULT( 4 ) = ULPINV
50 CONTINUE
IF( MNMIN.GE.1 ) THEN
IF( SSAV( MNMIN ).LT.ZERO )
$ RESULT( 4 ) = ULPINV
END IF
*
* Do partial SVDs, comparing to SSAV, USAV, and VTSAV
*
RESULT( 5 ) = ZERO
RESULT( 6 ) = ZERO
RESULT( 7 ) = ZERO
DO 80 IJU = 0, 3
DO 70 IJVT = 0, 3
IF( ( IJU.EQ.3 .AND. IJVT.EQ.3 ) .OR.
$ ( IJU.EQ.1 .AND. IJVT.EQ.1 ) )GO TO 70
JOBU = CJOB( IJU+1 )
JOBVT = CJOB( IJVT+1 )
CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
SRNAMT = 'DGESVD'
CALL DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU,
$ VT, LDVT, WORK, LSWORK, IINFO )
*
* Compare U
*
DIF = ZERO
IF( M.GT.0 .AND. N.GT.0 ) THEN
IF( IJU.EQ.1 ) THEN
CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
$ LDU, A, LDA, WORK, LWORK, DIF,
$ IINFO )
ELSE IF( IJU.EQ.2 ) THEN
CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
$ LDU, U, LDU, WORK, LWORK, DIF,
$ IINFO )
ELSE IF( IJU.EQ.3 ) THEN
CALL DORT03( 'C', M, M, M, MNMIN, USAV, LDU,
$ U, LDU, WORK, LWORK, DIF,
$ IINFO )
END IF
END IF
RESULT( 5 ) = MAX( RESULT( 5 ), DIF )
*
* Compare VT
*
DIF = ZERO
IF( M.GT.0 .AND. N.GT.0 ) THEN
IF( IJVT.EQ.1 ) THEN
CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
$ LDVT, A, LDA, WORK, LWORK, DIF,
$ IINFO )
ELSE IF( IJVT.EQ.2 ) THEN
CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
$ LDVT, VT, LDVT, WORK, LWORK,
$ DIF, IINFO )
ELSE IF( IJVT.EQ.3 ) THEN
CALL DORT03( 'R', N, N, N, MNMIN, VTSAV,
$ LDVT, VT, LDVT, WORK, LWORK,
$ DIF, IINFO )
END IF
END IF
RESULT( 6 ) = MAX( RESULT( 6 ), DIF )
*
* Compare S
*
DIF = ZERO
DIV = MAX( DBLE( MNMIN )*ULP*S( 1 ), UNFL )
DO 60 I = 1, MNMIN - 1
IF( SSAV( I ).LT.SSAV( I+1 ) )
$ DIF = ULPINV
IF( SSAV( I ).LT.ZERO )
$ DIF = ULPINV
DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
60 CONTINUE
RESULT( 7 ) = MAX( RESULT( 7 ), DIF )
70 CONTINUE
80 CONTINUE
*
* Test DGESDD: Factorize A
*
IWTMP = 5*MNMIN*MNMIN + 9*MNMIN + MAX( M, N )
LSWORK = IWTMP + ( IWS-1 )*( LWORK-IWTMP ) / 3
LSWORK = MIN( LSWORK, LWORK )
LSWORK = MAX( LSWORK, 1 )
IF( IWS.EQ.4 )
$ LSWORK = LWORK
*
CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
SRNAMT = 'DGESDD'
CALL DGESDD( 'A', M, N, A, LDA, SSAV, USAV, LDU, VTSAV,
$ LDVT, WORK, LSWORK, IWORK, IINFO )
IF( IINFO.NE.0 ) THEN
WRITE( NOUT, FMT = 9995 )'GESDD', IINFO, M, N, JTYPE,
$ LSWORK, IOLDSD
INFO = ABS( IINFO )
RETURN
END IF
*
* Do tests 8--11
*
CALL DBDT01( M, N, 0, ASAV, LDA, USAV, LDU, SSAV, E,
$ VTSAV, LDVT, WORK, RESULT( 8 ) )
IF( M.NE.0 .AND. N.NE.0 ) THEN
CALL DORT01( 'Columns', M, M, USAV, LDU, WORK, LWORK,
$ RESULT( 9 ) )
CALL DORT01( 'Rows', N, N, VTSAV, LDVT, WORK, LWORK,
$ RESULT( 10 ) )
END IF
RESULT( 11 ) = ZERO
DO 90 I = 1, MNMIN - 1
IF( SSAV( I ).LT.SSAV( I+1 ) )
$ RESULT( 11 ) = ULPINV
IF( SSAV( I ).LT.ZERO )
$ RESULT( 11 ) = ULPINV
90 CONTINUE
IF( MNMIN.GE.1 ) THEN
IF( SSAV( MNMIN ).LT.ZERO )
$ RESULT( 11 ) = ULPINV
END IF
*
* Do partial SVDs, comparing to SSAV, USAV, and VTSAV
*
RESULT( 12 ) = ZERO
RESULT( 13 ) = ZERO
RESULT( 14 ) = ZERO
DO 110 IJQ = 0, 2
JOBQ = CJOB( IJQ+1 )
CALL DLACPY( 'F', M, N, ASAV, LDA, A, LDA )
SRNAMT = 'DGESDD'
CALL DGESDD( JOBQ, M, N, A, LDA, S, U, LDU, VT, LDVT,
$ WORK, LSWORK, IWORK, IINFO )
*
* Compare U
*
DIF = ZERO
IF( M.GT.0 .AND. N.GT.0 ) THEN
IF( IJQ.EQ.1 ) THEN
IF( M.GE.N ) THEN
CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
$ LDU, A, LDA, WORK, LWORK, DIF,
$ INFO )
ELSE
CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV,
$ LDU, U, LDU, WORK, LWORK, DIF,
$ INFO )
END IF
ELSE IF( IJQ.EQ.2 ) THEN
CALL DORT03( 'C', M, MNMIN, M, MNMIN, USAV, LDU,
$ U, LDU, WORK, LWORK, DIF, INFO )
END IF
END IF
RESULT( 12 ) = MAX( RESULT( 12 ), DIF )
*
* Compare VT
*
DIF = ZERO
IF( M.GT.0 .AND. N.GT.0 ) THEN
IF( IJQ.EQ.1 ) THEN
IF( M.GE.N ) THEN
CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
$ LDVT, VT, LDVT, WORK, LWORK,
$ DIF, INFO )
ELSE
CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
$ LDVT, A, LDA, WORK, LWORK, DIF,
$ INFO )
END IF
ELSE IF( IJQ.EQ.2 ) THEN
CALL DORT03( 'R', N, MNMIN, N, MNMIN, VTSAV,
$ LDVT, VT, LDVT, WORK, LWORK, DIF,
$ INFO )
END IF
END IF
RESULT( 13 ) = MAX( RESULT( 13 ), DIF )
*
* Compare S
*
DIF = ZERO
DIV = MAX( DBLE( MNMIN )*ULP*S( 1 ), UNFL )
DO 100 I = 1, MNMIN - 1
IF( SSAV( I ).LT.SSAV( I+1 ) )
$ DIF = ULPINV
IF( SSAV( I ).LT.ZERO )
$ DIF = ULPINV
DIF = MAX( DIF, ABS( SSAV( I )-S( I ) ) / DIV )
100 CONTINUE
RESULT( 14 ) = MAX( RESULT( 14 ), DIF )
110 CONTINUE
*
* End of Loop -- Check for RESULT(j) > THRESH
*
DO 120 J = 1, 14
IF( RESULT( J ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 ) THEN
WRITE( NOUT, FMT = 9999 )
WRITE( NOUT, FMT = 9998 )
END IF
WRITE( NOUT, FMT = 9997 )M, N, JTYPE, IWS, IOLDSD,
$ J, RESULT( J )
NFAIL = NFAIL + 1
END IF
120 CONTINUE
NTEST = NTEST + 14
*
130 CONTINUE
140 CONTINUE
150 CONTINUE
*
* Summary
*
CALL ALASVM( PATH, NOUT, NFAIL, NTEST, 0 )
*
9999 FORMAT( ' SVD -- Real Singular Value Decomposition Driver ',
$ / ' Matrix types (see DDRVBD for details):',
$ / / ' 1 = Zero matrix', / ' 2 = Identity matrix',
$ / ' 3 = Evenly spaced singular values near 1',
$ / ' 4 = Evenly spaced singular values near underflow',
$ / ' 5 = Evenly spaced singular values near overflow', / /
$ ' Tests performed: ( A is dense, U and V are orthogonal,',
$ / 19X, ' S is an array, and Upartial, VTpartial, and',
$ / 19X, ' Spartial are partially computed U, VT and S),', / )
9998 FORMAT( ' 1 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
$ / ' 2 = | I - U**T U | / ( M ulp ) ',
$ / ' 3 = | I - VT VT**T | / ( N ulp ) ',
$ / ' 4 = 0 if S contains min(M,N) nonnegative values in',
$ ' decreasing order, else 1/ulp',
$ / ' 5 = | U - Upartial | / ( M ulp )',
$ / ' 6 = | VT - VTpartial | / ( N ulp )',
$ / ' 7 = | S - Spartial | / ( min(M,N) ulp |S| )',
$ / ' 8 = | A - U diag(S) VT | / ( |A| max(M,N) ulp ) ',
$ / ' 9 = | I - U**T U | / ( M ulp ) ',
$ / '10 = | I - VT VT**T | / ( N ulp ) ',
$ / '11 = 0 if S contains min(M,N) nonnegative values in',
$ ' decreasing order, else 1/ulp',
$ / '12 = | U - Upartial | / ( M ulp )',
$ / '13 = | VT - VTpartial | / ( N ulp )',
$ / '14 = | S - Spartial | / ( min(M,N) ulp |S| )', / / )
9997 FORMAT( ' M=', I5, ', N=', I5, ', type ', I1, ', IWS=', I1,
$ ', seed=', 4( I4, ',' ), ' test(', I2, ')=', G11.4 )
9996 FORMAT( ' DDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
$ I6, ', N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
$ I5, ')' )
9995 FORMAT( ' DDRVBD: ', A, ' returned INFO=', I6, '.', / 9X, 'M=',
$ I6, ', N=', I6, ', JTYPE=', I6, ', LSWORK=', I6, / 9X,
$ 'ISEED=(', 3( I5, ',' ), I5, ')' )
*
RETURN
*
* End of DDRVBD
*
END
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?