![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGBT01 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 ZGBT01( 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 * COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), WORK( * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> ZGBT01 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 COMPLEX*16 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 COMPLEX*16 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 *> ZGBTRF. 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 ZGBTRF 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 ZGBTRF. 00100 *> \endverbatim 00101 *> 00102 *> \param[out] WORK 00103 *> \verbatim 00104 *> WORK is COMPLEX*16 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 complex16_lin 00124 * 00125 * ===================================================================== 00126 SUBROUTINE ZGBT01( 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 COMPLEX*16 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 00152 COMPLEX*16 T 00153 * .. 00154 * .. External Functions .. 00155 DOUBLE PRECISION DLAMCH, DZASUM 00156 EXTERNAL DLAMCH, DZASUM 00157 * .. 00158 * .. External Subroutines .. 00159 EXTERNAL ZAXPY, ZCOPY 00160 * .. 00161 * .. Intrinsic Functions .. 00162 INTRINSIC DBLE, DCMPLX, MAX, MIN 00163 * .. 00164 * .. Executable Statements .. 00165 * 00166 * Quick exit if M = 0 or N = 0. 00167 * 00168 RESID = ZERO 00169 IF( M.LE.0 .OR. N.LE.0 ) 00170 $ RETURN 00171 * 00172 * Determine EPS and the norm of A. 00173 * 00174 EPS = DLAMCH( 'Epsilon' ) 00175 KD = KU + 1 00176 ANORM = ZERO 00177 DO 10 J = 1, N 00178 I1 = MAX( KD+1-J, 1 ) 00179 I2 = MIN( KD+M-J, KL+KD ) 00180 IF( I2.GE.I1 ) 00181 $ ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) ) 00182 10 CONTINUE 00183 * 00184 * Compute one column at a time of L*U - A. 00185 * 00186 KD = KL + KU + 1 00187 DO 40 J = 1, N 00188 * 00189 * Copy the J-th column of U to WORK. 00190 * 00191 JU = MIN( KL+KU, J-1 ) 00192 JL = MIN( KL, M-J ) 00193 LENJ = MIN( M, J ) - J + JU + 1 00194 IF( LENJ.GT.0 ) THEN 00195 CALL ZCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 ) 00196 DO 20 I = LENJ + 1, JU + JL + 1 00197 WORK( I ) = ZERO 00198 20 CONTINUE 00199 * 00200 * Multiply by the unit lower triangular matrix L. Note that L 00201 * is stored as a product of transformations and permutations. 00202 * 00203 DO 30 I = MIN( M-1, J ), J - JU, -1 00204 IL = MIN( KL, M-I ) 00205 IF( IL.GT.0 ) THEN 00206 IW = I - J + JU + 1 00207 T = WORK( IW ) 00208 CALL ZAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ), 00209 $ 1 ) 00210 IP = IPIV( I ) 00211 IF( I.NE.IP ) THEN 00212 IP = IP - J + JU + 1 00213 WORK( IW ) = WORK( IP ) 00214 WORK( IP ) = T 00215 END IF 00216 END IF 00217 30 CONTINUE 00218 * 00219 * Subtract the corresponding column of A. 00220 * 00221 JUA = MIN( JU, KU ) 00222 IF( JUA+JL+1.GT.0 ) 00223 $ CALL ZAXPY( JUA+JL+1, -DCMPLX( ONE ), A( KU+1-JUA, J ), 00224 $ 1, WORK( JU+1-JUA ), 1 ) 00225 * 00226 * Compute the 1-norm of the column. 00227 * 00228 RESID = MAX( RESID, DZASUM( JU+JL+1, WORK, 1 ) ) 00229 END IF 00230 40 CONTINUE 00231 * 00232 * Compute norm( L*U - A ) / ( N * norm(A) * EPS ) 00233 * 00234 IF( ANORM.LE.ZERO ) THEN 00235 IF( RESID.NE.ZERO ) 00236 $ RESID = ONE / EPS 00237 ELSE 00238 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS 00239 END IF 00240 * 00241 RETURN 00242 * 00243 * End of ZGBT01 00244 * 00245 END