sdrvpt.f

来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 491 行 · 第 1/2 页

                  Z( 2 ) = D( 1 )
                  D( 1 ) = ZERO
                  IF( N.GT.1 ) THEN
                     Z( 3 ) = E( 1 )
                     E( 1 ) = ZERO
                  END IF
               ELSE IF( IMAT.EQ.9 ) THEN
                  IZERO = N
                  IF( N.GT.1 ) THEN
                     Z( 1 ) = E( N-1 )
                     E( N-1 ) = ZERO
                  END IF
                  Z( 2 ) = D( N )
                  D( N ) = ZERO
               ELSE IF( IMAT.EQ.10 ) THEN
                  IZERO = ( N+1 ) / 2
                  IF( IZERO.GT.1 ) THEN
                     Z( 1 ) = E( IZERO-1 )
                     Z( 3 ) = E( IZERO )
                     E( IZERO-1 ) = ZERO
                     E( IZERO ) = ZERO
                  END IF
                  Z( 2 ) = D( IZERO )
                  D( IZERO ) = ZERO
               END IF
            END IF
*
*           Generate NRHS random solution vectors.
*
            IX = 1
            DO 40 J = 1, NRHS
               CALL SLARNV( 2, ISEED, N, XACT( IX ) )
               IX = IX + LDA
   40       CONTINUE
*
*           Set the right hand side.
*
            CALL SLAPTM( N, NRHS, ONE, D, E, XACT, LDA, ZERO, B, LDA )
*
            DO 100 IFACT = 1, 2
               IF( IFACT.EQ.1 ) THEN
                  FACT = 'F'
               ELSE
                  FACT = 'N'
               END IF
*
*              Compute the condition number for comparison with
*              the value returned by SPTSVX.
*
               IF( ZEROT ) THEN
                  IF( IFACT.EQ.1 )
     $               GO TO 100
                  RCONDC = ZERO
*
               ELSE IF( IFACT.EQ.1 ) THEN
*
*                 Compute the 1-norm of A.
*
                  ANORM = SLANST( '1', N, D, E )
*
                  CALL SCOPY( N, D, 1, D( N+1 ), 1 )
                  IF( N.GT.1 )
     $               CALL SCOPY( N-1, E, 1, E( N+1 ), 1 )
*
*                 Factor the matrix A.
*
                  CALL SPTTRF( N, D( N+1 ), E( N+1 ), INFO )
*
*                 Use SPTTRS to solve for one column at a time of
*                 inv(A), computing the maximum column sum as we go.
*
                  AINVNM = ZERO
                  DO 60 I = 1, N
                     DO 50 J = 1, N
                        X( J ) = ZERO
   50                CONTINUE
                     X( I ) = ONE
                     CALL SPTTRS( N, 1, D( N+1 ), E( N+1 ), X, LDA,
     $                            INFO )
                     AINVNM = MAX( AINVNM, SASUM( N, X, 1 ) )
   60             CONTINUE
*
*                 Compute the 1-norm condition number of A.
*
                  IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
                     RCONDC = ONE
                  ELSE
                     RCONDC = ( ONE / ANORM ) / AINVNM
                  END IF
               END IF
*
               IF( IFACT.EQ.2 ) THEN
*
*                 --- Test SPTSV --
*
                  CALL SCOPY( N, D, 1, D( N+1 ), 1 )
                  IF( N.GT.1 )
     $               CALL SCOPY( N-1, E, 1, E( N+1 ), 1 )
                  CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
*
*                 Factor A as L*D*L' and solve the system A*X = B.
*
                  SRNAMT = 'SPTSV '
                  CALL SPTSV( N, NRHS, D( N+1 ), E( N+1 ), X, LDA,
     $                        INFO )
*
*                 Check error code from SPTSV .
*
                  IF( INFO.NE.IZERO )
     $               CALL ALAERH( PATH, 'SPTSV ', INFO, IZERO, ' ', N,
     $                            N, 1, 1, NRHS, IMAT, NFAIL, NERRS,
     $                            NOUT )
                  NT = 0
                  IF( IZERO.EQ.0 ) THEN
*
*                    Check the factorization by computing the ratio
*                       norm(L*D*L' - A) / (n * norm(A) * EPS )
*
                     CALL SPTT01( N, D, E, D( N+1 ), E( N+1 ), WORK,
     $                            RESULT( 1 ) )
*
*                    Compute the residual in the solution.
*
                     CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
                     CALL SPTT02( N, NRHS, D, E, X, LDA, WORK, LDA,
     $                            RESULT( 2 ) )
*
*                    Check solution from generated exact solution.
*
                     CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
     $                            RESULT( 3 ) )
                     NT = 3
                  END IF
*
*                 Print information about the tests that did not pass
*                 the threshold.
*
                  DO 70 K = 1, NT
                     IF( RESULT( K ).GE.THRESH ) THEN
                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                     CALL ALADHD( NOUT, PATH )
                        WRITE( NOUT, FMT = 9999 )'SPTSV ', N, IMAT, K,
     $                     RESULT( K )
                        NFAIL = NFAIL + 1
                     END IF
   70             CONTINUE
                  NRUN = NRUN + NT
               END IF
*
*              --- Test SPTSVX ---
*
               IF( IFACT.GT.1 ) THEN
*
*                 Initialize D( N+1:2*N ) and E( N+1:2*N ) to zero.
*
                  DO 80 I = 1, N - 1
                     D( N+I ) = ZERO
                     E( N+I ) = ZERO
   80             CONTINUE
                  IF( N.GT.0 )
     $               D( N+N ) = ZERO
               END IF
*
               CALL SLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDA )
*
*              Solve the system and compute the condition number and
*              error bounds using SPTSVX.
*
               SRNAMT = 'SPTSVX'
               CALL SPTSVX( FACT, N, NRHS, D, E, D( N+1 ), E( N+1 ), B,
     $                      LDA, X, LDA, RCOND, RWORK, RWORK( NRHS+1 ),
     $                      WORK, INFO )
*
*              Check the error code from SPTSVX.
*
               IF( INFO.NE.IZERO )
     $            CALL ALAERH( PATH, 'SPTSVX', INFO, IZERO, FACT, N, N,
     $                         1, 1, NRHS, IMAT, NFAIL, NERRS, NOUT )
               IF( IZERO.EQ.0 ) THEN
                  IF( IFACT.EQ.2 ) THEN
*
*                    Check the factorization by computing the ratio
*                       norm(L*D*L' - A) / (n * norm(A) * EPS )
*
                     K1 = 1
                     CALL SPTT01( N, D, E, D( N+1 ), E( N+1 ), WORK,
     $                            RESULT( 1 ) )
                  ELSE
                     K1 = 2
                  END IF
*
*                 Compute the residual in the solution.
*
                  CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
                  CALL SPTT02( N, NRHS, D, E, X, LDA, WORK, LDA,
     $                         RESULT( 2 ) )
*
*                 Check solution from generated exact solution.
*
                  CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
     $                         RESULT( 3 ) )
*
*                 Check error bounds from iterative refinement.
*
                  CALL SPTT05( N, NRHS, D, E, B, LDA, X, LDA, XACT, LDA,
     $                         RWORK, RWORK( NRHS+1 ), RESULT( 4 ) )
               ELSE
                  K1 = 6
               END IF
*
*              Check the reciprocal of the condition number.
*
               RESULT( 6 ) = SGET06( RCOND, RCONDC )
*
*              Print information about the tests that did not pass
*              the threshold.
*
               DO 90 K = K1, 6
                  IF( RESULT( K ).GE.THRESH ) THEN
                     IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                  CALL ALADHD( NOUT, PATH )
                     WRITE( NOUT, FMT = 9998 )'SPTSVX', FACT, N, IMAT,
     $                  K, RESULT( K )
                     NFAIL = NFAIL + 1
                  END IF
   90          CONTINUE
               NRUN = NRUN + 7 - K1
  100       CONTINUE
  110    CONTINUE
  120 CONTINUE
*
*     Print a summary of the results.
*
      CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
 9999 FORMAT( 1X, A6, ', N =', I5, ', type ', I2, ', test ', I2,
     $      ', ratio = ', G12.5 )
 9998 FORMAT( 1X, A6, ', FACT=''', A1, ''', N =', I5, ', type ', I2,
     $      ', test ', I2, ', ratio = ', G12.5 )
      RETURN
*
*     End of SDRVPT
*
      END