00001 SUBROUTINE SPPT02( UPLO, N, NRHS, A, X, LDX, B, LDB, RWORK, 00002 $ RESID ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER LDB, LDX, N, NRHS 00011 REAL RESID 00012 * .. 00013 * .. Array Arguments .. 00014 REAL A( * ), B( LDB, * ), RWORK( * ), X( LDX, * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * SPPT02 computes the residual in the solution of a symmetric system 00021 * of linear equations A*x = b when packed storage is used for the 00022 * coefficient matrix. The ratio computed is 00023 * 00024 * RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS), 00025 * 00026 * where EPS is the machine precision. 00027 * 00028 * Arguments 00029 * ========= 00030 * 00031 * UPLO (input) CHARACTER*1 00032 * Specifies whether the upper or lower triangular part of the 00033 * symmetric matrix A is stored: 00034 * = 'U': Upper triangular 00035 * = 'L': Lower triangular 00036 * 00037 * N (input) INTEGER 00038 * The number of rows and columns of the matrix A. N >= 0. 00039 * 00040 * NRHS (input) INTEGER 00041 * The number of columns of B, the matrix of right hand sides. 00042 * NRHS >= 0. 00043 * 00044 * A (input) REAL array, dimension (N*(N+1)/2) 00045 * The original symmetric matrix A, stored as a packed 00046 * triangular matrix. 00047 * 00048 * X (input) REAL array, dimension (LDX,NRHS) 00049 * The computed solution vectors for the system of linear 00050 * equations. 00051 * 00052 * LDX (input) INTEGER 00053 * The leading dimension of the array X. LDX >= max(1,N). 00054 * 00055 * B (input/output) REAL array, dimension (LDB,NRHS) 00056 * On entry, the right hand side vectors for the system of 00057 * linear equations. 00058 * On exit, B is overwritten with the difference B - A*X. 00059 * 00060 * LDB (input) INTEGER 00061 * The leading dimension of the array B. LDB >= max(1,N). 00062 * 00063 * RWORK (workspace) REAL array, dimension (N) 00064 * 00065 * RESID (output) REAL 00066 * The maximum over the number of right hand sides of 00067 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ). 00068 * 00069 * ===================================================================== 00070 * 00071 * .. Parameters .. 00072 REAL ZERO, ONE 00073 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00074 * .. 00075 * .. Local Scalars .. 00076 INTEGER J 00077 REAL ANORM, BNORM, EPS, XNORM 00078 * .. 00079 * .. External Functions .. 00080 REAL SASUM, SLAMCH, SLANSP 00081 EXTERNAL SASUM, SLAMCH, SLANSP 00082 * .. 00083 * .. External Subroutines .. 00084 EXTERNAL SSPMV 00085 * .. 00086 * .. Intrinsic Functions .. 00087 INTRINSIC MAX 00088 * .. 00089 * .. Executable Statements .. 00090 * 00091 * Quick exit if N = 0 or NRHS = 0. 00092 * 00093 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN 00094 RESID = ZERO 00095 RETURN 00096 END IF 00097 * 00098 * Exit with RESID = 1/EPS if ANORM = 0. 00099 * 00100 EPS = SLAMCH( 'Epsilon' ) 00101 ANORM = SLANSP( '1', UPLO, N, A, RWORK ) 00102 IF( ANORM.LE.ZERO ) THEN 00103 RESID = ONE / EPS 00104 RETURN 00105 END IF 00106 * 00107 * Compute B - A*X for the matrix of right hand sides B. 00108 * 00109 DO 10 J = 1, NRHS 00110 CALL SSPMV( UPLO, N, -ONE, A, X( 1, J ), 1, ONE, B( 1, J ), 1 ) 00111 10 CONTINUE 00112 * 00113 * Compute the maximum over the number of right hand sides of 00114 * norm( B - A*X ) / ( norm(A) * norm(X) * EPS ) . 00115 * 00116 RESID = ZERO 00117 DO 20 J = 1, NRHS 00118 BNORM = SASUM( N, B( 1, J ), 1 ) 00119 XNORM = SASUM( N, X( 1, J ), 1 ) 00120 IF( XNORM.LE.ZERO ) THEN 00121 RESID = ONE / EPS 00122 ELSE 00123 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) 00124 END IF 00125 20 CONTINUE 00126 * 00127 RETURN 00128 * 00129 * End of SPPT02 00130 * 00131 END
1.7.2