LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dgbt01()

subroutine dgbt01 ( integer  M,
integer  N,
integer  KL,
integer  KU,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldafac, * )  AFAC,
integer  LDAFAC,
integer, dimension( * )  IPIV,
double precision, dimension( * )  WORK,
double precision  RESID 
)

DGBT01

Purpose:
 DGBT01 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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
          DGBTRF.  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 DGBTRF 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 DGBTRF.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (2*KL+KU+1)
[out]RESID
          RESID is DOUBLE PRECISION
          norm(L*U - A) / ( N * norm(A) * EPS )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 128 of file dgbt01.f.

128 *
129 * -- LAPACK test routine (version 3.7.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 * December 2016
133 *
134 * .. Scalar Arguments ..
135  INTEGER kl, ku, lda, ldafac, m, n
136  DOUBLE PRECISION resid
137 * ..
138 * .. Array Arguments ..
139  INTEGER ipiv( * )
140  DOUBLE PRECISION 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, t
152 * ..
153 * .. External Functions ..
154  DOUBLE PRECISION dasum, dlamch
155  EXTERNAL dasum, dlamch
156 * ..
157 * .. External Subroutines ..
158  EXTERNAL daxpy, dcopy
159 * ..
160 * .. Intrinsic Functions ..
161  INTRINSIC dble, max, min
162 * ..
163 * .. Executable Statements ..
164 *
165 * Quick exit if M = 0 or N = 0.
166 *
167  resid = zero
168  IF( m.LE.0 .OR. n.LE.0 )
169  $ RETURN
170 *
171 * Determine EPS and the norm of A.
172 *
173  eps = dlamch( 'Epsilon' )
174  kd = ku + 1
175  anorm = zero
176  DO 10 j = 1, n
177  i1 = max( kd+1-j, 1 )
178  i2 = min( kd+m-j, kl+kd )
179  IF( i2.GE.i1 )
180  $ anorm = max( anorm, dasum( i2-i1+1, a( i1, j ), 1 ) )
181  10 CONTINUE
182 *
183 * Compute one column at a time of L*U - A.
184 *
185  kd = kl + ku + 1
186  DO 40 j = 1, n
187 *
188 * Copy the J-th column of U to WORK.
189 *
190  ju = min( kl+ku, j-1 )
191  jl = min( kl, m-j )
192  lenj = min( m, j ) - j + ju + 1
193  IF( lenj.GT.0 ) THEN
194  CALL dcopy( lenj, afac( kd-ju, j ), 1, work, 1 )
195  DO 20 i = lenj + 1, ju + jl + 1
196  work( i ) = zero
197  20 CONTINUE
198 *
199 * Multiply by the unit lower triangular matrix L. Note that L
200 * is stored as a product of transformations and permutations.
201 *
202  DO 30 i = min( m-1, j ), j - ju, -1
203  il = min( kl, m-i )
204  IF( il.GT.0 ) THEN
205  iw = i - j + ju + 1
206  t = work( iw )
207  CALL daxpy( il, t, afac( kd+1, i ), 1, work( iw+1 ),
208  $ 1 )
209  ip = ipiv( i )
210  IF( i.NE.ip ) THEN
211  ip = ip - j + ju + 1
212  work( iw ) = work( ip )
213  work( ip ) = t
214  END IF
215  END IF
216  30 CONTINUE
217 *
218 * Subtract the corresponding column of A.
219 *
220  jua = min( ju, ku )
221  IF( jua+jl+1.GT.0 )
222  $ CALL daxpy( jua+jl+1, -one, a( ku+1-jua, j ), 1,
223  $ work( ju+1-jua ), 1 )
224 *
225 * Compute the 1-norm of the column.
226 *
227  resid = max( resid, dasum( ju+jl+1, work, 1 ) )
228  END IF
229  40 CONTINUE
230 *
231 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
232 *
233  IF( anorm.LE.zero ) THEN
234  IF( resid.NE.zero )
235  $ resid = one / eps
236  ELSE
237  resid = ( ( resid / dble( n ) ) / anorm ) / eps
238  END IF
239 *
240  RETURN
241 *
242 * End of DGBT01
243 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:91
double precision function dasum(N, DX, INCX)
DASUM
Definition: dasum.f:73
Here is the call graph for this function:
Here is the caller graph for this function: