sdrvls.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 621 行 · 第 1/2 页
F
621 行
20 CONTINUE
NRUN = NRUN + 2
30 CONTINUE
40 CONTINUE
END IF
*
* Generate a matrix of scaling type ISCALE and rank
* type IRANK.
*
CALL SQRT15( ISCALE, IRANK, M, N, NRHS, COPYA, LDA,
$ COPYB, LDB, COPYS, RANK, NORMA, NORMB,
$ ISEED, WORK, LWORK )
*
* workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
*
* Initialize vector IWORK.
*
DO 50 J = 1, N
IWORK( J ) = 0
50 CONTINUE
LDWORK = MAX( 1, M )
*
* Test SGELSX
*
* SGELSX: Compute the minimum-norm solution X
* to min( norm( A * X - B ) ) using a complete
* orthogonal factorization.
*
CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B, LDB )
*
SRNAMT = 'SGELSX'
CALL SGELSX( M, N, NRHS, A, LDA, B, LDB, IWORK,
$ RCOND, CRANK, WORK, INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'SGELSX', INFO, 0, ' ', M, N,
$ NRHS, -1, NB, ITYPE, NFAIL, NERRS,
$ NOUT )
*
* workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS )
*
* Test 3: Compute relative error in svd
* workspace: M*N + 4*MIN(M,N) + MAX(M,N)
*
RESULT( 3 ) = SQRT12( CRANK, CRANK, A, LDA, COPYS,
$ WORK, LWORK )
*
* Test 4: Compute error in solution
* workspace: M*NRHS + M
*
CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK )
CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK,
$ WORK( M*NRHS+1 ), RESULT( 4 ) )
*
* Test 5: Check norm of r'*A
* workspace: NRHS*(M+N)
*
RESULT( 5 ) = ZERO
IF( M.GT.CRANK )
$ RESULT( 5 ) = SQRT17( 'No transpose', 1, M, N,
$ NRHS, COPYA, LDA, B, LDB, COPYB,
$ LDB, C, WORK, LWORK )
*
* Test 6: Check if x is in the rowspace of A
* workspace: (M+NRHS)*(N+2)
*
RESULT( 6 ) = ZERO
*
IF( N.GT.CRANK )
$ RESULT( 6 ) = SQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB, WORK,
$ LWORK )
*
* Print information about the tests that did not
* pass the threshold.
*
DO 60 K = 3, 6
IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
$ ITYPE, K, RESULT( K )
NFAIL = NFAIL + 1
END IF
60 CONTINUE
NRUN = NRUN + 4
*
* Loop for testing different block sizes.
*
DO 100 INB = 1, NNB
NB = NBVAL( INB )
CALL XLAENV( 1, NB )
CALL XLAENV( 3, NXVAL( INB ) )
*
* Test SGELSY
*
* SGELSY: Compute the minimum-norm solution X
* to min( norm( A * X - B ) )
* using the rank-revealing orthogonal
* factorization.
*
* Initialize vector IWORK.
*
DO 70 J = 1, N
IWORK( J ) = 0
70 CONTINUE
*
* Set LWLSY to the adequate value.
*
LWLSY = MAX( 1, MNMIN+2*N+NB*( N+1 ),
$ 2*MNMIN+NB*NRHS )
*
CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
$ LDB )
*
SRNAMT = 'SGELSY'
CALL SGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
$ RCOND, CRANK, WORK, LWLSY, INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'SGELSY', INFO, 0, ' ', M,
$ N, NRHS, -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT )
*
* Test 7: Compute relative error in svd
* workspace: M*N + 4*MIN(M,N) + MAX(M,N)
*
RESULT( 7 ) = SQRT12( CRANK, CRANK, A, LDA,
$ COPYS, WORK, LWORK )
*
* Test 8: Compute error in solution
* workspace: M*NRHS + M
*
CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK )
CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK,
$ WORK( M*NRHS+1 ), RESULT( 8 ) )
*
* Test 9: Check norm of r'*A
* workspace: NRHS*(M+N)
*
RESULT( 9 ) = ZERO
IF( M.GT.CRANK )
$ RESULT( 9 ) = SQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK )
*
* Test 10: Check if x is in the rowspace of A
* workspace: (M+NRHS)*(N+2)
*
RESULT( 10 ) = ZERO
*
IF( N.GT.CRANK )
$ RESULT( 10 ) = SQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK )
*
* Test SGELSS
*
* SGELSS: Compute the minimum-norm solution X
* to min( norm( A * X - B ) )
* using the SVD.
*
CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
$ LDB )
SRNAMT = 'SGELSS'
CALL SGELSS( M, N, NRHS, A, LDA, B, LDB, S,
$ RCOND, CRANK, WORK, LWORK, INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'SGELSS', INFO, 0, ' ', M,
$ N, NRHS, -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT )
*
* workspace used: 3*min(m,n) +
* max(2*min(m,n),nrhs,max(m,n))
*
* Test 11: Compute relative error in svd
*
IF( RANK.GT.0 ) THEN
CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
RESULT( 11 ) = SASUM( MNMIN, S, 1 ) /
$ SASUM( MNMIN, COPYS, 1 ) /
$ ( EPS*REAL( MNMIN ) )
ELSE
RESULT( 11 ) = ZERO
END IF
*
* Test 12: Compute error in solution
*
CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK )
CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK,
$ WORK( M*NRHS+1 ), RESULT( 12 ) )
*
* Test 13: Check norm of r'*A
*
RESULT( 13 ) = ZERO
IF( M.GT.CRANK )
$ RESULT( 13 ) = SQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK )
*
* Test 14: Check if x is in the rowspace of A
*
RESULT( 14 ) = ZERO
IF( N.GT.CRANK )
$ RESULT( 14 ) = SQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK )
*
* Test SGELSD
*
* SGELSD: Compute the minimum-norm solution X
* to min( norm( A * X - B ) ) using a
* divide and conquer SVD.
*
* Initialize vector IWORK.
*
DO 80 J = 1, N
IWORK( J ) = 0
80 CONTINUE
*
CALL SLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, B,
$ LDB )
*
SRNAMT = 'SGELSD'
CALL SGELSD( M, N, NRHS, A, LDA, B, LDB, S,
$ RCOND, CRANK, WORK, LWORK, IWORK,
$ INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'SGELSD', INFO, 0, ' ', M,
$ N, NRHS, -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT )
*
* Test 15: Compute relative error in svd
*
IF( RANK.GT.0 ) THEN
CALL SAXPY( MNMIN, -ONE, COPYS, 1, S, 1 )
RESULT( 15 ) = SASUM( MNMIN, S, 1 ) /
$ SASUM( MNMIN, COPYS, 1 ) /
$ ( EPS*REAL( MNMIN ) )
ELSE
RESULT( 15 ) = ZERO
END IF
*
* Test 16: Compute error in solution
*
CALL SLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK )
CALL SQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK,
$ WORK( M*NRHS+1 ), RESULT( 16 ) )
*
* Test 17: Check norm of r'*A
*
RESULT( 17 ) = ZERO
IF( M.GT.CRANK )
$ RESULT( 17 ) = SQRT17( 'No transpose', 1, M,
$ N, NRHS, COPYA, LDA, B, LDB,
$ COPYB, LDB, C, WORK, LWORK )
*
* Test 18: Check if x is in the rowspace of A
*
RESULT( 18 ) = ZERO
IF( N.GT.CRANK )
$ RESULT( 18 ) = SQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK )
*
* Print information about the tests that did not
* pass the threshold.
*
DO 90 K = 7, NTESTS
IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
WRITE( NOUT, FMT = 9998 )M, N, NRHS, NB,
$ ITYPE, K, RESULT( K )
NFAIL = NFAIL + 1
END IF
90 CONTINUE
NRUN = NRUN + 12
*
100 CONTINUE
110 CONTINUE
120 CONTINUE
130 CONTINUE
140 CONTINUE
150 CONTINUE
*
* Print a summary of the results.
*
CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
9999 FORMAT( ' TRANS=''', A1, ''', M=', I5, ', N=', I5, ', NRHS=', I4,
$ ', NB=', I4, ', type', I2, ', test(', I2, ')=', G12.5 )
9998 FORMAT( ' M=', I5, ', N=', I5, ', NRHS=', I4, ', NB=', I4,
$ ', type', I2, ', test(', I2, ')=', G12.5 )
RETURN
*
* End of SDRVLS
*
END
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?