00001 SUBROUTINE SSPT01( UPLO, N, A, AFAC, IPIV, C, LDC, RWORK, RESID ) 00002 * 00003 * -- LAPACK test routine (version 3.1) -- 00004 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00005 * November 2006 00006 * 00007 * .. Scalar Arguments .. 00008 CHARACTER UPLO 00009 INTEGER LDC, N 00010 REAL RESID 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IPIV( * ) 00014 REAL A( * ), AFAC( * ), C( LDC, * ), RWORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * SSPT01 reconstructs a symmetric indefinite packed matrix A from its 00021 * block L*D*L' or U*D*U' factorization and computes the residual 00022 * norm( C - A ) / ( N * norm(A) * EPS ), 00023 * where C is the reconstructed matrix and EPS is the machine epsilon. 00024 * 00025 * Arguments 00026 * ========== 00027 * 00028 * UPLO (input) CHARACTER*1 00029 * Specifies whether the upper or lower triangular part of the 00030 * symmetric matrix A is stored: 00031 * = 'U': Upper triangular 00032 * = 'L': Lower triangular 00033 * 00034 * N (input) INTEGER 00035 * The number of rows and columns of the matrix A. N >= 0. 00036 * 00037 * A (input) REAL array, dimension (N*(N+1)/2) 00038 * The original symmetric matrix A, stored as a packed 00039 * triangular matrix. 00040 * 00041 * AFAC (input) REAL array, dimension (N*(N+1)/2) 00042 * The factored form of the matrix A, stored as a packed 00043 * triangular matrix. AFAC contains the block diagonal matrix D 00044 * and the multipliers used to obtain the factor L or U from the 00045 * block L*D*L' or U*D*U' factorization as computed by SSPTRF. 00046 * 00047 * IPIV (input) INTEGER array, dimension (N) 00048 * The pivot indices from SSPTRF. 00049 * 00050 * C (workspace) REAL array, dimension (LDC,N) 00051 * 00052 * LDC (integer) INTEGER 00053 * The leading dimension of the array C. LDC >= max(1,N). 00054 * 00055 * RWORK (workspace) REAL array, dimension (N) 00056 * 00057 * RESID (output) REAL 00058 * If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) 00059 * If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) 00060 * 00061 * ===================================================================== 00062 * 00063 * .. Parameters .. 00064 REAL ZERO, ONE 00065 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00066 * .. 00067 * .. Local Scalars .. 00068 INTEGER I, INFO, J, JC 00069 REAL ANORM, EPS 00070 * .. 00071 * .. External Functions .. 00072 LOGICAL LSAME 00073 REAL SLAMCH, SLANSP, SLANSY 00074 EXTERNAL LSAME, SLAMCH, SLANSP, SLANSY 00075 * .. 00076 * .. External Subroutines .. 00077 EXTERNAL SLAVSP, SLASET 00078 * .. 00079 * .. Intrinsic Functions .. 00080 INTRINSIC REAL 00081 * .. 00082 * .. Executable Statements .. 00083 * 00084 * Quick exit if N = 0. 00085 * 00086 IF( N.LE.0 ) THEN 00087 RESID = ZERO 00088 RETURN 00089 END IF 00090 * 00091 * Determine EPS and the norm of A. 00092 * 00093 EPS = SLAMCH( 'Epsilon' ) 00094 ANORM = SLANSP( '1', UPLO, N, A, RWORK ) 00095 * 00096 * Initialize C to the identity matrix. 00097 * 00098 CALL SLASET( 'Full', N, N, ZERO, ONE, C, LDC ) 00099 * 00100 * Call SLAVSP to form the product D * U' (or D * L' ). 00101 * 00102 CALL SLAVSP( UPLO, 'Transpose', 'Non-unit', N, N, AFAC, IPIV, C, 00103 $ LDC, INFO ) 00104 * 00105 * Call SLAVSP again to multiply by U ( or L ). 00106 * 00107 CALL SLAVSP( UPLO, 'No transpose', 'Unit', N, N, AFAC, IPIV, C, 00108 $ LDC, INFO ) 00109 * 00110 * Compute the difference C - A . 00111 * 00112 IF( LSAME( UPLO, 'U' ) ) THEN 00113 JC = 0 00114 DO 20 J = 1, N 00115 DO 10 I = 1, J 00116 C( I, J ) = C( I, J ) - A( JC+I ) 00117 10 CONTINUE 00118 JC = JC + J 00119 20 CONTINUE 00120 ELSE 00121 JC = 1 00122 DO 40 J = 1, N 00123 DO 30 I = J, N 00124 C( I, J ) = C( I, J ) - A( JC+I-J ) 00125 30 CONTINUE 00126 JC = JC + N - J + 1 00127 40 CONTINUE 00128 END IF 00129 * 00130 * Compute norm( C - A ) / ( N * norm(A) * EPS ) 00131 * 00132 RESID = SLANSY( '1', UPLO, N, C, LDC, RWORK ) 00133 * 00134 IF( ANORM.LE.ZERO ) THEN 00135 IF( RESID.NE.ZERO ) 00136 $ RESID = ONE / EPS 00137 ELSE 00138 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 00139 END IF 00140 * 00141 RETURN 00142 * 00143 * End of SSPT01 00144 * 00145 END
1.7.2