3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
11 * SUBROUTINE ZGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
14 * .. Scalar Arguments ..
15 * INTEGER KL, KU, LDA, LDAFAC, M, N
16 * DOUBLE PRECISION RESID
18 * .. Array Arguments ..
20 * COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
29 *> ZGBT01 reconstructs a band matrix A from its L*U factorization and
30 *> computes the residual:
31 *> norm(L*U - A) / ( N * norm(A) * EPS ),
32 *> where EPS is the machine epsilon.
34 *> The expression L*U - A is computed one column at a time, so A and
35 *> AFAC are not modified.
44 *> The number of rows of the matrix A. M >= 0.
50 *> The number of columns of the matrix A. N >= 0.
56 *> The number of subdiagonals within the band of A. KL >= 0.
62 *> The number of superdiagonals within the band of A. KU >= 0.
67 *> A is COMPLEX*16 array, dimension (LDA,N)
68 *> The original matrix A in band storage, stored in rows 1 to
75 *> The leading dimension of the array A. LDA >= max(1,KL+KU+1).
80 *> AFAC is COMPLEX*16 array, dimension (LDAFAC,N)
81 *> The factored form of the matrix A. AFAC contains the banded
82 *> factors L and U from the L*U factorization, as computed by
83 *> ZGBTRF. U is stored as an upper triangular band matrix with
84 *> KL+KU superdiagonals in rows 1 to KL+KU+1, and the
85 *> multipliers used during the factorization are stored in rows
86 *> KL+KU+2 to 2*KL+KU+1. See ZGBTRF for further details.
92 *> The leading dimension of the array AFAC.
93 *> LDAFAC >= max(1,2*KL*KU+1).
98 *> IPIV is INTEGER array, dimension (min(M,N))
99 *> The pivot indices from ZGBTRF.
104 *> WORK is COMPLEX*16 array, dimension (2*KL+KU+1)
109 *> RESID is DOUBLE PRECISION
110 *> norm(L*U - A) / ( N * norm(A) * EPS )
116 *> \author Univ. of Tennessee
117 *> \author Univ. of California Berkeley
118 *> \author Univ. of Colorado Denver
121 *> \date November 2011
123 *> \ingroup complex16_lin
125 * =====================================================================
126 SUBROUTINE ZGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
129 * -- LAPACK test routine (version 3.4.0) --
130 * -- LAPACK is a software package provided by Univ. of Tennessee, --
131 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
134 * .. Scalar Arguments ..
135 INTEGER KL, KU, LDA, LDAFAC, M, N
136 DOUBLE PRECISION RESID
138 * .. Array Arguments ..
140 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
143 * =====================================================================
146 DOUBLE PRECISION ZERO, ONE
147 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
149 * .. Local Scalars ..
150 INTEGER I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ
151 DOUBLE PRECISION ANORM, EPS
154 * .. External Functions ..
155 DOUBLE PRECISION DLAMCH, DZASUM
156 EXTERNAL DLAMCH, DZASUM
158 * .. External Subroutines ..
159 EXTERNAL ZAXPY, ZCOPY
161 * .. Intrinsic Functions ..
162 INTRINSIC DBLE, DCMPLX, MAX, MIN
164 * .. Executable Statements ..
166 * Quick exit if M = 0 or N = 0.
169 IF( M.LE.0 .OR. N.LE.0 )
172 * Determine EPS and the norm of A.
174 EPS = DLAMCH( 'Epsilon' )
178 I1 = MAX( KD+1-J, 1 )
179 I2 = MIN( KD+M-J, KL+KD )
181 $ ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) )
184 * Compute one column at a time of L*U - A.
189 * Copy the J-th column of U to WORK.
191 JU = MIN( KL+KU, J-1 )
193 LENJ = MIN( M, J ) - J + JU + 1
195 CALL ZCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 )
196 DO 20 I = LENJ + 1, JU + JL + 1
200 * Multiply by the unit lower triangular matrix L. Note that L
201 * is stored as a product of transformations and permutations.
203 DO 30 I = MIN( M-1, J ), J - JU, -1
208 CALL ZAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ),
213 WORK( IW ) = WORK( IP )
219 * Subtract the corresponding column of A.
223 $ CALL ZAXPY( JUA+JL+1, -DCMPLX( ONE ), A( KU+1-JUA, J ),
224 $ 1, WORK( JU+1-JUA ), 1 )
226 * Compute the 1-norm of the column.
228 RESID = MAX( RESID, DZASUM( JU+JL+1, WORK, 1 ) )
232 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
234 IF( ANORM.LE.ZERO ) THEN
238 RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS