LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cgbt01()

subroutine cgbt01 ( integer  M,
integer  N,
integer  KL,
integer  KU,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( ldafac, * )  AFAC,
integer  LDAFAC,
integer, dimension( * )  IPIV,
complex, dimension( * )  WORK,
real  RESID 
)

CGBT01

Purpose:
 CGBT01 reconstructs a band matrix A from its L*U factorization and
 computes the residual:
    norm(L*U - A) / ( N * norm(A) * EPS ),
 where EPS is the machine epsilon.

 The expression L*U - A is computed one column at a time, so A and
 AFAC are not modified.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in]KL
          KL is INTEGER
          The number of subdiagonals within the band of A.  KL >= 0.
[in]KU
          KU is INTEGER
          The number of superdiagonals within the band of A.  KU >= 0.
[in,out]A
          A is COMPLEX array, dimension (LDA,N)
          The original matrix A in band storage, stored in rows 1 to
          KL+KU+1.
[in]LDA
          LDA is INTEGER.
          The leading dimension of the array A.  LDA >= max(1,KL+KU+1).
[in]AFAC
          AFAC is COMPLEX array, dimension (LDAFAC,N)
          The factored form of the matrix A.  AFAC contains the banded
          factors L and U from the L*U factorization, as computed by
          CGBTRF.  U is stored as an upper triangular band matrix with
          KL+KU superdiagonals in rows 1 to KL+KU+1, and the
          multipliers used during the factorization are stored in rows
          KL+KU+2 to 2*KL+KU+1.  See CGBTRF for further details.
[in]LDAFAC
          LDAFAC is INTEGER
          The leading dimension of the array AFAC.
          LDAFAC >= max(1,2*KL*KU+1).
[in]IPIV
          IPIV is INTEGER array, dimension (min(M,N))
          The pivot indices from CGBTRF.
[out]WORK
          WORK is COMPLEX array, dimension (2*KL+KU+1)
[out]RESID
          RESID is REAL
          norm(L*U - A) / ( N * norm(A) * EPS )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 124 of file cgbt01.f.

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  REAL RESID
134 * ..
135 * .. Array Arguments ..
136  INTEGER IPIV( * )
137  COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
138 * ..
139 *
140 * =====================================================================
141 *
142 * .. Parameters ..
143  REAL ZERO, ONE
144  parameter( zero = 0.0e+0, one = 1.0e+0 )
145 * ..
146 * .. Local Scalars ..
147  INTEGER I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ
148  REAL ANORM, EPS
149  COMPLEX T
150 * ..
151 * .. External Functions ..
152  REAL SCASUM, SLAMCH
153  EXTERNAL scasum, slamch
154 * ..
155 * .. External Subroutines ..
156  EXTERNAL caxpy, ccopy
157 * ..
158 * .. Intrinsic Functions ..
159  INTRINSIC cmplx, max, min, real
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 = slamch( '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, scasum( 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 ccopy( 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 caxpy( 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 caxpy( jua+jl+1, -cmplx( one ), a( ku+1-jua, j ), 1,
221  $ work( ju+1-jua ), 1 )
222 *
223 * Compute the 1-norm of the column.
224 *
225  resid = max( resid, scasum( 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 / real( n ) ) / anorm ) / eps
236  END IF
237 *
238  RETURN
239 *
240 * End of CGBT01
241 *
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
real function scasum(N, CX, INCX)
SCASUM
Definition: scasum.f:72
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: