LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zgbt01.f
Go to the documentation of this file.
1*> \brief \b ZGBT01
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE ZGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
12* RESID )
13*
14* .. Scalar Arguments ..
15* INTEGER KL, KU, LDA, LDAFAC, M, N
16* DOUBLE PRECISION RESID
17* ..
18* .. Array Arguments ..
19* INTEGER IPIV( * )
20* COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
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.
33*>
34*> The expression L*U - A is computed one column at a time, so A and
35*> AFAC are not modified.
36*> \endverbatim
37*
38* Arguments:
39* ==========
40*
41*> \param[in] M
42*> \verbatim
43*> M is INTEGER
44*> The number of rows of the matrix A. M >= 0.
45*> \endverbatim
46*>
47*> \param[in] N
48*> \verbatim
49*> N is INTEGER
50*> The number of columns of the matrix A. N >= 0.
51*> \endverbatim
52*>
53*> \param[in] KL
54*> \verbatim
55*> KL is INTEGER
56*> The number of subdiagonals within the band of A. KL >= 0.
57*> \endverbatim
58*>
59*> \param[in] KU
60*> \verbatim
61*> KU is INTEGER
62*> The number of superdiagonals within the band of A. KU >= 0.
63*> \endverbatim
64*>
65*> \param[in,out] A
66*> \verbatim
67*> A is COMPLEX*16 array, dimension (LDA,N)
68*> The original matrix A in band storage, stored in rows 1 to
69*> KL+KU+1.
70*> \endverbatim
71*>
72*> \param[in] LDA
73*> \verbatim
74*> LDA is INTEGER.
75*> The leading dimension of the array A. LDA >= max(1,KL+KU+1).
76*> \endverbatim
77*>
78*> \param[in] AFAC
79*> \verbatim
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.
87*> \endverbatim
88*>
89*> \param[in] LDAFAC
90*> \verbatim
91*> LDAFAC is INTEGER
92*> The leading dimension of the array AFAC.
93*> LDAFAC >= max(1,2*KL*KU+1).
94*> \endverbatim
95*>
96*> \param[in] IPIV
97*> \verbatim
98*> IPIV is INTEGER array, dimension (min(M,N))
99*> The pivot indices from ZGBTRF.
100*> \endverbatim
101*>
102*> \param[out] WORK
103*> \verbatim
104*> WORK is COMPLEX*16 array, dimension (2*KL+KU+1)
105*> \endverbatim
106*>
107*> \param[out] RESID
108*> \verbatim
109*> RESID is DOUBLE PRECISION
110*> norm(L*U - A) / ( N * norm(A) * EPS )
111*> \endverbatim
112*
113* Authors:
114* ========
115*
116*> \author Univ. of Tennessee
117*> \author Univ. of California Berkeley
118*> \author Univ. of Colorado Denver
119*> \author NAG Ltd.
120*
121*> \ingroup complex16_lin
122*
123* =====================================================================
124 SUBROUTINE zgbt01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
125 \$ RESID )
126*
127* -- LAPACK test routine --
128* -- LAPACK is a software package provided by Univ. of Tennessee, --
129* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
130*
131* .. Scalar Arguments ..
132 INTEGER KL, KU, LDA, LDAFAC, M, N
133 DOUBLE PRECISION RESID
134* ..
135* .. Array Arguments ..
136 INTEGER IPIV( * )
137 COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
138* ..
139*
140* =====================================================================
141*
142* .. Parameters ..
143 DOUBLE PRECISION ZERO, ONE
144 parameter( zero = 0.0d+0, one = 1.0d+0 )
145* ..
146* .. Local Scalars ..
147 INTEGER I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ
148 DOUBLE PRECISION ANORM, EPS
149 COMPLEX*16 T
150* ..
151* .. External Functions ..
152 DOUBLE PRECISION DLAMCH, DZASUM
153 EXTERNAL dlamch, dzasum
154* ..
155* .. External Subroutines ..
156 EXTERNAL zaxpy, zcopy
157* ..
158* .. Intrinsic Functions ..
159 INTRINSIC dble, dcmplx, max, min
160* ..
161* .. Executable Statements ..
162*
163* Quick exit if M = 0 or N = 0.
164*
165 resid = zero
166 IF( m.LE.0 .OR. n.LE.0 )
167 \$ RETURN
168*
169* Determine EPS and the norm of A.
170*
171 eps = dlamch( 'Epsilon' )
172 kd = ku + 1
173 anorm = zero
174 DO 10 j = 1, n
175 i1 = max( kd+1-j, 1 )
176 i2 = min( kd+m-j, kl+kd )
177 IF( i2.GE.i1 )
178 \$ anorm = max( anorm, dzasum( i2-i1+1, a( i1, j ), 1 ) )
179 10 CONTINUE
180*
181* Compute one column at a time of L*U - A.
182*
183 kd = kl + ku + 1
184 DO 40 j = 1, n
185*
186* Copy the J-th column of U to WORK.
187*
188 ju = min( kl+ku, j-1 )
189 jl = min( kl, m-j )
190 lenj = min( m, j ) - j + ju + 1
191 IF( lenj.GT.0 ) THEN
192 CALL zcopy( lenj, afac( kd-ju, j ), 1, work, 1 )
193 DO 20 i = lenj + 1, ju + jl + 1
194 work( i ) = zero
195 20 CONTINUE
196*
197* Multiply by the unit lower triangular matrix L. Note that L
198* is stored as a product of transformations and permutations.
199*
200 DO 30 i = min( m-1, j ), j - ju, -1
201 il = min( kl, m-i )
202 IF( il.GT.0 ) THEN
203 iw = i - j + ju + 1
204 t = work( iw )
205 CALL zaxpy( il, t, afac( kd+1, i ), 1, work( iw+1 ),
206 \$ 1 )
207 ip = ipiv( i )
208 IF( i.NE.ip ) THEN
209 ip = ip - j + ju + 1
210 work( iw ) = work( ip )
211 work( ip ) = t
212 END IF
213 END IF
214 30 CONTINUE
215*
216* Subtract the corresponding column of A.
217*
218 jua = min( ju, ku )
219 IF( jua+jl+1.GT.0 )
220 \$ CALL zaxpy( jua+jl+1, -dcmplx( one ), a( ku+1-jua, j ),
221 \$ 1, work( ju+1-jua ), 1 )
222*
223* Compute the 1-norm of the column.
224*
225 resid = max( resid, dzasum( ju+jl+1, work, 1 ) )
226 END IF
227 40 CONTINUE
228*
229* Compute norm(L*U - A) / ( N * norm(A) * EPS )
230*
231 IF( anorm.LE.zero ) THEN
232 IF( resid.NE.zero )
233 \$ resid = one / eps
234 ELSE
235 resid = ( ( resid / dble( n ) ) / anorm ) / eps
236 END IF
237*
238 RETURN
239*
240* End of ZGBT01
241*
242 END
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
subroutine zgbt01(m, n, kl, ku, a, lda, afac, ldafac, ipiv, work, resid)
ZGBT01
Definition zgbt01.f:126