![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGBT01 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE DGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK, 00012 * RESID ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER KL, KU, LDA, LDAFAC, M, N 00016 * DOUBLE PRECISION RESID 00017 * .. 00018 * .. Array Arguments .. 00019 * INTEGER IPIV( * ) 00020 * DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), WORK( * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> DGBT01 reconstructs a band matrix A from its L*U factorization and 00030 *> computes the residual: 00031 *> norm(L*U - A) / ( N * norm(A) * EPS ), 00032 *> where EPS is the machine epsilon. 00033 *> 00034 *> The expression L*U - A is computed one column at a time, so A and 00035 *> AFAC are not modified. 00036 *> \endverbatim 00037 * 00038 * Arguments: 00039 * ========== 00040 * 00041 *> \param[in] M 00042 *> \verbatim 00043 *> M is INTEGER 00044 *> The number of rows of the matrix A. M >= 0. 00045 *> \endverbatim 00046 *> 00047 *> \param[in] N 00048 *> \verbatim 00049 *> N is INTEGER 00050 *> The number of columns of the matrix A. N >= 0. 00051 *> \endverbatim 00052 *> 00053 *> \param[in] KL 00054 *> \verbatim 00055 *> KL is INTEGER 00056 *> The number of subdiagonals within the band of A. KL >= 0. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] KU 00060 *> \verbatim 00061 *> KU is INTEGER 00062 *> The number of superdiagonals within the band of A. KU >= 0. 00063 *> \endverbatim 00064 *> 00065 *> \param[in,out] A 00066 *> \verbatim 00067 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00068 *> The original matrix A in band storage, stored in rows 1 to 00069 *> KL+KU+1. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] LDA 00073 *> \verbatim 00074 *> LDA is INTEGER. 00075 *> The leading dimension of the array A. LDA >= max(1,KL+KU+1). 00076 *> \endverbatim 00077 *> 00078 *> \param[in] AFAC 00079 *> \verbatim 00080 *> AFAC is DOUBLE PRECISION array, dimension (LDAFAC,N) 00081 *> The factored form of the matrix A. AFAC contains the banded 00082 *> factors L and U from the L*U factorization, as computed by 00083 *> DGBTRF. U is stored as an upper triangular band matrix with 00084 *> KL+KU superdiagonals in rows 1 to KL+KU+1, and the 00085 *> multipliers used during the factorization are stored in rows 00086 *> KL+KU+2 to 2*KL+KU+1. See DGBTRF for further details. 00087 *> \endverbatim 00088 *> 00089 *> \param[in] LDAFAC 00090 *> \verbatim 00091 *> LDAFAC is INTEGER 00092 *> The leading dimension of the array AFAC. 00093 *> LDAFAC >= max(1,2*KL*KU+1). 00094 *> \endverbatim 00095 *> 00096 *> \param[in] IPIV 00097 *> \verbatim 00098 *> IPIV is INTEGER array, dimension (min(M,N)) 00099 *> The pivot indices from DGBTRF. 00100 *> \endverbatim 00101 *> 00102 *> \param[out] WORK 00103 *> \verbatim 00104 *> WORK is DOUBLE PRECISION array, dimension (2*KL+KU+1) 00105 *> \endverbatim 00106 *> 00107 *> \param[out] RESID 00108 *> \verbatim 00109 *> RESID is DOUBLE PRECISION 00110 *> norm(L*U - A) / ( N * norm(A) * EPS ) 00111 *> \endverbatim 00112 * 00113 * Authors: 00114 * ======== 00115 * 00116 *> \author Univ. of Tennessee 00117 *> \author Univ. of California Berkeley 00118 *> \author Univ. of Colorado Denver 00119 *> \author NAG Ltd. 00120 * 00121 *> \date November 2011 00122 * 00123 *> \ingroup double_lin 00124 * 00125 * ===================================================================== 00126 SUBROUTINE DGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK, 00127 $ RESID ) 00128 * 00129 * -- LAPACK test routine (version 3.4.0) -- 00130 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00131 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00132 * November 2011 00133 * 00134 * .. Scalar Arguments .. 00135 INTEGER KL, KU, LDA, LDAFAC, M, N 00136 DOUBLE PRECISION RESID 00137 * .. 00138 * .. Array Arguments .. 00139 INTEGER IPIV( * ) 00140 DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), WORK( * ) 00141 * .. 00142 * 00143 * ===================================================================== 00144 * 00145 * .. Parameters .. 00146 DOUBLE PRECISION ZERO, ONE 00147 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00148 * .. 00149 * .. Local Scalars .. 00150 INTEGER I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ 00151 DOUBLE PRECISION ANORM, EPS, T 00152 * .. 00153 * .. External Functions .. 00154 DOUBLE PRECISION DASUM, DLAMCH 00155 EXTERNAL DASUM, DLAMCH 00156 * .. 00157 * .. External Subroutines .. 00158 EXTERNAL DAXPY, DCOPY 00159 * .. 00160 * .. Intrinsic Functions .. 00161 INTRINSIC DBLE, MAX, MIN 00162 * .. 00163 * .. Executable Statements .. 00164 * 00165 * Quick exit if M = 0 or N = 0. 00166 * 00167 RESID = ZERO 00168 IF( M.LE.0 .OR. N.LE.0 ) 00169 $ RETURN 00170 * 00171 * Determine EPS and the norm of A. 00172 * 00173 EPS = DLAMCH( 'Epsilon' ) 00174 KD = KU + 1 00175 ANORM = ZERO 00176 DO 10 J = 1, N 00177 I1 = MAX( KD+1-J, 1 ) 00178 I2 = MIN( KD+M-J, KL+KD ) 00179 IF( I2.GE.I1 ) 00180 $ ANORM = MAX( ANORM, DASUM( I2-I1+1, A( I1, J ), 1 ) ) 00181 10 CONTINUE 00182 * 00183 * Compute one column at a time of L*U - A. 00184 * 00185 KD = KL + KU + 1 00186 DO 40 J = 1, N 00187 * 00188 * Copy the J-th column of U to WORK. 00189 * 00190 JU = MIN( KL+KU, J-1 ) 00191 JL = MIN( KL, M-J ) 00192 LENJ = MIN( M, J ) - J + JU + 1 00193 IF( LENJ.GT.0 ) THEN 00194 CALL DCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 ) 00195 DO 20 I = LENJ + 1, JU + JL + 1 00196 WORK( I ) = ZERO 00197 20 CONTINUE 00198 * 00199 * Multiply by the unit lower triangular matrix L. Note that L 00200 * is stored as a product of transformations and permutations. 00201 * 00202 DO 30 I = MIN( M-1, J ), J - JU, -1 00203 IL = MIN( KL, M-I ) 00204 IF( IL.GT.0 ) THEN 00205 IW = I - J + JU + 1 00206 T = WORK( IW ) 00207 CALL DAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ), 00208 $ 1 ) 00209 IP = IPIV( I ) 00210 IF( I.NE.IP ) THEN 00211 IP = IP - J + JU + 1 00212 WORK( IW ) = WORK( IP ) 00213 WORK( IP ) = T 00214 END IF 00215 END IF 00216 30 CONTINUE 00217 * 00218 * Subtract the corresponding column of A. 00219 * 00220 JUA = MIN( JU, KU ) 00221 IF( JUA+JL+1.GT.0 ) 00222 $ CALL DAXPY( JUA+JL+1, -ONE, A( KU+1-JUA, J ), 1, 00223 $ WORK( JU+1-JUA ), 1 ) 00224 * 00225 * Compute the 1-norm of the column. 00226 * 00227 RESID = MAX( RESID, DASUM( JU+JL+1, WORK, 1 ) ) 00228 END IF 00229 40 CONTINUE 00230 * 00231 * Compute norm( L*U - A ) / ( N * norm(A) * EPS ) 00232 * 00233 IF( ANORM.LE.ZERO ) THEN 00234 IF( RESID.NE.ZERO ) 00235 $ RESID = ONE / EPS 00236 ELSE 00237 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 00238 END IF 00239 * 00240 RETURN 00241 * 00242 * End of DGBT01 00243 * 00244 END