cdrvls.f
来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 626 行 · 第 1/2 页
F
626 行
NFAIL = NFAIL + 1
END IF
20 CONTINUE
NRUN = NRUN + 2
30 CONTINUE
40 CONTINUE
END IF
*
* Generate a matrix of scaling type ISCALE and rank
* type IRANK.
*
CALL CQRT15( 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)
*
DO 50 J = 1, N
IWORK( J ) = 0
50 CONTINUE
LDWORK = MAX( 1, M )
*
* Test CGELSX
*
* CGELSX: Compute the minimum-norm solution X
* to min( norm( A * X - B ) )
* using a complete orthogonal factorization.
*
CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B, LDB )
*
SRNAMT = 'CGELSX'
CALL CGELSX( M, N, NRHS, A, LDA, B, LDB, IWORK,
$ RCOND, CRANK, WORK, RWORK, INFO )
*
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'CGELSX', 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 ) = CQRT12( CRANK, CRANK, A, LDA, COPYS,
$ WORK, LWORK, RWORK )
*
* Test 4: Compute error in solution
* workspace: M*NRHS + M
*
CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK )
CALL CQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, RWORK,
$ RESULT( 4 ) )
*
* Test 5: Check norm of r'*A
* workspace: NRHS*(M+N)
*
RESULT( 5 ) = ZERO
IF( M.GT.CRANK )
$ RESULT( 5 ) = CQRT17( '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 ) = CQRT14( '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, 0,
$ ITYPE, K, RESULT( K )
NFAIL = NFAIL + 1
END IF
60 CONTINUE
NRUN = NRUN + 4
*
* Loop for testing different block sizes.
*
DO 90 INB = 1, NNB
NB = NBVAL( INB )
CALL XLAENV( 1, NB )
CALL XLAENV( 3, NXVAL( INB ) )
*
* Test CGELSY
*
* CGELSY: Compute the minimum-norm solution
* X to min( norm( A * X - B ) )
* using the rank-revealing orthogonal
* factorization.
*
CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B,
$ LDB )
*
* Initialize vector IWORK.
*
DO 70 J = 1, N
IWORK( J ) = 0
70 CONTINUE
*
* Set LWLSY to the adequate value.
*
LWLSY = MNMIN + MAX( 2*MNMIN, NB*( N+1 ),
$ MNMIN+NB*NRHS )
LWLSY = MAX( 1, LWLSY )
*
SRNAMT = 'CGELSY'
CALL CGELSY( M, N, NRHS, A, LDA, B, LDB, IWORK,
$ RCOND, CRANK, WORK, LWLSY, RWORK,
$ INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'CGELSY', INFO, 0, ' ', M,
$ N, NRHS, -1, NB, ITYPE, NFAIL,
$ NERRS, NOUT )
*
* workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS)
*
* Test 7: Compute relative error in svd
* workspace: M*N + 4*MIN(M,N) + MAX(M,N)
*
RESULT( 7 ) = CQRT12( CRANK, CRANK, A, LDA,
$ COPYS, WORK, LWORK, RWORK )
*
* Test 8: Compute error in solution
* workspace: M*NRHS + M
*
CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK )
CALL CQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, RWORK,
$ RESULT( 8 ) )
*
* Test 9: Check norm of r'*A
* workspace: NRHS*(M+N)
*
RESULT( 9 ) = ZERO
IF( M.GT.CRANK )
$ RESULT( 9 ) = CQRT17( '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 ) = CQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK )
*
* Test CGELSS
*
* CGELSS: Compute the minimum-norm solution
* X to min( norm( A * X - B ) )
* using the SVD.
*
CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B,
$ LDB )
SRNAMT = 'CGELSS'
CALL CGELSS( M, N, NRHS, A, LDA, B, LDB, S,
$ RCOND, CRANK, WORK, LWORK, RWORK,
$ INFO )
*
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'CGELSS', 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 CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK )
CALL CQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, RWORK,
$ RESULT( 12 ) )
*
* Test 13: Check norm of r'*A
*
RESULT( 13 ) = ZERO
IF( M.GT.CRANK )
$ RESULT( 13 ) = CQRT17( '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 ) = CQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK )
*
* Test CGELSD
*
* CGELSD: Compute the minimum-norm solution X
* to min( norm( A * X - B ) ) using a
* divide and conquer SVD.
*
CALL XLAENV( 9, 25 )
*
CALL CLACPY( 'Full', M, N, COPYA, LDA, A, LDA )
CALL CLACPY( 'Full', M, NRHS, COPYB, LDB, B,
$ LDB )
*
SRNAMT = 'CGELSD'
CALL CGELSD( M, N, NRHS, A, LDA, B, LDB, S,
$ RCOND, CRANK, WORK, LWORK, RWORK,
$ IWORK, INFO )
IF( INFO.NE.0 )
$ CALL ALAERH( PATH, 'CGELSD', 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 CLACPY( 'Full', M, NRHS, COPYB, LDB, WORK,
$ LDWORK )
CALL CQRT16( 'No transpose', M, N, NRHS, COPYA,
$ LDA, B, LDB, WORK, LDWORK, RWORK,
$ RESULT( 16 ) )
*
* Test 17: Check norm of r'*A
*
RESULT( 17 ) = ZERO
IF( M.GT.CRANK )
$ RESULT( 17 ) = CQRT17( '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 ) = CQRT14( 'No transpose', M, N,
$ NRHS, COPYA, LDA, B, LDB,
$ WORK, LWORK )
*
* Print information about the tests that did not
* pass the threshold.
*
DO 80 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
80 CONTINUE
NRUN = NRUN + 12
*
90 CONTINUE
100 CONTINUE
110 CONTINUE
120 CONTINUE
130 CONTINUE
140 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 CDRVLS
*
END
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?