LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 *> \date November 2011
122 *
123 *> \ingroup complex16_lin
124 *
125 * =====================================================================
126  SUBROUTINE zgbt01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
127  $ resid )
128 *
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..--
132 * November 2011
133 *
134 * .. Scalar Arguments ..
135  INTEGER KL, KU, LDA, LDAFAC, M, N
136  DOUBLE PRECISION RESID
137 * ..
138 * .. Array Arguments ..
139  INTEGER IPIV( * )
140  COMPLEX*16 A( lda, * ), AFAC( ldafac, * ), WORK( * )
141 * ..
142 *
143 * =====================================================================
144 *
145 * .. Parameters ..
146  DOUBLE PRECISION ZERO, ONE
147  parameter ( zero = 0.0d+0, one = 1.0d+0 )
148 * ..
149 * .. Local Scalars ..
150  INTEGER I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ
151  DOUBLE PRECISION ANORM, EPS
152  COMPLEX*16 T
153 * ..
154 * .. External Functions ..
155  DOUBLE PRECISION DLAMCH, DZASUM
156  EXTERNAL dlamch, dzasum
157 * ..
158 * .. External Subroutines ..
159  EXTERNAL zaxpy, zcopy
160 * ..
161 * .. Intrinsic Functions ..
162  INTRINSIC dble, dcmplx, max, min
163 * ..
164 * .. Executable Statements ..
165 *
166 * Quick exit if M = 0 or N = 0.
167 *
168  resid = zero
169  IF( m.LE.0 .OR. n.LE.0 )
170  $ RETURN
171 *
172 * Determine EPS and the norm of A.
173 *
174  eps = dlamch( 'Epsilon' )
175  kd = ku + 1
176  anorm = zero
177  DO 10 j = 1, n
178  i1 = max( kd+1-j, 1 )
179  i2 = min( kd+m-j, kl+kd )
180  IF( i2.GE.i1 )
181  $ anorm = max( anorm, dzasum( i2-i1+1, a( i1, j ), 1 ) )
182  10 CONTINUE
183 *
184 * Compute one column at a time of L*U - A.
185 *
186  kd = kl + ku + 1
187  DO 40 j = 1, n
188 *
189 * Copy the J-th column of U to WORK.
190 *
191  ju = min( kl+ku, j-1 )
192  jl = min( kl, m-j )
193  lenj = min( m, j ) - j + ju + 1
194  IF( lenj.GT.0 ) THEN
195  CALL zcopy( lenj, afac( kd-ju, j ), 1, work, 1 )
196  DO 20 i = lenj + 1, ju + jl + 1
197  work( i ) = zero
198  20 CONTINUE
199 *
200 * Multiply by the unit lower triangular matrix L. Note that L
201 * is stored as a product of transformations and permutations.
202 *
203  DO 30 i = min( m-1, j ), j - ju, -1
204  il = min( kl, m-i )
205  IF( il.GT.0 ) THEN
206  iw = i - j + ju + 1
207  t = work( iw )
208  CALL zaxpy( il, t, afac( kd+1, i ), 1, work( iw+1 ),
209  $ 1 )
210  ip = ipiv( i )
211  IF( i.NE.ip ) THEN
212  ip = ip - j + ju + 1
213  work( iw ) = work( ip )
214  work( ip ) = t
215  END IF
216  END IF
217  30 CONTINUE
218 *
219 * Subtract the corresponding column of A.
220 *
221  jua = min( ju, ku )
222  IF( jua+jl+1.GT.0 )
223  $ CALL zaxpy( jua+jl+1, -dcmplx( one ), a( ku+1-jua, j ),
224  $ 1, work( ju+1-jua ), 1 )
225 *
226 * Compute the 1-norm of the column.
227 *
228  resid = max( resid, dzasum( ju+jl+1, work, 1 ) )
229  END IF
230  40 CONTINUE
231 *
232 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
233 *
234  IF( anorm.LE.zero ) THEN
235  IF( resid.NE.zero )
236  $ resid = one / eps
237  ELSE
238  resid = ( ( resid / dble( n ) ) / anorm ) / eps
239  END IF
240 *
241  RETURN
242 *
243 * End of ZGBT01
244 *
245  END
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zgbt01(M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK, RESID)
ZGBT01
Definition: zgbt01.f:128
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53