00001 SUBROUTINE SPBT01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, 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 KD, LDA, LDAFAC, N 00011 REAL RESID 00012 * .. 00013 * .. Array Arguments .. 00014 REAL A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * SPBT01 reconstructs a symmetric positive definite band matrix A from 00021 * its L*L' or U'*U factorization and computes the residual 00022 * norm( L*L' - A ) / ( N * norm(A) * EPS ) or 00023 * norm( U'*U - A ) / ( N * norm(A) * EPS ), 00024 * where EPS is the machine epsilon, L' is the conjugate transpose of 00025 * L, and U' is the conjugate transpose of U. 00026 * 00027 * Arguments 00028 * ========= 00029 * 00030 * UPLO (input) CHARACTER*1 00031 * Specifies whether the upper or lower triangular part of the 00032 * symmetric matrix A is stored: 00033 * = 'U': Upper triangular 00034 * = 'L': Lower triangular 00035 * 00036 * N (input) INTEGER 00037 * The number of rows and columns of the matrix A. N >= 0. 00038 * 00039 * KD (input) INTEGER 00040 * The number of super-diagonals of the matrix A if UPLO = 'U', 00041 * or the number of sub-diagonals if UPLO = 'L'. KD >= 0. 00042 * 00043 * A (input) REAL array, dimension (LDA,N) 00044 * The original symmetric band matrix A. If UPLO = 'U', the 00045 * upper triangular part of A is stored as a band matrix; if 00046 * UPLO = 'L', the lower triangular part of A is stored. The 00047 * columns of the appropriate triangle are stored in the columns 00048 * of A and the diagonals of the triangle are stored in the rows 00049 * of A. See SPBTRF for further details. 00050 * 00051 * LDA (input) INTEGER. 00052 * The leading dimension of the array A. LDA >= max(1,KD+1). 00053 * 00054 * AFAC (input) REAL array, dimension (LDAFAC,N) 00055 * The factored form of the matrix A. AFAC contains the factor 00056 * L or U from the L*L' or U'*U factorization in band storage 00057 * format, as computed by SPBTRF. 00058 * 00059 * LDAFAC (input) INTEGER 00060 * The leading dimension of the array AFAC. 00061 * LDAFAC >= max(1,KD+1). 00062 * 00063 * RWORK (workspace) REAL array, dimension (N) 00064 * 00065 * RESID (output) REAL 00066 * If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS ) 00067 * If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS ) 00068 * 00069 * ===================================================================== 00070 * 00071 * 00072 * .. Parameters .. 00073 REAL ZERO, ONE 00074 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00075 * .. 00076 * .. Local Scalars .. 00077 INTEGER I, J, K, KC, KLEN, ML, MU 00078 REAL ANORM, EPS, T 00079 * .. 00080 * .. External Functions .. 00081 LOGICAL LSAME 00082 REAL SDOT, SLAMCH, SLANSB 00083 EXTERNAL LSAME, SDOT, SLAMCH, SLANSB 00084 * .. 00085 * .. External Subroutines .. 00086 EXTERNAL SSCAL, SSYR, STRMV 00087 * .. 00088 * .. Intrinsic Functions .. 00089 INTRINSIC MAX, MIN, REAL 00090 * .. 00091 * .. Executable Statements .. 00092 * 00093 * Quick exit if N = 0. 00094 * 00095 IF( N.LE.0 ) THEN 00096 RESID = ZERO 00097 RETURN 00098 END IF 00099 * 00100 * Exit with RESID = 1/EPS if ANORM = 0. 00101 * 00102 EPS = SLAMCH( 'Epsilon' ) 00103 ANORM = SLANSB( '1', UPLO, N, KD, A, LDA, RWORK ) 00104 IF( ANORM.LE.ZERO ) THEN 00105 RESID = ONE / EPS 00106 RETURN 00107 END IF 00108 * 00109 * Compute the product U'*U, overwriting U. 00110 * 00111 IF( LSAME( UPLO, 'U' ) ) THEN 00112 DO 10 K = N, 1, -1 00113 KC = MAX( 1, KD+2-K ) 00114 KLEN = KD + 1 - KC 00115 * 00116 * Compute the (K,K) element of the result. 00117 * 00118 T = SDOT( KLEN+1, AFAC( KC, K ), 1, AFAC( KC, K ), 1 ) 00119 AFAC( KD+1, K ) = T 00120 * 00121 * Compute the rest of column K. 00122 * 00123 IF( KLEN.GT.0 ) 00124 $ CALL STRMV( 'Upper', 'Transpose', 'Non-unit', KLEN, 00125 $ AFAC( KD+1, K-KLEN ), LDAFAC-1, 00126 $ AFAC( KC, K ), 1 ) 00127 * 00128 10 CONTINUE 00129 * 00130 * UPLO = 'L': Compute the product L*L', overwriting L. 00131 * 00132 ELSE 00133 DO 20 K = N, 1, -1 00134 KLEN = MIN( KD, N-K ) 00135 * 00136 * Add a multiple of column K of the factor L to each of 00137 * columns K+1 through N. 00138 * 00139 IF( KLEN.GT.0 ) 00140 $ CALL SSYR( 'Lower', KLEN, ONE, AFAC( 2, K ), 1, 00141 $ AFAC( 1, K+1 ), LDAFAC-1 ) 00142 * 00143 * Scale column K by the diagonal element. 00144 * 00145 T = AFAC( 1, K ) 00146 CALL SSCAL( KLEN+1, T, AFAC( 1, K ), 1 ) 00147 * 00148 20 CONTINUE 00149 END IF 00150 * 00151 * Compute the difference L*L' - A or U'*U - A. 00152 * 00153 IF( LSAME( UPLO, 'U' ) ) THEN 00154 DO 40 J = 1, N 00155 MU = MAX( 1, KD+2-J ) 00156 DO 30 I = MU, KD + 1 00157 AFAC( I, J ) = AFAC( I, J ) - A( I, J ) 00158 30 CONTINUE 00159 40 CONTINUE 00160 ELSE 00161 DO 60 J = 1, N 00162 ML = MIN( KD+1, N-J+1 ) 00163 DO 50 I = 1, ML 00164 AFAC( I, J ) = AFAC( I, J ) - A( I, J ) 00165 50 CONTINUE 00166 60 CONTINUE 00167 END IF 00168 * 00169 * Compute norm( L*L' - A ) / ( N * norm(A) * EPS ) 00170 * 00171 RESID = SLANSB( 'I', UPLO, N, KD, AFAC, LDAFAC, RWORK ) 00172 * 00173 RESID = ( ( RESID / REAL( N ) ) / ANORM ) / EPS 00174 * 00175 RETURN 00176 * 00177 * End of SPBT01 00178 * 00179 END
1.7.2