OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
double
LAPACK » General Band Matrix » Computational routines

Functions

subroutine dgbbrd (vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, c, ldc, work, info)
 DGBBRD
subroutine dgbcon (norm, n, kl, ku, ab, ldab, ipiv, anorm, rcond, work, iwork, info)
 DGBCON
subroutine dgbequ (m, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd, amax, info)
 DGBEQU
subroutine dgbequb (m, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd, amax, info)
 DGBEQUB
subroutine dgbrfs (trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, b, ldb, x, ldx, ferr, berr, work, iwork, info)
 DGBRFS
subroutine dgbrfsx (trans, equed, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, r, c, b, ldb, x, ldx, rcond, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info)
 DGBRFSX
subroutine dgbtf2 (m, n, kl, ku, ab, ldab, ipiv, info)
 DGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algorithm.
subroutine dgbtrf (m, n, kl, ku, ab, ldab, ipiv, info)
 DGBTRF
subroutine dgbtrs (trans, n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
 DGBTRS
subroutine dggbak (job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
 DGGBAK
subroutine dggbal (job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
 DGGBAL
subroutine dla_gbamv (trans, m, n, kl, ku, alpha, ab, ldab, x, incx, beta, y, incy)
 DLA_GBAMV performs a matrix-vector operation to calculate error bounds.
double precision function dla_gbrcond (trans, n, kl, ku, ab, ldab, afb, ldafb, ipiv, cmode, c, info, work, iwork)
 DLA_GBRCOND estimates the Skeel condition number for a general banded matrix.
subroutine dla_gbrfsx_extended (prec_type, trans_type, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, colequ, c, b, ldb, y, ldy, berr_out, n_norms, err_bnds_norm, err_bnds_comp, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)
 DLA_GBRFSX_EXTENDED improves the computed solution to a system of linear equations for general banded matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.
double precision function dla_gbrpvgrw (n, kl, ku, ncols, ab, ldab, afb, ldafb)
 DLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix.
subroutine dorgbr (vect, m, n, k, a, lda, tau, work, lwork, info)
 DORGBR

Detailed Description

This is the group of double computational functions for GB matrices

Function Documentation

◆ dgbbrd()

subroutine dgbbrd ( character vect,
integer m,
integer n,
integer ncc,
integer kl,
integer ku,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
double precision, dimension( ldq, * ) q,
integer ldq,
double precision, dimension( ldpt, * ) pt,
integer ldpt,
double precision, dimension( ldc, * ) c,
integer ldc,
double precision, dimension( * ) work,
integer info )

DGBBRD

Download DGBBRD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGBBRD reduces a real general m-by-n band matrix A to upper
!> bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
!>
!> The routine computes B, and optionally forms Q or P**T, or computes
!> Q**T*C for a given matrix C.
!>
Parameters
[in]VECT
!>          VECT is CHARACTER*1
!>          Specifies whether or not the matrices Q and P**T are to be
!>          formed.
!>          = 'N': do not form Q or P**T;
!>          = 'Q': form Q only;
!>          = 'P': form P**T only;
!>          = 'B': form both.
!>
[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]NCC
!>          NCC is INTEGER
!>          The number of columns of the matrix C.  NCC >= 0.
!>
[in]KL
!>          KL is INTEGER
!>          The number of subdiagonals of the matrix A. KL >= 0.
!>
[in]KU
!>          KU is INTEGER
!>          The number of superdiagonals of the matrix A. KU >= 0.
!>
[in,out]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>          On entry, the m-by-n band matrix A, stored in rows 1 to
!>          KL+KU+1. The j-th column of A is stored in the j-th column of
!>          the array AB as follows:
!>          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
!>          On exit, A is overwritten by values generated during the
!>          reduction.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array A. LDAB >= KL+KU+1.
!>
[out]D
!>          D is DOUBLE PRECISION array, dimension (min(M,N))
!>          The diagonal elements of the bidiagonal matrix B.
!>
[out]E
!>          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
!>          The superdiagonal elements of the bidiagonal matrix B.
!>
[out]Q
!>          Q is DOUBLE PRECISION array, dimension (LDQ,M)
!>          If VECT = 'Q' or 'B', the m-by-m orthogonal matrix Q.
!>          If VECT = 'N' or 'P', the array Q is not referenced.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q.
!>          LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise.
!>
[out]PT
!>          PT is DOUBLE PRECISION array, dimension (LDPT,N)
!>          If VECT = 'P' or 'B', the n-by-n orthogonal matrix P'.
!>          If VECT = 'N' or 'Q', the array PT is not referenced.
!>
[in]LDPT
!>          LDPT is INTEGER
!>          The leading dimension of the array PT.
!>          LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise.
!>
[in,out]C
!>          C is DOUBLE PRECISION array, dimension (LDC,NCC)
!>          On entry, an m-by-ncc matrix C.
!>          On exit, C is overwritten by Q**T*C.
!>          C is not referenced if NCC = 0.
!>
[in]LDC
!>          LDC is INTEGER
!>          The leading dimension of the array C.
!>          LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (2*max(M,N))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 185 of file dgbbrd.f.

187*
188* -- LAPACK computational routine --
189* -- LAPACK is a software package provided by Univ. of Tennessee, --
190* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
191*
192* .. Scalar Arguments ..
193 CHARACTER VECT
194 INTEGER INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
195* ..
196* .. Array Arguments ..
197 DOUBLE PRECISION AB( LDAB, * ), C( LDC, * ), D( * ), E( * ),
198 $ PT( LDPT, * ), Q( LDQ, * ), WORK( * )
199* ..
200*
201* =====================================================================
202*
203* .. Parameters ..
204 DOUBLE PRECISION ZERO, ONE
205 parameter( zero = 0.0d+0, one = 1.0d+0 )
206* ..
207* .. Local Scalars ..
208 LOGICAL WANTB, WANTC, WANTPT, WANTQ
209 INTEGER I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1,
210 $ KUN, L, MINMN, ML, ML0, MN, MU, MU0, NR, NRT
211 DOUBLE PRECISION RA, RB, RC, RS
212* ..
213* .. External Subroutines ..
214 EXTERNAL dlargv, dlartg, dlartv, dlaset, drot, xerbla
215* ..
216* .. Intrinsic Functions ..
217 INTRINSIC max, min
218* ..
219* .. External Functions ..
220 LOGICAL LSAME
221 EXTERNAL lsame
222* ..
223* .. Executable Statements ..
224*
225* Test the input parameters
226*
227 wantb = lsame( vect, 'B' )
228 wantq = lsame( vect, 'Q' ) .OR. wantb
229 wantpt = lsame( vect, 'P' ) .OR. wantb
230 wantc = ncc.GT.0
231 klu1 = kl + ku + 1
232 info = 0
233 IF( .NOT.wantq .AND. .NOT.wantpt .AND. .NOT.lsame( vect, 'N' ) )
234 $ THEN
235 info = -1
236 ELSE IF( m.LT.0 ) THEN
237 info = -2
238 ELSE IF( n.LT.0 ) THEN
239 info = -3
240 ELSE IF( ncc.LT.0 ) THEN
241 info = -4
242 ELSE IF( kl.LT.0 ) THEN
243 info = -5
244 ELSE IF( ku.LT.0 ) THEN
245 info = -6
246 ELSE IF( ldab.LT.klu1 ) THEN
247 info = -8
248 ELSE IF( ldq.LT.1 .OR. wantq .AND. ldq.LT.max( 1, m ) ) THEN
249 info = -12
250 ELSE IF( ldpt.LT.1 .OR. wantpt .AND. ldpt.LT.max( 1, n ) ) THEN
251 info = -14
252 ELSE IF( ldc.LT.1 .OR. wantc .AND. ldc.LT.max( 1, m ) ) THEN
253 info = -16
254 END IF
255 IF( info.NE.0 ) THEN
256 CALL xerbla( 'DGBBRD', -info )
257 RETURN
258 END IF
259*
260* Initialize Q and P**T to the unit matrix, if needed
261*
262 IF( wantq )
263 $ CALL dlaset( 'Full', m, m, zero, one, q, ldq )
264 IF( wantpt )
265 $ CALL dlaset( 'Full', n, n, zero, one, pt, ldpt )
266*
267* Quick return if possible.
268*
269 IF( m.EQ.0 .OR. n.EQ.0 )
270 $ RETURN
271*
272 minmn = min( m, n )
273*
274 IF( kl+ku.GT.1 ) THEN
275*
276* Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
277* first to lower bidiagonal form and then transform to upper
278* bidiagonal
279*
280 IF( ku.GT.0 ) THEN
281 ml0 = 1
282 mu0 = 2
283 ELSE
284 ml0 = 2
285 mu0 = 1
286 END IF
287*
288* Wherever possible, plane rotations are generated and applied in
289* vector operations of length NR over the index set J1:J2:KLU1.
290*
291* The sines of the plane rotations are stored in WORK(1:max(m,n))
292* and the cosines in WORK(max(m,n)+1:2*max(m,n)).
293*
294 mn = max( m, n )
295 klm = min( m-1, kl )
296 kun = min( n-1, ku )
297 kb = klm + kun
298 kb1 = kb + 1
299 inca = kb1*ldab
300 nr = 0
301 j1 = klm + 2
302 j2 = 1 - kun
303*
304 DO 90 i = 1, minmn
305*
306* Reduce i-th column and i-th row of matrix to bidiagonal form
307*
308 ml = klm + 1
309 mu = kun + 1
310 DO 80 kk = 1, kb
311 j1 = j1 + kb
312 j2 = j2 + kb
313*
314* generate plane rotations to annihilate nonzero elements
315* which have been created below the band
316*
317 IF( nr.GT.0 )
318 $ CALL dlargv( nr, ab( klu1, j1-klm-1 ), inca,
319 $ work( j1 ), kb1, work( mn+j1 ), kb1 )
320*
321* apply plane rotations from the left
322*
323 DO 10 l = 1, kb
324 IF( j2-klm+l-1.GT.n ) THEN
325 nrt = nr - 1
326 ELSE
327 nrt = nr
328 END IF
329 IF( nrt.GT.0 )
330 $ CALL dlartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,
331 $ ab( klu1-l+1, j1-klm+l-1 ), inca,
332 $ work( mn+j1 ), work( j1 ), kb1 )
333 10 CONTINUE
334*
335 IF( ml.GT.ml0 ) THEN
336 IF( ml.LE.m-i+1 ) THEN
337*
338* generate plane rotation to annihilate a(i+ml-1,i)
339* within the band, and apply rotation from the left
340*
341 CALL dlartg( ab( ku+ml-1, i ), ab( ku+ml, i ),
342 $ work( mn+i+ml-1 ), work( i+ml-1 ),
343 $ ra )
344 ab( ku+ml-1, i ) = ra
345 IF( i.LT.n )
346 $ CALL drot( min( ku+ml-2, n-i ),
347 $ ab( ku+ml-2, i+1 ), ldab-1,
348 $ ab( ku+ml-1, i+1 ), ldab-1,
349 $ work( mn+i+ml-1 ), work( i+ml-1 ) )
350 END IF
351 nr = nr + 1
352 j1 = j1 - kb1
353 END IF
354*
355 IF( wantq ) THEN
356*
357* accumulate product of plane rotations in Q
358*
359 DO 20 j = j1, j2, kb1
360 CALL drot( m, q( 1, j-1 ), 1, q( 1, j ), 1,
361 $ work( mn+j ), work( j ) )
362 20 CONTINUE
363 END IF
364*
365 IF( wantc ) THEN
366*
367* apply plane rotations to C
368*
369 DO 30 j = j1, j2, kb1
370 CALL drot( ncc, c( j-1, 1 ), ldc, c( j, 1 ), ldc,
371 $ work( mn+j ), work( j ) )
372 30 CONTINUE
373 END IF
374*
375 IF( j2+kun.GT.n ) THEN
376*
377* adjust J2 to keep within the bounds of the matrix
378*
379 nr = nr - 1
380 j2 = j2 - kb1
381 END IF
382*
383 DO 40 j = j1, j2, kb1
384*
385* create nonzero element a(j-1,j+ku) above the band
386* and store it in WORK(n+1:2*n)
387*
388 work( j+kun ) = work( j )*ab( 1, j+kun )
389 ab( 1, j+kun ) = work( mn+j )*ab( 1, j+kun )
390 40 CONTINUE
391*
392* generate plane rotations to annihilate nonzero elements
393* which have been generated above the band
394*
395 IF( nr.GT.0 )
396 $ CALL dlargv( nr, ab( 1, j1+kun-1 ), inca,
397 $ work( j1+kun ), kb1, work( mn+j1+kun ),
398 $ kb1 )
399*
400* apply plane rotations from the right
401*
402 DO 50 l = 1, kb
403 IF( j2+l-1.GT.m ) THEN
404 nrt = nr - 1
405 ELSE
406 nrt = nr
407 END IF
408 IF( nrt.GT.0 )
409 $ CALL dlartv( nrt, ab( l+1, j1+kun-1 ), inca,
410 $ ab( l, j1+kun ), inca,
411 $ work( mn+j1+kun ), work( j1+kun ),
412 $ kb1 )
413 50 CONTINUE
414*
415 IF( ml.EQ.ml0 .AND. mu.GT.mu0 ) THEN
416 IF( mu.LE.n-i+1 ) THEN
417*
418* generate plane rotation to annihilate a(i,i+mu-1)
419* within the band, and apply rotation from the right
420*
421 CALL dlartg( ab( ku-mu+3, i+mu-2 ),
422 $ ab( ku-mu+2, i+mu-1 ),
423 $ work( mn+i+mu-1 ), work( i+mu-1 ),
424 $ ra )
425 ab( ku-mu+3, i+mu-2 ) = ra
426 CALL drot( min( kl+mu-2, m-i ),
427 $ ab( ku-mu+4, i+mu-2 ), 1,
428 $ ab( ku-mu+3, i+mu-1 ), 1,
429 $ work( mn+i+mu-1 ), work( i+mu-1 ) )
430 END IF
431 nr = nr + 1
432 j1 = j1 - kb1
433 END IF
434*
435 IF( wantpt ) THEN
436*
437* accumulate product of plane rotations in P**T
438*
439 DO 60 j = j1, j2, kb1
440 CALL drot( n, pt( j+kun-1, 1 ), ldpt,
441 $ pt( j+kun, 1 ), ldpt, work( mn+j+kun ),
442 $ work( j+kun ) )
443 60 CONTINUE
444 END IF
445*
446 IF( j2+kb.GT.m ) THEN
447*
448* adjust J2 to keep within the bounds of the matrix
449*
450 nr = nr - 1
451 j2 = j2 - kb1
452 END IF
453*
454 DO 70 j = j1, j2, kb1
455*
456* create nonzero element a(j+kl+ku,j+ku-1) below the
457* band and store it in WORK(1:n)
458*
459 work( j+kb ) = work( j+kun )*ab( klu1, j+kun )
460 ab( klu1, j+kun ) = work( mn+j+kun )*ab( klu1, j+kun )
461 70 CONTINUE
462*
463 IF( ml.GT.ml0 ) THEN
464 ml = ml - 1
465 ELSE
466 mu = mu - 1
467 END IF
468 80 CONTINUE
469 90 CONTINUE
470 END IF
471*
472 IF( ku.EQ.0 .AND. kl.GT.0 ) THEN
473*
474* A has been reduced to lower bidiagonal form
475*
476* Transform lower bidiagonal form to upper bidiagonal by applying
477* plane rotations from the left, storing diagonal elements in D
478* and off-diagonal elements in E
479*
480 DO 100 i = 1, min( m-1, n )
481 CALL dlartg( ab( 1, i ), ab( 2, i ), rc, rs, ra )
482 d( i ) = ra
483 IF( i.LT.n ) THEN
484 e( i ) = rs*ab( 1, i+1 )
485 ab( 1, i+1 ) = rc*ab( 1, i+1 )
486 END IF
487 IF( wantq )
488 $ CALL drot( m, q( 1, i ), 1, q( 1, i+1 ), 1, rc, rs )
489 IF( wantc )
490 $ CALL drot( ncc, c( i, 1 ), ldc, c( i+1, 1 ), ldc, rc,
491 $ rs )
492 100 CONTINUE
493 IF( m.LE.n )
494 $ d( m ) = ab( 1, m )
495 ELSE IF( ku.GT.0 ) THEN
496*
497* A has been reduced to upper bidiagonal form
498*
499 IF( m.LT.n ) THEN
500*
501* Annihilate a(m,m+1) by applying plane rotations from the
502* right, storing diagonal elements in D and off-diagonal
503* elements in E
504*
505 rb = ab( ku, m+1 )
506 DO 110 i = m, 1, -1
507 CALL dlartg( ab( ku+1, i ), rb, rc, rs, ra )
508 d( i ) = ra
509 IF( i.GT.1 ) THEN
510 rb = -rs*ab( ku, i )
511 e( i-1 ) = rc*ab( ku, i )
512 END IF
513 IF( wantpt )
514 $ CALL drot( n, pt( i, 1 ), ldpt, pt( m+1, 1 ), ldpt,
515 $ rc, rs )
516 110 CONTINUE
517 ELSE
518*
519* Copy off-diagonal elements to E and diagonal elements to D
520*
521 DO 120 i = 1, minmn - 1
522 e( i ) = ab( ku, i+1 )
523 120 CONTINUE
524 DO 130 i = 1, minmn
525 d( i ) = ab( ku+1, i )
526 130 CONTINUE
527 END IF
528 ELSE
529*
530* A is diagonal. Set elements of E to zero and copy diagonal
531* elements to D.
532*
533 DO 140 i = 1, minmn - 1
534 e( i ) = zero
535 140 CONTINUE
536 DO 150 i = 1, minmn
537 d( i ) = ab( 1, i )
538 150 CONTINUE
539 END IF
540 RETURN
541*
542* End of DGBBRD
543*
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:113
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
dlartv
subroutine dlartv(n, x, incx, y, incy, c, s, incc)
DLARTV applies a vector of plane rotations with real cosines and real sines to the elements of a pair...
Definition dlartv.f:108
dlargv
subroutine dlargv(n, x, incx, y, incy, c, incc)
DLARGV generates a vector of plane rotations with real cosines and real sines.
Definition dlargv.f:104
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ dgbcon()

subroutine dgbcon ( character norm,
integer n,
integer kl,
integer ku,
double precision, dimension( ldab, * ) ab,
integer ldab,
integer, dimension( * ) ipiv,
double precision anorm,
double precision rcond,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

DGBCON

Download DGBCON + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGBCON estimates the reciprocal of the condition number of a real
!> general band matrix A, in either the 1-norm or the infinity-norm,
!> using the LU factorization computed by DGBTRF.
!>
!> An estimate is obtained for norm(inv(A)), and the reciprocal of the
!> condition number is computed as
!>    RCOND = 1 / ( norm(A) * norm(inv(A)) ).
!>
Parameters
[in]NORM
!>          NORM is CHARACTER*1
!>          Specifies whether the 1-norm condition number or the
!>          infinity-norm condition number is required:
!>          = '1' or 'O':  1-norm;
!>          = 'I':         Infinity-norm.
!>
[in]N
!>          N is INTEGER
!>          The order 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]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>          Details of the LU factorization of the band matrix A, 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.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices; for 1 <= i <= N, row i of the matrix was
!>          interchanged with row IPIV(i).
!>
[in]ANORM
!>          ANORM is DOUBLE PRECISION
!>          If NORM = '1' or 'O', the 1-norm of the original matrix A.
!>          If NORM = 'I', the infinity-norm of the original matrix A.
!>
[out]RCOND
!>          RCOND is DOUBLE PRECISION
!>          The reciprocal of the condition number of the matrix A,
!>          computed as RCOND = 1/(norm(A) * norm(inv(A))).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (3*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 144 of file dgbcon.f.

146*
147* -- LAPACK computational routine --
148* -- LAPACK is a software package provided by Univ. of Tennessee, --
149* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150*
151* .. Scalar Arguments ..
152 CHARACTER NORM
153 INTEGER INFO, KL, KU, LDAB, N
154 DOUBLE PRECISION ANORM, RCOND
155* ..
156* .. Array Arguments ..
157 INTEGER IPIV( * ), IWORK( * )
158 DOUBLE PRECISION AB( LDAB, * ), WORK( * )
159* ..
160*
161* =====================================================================
162*
163* .. Parameters ..
164 DOUBLE PRECISION ONE, ZERO
165 parameter( one = 1.0d+0, zero = 0.0d+0 )
166* ..
167* .. Local Scalars ..
168 LOGICAL LNOTI, ONENRM
169 CHARACTER NORMIN
170 INTEGER IX, J, JP, KASE, KASE1, KD, LM
171 DOUBLE PRECISION AINVNM, SCALE, SMLNUM, T
172* ..
173* .. Local Arrays ..
174 INTEGER ISAVE( 3 )
175* ..
176* .. External Functions ..
177 LOGICAL LSAME
178 INTEGER IDAMAX
179 DOUBLE PRECISION DDOT, DLAMCH
180 EXTERNAL lsame, idamax, ddot, dlamch
181* ..
182* .. External Subroutines ..
183 EXTERNAL daxpy, dlacn2, dlatbs, drscl, xerbla
184* ..
185* .. Intrinsic Functions ..
186 INTRINSIC abs, min
187* ..
188* .. Executable Statements ..
189*
190* Test the input parameters.
191*
192 info = 0
193 onenrm = norm.EQ.'1' .OR. lsame( norm, 'O' )
194 IF( .NOT.onenrm .AND. .NOT.lsame( norm, 'I' ) ) THEN
195 info = -1
196 ELSE IF( n.LT.0 ) THEN
197 info = -2
198 ELSE IF( kl.LT.0 ) THEN
199 info = -3
200 ELSE IF( ku.LT.0 ) THEN
201 info = -4
202 ELSE IF( ldab.LT.2*kl+ku+1 ) THEN
203 info = -6
204 ELSE IF( anorm.LT.zero ) THEN
205 info = -8
206 END IF
207 IF( info.NE.0 ) THEN
208 CALL xerbla( 'DGBCON', -info )
209 RETURN
210 END IF
211*
212* Quick return if possible
213*
214 rcond = zero
215 IF( n.EQ.0 ) THEN
216 rcond = one
217 RETURN
218 ELSE IF( anorm.EQ.zero ) THEN
219 RETURN
220 END IF
221*
222 smlnum = dlamch( 'Safe minimum' )
223*
224* Estimate the norm of inv(A).
225*
226 ainvnm = zero
227 normin = 'N'
228 IF( onenrm ) THEN
229 kase1 = 1
230 ELSE
231 kase1 = 2
232 END IF
233 kd = kl + ku + 1
234 lnoti = kl.GT.0
235 kase = 0
236 10 CONTINUE
237 CALL dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
238 IF( kase.NE.0 ) THEN
239 IF( kase.EQ.kase1 ) THEN
240*
241* Multiply by inv(L).
242*
243 IF( lnoti ) THEN
244 DO 20 j = 1, n - 1
245 lm = min( kl, n-j )
246 jp = ipiv( j )
247 t = work( jp )
248 IF( jp.NE.j ) THEN
249 work( jp ) = work( j )
250 work( j ) = t
251 END IF
252 CALL daxpy( lm, -t, ab( kd+1, j ), 1, work( j+1 ), 1 )
253 20 CONTINUE
254 END IF
255*
256* Multiply by inv(U).
257*
258 CALL dlatbs( 'Upper', 'No transpose', 'Non-unit', normin, n,
259 $ kl+ku, ab, ldab, work, scale, work( 2*n+1 ),
260 $ info )
261 ELSE
262*
263* Multiply by inv(U**T).
264*
265 CALL dlatbs( 'Upper', 'Transpose', 'Non-unit', normin, n,
266 $ kl+ku, ab, ldab, work, scale, work( 2*n+1 ),
267 $ info )
268*
269* Multiply by inv(L**T).
270*
271 IF( lnoti ) THEN
272 DO 30 j = n - 1, 1, -1
273 lm = min( kl, n-j )
274 work( j ) = work( j ) - ddot( lm, ab( kd+1, j ), 1,
275 $ work( j+1 ), 1 )
276 jp = ipiv( j )
277 IF( jp.NE.j ) THEN
278 t = work( jp )
279 work( jp ) = work( j )
280 work( j ) = t
281 END IF
282 30 CONTINUE
283 END IF
284 END IF
285*
286* Divide X by 1/SCALE if doing so will not cause overflow.
287*
288 normin = 'Y'
289 IF( scale.NE.one ) THEN
290 ix = idamax( n, work, 1 )
291 IF( scale.LT.abs( work( ix ) )*smlnum .OR. scale.EQ.zero )
292 $ GO TO 40
293 CALL drscl( n, scale, work, 1 )
294 END IF
295 GO TO 10
296 END IF
297*
298* Compute the estimate of the reciprocal condition number.
299*
300 IF( ainvnm.NE.zero )
301 $ rcond = ( one / ainvnm ) / anorm
302*
303 40 CONTINUE
304 RETURN
305*
306* End of DGBCON
307*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
drscl
subroutine drscl(n, sa, sx, incx)
DRSCL multiplies a vector by the reciprocal of a real scalar.
Definition drscl.f:84
dlatbs
subroutine dlatbs(uplo, trans, diag, normin, n, kd, ab, ldab, x, scale, cnorm, info)
DLATBS solves a triangular banded system of equations.
Definition dlatbs.f:242
dlacn2
subroutine dlacn2(n, v, x, isgn, est, kase, isave)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition dlacn2.f:136
ddot
double precision function ddot(n, dx, incx, dy, incy)
DDOT
Definition ddot.f:82
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69

◆ dgbequ()

subroutine dgbequ ( integer m,
integer n,
integer kl,
integer ku,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( * ) r,
double precision, dimension( * ) c,
double precision rowcnd,
double precision colcnd,
double precision amax,
integer info )

DGBEQU

Download DGBEQU + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGBEQU computes row and column scalings intended to equilibrate an
!> M-by-N band matrix A and reduce its condition number.  R returns the
!> row scale factors and C the column scale factors, chosen to try to
!> make the largest element in each row and column of the matrix B with
!> elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.
!>
!> R(i) and C(j) are restricted to be between SMLNUM = smallest safe
!> number and BIGNUM = largest safe number.  Use of these scaling
!> factors is not guaranteed to reduce the condition number of A but
!> works well in practice.
!>
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]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>          The band matrix A, stored in rows 1 to KL+KU+1.  The j-th
!>          column of A is stored in the j-th column of the array AB as
!>          follows:
!>          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KL+KU+1.
!>
[out]R
!>          R is DOUBLE PRECISION array, dimension (M)
!>          If INFO = 0, or INFO > M, R contains the row scale factors
!>          for A.
!>
[out]C
!>          C is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0, C contains the column scale factors for A.
!>
[out]ROWCND
!>          ROWCND is DOUBLE PRECISION
!>          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
!>          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
!>          AMAX is neither too large nor too small, it is not worth
!>          scaling by R.
!>
[out]COLCND
!>          COLCND is DOUBLE PRECISION
!>          If INFO = 0, COLCND contains the ratio of the smallest
!>          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
!>          worth scaling by C.
!>
[out]AMAX
!>          AMAX is DOUBLE PRECISION
!>          Absolute value of largest matrix element.  If AMAX is very
!>          close to overflow or very close to underflow, the matrix
!>          should be scaled.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, and i is
!>                <= M:  the i-th row of A is exactly zero
!>                >  M:  the (i-M)-th column of A is exactly zero
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 151 of file dgbequ.f.

153*
154* -- LAPACK computational routine --
155* -- LAPACK is a software package provided by Univ. of Tennessee, --
156* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
157*
158* .. Scalar Arguments ..
159 INTEGER INFO, KL, KU, LDAB, M, N
160 DOUBLE PRECISION AMAX, COLCND, ROWCND
161* ..
162* .. Array Arguments ..
163 DOUBLE PRECISION AB( LDAB, * ), C( * ), R( * )
164* ..
165*
166* =====================================================================
167*
168* .. Parameters ..
169 DOUBLE PRECISION ONE, ZERO
170 parameter( one = 1.0d+0, zero = 0.0d+0 )
171* ..
172* .. Local Scalars ..
173 INTEGER I, J, KD
174 DOUBLE PRECISION BIGNUM, RCMAX, RCMIN, SMLNUM
175* ..
176* .. External Functions ..
177 DOUBLE PRECISION DLAMCH
178 EXTERNAL dlamch
179* ..
180* .. External Subroutines ..
181 EXTERNAL xerbla
182* ..
183* .. Intrinsic Functions ..
184 INTRINSIC abs, max, min
185* ..
186* .. Executable Statements ..
187*
188* Test the input parameters
189*
190 info = 0
191 IF( m.LT.0 ) THEN
192 info = -1
193 ELSE IF( n.LT.0 ) THEN
194 info = -2
195 ELSE IF( kl.LT.0 ) THEN
196 info = -3
197 ELSE IF( ku.LT.0 ) THEN
198 info = -4
199 ELSE IF( ldab.LT.kl+ku+1 ) THEN
200 info = -6
201 END IF
202 IF( info.NE.0 ) THEN
203 CALL xerbla( 'DGBEQU', -info )
204 RETURN
205 END IF
206*
207* Quick return if possible
208*
209 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
210 rowcnd = one
211 colcnd = one
212 amax = zero
213 RETURN
214 END IF
215*
216* Get machine constants.
217*
218 smlnum = dlamch( 'S' )
219 bignum = one / smlnum
220*
221* Compute row scale factors.
222*
223 DO 10 i = 1, m
224 r( i ) = zero
225 10 CONTINUE
226*
227* Find the maximum element in each row.
228*
229 kd = ku + 1
230 DO 30 j = 1, n
231 DO 20 i = max( j-ku, 1 ), min( j+kl, m )
232 r( i ) = max( r( i ), abs( ab( kd+i-j, j ) ) )
233 20 CONTINUE
234 30 CONTINUE
235*
236* Find the maximum and minimum scale factors.
237*
238 rcmin = bignum
239 rcmax = zero
240 DO 40 i = 1, m
241 rcmax = max( rcmax, r( i ) )
242 rcmin = min( rcmin, r( i ) )
243 40 CONTINUE
244 amax = rcmax
245*
246 IF( rcmin.EQ.zero ) THEN
247*
248* Find the first zero scale factor and return an error code.
249*
250 DO 50 i = 1, m
251 IF( r( i ).EQ.zero ) THEN
252 info = i
253 RETURN
254 END IF
255 50 CONTINUE
256 ELSE
257*
258* Invert the scale factors.
259*
260 DO 60 i = 1, m
261 r( i ) = one / min( max( r( i ), smlnum ), bignum )
262 60 CONTINUE
263*
264* Compute ROWCND = min(R(I)) / max(R(I))
265*
266 rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
267 END IF
268*
269* Compute column scale factors
270*
271 DO 70 j = 1, n
272 c( j ) = zero
273 70 CONTINUE
274*
275* Find the maximum element in each column,
276* assuming the row scaling computed above.
277*
278 kd = ku + 1
279 DO 90 j = 1, n
280 DO 80 i = max( j-ku, 1 ), min( j+kl, m )
281 c( j ) = max( c( j ), abs( ab( kd+i-j, j ) )*r( i ) )
282 80 CONTINUE
283 90 CONTINUE
284*
285* Find the maximum and minimum scale factors.
286*
287 rcmin = bignum
288 rcmax = zero
289 DO 100 j = 1, n
290 rcmin = min( rcmin, c( j ) )
291 rcmax = max( rcmax, c( j ) )
292 100 CONTINUE
293*
294 IF( rcmin.EQ.zero ) THEN
295*
296* Find the first zero scale factor and return an error code.
297*
298 DO 110 j = 1, n
299 IF( c( j ).EQ.zero ) THEN
300 info = m + j
301 RETURN
302 END IF
303 110 CONTINUE
304 ELSE
305*
306* Invert the scale factors.
307*
308 DO 120 j = 1, n
309 c( j ) = one / min( max( c( j ), smlnum ), bignum )
310 120 CONTINUE
311*
312* Compute COLCND = min(C(J)) / max(C(J))
313*
314 colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
315 END IF
316*
317 RETURN
318*
319* End of DGBEQU
320*

◆ dgbequb()

subroutine dgbequb ( integer m,
integer n,
integer kl,
integer ku,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( * ) r,
double precision, dimension( * ) c,
double precision rowcnd,
double precision colcnd,
double precision amax,
integer info )

DGBEQUB

Download DGBEQUB + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGBEQUB computes row and column scalings intended to equilibrate an
!> M-by-N matrix A and reduce its condition number.  R returns the row
!> scale factors and C the column scale factors, chosen to try to make
!> the largest element in each row and column of the matrix B with
!> elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most
!> the radix.
!>
!> R(i) and C(j) are restricted to be a power of the radix between
!> SMLNUM = smallest safe number and BIGNUM = largest safe number.  Use
!> of these scaling factors is not guaranteed to reduce the condition
!> number of A but works well in practice.
!>
!> This routine differs from DGEEQU by restricting the scaling factors
!> to a power of the radix.  Barring over- and underflow, scaling by
!> these factors introduces no additional rounding errors.  However, the
!> scaled entries' magnitudes are no longer approximately 1 but lie
!> between sqrt(radix) and 1/sqrt(radix).
!>
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]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>          On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
!>          The j-th column of A is stored in the j-th column of the
!>          array AB as follows:
!>          AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array A.  LDAB >= max(1,M).
!>
[out]R
!>          R is DOUBLE PRECISION array, dimension (M)
!>          If INFO = 0 or INFO > M, R contains the row scale factors
!>          for A.
!>
[out]C
!>          C is DOUBLE PRECISION array, dimension (N)
!>          If INFO = 0,  C contains the column scale factors for A.
!>
[out]ROWCND
!>          ROWCND is DOUBLE PRECISION
!>          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
!>          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
!>          AMAX is neither too large nor too small, it is not worth
!>          scaling by R.
!>
[out]COLCND
!>          COLCND is DOUBLE PRECISION
!>          If INFO = 0, COLCND contains the ratio of the smallest
!>          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
!>          worth scaling by C.
!>
[out]AMAX
!>          AMAX is DOUBLE PRECISION
!>          Absolute value of largest matrix element.  If AMAX is very
!>          close to overflow or very close to underflow, the matrix
!>          should be scaled.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i,  and i is
!>                <= M:  the i-th row of A is exactly zero
!>                >  M:  the (i-M)-th column of A is exactly zero
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 158 of file dgbequb.f.

160*
161* -- LAPACK computational routine --
162* -- LAPACK is a software package provided by Univ. of Tennessee, --
163* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
164*
165* .. Scalar Arguments ..
166 INTEGER INFO, KL, KU, LDAB, M, N
167 DOUBLE PRECISION AMAX, COLCND, ROWCND
168* ..
169* .. Array Arguments ..
170 DOUBLE PRECISION AB( LDAB, * ), C( * ), R( * )
171* ..
172*
173* =====================================================================
174*
175* .. Parameters ..
176 DOUBLE PRECISION ONE, ZERO
177 parameter( one = 1.0d+0, zero = 0.0d+0 )
178* ..
179* .. Local Scalars ..
180 INTEGER I, J, KD
181 DOUBLE PRECISION BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX
182* ..
183* .. External Functions ..
184 DOUBLE PRECISION DLAMCH
185 EXTERNAL dlamch
186* ..
187* .. External Subroutines ..
188 EXTERNAL xerbla
189* ..
190* .. Intrinsic Functions ..
191 INTRINSIC abs, max, min, log
192* ..
193* .. Executable Statements ..
194*
195* Test the input parameters.
196*
197 info = 0
198 IF( m.LT.0 ) THEN
199 info = -1
200 ELSE IF( n.LT.0 ) THEN
201 info = -2
202 ELSE IF( kl.LT.0 ) THEN
203 info = -3
204 ELSE IF( ku.LT.0 ) THEN
205 info = -4
206 ELSE IF( ldab.LT.kl+ku+1 ) THEN
207 info = -6
208 END IF
209 IF( info.NE.0 ) THEN
210 CALL xerbla( 'DGBEQUB', -info )
211 RETURN
212 END IF
213*
214* Quick return if possible.
215*
216 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
217 rowcnd = one
218 colcnd = one
219 amax = zero
220 RETURN
221 END IF
222*
223* Get machine constants. Assume SMLNUM is a power of the radix.
224*
225 smlnum = dlamch( 'S' )
226 bignum = one / smlnum
227 radix = dlamch( 'B' )
228 logrdx = log(radix)
229*
230* Compute row scale factors.
231*
232 DO 10 i = 1, m
233 r( i ) = zero
234 10 CONTINUE
235*
236* Find the maximum element in each row.
237*
238 kd = ku + 1
239 DO 30 j = 1, n
240 DO 20 i = max( j-ku, 1 ), min( j+kl, m )
241 r( i ) = max( r( i ), abs( ab( kd+i-j, j ) ) )
242 20 CONTINUE
243 30 CONTINUE
244 DO i = 1, m
245 IF( r( i ).GT.zero ) THEN
246 r( i ) = radix**int( log( r( i ) ) / logrdx )
247 END IF
248 END DO
249*
250* Find the maximum and minimum scale factors.
251*
252 rcmin = bignum
253 rcmax = zero
254 DO 40 i = 1, m
255 rcmax = max( rcmax, r( i ) )
256 rcmin = min( rcmin, r( i ) )
257 40 CONTINUE
258 amax = rcmax
259*
260 IF( rcmin.EQ.zero ) THEN
261*
262* Find the first zero scale factor and return an error code.
263*
264 DO 50 i = 1, m
265 IF( r( i ).EQ.zero ) THEN
266 info = i
267 RETURN
268 END IF
269 50 CONTINUE
270 ELSE
271*
272* Invert the scale factors.
273*
274 DO 60 i = 1, m
275 r( i ) = one / min( max( r( i ), smlnum ), bignum )
276 60 CONTINUE
277*
278* Compute ROWCND = min(R(I)) / max(R(I)).
279*
280 rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
281 END IF
282*
283* Compute column scale factors.
284*
285 DO 70 j = 1, n
286 c( j ) = zero
287 70 CONTINUE
288*
289* Find the maximum element in each column,
290* assuming the row scaling computed above.
291*
292 DO 90 j = 1, n
293 DO 80 i = max( j-ku, 1 ), min( j+kl, m )
294 c( j ) = max( c( j ), abs( ab( kd+i-j, j ) )*r( i ) )
295 80 CONTINUE
296 IF( c( j ).GT.zero ) THEN
297 c( j ) = radix**int( log( c( j ) ) / logrdx )
298 END IF
299 90 CONTINUE
300*
301* Find the maximum and minimum scale factors.
302*
303 rcmin = bignum
304 rcmax = zero
305 DO 100 j = 1, n
306 rcmin = min( rcmin, c( j ) )
307 rcmax = max( rcmax, c( j ) )
308 100 CONTINUE
309*
310 IF( rcmin.EQ.zero ) THEN
311*
312* Find the first zero scale factor and return an error code.
313*
314 DO 110 j = 1, n
315 IF( c( j ).EQ.zero ) THEN
316 info = m + j
317 RETURN
318 END IF
319 110 CONTINUE
320 ELSE
321*
322* Invert the scale factors.
323*
324 DO 120 j = 1, n
325 c( j ) = one / min( max( c( j ), smlnum ), bignum )
326 120 CONTINUE
327*
328* Compute COLCND = min(C(J)) / max(C(J)).
329*
330 colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
331 END IF
332*
333 RETURN
334*
335* End of DGBEQUB
336*

◆ dgbrfs()

subroutine dgbrfs ( character trans,
integer n,
integer kl,
integer ku,
integer nrhs,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( ldafb, * ) afb,
integer ldafb,
integer, dimension( * ) ipiv,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( ldx, * ) x,
integer ldx,
double precision, dimension( * ) ferr,
double precision, dimension( * ) berr,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

DGBRFS

Download DGBRFS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGBRFS improves the computed solution to a system of linear
!> equations when the coefficient matrix is banded, and provides
!> error bounds and backward error estimates for the solution.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>          Specifies the form of the system of equations:
!>          = 'N':  A * X = B     (No transpose)
!>          = 'T':  A**T * X = B  (Transpose)
!>          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
!>
[in]N
!>          N is INTEGER
!>          The order 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]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrices B and X.  NRHS >= 0.
!>
[in]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>          The original band matrix A, stored in rows 1 to KL+KU+1.
!>          The j-th column of A is stored in the j-th column of the
!>          array AB as follows:
!>          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl).
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KL+KU+1.
!>
[in]AFB
!>          AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
!>          Details of the LU factorization of the band matrix A, 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.
!>
[in]LDAFB
!>          LDAFB is INTEGER
!>          The leading dimension of the array AFB.  LDAFB >= 2*KL*KU+1.
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices from DGBTRF; for 1<=i<=N, row i of the
!>          matrix was interchanged with row IPIV(i).
!>
[in]B
!>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
!>          The right hand side matrix B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[in,out]X
!>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
!>          On entry, the solution matrix X, as computed by DGBTRS.
!>          On exit, the improved solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]FERR
!>          FERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (3*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>
Internal Parameters:
!>  ITMAX is the maximum number of steps of iterative refinement.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 202 of file dgbrfs.f.

205*
206* -- LAPACK computational routine --
207* -- LAPACK is a software package provided by Univ. of Tennessee, --
208* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
209*
210* .. Scalar Arguments ..
211 CHARACTER TRANS
212 INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
213* ..
214* .. Array Arguments ..
215 INTEGER IPIV( * ), IWORK( * )
216 DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
217 $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
218* ..
219*
220* =====================================================================
221*
222* .. Parameters ..
223 INTEGER ITMAX
224 parameter( itmax = 5 )
225 DOUBLE PRECISION ZERO
226 parameter( zero = 0.0d+0 )
227 DOUBLE PRECISION ONE
228 parameter( one = 1.0d+0 )
229 DOUBLE PRECISION TWO
230 parameter( two = 2.0d+0 )
231 DOUBLE PRECISION THREE
232 parameter( three = 3.0d+0 )
233* ..
234* .. Local Scalars ..
235 LOGICAL NOTRAN
236 CHARACTER TRANST
237 INTEGER COUNT, I, J, K, KASE, KK, NZ
238 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
239* ..
240* .. Local Arrays ..
241 INTEGER ISAVE( 3 )
242* ..
243* .. External Subroutines ..
244 EXTERNAL daxpy, dcopy, dgbmv, dgbtrs, dlacn2, xerbla
245* ..
246* .. Intrinsic Functions ..
247 INTRINSIC abs, max, min
248* ..
249* .. External Functions ..
250 LOGICAL LSAME
251 DOUBLE PRECISION DLAMCH
252 EXTERNAL lsame, dlamch
253* ..
254* .. Executable Statements ..
255*
256* Test the input parameters.
257*
258 info = 0
259 notran = lsame( trans, 'N' )
260 IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
261 $ lsame( trans, 'C' ) ) THEN
262 info = -1
263 ELSE IF( n.LT.0 ) THEN
264 info = -2
265 ELSE IF( kl.LT.0 ) THEN
266 info = -3
267 ELSE IF( ku.LT.0 ) THEN
268 info = -4
269 ELSE IF( nrhs.LT.0 ) THEN
270 info = -5
271 ELSE IF( ldab.LT.kl+ku+1 ) THEN
272 info = -7
273 ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
274 info = -9
275 ELSE IF( ldb.LT.max( 1, n ) ) THEN
276 info = -12
277 ELSE IF( ldx.LT.max( 1, n ) ) THEN
278 info = -14
279 END IF
280 IF( info.NE.0 ) THEN
281 CALL xerbla( 'DGBRFS', -info )
282 RETURN
283 END IF
284*
285* Quick return if possible
286*
287 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
288 DO 10 j = 1, nrhs
289 ferr( j ) = zero
290 berr( j ) = zero
291 10 CONTINUE
292 RETURN
293 END IF
294*
295 IF( notran ) THEN
296 transt = 'T'
297 ELSE
298 transt = 'N'
299 END IF
300*
301* NZ = maximum number of nonzero elements in each row of A, plus 1
302*
303 nz = min( kl+ku+2, n+1 )
304 eps = dlamch( 'Epsilon' )
305 safmin = dlamch( 'Safe minimum' )
306 safe1 = nz*safmin
307 safe2 = safe1 / eps
308*
309* Do for each right hand side
310*
311 DO 140 j = 1, nrhs
312*
313 count = 1
314 lstres = three
315 20 CONTINUE
316*
317* Loop until stopping criterion is satisfied.
318*
319* Compute residual R = B - op(A) * X,
320* where op(A) = A, A**T, or A**H, depending on TRANS.
321*
322 CALL dcopy( n, b( 1, j ), 1, work( n+1 ), 1 )
323 CALL dgbmv( trans, n, n, kl, ku, -one, ab, ldab, x( 1, j ), 1,
324 $ one, work( n+1 ), 1 )
325*
326* Compute componentwise relative backward error from formula
327*
328* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
329*
330* where abs(Z) is the componentwise absolute value of the matrix
331* or vector Z. If the i-th component of the denominator is less
332* than SAFE2, then SAFE1 is added to the i-th components of the
333* numerator and denominator before dividing.
334*
335 DO 30 i = 1, n
336 work( i ) = abs( b( i, j ) )
337 30 CONTINUE
338*
339* Compute abs(op(A))*abs(X) + abs(B).
340*
341 IF( notran ) THEN
342 DO 50 k = 1, n
343 kk = ku + 1 - k
344 xk = abs( x( k, j ) )
345 DO 40 i = max( 1, k-ku ), min( n, k+kl )
346 work( i ) = work( i ) + abs( ab( kk+i, k ) )*xk
347 40 CONTINUE
348 50 CONTINUE
349 ELSE
350 DO 70 k = 1, n
351 s = zero
352 kk = ku + 1 - k
353 DO 60 i = max( 1, k-ku ), min( n, k+kl )
354 s = s + abs( ab( kk+i, k ) )*abs( x( i, j ) )
355 60 CONTINUE
356 work( k ) = work( k ) + s
357 70 CONTINUE
358 END IF
359 s = zero
360 DO 80 i = 1, n
361 IF( work( i ).GT.safe2 ) THEN
362 s = max( s, abs( work( n+i ) ) / work( i ) )
363 ELSE
364 s = max( s, ( abs( work( n+i ) )+safe1 ) /
365 $ ( work( i )+safe1 ) )
366 END IF
367 80 CONTINUE
368 berr( j ) = s
369*
370* Test stopping criterion. Continue iterating if
371* 1) The residual BERR(J) is larger than machine epsilon, and
372* 2) BERR(J) decreased by at least a factor of 2 during the
373* last iteration, and
374* 3) At most ITMAX iterations tried.
375*
376 IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
377 $ count.LE.itmax ) THEN
378*
379* Update solution and try again.
380*
381 CALL dgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv,
382 $ work( n+1 ), n, info )
383 CALL daxpy( n, one, work( n+1 ), 1, x( 1, j ), 1 )
384 lstres = berr( j )
385 count = count + 1
386 GO TO 20
387 END IF
388*
389* Bound error from formula
390*
391* norm(X - XTRUE) / norm(X) .le. FERR =
392* norm( abs(inv(op(A)))*
393* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
394*
395* where
396* norm(Z) is the magnitude of the largest component of Z
397* inv(op(A)) is the inverse of op(A)
398* abs(Z) is the componentwise absolute value of the matrix or
399* vector Z
400* NZ is the maximum number of nonzeros in any row of A, plus 1
401* EPS is machine epsilon
402*
403* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
404* is incremented by SAFE1 if the i-th component of
405* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
406*
407* Use DLACN2 to estimate the infinity-norm of the matrix
408* inv(op(A)) * diag(W),
409* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
410*
411 DO 90 i = 1, n
412 IF( work( i ).GT.safe2 ) THEN
413 work( i ) = abs( work( n+i ) ) + nz*eps*work( i )
414 ELSE
415 work( i ) = abs( work( n+i ) ) + nz*eps*work( i ) + safe1
416 END IF
417 90 CONTINUE
418*
419 kase = 0
420 100 CONTINUE
421 CALL dlacn2( n, work( 2*n+1 ), work( n+1 ), iwork, ferr( j ),
422 $ kase, isave )
423 IF( kase.NE.0 ) THEN
424 IF( kase.EQ.1 ) THEN
425*
426* Multiply by diag(W)*inv(op(A)**T).
427*
428 CALL dgbtrs( transt, n, kl, ku, 1, afb, ldafb, ipiv,
429 $ work( n+1 ), n, info )
430 DO 110 i = 1, n
431 work( n+i ) = work( n+i )*work( i )
432 110 CONTINUE
433 ELSE
434*
435* Multiply by inv(op(A))*diag(W).
436*
437 DO 120 i = 1, n
438 work( n+i ) = work( n+i )*work( i )
439 120 CONTINUE
440 CALL dgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv,
441 $ work( n+1 ), n, info )
442 END IF
443 GO TO 100
444 END IF
445*
446* Normalize error.
447*
448 lstres = zero
449 DO 130 i = 1, n
450 lstres = max( lstres, abs( x( i, j ) ) )
451 130 CONTINUE
452 IF( lstres.NE.zero )
453 $ ferr( j ) = ferr( j ) / lstres
454*
455 140 CONTINUE
456*
457 RETURN
458*
459* End of DGBRFS
460*
dgbtrs
subroutine dgbtrs(trans, n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
DGBTRS
Definition dgbtrs.f:138
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgbmv
subroutine dgbmv(trans, m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy)
DGBMV
Definition dgbmv.f:185

◆ dgbrfsx()

subroutine dgbrfsx ( character trans,
character equed,
integer n,
integer kl,
integer ku,
integer nrhs,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( ldafb, * ) afb,
integer ldafb,
integer, dimension( * ) ipiv,
double precision, dimension( * ) r,
double precision, dimension( * ) c,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( ldx , * ) x,
integer ldx,
double precision rcond,
double precision, dimension( * ) berr,
integer n_err_bnds,
double precision, dimension( nrhs, * ) err_bnds_norm,
double precision, dimension( nrhs, * ) err_bnds_comp,
integer nparams,
double precision, dimension( * ) params,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

DGBRFSX

Download DGBRFSX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!>    DGBRFSX improves the computed solution to a system of linear
!>    equations and provides error bounds and backward error estimates
!>    for the solution.  In addition to normwise error bound, the code
!>    provides maximum componentwise error bound if possible.  See
!>    comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the
!>    error bounds.
!>
!>    The original system of linear equations may have been equilibrated
!>    before calling this routine, as described by arguments EQUED, R
!>    and C below. In this case, the solution and error bounds returned
!>    are for the original unequilibrated system.
!> 
!>     Some optional parameters are bundled in the PARAMS array.  These
!>     settings determine how refinement is performed, but often the
!>     defaults are acceptable.  If the defaults are acceptable, users
!>     can pass NPARAMS = 0 which prevents the source code from accessing
!>     the PARAMS argument.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>     Specifies the form of the system of equations:
!>       = 'N':  A * X = B     (No transpose)
!>       = 'T':  A**T * X = B  (Transpose)
!>       = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
!>
[in]EQUED
!>          EQUED is CHARACTER*1
!>     Specifies the form of equilibration that was done to A
!>     before calling this routine. This is needed to compute
!>     the solution and error bounds correctly.
!>       = 'N':  No equilibration
!>       = 'R':  Row equilibration, i.e., A has been premultiplied by
!>               diag(R).
!>       = 'C':  Column equilibration, i.e., A has been postmultiplied
!>               by diag(C).
!>       = 'B':  Both row and column equilibration, i.e., A has been
!>               replaced by diag(R) * A * diag(C).
!>               The right hand side B has been changed accordingly.
!>
[in]N
!>          N is INTEGER
!>     The order 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]NRHS
!>          NRHS is INTEGER
!>     The number of right hand sides, i.e., the number of columns
!>     of the matrices B and X.  NRHS >= 0.
!>
[in]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>     The original band matrix A, stored in rows 1 to KL+KU+1.
!>     The j-th column of A is stored in the j-th column of the
!>     array AB as follows:
!>     AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(n,j+kl).
!>
[in]LDAB
!>          LDAB is INTEGER
!>     The leading dimension of the array AB.  LDAB >= KL+KU+1.
!>
[in]AFB
!>          AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
!>     Details of the LU factorization of the band matrix A, 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.
!>
[in]LDAFB
!>          LDAFB is INTEGER
!>     The leading dimension of the array AFB.  LDAFB >= 2*KL*KU+1.
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>     The pivot indices from DGETRF; for 1<=i<=N, row i of the
!>     matrix was interchanged with row IPIV(i).
!>
[in,out]R
!>          R is DOUBLE PRECISION array, dimension (N)
!>     The row scale factors for A.  If EQUED = 'R' or 'B', A is
!>     multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
!>     is not accessed.  R is an input argument if FACT = 'F';
!>     otherwise, R is an output argument.  If FACT = 'F' and
!>     EQUED = 'R' or 'B', each element of R must be positive.
!>     If R is output, each element of R is a power of the radix.
!>     If R is input, each element of R should be a power of the radix
!>     to ensure a reliable solution and error estimates. Scaling by
!>     powers of the radix does not cause rounding errors unless the
!>     result underflows or overflows. Rounding errors during scaling
!>     lead to refining with a matrix that is not equivalent to the
!>     input matrix, producing error estimates that may not be
!>     reliable.
!>
[in,out]C
!>          C is DOUBLE PRECISION array, dimension (N)
!>     The column scale factors for A.  If EQUED = 'C' or 'B', A is
!>     multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
!>     is not accessed.  C is an input argument if FACT = 'F';
!>     otherwise, C is an output argument.  If FACT = 'F' and
!>     EQUED = 'C' or 'B', each element of C must be positive.
!>     If C is output, each element of C is a power of the radix.
!>     If C is input, each element of C should be a power of the radix
!>     to ensure a reliable solution and error estimates. Scaling by
!>     powers of the radix does not cause rounding errors unless the
!>     result underflows or overflows. Rounding errors during scaling
!>     lead to refining with a matrix that is not equivalent to the
!>     input matrix, producing error estimates that may not be
!>     reliable.
!>
[in]B
!>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
!>     The right hand side matrix B.
!>
[in]LDB
!>          LDB is INTEGER
!>     The leading dimension of the array B.  LDB >= max(1,N).
!>
[in,out]X
!>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
!>     On entry, the solution matrix X, as computed by DGETRS.
!>     On exit, the improved solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>     The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is DOUBLE PRECISION
!>     Reciprocal scaled condition number.  This is an estimate of the
!>     reciprocal Skeel condition number of the matrix A after
!>     equilibration (if done).  If this is less than the machine
!>     precision (in particular, if it is zero), the matrix is singular
!>     to working precision.  Note that the error may still be small even
!>     if this number is very small and the matrix appears ill-
!>     conditioned.
!>
[out]BERR
!>          BERR is DOUBLE PRECISION array, dimension (NRHS)
!>     Componentwise relative backward error.  This is the
!>     componentwise relative backward error of each solution vector X(j)
!>     (i.e., the smallest relative change in any element of A or B that
!>     makes X(j) an exact solution).
!>
[in]N_ERR_BNDS
!>          N_ERR_BNDS is INTEGER
!>     Number of error bounds to return for each right hand side
!>     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
!>     ERR_BNDS_COMP below.
!>
[out]ERR_BNDS_NORM
!>          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     normwise relative error, which is defined as follows:
!>
!>     Normwise relative error in the ith solution vector:
!>             max_j (abs(XTRUE(j,i) - X(j,i)))
!>            ------------------------------
!>                  max_j abs(X(j,i))
!>
!>     The array is indexed by the type of error information as described
!>     below. There currently are up to three pieces of information
!>     returned.
!>
!>     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_NORM(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * dlamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * dlamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated normwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * dlamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*A, where S scales each row by a power of the
!>              radix so all absolute row sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[out]ERR_BNDS_COMP
!>          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     componentwise relative error, which is defined as follows:
!>
!>     Componentwise relative error in the ith solution vector:
!>                    abs(XTRUE(j,i) - X(j,i))
!>             max_j ----------------------
!>                         abs(X(j,i))
!>
!>     The array is indexed by the right-hand side i (on which the
!>     componentwise relative error depends), and the type of error
!>     information as described below. There currently are up to three
!>     pieces of information returned for each right-hand side. If
!>     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
!>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS < 3, then at most
!>     the first (:,N_ERR_BNDS) entries are returned.
!>
!>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_COMP(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * dlamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * dlamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated componentwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * dlamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*(A*diag(x)), where x is the solution for the
!>              current right-hand side and S scales each row of
!>              A*diag(x) by a power of the radix so all absolute row
!>              sums of Z are approximately 1.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in]NPARAMS
!>          NPARAMS is INTEGER
!>     Specifies the number of parameters set in PARAMS.  If <= 0, the
!>     PARAMS array is never referenced and default values are used.
!>
[in,out]PARAMS
!>          PARAMS is DOUBLE PRECISION array, dimension (NPARAMS)
!>     Specifies algorithm parameters.  If an entry is < 0.0, then
!>     that entry will be filled with default value used for that
!>     parameter.  Only positions up to NPARAMS are accessed; defaults
!>     are used for higher-numbered parameters.
!>
!>       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
!>            refinement or not.
!>         Default: 1.0D+0
!>            = 0.0:  No refinement is performed, and no error bounds are
!>                    computed.
!>            = 1.0:  Use the double-precision refinement algorithm,
!>                    possibly with doubled-single computations if the
!>                    compilation environment does not support DOUBLE
!>                    PRECISION.
!>              (other values are reserved for future use)
!>
!>       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
!>            computations allowed for refinement.
!>         Default: 10
!>         Aggressive: Set to 100 to permit convergence using approximate
!>                     factorizations or factorizations other than LU. If
!>                     the factorization uses a technique other than
!>                     Gaussian elimination, the guarantees in
!>                     err_bnds_norm and err_bnds_comp may no longer be
!>                     trustworthy.
!>
!>       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
!>            will attempt to find a solution with small componentwise
!>            relative error in the double-precision algorithm.  Positive
!>            is true, 0.0 is false.
!>         Default: 1.0 (attempt componentwise convergence)
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (4*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit. The solution to every right-hand side is
!>         guaranteed.
!>       < 0:  If INFO = -i, the i-th argument had an illegal value
!>       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization
!>         has been completed, but the factor U is exactly singular, so
!>         the solution and error bounds could not be computed. RCOND = 0
!>         is returned.
!>       = N+J: The solution corresponding to the Jth right-hand side is
!>         not guaranteed. The solutions corresponding to other right-
!>         hand sides K with K > J may not be guaranteed as well, but
!>         only the first such right-hand side is reported. If a small
!>         componentwise error is not requested (PARAMS(3) = 0.0) then
!>         the Jth right-hand side is the first with a normwise error
!>         bound that is not guaranteed (the smallest J such
!>         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
!>         the Jth right-hand side is the first with either a normwise or
!>         componentwise error bound that is not guaranteed (the smallest
!>         J such that either ERR_BNDS_NORM(J,1) = 0.0 or
!>         ERR_BNDS_COMP(J,1) = 0.0). See the definition of
!>         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
!>         about all of the right-hand sides check ERR_BNDS_NORM or
!>         ERR_BNDS_COMP.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 435 of file dgbrfsx.f.

440*
441* -- LAPACK computational routine --
442* -- LAPACK is a software package provided by Univ. of Tennessee, --
443* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
444*
445* .. Scalar Arguments ..
446 CHARACTER TRANS, EQUED
447 INTEGER INFO, LDAB, LDAFB, LDB, LDX, N, KL, KU, NRHS,
448 $ NPARAMS, N_ERR_BNDS
449 DOUBLE PRECISION RCOND
450* ..
451* .. Array Arguments ..
452 INTEGER IPIV( * ), IWORK( * )
453 DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
454 $ X( LDX , * ),WORK( * )
455 DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ),
456 $ ERR_BNDS_NORM( NRHS, * ),
457 $ ERR_BNDS_COMP( NRHS, * )
458* ..
459*
460* ==================================================================
461*
462* .. Parameters ..
463 DOUBLE PRECISION ZERO, ONE
464 parameter( zero = 0.0d+0, one = 1.0d+0 )
465 DOUBLE PRECISION ITREF_DEFAULT, ITHRESH_DEFAULT
466 DOUBLE PRECISION COMPONENTWISE_DEFAULT, RTHRESH_DEFAULT
467 DOUBLE PRECISION DZTHRESH_DEFAULT
468 parameter( itref_default = 1.0d+0 )
469 parameter( ithresh_default = 10.0d+0 )
470 parameter( componentwise_default = 1.0d+0 )
471 parameter( rthresh_default = 0.5d+0 )
472 parameter( dzthresh_default = 0.25d+0 )
473 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
474 $ LA_LINRX_CWISE_I
475 parameter( la_linrx_itref_i = 1,
476 $ la_linrx_ithresh_i = 2 )
477 parameter( la_linrx_cwise_i = 3 )
478 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
479 $ LA_LINRX_RCOND_I
480 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
481 parameter( la_linrx_rcond_i = 3 )
482* ..
483* .. Local Scalars ..
484 CHARACTER(1) NORM
485 LOGICAL ROWEQU, COLEQU, NOTRAN
486 INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE
487 INTEGER N_NORMS
488 DOUBLE PRECISION ANORM, RCOND_TMP
489 DOUBLE PRECISION ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
490 LOGICAL IGNORE_CWISE
491 INTEGER ITHRESH
492 DOUBLE PRECISION RTHRESH, UNSTABLE_THRESH
493* ..
494* .. External Subroutines ..
495 EXTERNAL xerbla, dgbcon
496 EXTERNAL dla_gbrfsx_extended
497* ..
498* .. Intrinsic Functions ..
499 INTRINSIC max, sqrt
500* ..
501* .. External Functions ..
502 EXTERNAL lsame, ilatrans, ilaprec
503 EXTERNAL dlamch, dlangb, dla_gbrcond
504 DOUBLE PRECISION DLAMCH, DLANGB, DLA_GBRCOND
505 LOGICAL LSAME
506 INTEGER ILATRANS, ILAPREC
507* ..
508* .. Executable Statements ..
509*
510* Check the input parameters.
511*
512 info = 0
513 trans_type = ilatrans( trans )
514 ref_type = int( itref_default )
515 IF ( nparams .GE. la_linrx_itref_i ) THEN
516 IF ( params( la_linrx_itref_i ) .LT. 0.0d+0 ) THEN
517 params( la_linrx_itref_i ) = itref_default
518 ELSE
519 ref_type = params( la_linrx_itref_i )
520 END IF
521 END IF
522*
523* Set default parameters.
524*
525 illrcond_thresh = dble( n ) * dlamch( 'Epsilon' )
526 ithresh = int( ithresh_default )
527 rthresh = rthresh_default
528 unstable_thresh = dzthresh_default
529 ignore_cwise = componentwise_default .EQ. 0.0d+0
530*
531 IF ( nparams.GE.la_linrx_ithresh_i ) THEN
532 IF ( params( la_linrx_ithresh_i ).LT.0.0d+0 ) THEN
533 params( la_linrx_ithresh_i ) = ithresh
534 ELSE
535 ithresh = int( params( la_linrx_ithresh_i ) )
536 END IF
537 END IF
538 IF ( nparams.GE.la_linrx_cwise_i ) THEN
539 IF ( params( la_linrx_cwise_i ).LT.0.0d+0 ) THEN
540 IF ( ignore_cwise ) THEN
541 params( la_linrx_cwise_i ) = 0.0d+0
542 ELSE
543 params( la_linrx_cwise_i ) = 1.0d+0
544 END IF
545 ELSE
546 ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0d+0
547 END IF
548 END IF
549 IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
550 n_norms = 0
551 ELSE IF ( ignore_cwise ) THEN
552 n_norms = 1
553 ELSE
554 n_norms = 2
555 END IF
556*
557 notran = lsame( trans, 'N' )
558 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
559 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
560*
561* Test input parameters.
562*
563 IF( trans_type.EQ.-1 ) THEN
564 info = -1
565 ELSE IF( .NOT.rowequ .AND. .NOT.colequ .AND.
566 $ .NOT.lsame( equed, 'N' ) ) THEN
567 info = -2
568 ELSE IF( n.LT.0 ) THEN
569 info = -3
570 ELSE IF( kl.LT.0 ) THEN
571 info = -4
572 ELSE IF( ku.LT.0 ) THEN
573 info = -5
574 ELSE IF( nrhs.LT.0 ) THEN
575 info = -6
576 ELSE IF( ldab.LT.kl+ku+1 ) THEN
577 info = -8
578 ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
579 info = -10
580 ELSE IF( ldb.LT.max( 1, n ) ) THEN
581 info = -13
582 ELSE IF( ldx.LT.max( 1, n ) ) THEN
583 info = -15
584 END IF
585 IF( info.NE.0 ) THEN
586 CALL xerbla( 'DGBRFSX', -info )
587 RETURN
588 END IF
589*
590* Quick return if possible.
591*
592 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
593 rcond = 1.0d+0
594 DO j = 1, nrhs
595 berr( j ) = 0.0d+0
596 IF ( n_err_bnds .GE. 1 ) THEN
597 err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
598 err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
599 END IF
600 IF ( n_err_bnds .GE. 2 ) THEN
601 err_bnds_norm( j, la_linrx_err_i ) = 0.0d+0
602 err_bnds_comp( j, la_linrx_err_i ) = 0.0d+0
603 END IF
604 IF ( n_err_bnds .GE. 3 ) THEN
605 err_bnds_norm( j, la_linrx_rcond_i ) = 1.0d+0
606 err_bnds_comp( j, la_linrx_rcond_i ) = 1.0d+0
607 END IF
608 END DO
609 RETURN
610 END IF
611*
612* Default to failure.
613*
614 rcond = 0.0d+0
615 DO j = 1, nrhs
616 berr( j ) = 1.0d+0
617 IF ( n_err_bnds .GE. 1 ) THEN
618 err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
619 err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
620 END IF
621 IF ( n_err_bnds .GE. 2 ) THEN
622 err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
623 err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
624 END IF
625 IF ( n_err_bnds .GE. 3 ) THEN
626 err_bnds_norm( j, la_linrx_rcond_i ) = 0.0d+0
627 err_bnds_comp( j, la_linrx_rcond_i ) = 0.0d+0
628 END IF
629 END DO
630*
631* Compute the norm of A and the reciprocal of the condition
632* number of A.
633*
634 IF( notran ) THEN
635 norm = 'I'
636 ELSE
637 norm = '1'
638 END IF
639 anorm = dlangb( norm, n, kl, ku, ab, ldab, work )
640 CALL dgbcon( norm, n, kl, ku, afb, ldafb, ipiv, anorm, rcond,
641 $ work, iwork, info )
642*
643* Perform refinement on each right-hand side
644*
645 IF ( ref_type .NE. 0 .AND. info .EQ. 0 ) THEN
646
647 prec_type = ilaprec( 'E' )
648
649 IF ( notran ) THEN
650 CALL dla_gbrfsx_extended( prec_type, trans_type, n, kl, ku,
651 $ nrhs, ab, ldab, afb, ldafb, ipiv, colequ, c, b,
652 $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
653 $ err_bnds_comp, work( n+1 ), work( 1 ), work( 2*n+1 ),
654 $ work( 1 ), rcond, ithresh, rthresh, unstable_thresh,
655 $ ignore_cwise, info )
656 ELSE
657 CALL dla_gbrfsx_extended( prec_type, trans_type, n, kl, ku,
658 $ nrhs, ab, ldab, afb, ldafb, ipiv, rowequ, r, b,
659 $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
660 $ err_bnds_comp, work( n+1 ), work( 1 ), work( 2*n+1 ),
661 $ work( 1 ), rcond, ithresh, rthresh, unstable_thresh,
662 $ ignore_cwise, info )
663 END IF
664 END IF
665
666 err_lbnd = max( 10.0d+0, sqrt( dble( n ) ) ) * dlamch( 'Epsilon' )
667 IF ( n_err_bnds .GE. 1 .AND. n_norms .GE. 1 ) THEN
668*
669* Compute scaled normwise condition number cond(A*C).
670*
671 IF ( colequ .AND. notran ) THEN
672 rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
673 $ ldafb, ipiv, -1, c, info, work, iwork )
674 ELSE IF ( rowequ .AND. .NOT. notran ) THEN
675 rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
676 $ ldafb, ipiv, -1, r, info, work, iwork )
677 ELSE
678 rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
679 $ ldafb, ipiv, 0, r, info, work, iwork )
680 END IF
681 DO j = 1, nrhs
682*
683* Cap the error at 1.0.
684*
685 IF ( n_err_bnds .GE. la_linrx_err_i
686 $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0d+0 )
687 $ err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
688*
689* Threshold the error (see LAWN).
690*
691 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
692 err_bnds_norm( j, la_linrx_err_i ) = 1.0d+0
693 err_bnds_norm( j, la_linrx_trust_i ) = 0.0d+0
694 IF ( info .LE. n ) info = n + j
695 ELSE IF ( err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd )
696 $ THEN
697 err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
698 err_bnds_norm( j, la_linrx_trust_i ) = 1.0d+0
699 END IF
700*
701* Save the condition number.
702*
703 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
704 err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
705 END IF
706
707 END DO
708 END IF
709
710 IF (n_err_bnds .GE. 1 .AND. n_norms .GE. 2) THEN
711*
712* Compute componentwise condition number cond(A*diag(Y(:,J))) for
713* each right-hand side using the current solution as an estimate of
714* the true solution. If the componentwise error estimate is too
715* large, then the solution is a lousy estimate of truth and the
716* estimated RCOND may be too optimistic. To avoid misleading users,
717* the inverse condition number is set to 0.0 when the estimated
718* cwise error is at least CWISE_WRONG.
719*
720 cwise_wrong = sqrt( dlamch( 'Epsilon' ) )
721 DO j = 1, nrhs
722 IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
723 $ THEN
724 rcond_tmp = dla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
725 $ ldafb, ipiv, 1, x( 1, j ), info, work, iwork )
726 ELSE
727 rcond_tmp = 0.0d+0
728 END IF
729*
730* Cap the error at 1.0.
731*
732 IF ( n_err_bnds .GE. la_linrx_err_i
733 $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0d+0 )
734 $ err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
735*
736* Threshold the error (see LAWN).
737*
738 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
739 err_bnds_comp( j, la_linrx_err_i ) = 1.0d+0
740 err_bnds_comp( j, la_linrx_trust_i ) = 0.0d+0
741 IF ( params( la_linrx_cwise_i ) .EQ. 1.0d+0
742 $ .AND. info.LT.n + j ) info = n + j
743 ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
744 $ .LT. err_lbnd ) THEN
745 err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
746 err_bnds_comp( j, la_linrx_trust_i ) = 1.0d+0
747 END IF
748*
749* Save the condition number.
750*
751 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
752 err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
753 END IF
754
755 END DO
756 END IF
757*
758 RETURN
759*
760* End of DGBRFSX
761*
ilaprec
integer function ilaprec(prec)
ILAPREC
Definition ilaprec.f:58
ilatrans
integer function ilatrans(trans)
ILATRANS
Definition ilatrans.f:58
dlangb
double precision function dlangb(norm, n, kl, ku, ab, ldab, work)
DLANGB returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlangb.f:124
dla_gbrcond
double precision function dla_gbrcond(trans, n, kl, ku, ab, ldab, afb, ldafb, ipiv, cmode, c, info, work, iwork)
DLA_GBRCOND estimates the Skeel condition number for a general banded matrix.
Definition dla_gbrcond.f:170
dgbcon
subroutine dgbcon(norm, n, kl, ku, ab, ldab, ipiv, anorm, rcond, work, iwork, info)
DGBCON
Definition dgbcon.f:146
dla_gbrfsx_extended
subroutine dla_gbrfsx_extended(prec_type, trans_type, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, colequ, c, b, ldb, y, ldy, berr_out, n_norms, err_bnds_norm, err_bnds_comp, res, ayb, dy, y_tail, rcond, ithresh, rthresh, dz_ub, ignore_cwise, info)
DLA_GBRFSX_EXTENDED improves the computed solution to a system of linear equations for general banded...
Definition dla_gbrfsx_extended.f:411

◆ dgbtf2()

subroutine dgbtf2 ( integer m,
integer n,
integer kl,
integer ku,
double precision, dimension( ldab, * ) ab,
integer ldab,
integer, dimension( * ) ipiv,
integer info )

DGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algorithm.

Download DGBTF2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGBTF2 computes an LU factorization of a real m-by-n band matrix A
!> using partial pivoting with row interchanges.
!>
!> This is the unblocked version of the algorithm, calling Level 2 BLAS.
!>
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]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>          On entry, the matrix A in band storage, in rows KL+1 to
!>          2*KL+KU+1; rows 1 to KL of the array need not be set.
!>          The j-th column of A is stored in the j-th column of the
!>          array AB as follows:
!>          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
!>
!>          On exit, details of the factorization: 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 below for further details.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
!>
[out]IPIV
!>          IPIV is INTEGER array, dimension (min(M,N))
!>          The pivot indices; for 1 <= i <= min(M,N), row i of the
!>          matrix was interchanged with row IPIV(i).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>          > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
!>               has been completed, but the factor U is exactly
!>               singular, and division by zero will occur if it is used
!>               to solve a system of equations.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The band storage scheme is illustrated by the following example, when
!>  M = N = 6, KL = 2, KU = 1:
!>
!>  On entry:                       On exit:
!>
!>      *    *    *    +    +    +       *    *    *   u14  u25  u36
!>      *    *    +    +    +    +       *    *   u13  u24  u35  u46
!>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
!>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
!>     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
!>     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
!>
!>  Array elements marked * are not used by the routine; elements marked
!>  + need not be set on entry, but are required by the routine to store
!>  elements of U, because of fill-in resulting from the row
!>  interchanges.
!>

Definition at line 144 of file dgbtf2.f.

145*
146* -- LAPACK computational routine --
147* -- LAPACK is a software package provided by Univ. of Tennessee, --
148* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
149*
150* .. Scalar Arguments ..
151 INTEGER INFO, KL, KU, LDAB, M, N
152* ..
153* .. Array Arguments ..
154 INTEGER IPIV( * )
155 DOUBLE PRECISION AB( LDAB, * )
156* ..
157*
158* =====================================================================
159*
160* .. Parameters ..
161 DOUBLE PRECISION ONE, ZERO
162 parameter( one = 1.0d+0, zero = 0.0d+0 )
163* ..
164* .. Local Scalars ..
165 INTEGER I, J, JP, JU, KM, KV
166* ..
167* .. External Functions ..
168 INTEGER IDAMAX
169 EXTERNAL idamax
170* ..
171* .. External Subroutines ..
172 EXTERNAL dger, dscal, dswap, xerbla
173* ..
174* .. Intrinsic Functions ..
175 INTRINSIC max, min
176* ..
177* .. Executable Statements ..
178*
179* KV is the number of superdiagonals in the factor U, allowing for
180* fill-in.
181*
182 kv = ku + kl
183*
184* Test the input parameters.
185*
186 info = 0
187 IF( m.LT.0 ) THEN
188 info = -1
189 ELSE IF( n.LT.0 ) THEN
190 info = -2
191 ELSE IF( kl.LT.0 ) THEN
192 info = -3
193 ELSE IF( ku.LT.0 ) THEN
194 info = -4
195 ELSE IF( ldab.LT.kl+kv+1 ) THEN
196 info = -6
197 END IF
198 IF( info.NE.0 ) THEN
199 CALL xerbla( 'DGBTF2', -info )
200 RETURN
201 END IF
202*
203* Quick return if possible
204*
205 IF( m.EQ.0 .OR. n.EQ.0 )
206 $ RETURN
207*
208* Gaussian elimination with partial pivoting
209*
210* Set fill-in elements in columns KU+2 to KV to zero.
211*
212 DO 20 j = ku + 2, min( kv, n )
213 DO 10 i = kv - j + 2, kl
214 ab( i, j ) = zero
215 10 CONTINUE
216 20 CONTINUE
217*
218* JU is the index of the last column affected by the current stage
219* of the factorization.
220*
221 ju = 1
222*
223 DO 40 j = 1, min( m, n )
224*
225* Set fill-in elements in column J+KV to zero.
226*
227 IF( j+kv.LE.n ) THEN
228 DO 30 i = 1, kl
229 ab( i, j+kv ) = zero
230 30 CONTINUE
231 END IF
232*
233* Find pivot and test for singularity. KM is the number of
234* subdiagonal elements in the current column.
235*
236 km = min( kl, m-j )
237 jp = idamax( km+1, ab( kv+1, j ), 1 )
238 ipiv( j ) = jp + j - 1
239 IF( ab( kv+jp, j ).NE.zero ) THEN
240 ju = max( ju, min( j+ku+jp-1, n ) )
241*
242* Apply interchange to columns J to JU.
243*
244 IF( jp.NE.1 )
245 $ CALL dswap( ju-j+1, ab( kv+jp, j ), ldab-1,
246 $ ab( kv+1, j ), ldab-1 )
247*
248 IF( km.GT.0 ) THEN
249*
250* Compute multipliers.
251*
252 CALL dscal( km, one / ab( kv+1, j ), ab( kv+2, j ), 1 )
253*
254* Update trailing submatrix within the band.
255*
256 IF( ju.GT.j )
257 $ CALL dger( km, ju-j, -one, ab( kv+2, j ), 1,
258 $ ab( kv, j+1 ), ldab-1, ab( kv+1, j+1 ),
259 $ ldab-1 )
260 END IF
261 ELSE
262*
263* If pivot is zero, set INFO to the index of the pivot
264* unless a zero pivot has already been found.
265*
266 IF( info.EQ.0 )
267 $ info = j
268 END IF
269 40 CONTINUE
270 RETURN
271*
272* End of DGBTF2
273*
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
dger
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130

◆ dgbtrf()

subroutine dgbtrf ( integer m,
integer n,
integer kl,
integer ku,
double precision, dimension( ldab, * ) ab,
integer ldab,
integer, dimension( * ) ipiv,
integer info )

DGBTRF

Download DGBTRF + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGBTRF computes an LU factorization of a real m-by-n band matrix A
!> using partial pivoting with row interchanges.
!>
!> This is the blocked version of the algorithm, calling Level 3 BLAS.
!>
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]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>          On entry, the matrix A in band storage, in rows KL+1 to
!>          2*KL+KU+1; rows 1 to KL of the array need not be set.
!>          The j-th column of A is stored in the j-th column of the
!>          array AB as follows:
!>          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
!>
!>          On exit, details of the factorization: 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 below for further details.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
!>
[out]IPIV
!>          IPIV is INTEGER array, dimension (min(M,N))
!>          The pivot indices; for 1 <= i <= min(M,N), row i of the
!>          matrix was interchanged with row IPIV(i).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>          > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
!>               has been completed, but the factor U is exactly
!>               singular, and division by zero will occur if it is used
!>               to solve a system of equations.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The band storage scheme is illustrated by the following example, when
!>  M = N = 6, KL = 2, KU = 1:
!>
!>  On entry:                       On exit:
!>
!>      *    *    *    +    +    +       *    *    *   u14  u25  u36
!>      *    *    +    +    +    +       *    *   u13  u24  u35  u46
!>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
!>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
!>     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
!>     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
!>
!>  Array elements marked * are not used by the routine; elements marked
!>  + need not be set on entry, but are required by the routine to store
!>  elements of U because of fill-in resulting from the row interchanges.
!>

Definition at line 143 of file dgbtrf.f.

144*
145* -- LAPACK computational routine --
146* -- LAPACK is a software package provided by Univ. of Tennessee, --
147* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148*
149* .. Scalar Arguments ..
150 INTEGER INFO, KL, KU, LDAB, M, N
151* ..
152* .. Array Arguments ..
153 INTEGER IPIV( * )
154 DOUBLE PRECISION AB( LDAB, * )
155* ..
156*
157* =====================================================================
158*
159* .. Parameters ..
160 DOUBLE PRECISION ONE, ZERO
161 parameter( one = 1.0d+0, zero = 0.0d+0 )
162 INTEGER NBMAX, LDWORK
163 parameter( nbmax = 64, ldwork = nbmax+1 )
164* ..
165* .. Local Scalars ..
166 INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP,
167 $ JU, K2, KM, KV, NB, NW
168 DOUBLE PRECISION TEMP
169* ..
170* .. Local Arrays ..
171 DOUBLE PRECISION WORK13( LDWORK, NBMAX ),
172 $ WORK31( LDWORK, NBMAX )
173* ..
174* .. External Functions ..
175 INTEGER IDAMAX, ILAENV
176 EXTERNAL idamax, ilaenv
177* ..
178* .. External Subroutines ..
179 EXTERNAL dcopy, dgbtf2, dgemm, dger, dlaswp, dscal,
180 $ dswap, dtrsm, xerbla
181* ..
182* .. Intrinsic Functions ..
183 INTRINSIC max, min
184* ..
185* .. Executable Statements ..
186*
187* KV is the number of superdiagonals in the factor U, allowing for
188* fill-in
189*
190 kv = ku + kl
191*
192* Test the input parameters.
193*
194 info = 0
195 IF( m.LT.0 ) THEN
196 info = -1
197 ELSE IF( n.LT.0 ) THEN
198 info = -2
199 ELSE IF( kl.LT.0 ) THEN
200 info = -3
201 ELSE IF( ku.LT.0 ) THEN
202 info = -4
203 ELSE IF( ldab.LT.kl+kv+1 ) THEN
204 info = -6
205 END IF
206 IF( info.NE.0 ) THEN
207 CALL xerbla( 'DGBTRF', -info )
208 RETURN
209 END IF
210*
211* Quick return if possible
212*
213 IF( m.EQ.0 .OR. n.EQ.0 )
214 $ RETURN
215*
216* Determine the block size for this environment
217*
218 nb = ilaenv( 1, 'DGBTRF', ' ', m, n, kl, ku )
219*
220* The block size must not exceed the limit set by the size of the
221* local arrays WORK13 and WORK31.
222*
223 nb = min( nb, nbmax )
224*
225 IF( nb.LE.1 .OR. nb.GT.kl ) THEN
226*
227* Use unblocked code
228*
229 CALL dgbtf2( m, n, kl, ku, ab, ldab, ipiv, info )
230 ELSE
231*
232* Use blocked code
233*
234* Zero the superdiagonal elements of the work array WORK13
235*
236 DO 20 j = 1, nb
237 DO 10 i = 1, j - 1
238 work13( i, j ) = zero
239 10 CONTINUE
240 20 CONTINUE
241*
242* Zero the subdiagonal elements of the work array WORK31
243*
244 DO 40 j = 1, nb
245 DO 30 i = j + 1, nb
246 work31( i, j ) = zero
247 30 CONTINUE
248 40 CONTINUE
249*
250* Gaussian elimination with partial pivoting
251*
252* Set fill-in elements in columns KU+2 to KV to zero
253*
254 DO 60 j = ku + 2, min( kv, n )
255 DO 50 i = kv - j + 2, kl
256 ab( i, j ) = zero
257 50 CONTINUE
258 60 CONTINUE
259*
260* JU is the index of the last column affected by the current
261* stage of the factorization
262*
263 ju = 1
264*
265 DO 180 j = 1, min( m, n ), nb
266 jb = min( nb, min( m, n )-j+1 )
267*
268* The active part of the matrix is partitioned
269*
270* A11 A12 A13
271* A21 A22 A23
272* A31 A32 A33
273*
274* Here A11, A21 and A31 denote the current block of JB columns
275* which is about to be factorized. The number of rows in the
276* partitioning are JB, I2, I3 respectively, and the numbers
277* of columns are JB, J2, J3. The superdiagonal elements of A13
278* and the subdiagonal elements of A31 lie outside the band.
279*
280 i2 = min( kl-jb, m-j-jb+1 )
281 i3 = min( jb, m-j-kl+1 )
282*
283* J2 and J3 are computed after JU has been updated.
284*
285* Factorize the current block of JB columns
286*
287 DO 80 jj = j, j + jb - 1
288*
289* Set fill-in elements in column JJ+KV to zero
290*
291 IF( jj+kv.LE.n ) THEN
292 DO 70 i = 1, kl
293 ab( i, jj+kv ) = zero
294 70 CONTINUE
295 END IF
296*
297* Find pivot and test for singularity. KM is the number of
298* subdiagonal elements in the current column.
299*
300 km = min( kl, m-jj )
301 jp = idamax( km+1, ab( kv+1, jj ), 1 )
302 ipiv( jj ) = jp + jj - j
303 IF( ab( kv+jp, jj ).NE.zero ) THEN
304 ju = max( ju, min( jj+ku+jp-1, n ) )
305 IF( jp.NE.1 ) THEN
306*
307* Apply interchange to columns J to J+JB-1
308*
309 IF( jp+jj-1.LT.j+kl ) THEN
310*
311 CALL dswap( jb, ab( kv+1+jj-j, j ), ldab-1,
312 $ ab( kv+jp+jj-j, j ), ldab-1 )
313 ELSE
314*
315* The interchange affects columns J to JJ-1 of A31
316* which are stored in the work array WORK31
317*
318 CALL dswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
319 $ work31( jp+jj-j-kl, 1 ), ldwork )
320 CALL dswap( j+jb-jj, ab( kv+1, jj ), ldab-1,
321 $ ab( kv+jp, jj ), ldab-1 )
322 END IF
323 END IF
324*
325* Compute multipliers
326*
327 CALL dscal( km, one / ab( kv+1, jj ), ab( kv+2, jj ),
328 $ 1 )
329*
330* Update trailing submatrix within the band and within
331* the current block. JM is the index of the last column
332* which needs to be updated.
333*
334 jm = min( ju, j+jb-1 )
335 IF( jm.GT.jj )
336 $ CALL dger( km, jm-jj, -one, ab( kv+2, jj ), 1,
337 $ ab( kv, jj+1 ), ldab-1,
338 $ ab( kv+1, jj+1 ), ldab-1 )
339 ELSE
340*
341* If pivot is zero, set INFO to the index of the pivot
342* unless a zero pivot has already been found.
343*
344 IF( info.EQ.0 )
345 $ info = jj
346 END IF
347*
348* Copy current column of A31 into the work array WORK31
349*
350 nw = min( jj-j+1, i3 )
351 IF( nw.GT.0 )
352 $ CALL dcopy( nw, ab( kv+kl+1-jj+j, jj ), 1,
353 $ work31( 1, jj-j+1 ), 1 )
354 80 CONTINUE
355 IF( j+jb.LE.n ) THEN
356*
357* Apply the row interchanges to the other blocks.
358*
359 j2 = min( ju-j+1, kv ) - jb
360 j3 = max( 0, ju-j-kv+1 )
361*
362* Use DLASWP to apply the row interchanges to A12, A22, and
363* A32.
364*
365 CALL dlaswp( j2, ab( kv+1-jb, j+jb ), ldab-1, 1, jb,
366 $ ipiv( j ), 1 )
367*
368* Adjust the pivot indices.
369*
370 DO 90 i = j, j + jb - 1
371 ipiv( i ) = ipiv( i ) + j - 1
372 90 CONTINUE
373*
374* Apply the row interchanges to A13, A23, and A33
375* columnwise.
376*
377 k2 = j - 1 + jb + j2
378 DO 110 i = 1, j3
379 jj = k2 + i
380 DO 100 ii = j + i - 1, j + jb - 1
381 ip = ipiv( ii )
382 IF( ip.NE.ii ) THEN
383 temp = ab( kv+1+ii-jj, jj )
384 ab( kv+1+ii-jj, jj ) = ab( kv+1+ip-jj, jj )
385 ab( kv+1+ip-jj, jj ) = temp
386 END IF
387 100 CONTINUE
388 110 CONTINUE
389*
390* Update the relevant part of the trailing submatrix
391*
392 IF( j2.GT.0 ) THEN
393*
394* Update A12
395*
396 CALL dtrsm( 'Left', 'Lower', 'No transpose', 'Unit',
397 $ jb, j2, one, ab( kv+1, j ), ldab-1,
398 $ ab( kv+1-jb, j+jb ), ldab-1 )
399*
400 IF( i2.GT.0 ) THEN
401*
402* Update A22
403*
404 CALL dgemm( 'No transpose', 'No transpose', i2, j2,
405 $ jb, -one, ab( kv+1+jb, j ), ldab-1,
406 $ ab( kv+1-jb, j+jb ), ldab-1, one,
407 $ ab( kv+1, j+jb ), ldab-1 )
408 END IF
409*
410 IF( i3.GT.0 ) THEN
411*
412* Update A32
413*
414 CALL dgemm( 'No transpose', 'No transpose', i3, j2,
415 $ jb, -one, work31, ldwork,
416 $ ab( kv+1-jb, j+jb ), ldab-1, one,
417 $ ab( kv+kl+1-jb, j+jb ), ldab-1 )
418 END IF
419 END IF
420*
421 IF( j3.GT.0 ) THEN
422*
423* Copy the lower triangle of A13 into the work array
424* WORK13
425*
426 DO 130 jj = 1, j3
427 DO 120 ii = jj, jb
428 work13( ii, jj ) = ab( ii-jj+1, jj+j+kv-1 )
429 120 CONTINUE
430 130 CONTINUE
431*
432* Update A13 in the work array
433*
434 CALL dtrsm( 'Left', 'Lower', 'No transpose', 'Unit',
435 $ jb, j3, one, ab( kv+1, j ), ldab-1,
436 $ work13, ldwork )
437*
438 IF( i2.GT.0 ) THEN
439*
440* Update A23
441*
442 CALL dgemm( 'No transpose', 'No transpose', i2, j3,
443 $ jb, -one, ab( kv+1+jb, j ), ldab-1,
444 $ work13, ldwork, one, ab( 1+jb, j+kv ),
445 $ ldab-1 )
446 END IF
447*
448 IF( i3.GT.0 ) THEN
449*
450* Update A33
451*
452 CALL dgemm( 'No transpose', 'No transpose', i3, j3,
453 $ jb, -one, work31, ldwork, work13,
454 $ ldwork, one, ab( 1+kl, j+kv ), ldab-1 )
455 END IF
456*
457* Copy the lower triangle of A13 back into place
458*
459 DO 150 jj = 1, j3
460 DO 140 ii = jj, jb
461 ab( ii-jj+1, jj+j+kv-1 ) = work13( ii, jj )
462 140 CONTINUE
463 150 CONTINUE
464 END IF
465 ELSE
466*
467* Adjust the pivot indices.
468*
469 DO 160 i = j, j + jb - 1
470 ipiv( i ) = ipiv( i ) + j - 1
471 160 CONTINUE
472 END IF
473*
474* Partially undo the interchanges in the current block to
475* restore the upper triangular form of A31 and copy the upper
476* triangle of A31 back into place
477*
478 DO 170 jj = j + jb - 1, j, -1
479 jp = ipiv( jj ) - jj + 1
480 IF( jp.NE.1 ) THEN
481*
482* Apply interchange to columns J to JJ-1
483*
484 IF( jp+jj-1.LT.j+kl ) THEN
485*
486* The interchange does not affect A31
487*
488 CALL dswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
489 $ ab( kv+jp+jj-j, j ), ldab-1 )
490 ELSE
491*
492* The interchange does affect A31
493*
494 CALL dswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
495 $ work31( jp+jj-j-kl, 1 ), ldwork )
496 END IF
497 END IF
498*
499* Copy the current column of A31 back into place
500*
501 nw = min( i3, jj-j+1 )
502 IF( nw.GT.0 )
503 $ CALL dcopy( nw, work31( 1, jj-j+1 ), 1,
504 $ ab( kv+kl+1-jj+j, jj ), 1 )
505 170 CONTINUE
506 180 CONTINUE
507 END IF
508*
509 RETURN
510*
511* End of DGBTRF
512*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
dgbtf2
subroutine dgbtf2(m, n, kl, ku, ab, ldab, ipiv, info)
DGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algo...
Definition dgbtf2.f:145
dlaswp
subroutine dlaswp(n, a, lda, k1, k2, ipiv, incx)
DLASWP performs a series of row interchanges on a general rectangular matrix.
Definition dlaswp.f:115
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187
dtrsm
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181

◆ dgbtrs()

subroutine dgbtrs ( character trans,
integer n,
integer kl,
integer ku,
integer nrhs,
double precision, dimension( ldab, * ) ab,
integer ldab,
integer, dimension( * ) ipiv,
double precision, dimension( ldb, * ) b,
integer ldb,
integer info )

DGBTRS

Download DGBTRS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGBTRS solves a system of linear equations
!>    A * X = B  or  A**T * X = B
!> with a general band matrix A using the LU factorization computed
!> by DGBTRF.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>          Specifies the form of the system of equations.
!>          = 'N':  A * X = B  (No transpose)
!>          = 'T':  A**T* X = B  (Transpose)
!>          = 'C':  A**T* X = B  (Conjugate transpose = Transpose)
!>
[in]N
!>          N is INTEGER
!>          The order 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]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>          Details of the LU factorization of the band matrix A, 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.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          The pivot indices; for 1 <= i <= N, row i of the matrix was
!>          interchanged with row IPIV(i).
!>
[in,out]B
!>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
!>          On entry, the right hand side matrix B.
!>          On exit, the solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 136 of file dgbtrs.f.

138*
139* -- LAPACK computational routine --
140* -- LAPACK is a software package provided by Univ. of Tennessee, --
141* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142*
143* .. Scalar Arguments ..
144 CHARACTER TRANS
145 INTEGER INFO, KL, KU, LDAB, LDB, N, NRHS
146* ..
147* .. Array Arguments ..
148 INTEGER IPIV( * )
149 DOUBLE PRECISION AB( LDAB, * ), B( LDB, * )
150* ..
151*
152* =====================================================================
153*
154* .. Parameters ..
155 DOUBLE PRECISION ONE
156 parameter( one = 1.0d+0 )
157* ..
158* .. Local Scalars ..
159 LOGICAL LNOTI, NOTRAN
160 INTEGER I, J, KD, L, LM
161* ..
162* .. External Functions ..
163 LOGICAL LSAME
164 EXTERNAL lsame
165* ..
166* .. External Subroutines ..
167 EXTERNAL dgemv, dger, dswap, dtbsv, xerbla
168* ..
169* .. Intrinsic Functions ..
170 INTRINSIC max, min
171* ..
172* .. Executable Statements ..
173*
174* Test the input parameters.
175*
176 info = 0
177 notran = lsame( trans, 'N' )
178 IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
179 $ lsame( trans, 'C' ) ) THEN
180 info = -1
181 ELSE IF( n.LT.0 ) THEN
182 info = -2
183 ELSE IF( kl.LT.0 ) THEN
184 info = -3
185 ELSE IF( ku.LT.0 ) THEN
186 info = -4
187 ELSE IF( nrhs.LT.0 ) THEN
188 info = -5
189 ELSE IF( ldab.LT.( 2*kl+ku+1 ) ) THEN
190 info = -7
191 ELSE IF( ldb.LT.max( 1, n ) ) THEN
192 info = -10
193 END IF
194 IF( info.NE.0 ) THEN
195 CALL xerbla( 'DGBTRS', -info )
196 RETURN
197 END IF
198*
199* Quick return if possible
200*
201 IF( n.EQ.0 .OR. nrhs.EQ.0 )
202 $ RETURN
203*
204 kd = ku + kl + 1
205 lnoti = kl.GT.0
206*
207 IF( notran ) THEN
208*
209* Solve A*X = B.
210*
211* Solve L*X = B, overwriting B with X.
212*
213* L is represented as a product of permutations and unit lower
214* triangular matrices L = P(1) * L(1) * ... * P(n-1) * L(n-1),
215* where each transformation L(i) is a rank-one modification of
216* the identity matrix.
217*
218 IF( lnoti ) THEN
219 DO 10 j = 1, n - 1
220 lm = min( kl, n-j )
221 l = ipiv( j )
222 IF( l.NE.j )
223 $ CALL dswap( nrhs, b( l, 1 ), ldb, b( j, 1 ), ldb )
224 CALL dger( lm, nrhs, -one, ab( kd+1, j ), 1, b( j, 1 ),
225 $ ldb, b( j+1, 1 ), ldb )
226 10 CONTINUE
227 END IF
228*
229 DO 20 i = 1, nrhs
230*
231* Solve U*X = B, overwriting B with X.
232*
233 CALL dtbsv( 'Upper', 'No transpose', 'Non-unit', n, kl+ku,
234 $ ab, ldab, b( 1, i ), 1 )
235 20 CONTINUE
236*
237 ELSE
238*
239* Solve A**T*X = B.
240*
241 DO 30 i = 1, nrhs
242*
243* Solve U**T*X = B, overwriting B with X.
244*
245 CALL dtbsv( 'Upper', 'Transpose', 'Non-unit', n, kl+ku, ab,
246 $ ldab, b( 1, i ), 1 )
247 30 CONTINUE
248*
249* Solve L**T*X = B, overwriting B with X.
250*
251 IF( lnoti ) THEN
252 DO 40 j = n - 1, 1, -1
253 lm = min( kl, n-j )
254 CALL dgemv( 'Transpose', lm, nrhs, -one, b( j+1, 1 ),
255 $ ldb, ab( kd+1, j ), 1, one, b( j, 1 ), ldb )
256 l = ipiv( j )
257 IF( l.NE.j )
258 $ CALL dswap( nrhs, b( l, 1 ), ldb, b( j, 1 ), ldb )
259 40 CONTINUE
260 END IF
261 END IF
262 RETURN
263*
264* End of DGBTRS
265*
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:156
dtbsv
subroutine dtbsv(uplo, trans, diag, n, k, a, lda, x, incx)
DTBSV
Definition dtbsv.f:189

◆ dggbak()

subroutine dggbak ( character job,
character side,
integer n,
integer ilo,
integer ihi,
double precision, dimension( * ) lscale,
double precision, dimension( * ) rscale,
integer m,
double precision, dimension( ldv, * ) v,
integer ldv,
integer info )

DGGBAK

Download DGGBAK + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGGBAK forms the right or left eigenvectors of a real generalized
!> eigenvalue problem A*x = lambda*B*x, by backward transformation on
!> the computed eigenvectors of the balanced pair of matrices output by
!> DGGBAL.
!>
Parameters
[in]JOB
!>          JOB is CHARACTER*1
!>          Specifies the type of backward transformation required:
!>          = 'N':  do nothing, return immediately;
!>          = 'P':  do backward transformation for permutation only;
!>          = 'S':  do backward transformation for scaling only;
!>          = 'B':  do backward transformations for both permutation and
!>                  scaling.
!>          JOB must be the same as the argument JOB supplied to DGGBAL.
!>
[in]SIDE
!>          SIDE is CHARACTER*1
!>          = 'R':  V contains right eigenvectors;
!>          = 'L':  V contains left eigenvectors.
!>
[in]N
!>          N is INTEGER
!>          The number of rows of the matrix V.  N >= 0.
!>
[in]ILO
!>          ILO is INTEGER
!>
[in]IHI
!>          IHI is INTEGER
!>          The integers ILO and IHI determined by DGGBAL.
!>          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
!>
[in]LSCALE
!>          LSCALE is DOUBLE PRECISION array, dimension (N)
!>          Details of the permutations and/or scaling factors applied
!>          to the left side of A and B, as returned by DGGBAL.
!>
[in]RSCALE
!>          RSCALE is DOUBLE PRECISION array, dimension (N)
!>          Details of the permutations and/or scaling factors applied
!>          to the right side of A and B, as returned by DGGBAL.
!>
[in]M
!>          M is INTEGER
!>          The number of columns of the matrix V.  M >= 0.
!>
[in,out]V
!>          V is DOUBLE PRECISION array, dimension (LDV,M)
!>          On entry, the matrix of right or left eigenvectors to be
!>          transformed, as returned by DTGEVC.
!>          On exit, V is overwritten by the transformed eigenvectors.
!>
[in]LDV
!>          LDV is INTEGER
!>          The leading dimension of the matrix V. LDV >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  See R.C. Ward, Balancing the generalized eigenvalue problem,
!>                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
!>

Definition at line 145 of file dggbak.f.

147*
148* -- LAPACK computational routine --
149* -- LAPACK is a software package provided by Univ. of Tennessee, --
150* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151*
152* .. Scalar Arguments ..
153 CHARACTER JOB, SIDE
154 INTEGER IHI, ILO, INFO, LDV, M, N
155* ..
156* .. Array Arguments ..
157 DOUBLE PRECISION LSCALE( * ), RSCALE( * ), V( LDV, * )
158* ..
159*
160* =====================================================================
161*
162* .. Local Scalars ..
163 LOGICAL LEFTV, RIGHTV
164 INTEGER I, K
165* ..
166* .. External Functions ..
167 LOGICAL LSAME
168 EXTERNAL lsame
169* ..
170* .. External Subroutines ..
171 EXTERNAL dscal, dswap, xerbla
172* ..
173* .. Intrinsic Functions ..
174 INTRINSIC max, int
175* ..
176* .. Executable Statements ..
177*
178* Test the input parameters
179*
180 rightv = lsame( side, 'R' )
181 leftv = lsame( side, 'L' )
182*
183 info = 0
184 IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
185 $ .NOT.lsame( job, 'S' ) .AND. .NOT.lsame( job, 'B' ) ) THEN
186 info = -1
187 ELSE IF( .NOT.rightv .AND. .NOT.leftv ) THEN
188 info = -2
189 ELSE IF( n.LT.0 ) THEN
190 info = -3
191 ELSE IF( ilo.LT.1 ) THEN
192 info = -4
193 ELSE IF( n.EQ.0 .AND. ihi.EQ.0 .AND. ilo.NE.1 ) THEN
194 info = -4
195 ELSE IF( n.GT.0 .AND. ( ihi.LT.ilo .OR. ihi.GT.max( 1, n ) ) )
196 $ THEN
197 info = -5
198 ELSE IF( n.EQ.0 .AND. ilo.EQ.1 .AND. ihi.NE.0 ) THEN
199 info = -5
200 ELSE IF( m.LT.0 ) THEN
201 info = -8
202 ELSE IF( ldv.LT.max( 1, n ) ) THEN
203 info = -10
204 END IF
205 IF( info.NE.0 ) THEN
206 CALL xerbla( 'DGGBAK', -info )
207 RETURN
208 END IF
209*
210* Quick return if possible
211*
212 IF( n.EQ.0 )
213 $ RETURN
214 IF( m.EQ.0 )
215 $ RETURN
216 IF( lsame( job, 'N' ) )
217 $ RETURN
218*
219 IF( ilo.EQ.ihi )
220 $ GO TO 30
221*
222* Backward balance
223*
224 IF( lsame( job, 'S' ) .OR. lsame( job, 'B' ) ) THEN
225*
226* Backward transformation on right eigenvectors
227*
228 IF( rightv ) THEN
229 DO 10 i = ilo, ihi
230 CALL dscal( m, rscale( i ), v( i, 1 ), ldv )
231 10 CONTINUE
232 END IF
233*
234* Backward transformation on left eigenvectors
235*
236 IF( leftv ) THEN
237 DO 20 i = ilo, ihi
238 CALL dscal( m, lscale( i ), v( i, 1 ), ldv )
239 20 CONTINUE
240 END IF
241 END IF
242*
243* Backward permutation
244*
245 30 CONTINUE
246 IF( lsame( job, 'P' ) .OR. lsame( job, 'B' ) ) THEN
247*
248* Backward permutation on right eigenvectors
249*
250 IF( rightv ) THEN
251 IF( ilo.EQ.1 )
252 $ GO TO 50
253*
254 DO 40 i = ilo - 1, 1, -1
255 k = int(rscale( i ))
256 IF( k.EQ.i )
257 $ GO TO 40
258 CALL dswap( m, v( i, 1 ), ldv, v( k, 1 ), ldv )
259 40 CONTINUE
260*
261 50 CONTINUE
262 IF( ihi.EQ.n )
263 $ GO TO 70
264 DO 60 i = ihi + 1, n
265 k = int(rscale( i ))
266 IF( k.EQ.i )
267 $ GO TO 60
268 CALL dswap( m, v( i, 1 ), ldv, v( k, 1 ), ldv )
269 60 CONTINUE
270 END IF
271*
272* Backward permutation on left eigenvectors
273*
274 70 CONTINUE
275 IF( leftv ) THEN
276 IF( ilo.EQ.1 )
277 $ GO TO 90
278 DO 80 i = ilo - 1, 1, -1
279 k = int(lscale( i ))
280 IF( k.EQ.i )
281 $ GO TO 80
282 CALL dswap( m, v( i, 1 ), ldv, v( k, 1 ), ldv )
283 80 CONTINUE
284*
285 90 CONTINUE
286 IF( ihi.EQ.n )
287 $ GO TO 110
288 DO 100 i = ihi + 1, n
289 k = int(lscale( i ))
290 IF( k.EQ.i )
291 $ GO TO 100
292 CALL dswap( m, v( i, 1 ), ldv, v( k, 1 ), ldv )
293 100 CONTINUE
294 END IF
295 END IF
296*
297 110 CONTINUE
298*
299 RETURN
300*
301* End of DGGBAK
302*

◆ dggbal()

subroutine dggbal ( character job,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
integer ilo,
integer ihi,
double precision, dimension( * ) lscale,
double precision, dimension( * ) rscale,
double precision, dimension( * ) work,
integer info )

DGGBAL

Download DGGBAL + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGGBAL balances a pair of general real matrices (A,B).  This
!> involves, first, permuting A and B by similarity transformations to
!> isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
!> elements on the diagonal; and second, applying a diagonal similarity
!> transformation to rows and columns ILO to IHI to make the rows
!> and columns as close in norm as possible. Both steps are optional.
!>
!> Balancing may reduce the 1-norm of the matrices, and improve the
!> accuracy of the computed eigenvalues and/or eigenvectors in the
!> generalized eigenvalue problem A*x = lambda*B*x.
!>
Parameters
[in]JOB
!>          JOB is CHARACTER*1
!>          Specifies the operations to be performed on A and B:
!>          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
!>                  and RSCALE(I) = 1.0 for i = 1,...,N.
!>          = 'P':  permute only;
!>          = 'S':  scale only;
!>          = 'B':  both permute and scale.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  N >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>          On entry, the input matrix A.
!>          On exit,  A is overwritten by the balanced matrix.
!>          If JOB = 'N', A is not referenced.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,N).
!>
[in,out]B
!>          B is DOUBLE PRECISION array, dimension (LDB,N)
!>          On entry, the input matrix B.
!>          On exit,  B is overwritten by the balanced matrix.
!>          If JOB = 'N', B is not referenced.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,N).
!>
[out]ILO
!>          ILO is INTEGER
!>
[out]IHI
!>          IHI is INTEGER
!>          ILO and IHI are set to integers such that on exit
!>          A(i,j) = 0 and B(i,j) = 0 if i > j and
!>          j = 1,...,ILO-1 or i = IHI+1,...,N.
!>          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
!>
[out]LSCALE
!>          LSCALE is DOUBLE PRECISION array, dimension (N)
!>          Details of the permutations and scaling factors applied
!>          to the left side of A and B.  If P(j) is the index of the
!>          row interchanged with row j, and D(j)
!>          is the scaling factor applied to row j, then
!>            LSCALE(j) = P(j)    for J = 1,...,ILO-1
!>                      = D(j)    for J = ILO,...,IHI
!>                      = P(j)    for J = IHI+1,...,N.
!>          The order in which the interchanges are made is N to IHI+1,
!>          then 1 to ILO-1.
!>
[out]RSCALE
!>          RSCALE is DOUBLE PRECISION array, dimension (N)
!>          Details of the permutations and scaling factors applied
!>          to the right side of A and B.  If P(j) is the index of the
!>          column interchanged with column j, and D(j)
!>          is the scaling factor applied to column j, then
!>            LSCALE(j) = P(j)    for J = 1,...,ILO-1
!>                      = D(j)    for J = ILO,...,IHI
!>                      = P(j)    for J = IHI+1,...,N.
!>          The order in which the interchanges are made is N to IHI+1,
!>          then 1 to ILO-1.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (lwork)
!>          lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
!>          at least 1 when JOB = 'N' or 'P'.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  See R.C. WARD, Balancing the generalized eigenvalue problem,
!>                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
!>

Definition at line 175 of file dggbal.f.

177*
178* -- LAPACK computational routine --
179* -- LAPACK is a software package provided by Univ. of Tennessee, --
180* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
181*
182* .. Scalar Arguments ..
183 CHARACTER JOB
184 INTEGER IHI, ILO, INFO, LDA, LDB, N
185* ..
186* .. Array Arguments ..
187 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ),
188 $ RSCALE( * ), WORK( * )
189* ..
190*
191* =====================================================================
192*
193* .. Parameters ..
194 DOUBLE PRECISION ZERO, HALF, ONE
195 parameter( zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0 )
196 DOUBLE PRECISION THREE, SCLFAC
197 parameter( three = 3.0d+0, sclfac = 1.0d+1 )
198* ..
199* .. Local Scalars ..
200 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
201 $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
202 $ M, NR, NRP2
203 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
204 $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
205 $ SFMIN, SUM, T, TA, TB, TC
206* ..
207* .. External Functions ..
208 LOGICAL LSAME
209 INTEGER IDAMAX
210 DOUBLE PRECISION DDOT, DLAMCH
211 EXTERNAL lsame, idamax, ddot, dlamch
212* ..
213* .. External Subroutines ..
214 EXTERNAL daxpy, dscal, dswap, xerbla
215* ..
216* .. Intrinsic Functions ..
217 INTRINSIC abs, dble, int, log10, max, min, sign
218* ..
219* .. Executable Statements ..
220*
221* Test the input parameters
222*
223 info = 0
224 IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
225 $ .NOT.lsame( job, 'S' ) .AND. .NOT.lsame( job, 'B' ) ) THEN
226 info = -1
227 ELSE IF( n.LT.0 ) THEN
228 info = -2
229 ELSE IF( lda.LT.max( 1, n ) ) THEN
230 info = -4
231 ELSE IF( ldb.LT.max( 1, n ) ) THEN
232 info = -6
233 END IF
234 IF( info.NE.0 ) THEN
235 CALL xerbla( 'DGGBAL', -info )
236 RETURN
237 END IF
238*
239* Quick return if possible
240*
241 IF( n.EQ.0 ) THEN
242 ilo = 1
243 ihi = n
244 RETURN
245 END IF
246*
247 IF( n.EQ.1 ) THEN
248 ilo = 1
249 ihi = n
250 lscale( 1 ) = one
251 rscale( 1 ) = one
252 RETURN
253 END IF
254*
255 IF( lsame( job, 'N' ) ) THEN
256 ilo = 1
257 ihi = n
258 DO 10 i = 1, n
259 lscale( i ) = one
260 rscale( i ) = one
261 10 CONTINUE
262 RETURN
263 END IF
264*
265 k = 1
266 l = n
267 IF( lsame( job, 'S' ) )
268 $ GO TO 190
269*
270 GO TO 30
271*
272* Permute the matrices A and B to isolate the eigenvalues.
273*
274* Find row with one nonzero in columns 1 through L
275*
276 20 CONTINUE
277 l = lm1
278 IF( l.NE.1 )
279 $ GO TO 30
280*
281 rscale( 1 ) = one
282 lscale( 1 ) = one
283 GO TO 190
284*
285 30 CONTINUE
286 lm1 = l - 1
287 DO 80 i = l, 1, -1
288 DO 40 j = 1, lm1
289 jp1 = j + 1
290 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
291 $ GO TO 50
292 40 CONTINUE
293 j = l
294 GO TO 70
295*
296 50 CONTINUE
297 DO 60 j = jp1, l
298 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
299 $ GO TO 80
300 60 CONTINUE
301 j = jp1 - 1
302*
303 70 CONTINUE
304 m = l
305 iflow = 1
306 GO TO 160
307 80 CONTINUE
308 GO TO 100
309*
310* Find column with one nonzero in rows K through N
311*
312 90 CONTINUE
313 k = k + 1
314*
315 100 CONTINUE
316 DO 150 j = k, l
317 DO 110 i = k, lm1
318 ip1 = i + 1
319 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
320 $ GO TO 120
321 110 CONTINUE
322 i = l
323 GO TO 140
324 120 CONTINUE
325 DO 130 i = ip1, l
326 IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
327 $ GO TO 150
328 130 CONTINUE
329 i = ip1 - 1
330 140 CONTINUE
331 m = k
332 iflow = 2
333 GO TO 160
334 150 CONTINUE
335 GO TO 190
336*
337* Permute rows M and I
338*
339 160 CONTINUE
340 lscale( m ) = i
341 IF( i.EQ.m )
342 $ GO TO 170
343 CALL dswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
344 CALL dswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
345*
346* Permute columns M and J
347*
348 170 CONTINUE
349 rscale( m ) = j
350 IF( j.EQ.m )
351 $ GO TO 180
352 CALL dswap( l, a( 1, j ), 1, a( 1, m ), 1 )
353 CALL dswap( l, b( 1, j ), 1, b( 1, m ), 1 )
354*
355 180 CONTINUE
356 GO TO ( 20, 90 )iflow
357*
358 190 CONTINUE
359 ilo = k
360 ihi = l
361*
362 IF( lsame( job, 'P' ) ) THEN
363 DO 195 i = ilo, ihi
364 lscale( i ) = one
365 rscale( i ) = one
366 195 CONTINUE
367 RETURN
368 END IF
369*
370 IF( ilo.EQ.ihi )
371 $ RETURN
372*
373* Balance the submatrix in rows ILO to IHI.
374*
375 nr = ihi - ilo + 1
376 DO 200 i = ilo, ihi
377 rscale( i ) = zero
378 lscale( i ) = zero
379*
380 work( i ) = zero
381 work( i+n ) = zero
382 work( i+2*n ) = zero
383 work( i+3*n ) = zero
384 work( i+4*n ) = zero
385 work( i+5*n ) = zero
386 200 CONTINUE
387*
388* Compute right side vector in resulting linear equations
389*
390 basl = log10( sclfac )
391 DO 240 i = ilo, ihi
392 DO 230 j = ilo, ihi
393 tb = b( i, j )
394 ta = a( i, j )
395 IF( ta.EQ.zero )
396 $ GO TO 210
397 ta = log10( abs( ta ) ) / basl
398 210 CONTINUE
399 IF( tb.EQ.zero )
400 $ GO TO 220
401 tb = log10( abs( tb ) ) / basl
402 220 CONTINUE
403 work( i+4*n ) = work( i+4*n ) - ta - tb
404 work( j+5*n ) = work( j+5*n ) - ta - tb
405 230 CONTINUE
406 240 CONTINUE
407*
408 coef = one / dble( 2*nr )
409 coef2 = coef*coef
410 coef5 = half*coef2
411 nrp2 = nr + 2
412 beta = zero
413 it = 1
414*
415* Start generalized conjugate gradient iteration
416*
417 250 CONTINUE
418*
419 gamma = ddot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
420 $ ddot( nr, work( ilo+5*n ), 1, work( ilo+5*n ), 1 )
421*
422 ew = zero
423 ewc = zero
424 DO 260 i = ilo, ihi
425 ew = ew + work( i+4*n )
426 ewc = ewc + work( i+5*n )
427 260 CONTINUE
428*
429 gamma = coef*gamma - coef2*( ew**2+ewc**2 ) - coef5*( ew-ewc )**2
430 IF( gamma.EQ.zero )
431 $ GO TO 350
432 IF( it.NE.1 )
433 $ beta = gamma / pgamma
434 t = coef5*( ewc-three*ew )
435 tc = coef5*( ew-three*ewc )
436*
437 CALL dscal( nr, beta, work( ilo ), 1 )
438 CALL dscal( nr, beta, work( ilo+n ), 1 )
439*
440 CALL daxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
441 CALL daxpy( nr, coef, work( ilo+5*n ), 1, work( ilo ), 1 )
442*
443 DO 270 i = ilo, ihi
444 work( i ) = work( i ) + tc
445 work( i+n ) = work( i+n ) + t
446 270 CONTINUE
447*
448* Apply matrix to vector
449*
450 DO 300 i = ilo, ihi
451 kount = 0
452 sum = zero
453 DO 290 j = ilo, ihi
454 IF( a( i, j ).EQ.zero )
455 $ GO TO 280
456 kount = kount + 1
457 sum = sum + work( j )
458 280 CONTINUE
459 IF( b( i, j ).EQ.zero )
460 $ GO TO 290
461 kount = kount + 1
462 sum = sum + work( j )
463 290 CONTINUE
464 work( i+2*n ) = dble( kount )*work( i+n ) + sum
465 300 CONTINUE
466*
467 DO 330 j = ilo, ihi
468 kount = 0
469 sum = zero
470 DO 320 i = ilo, ihi
471 IF( a( i, j ).EQ.zero )
472 $ GO TO 310
473 kount = kount + 1
474 sum = sum + work( i+n )
475 310 CONTINUE
476 IF( b( i, j ).EQ.zero )
477 $ GO TO 320
478 kount = kount + 1
479 sum = sum + work( i+n )
480 320 CONTINUE
481 work( j+3*n ) = dble( kount )*work( j ) + sum
482 330 CONTINUE
483*
484 sum = ddot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
485 $ ddot( nr, work( ilo ), 1, work( ilo+3*n ), 1 )
486 alpha = gamma / sum
487*
488* Determine correction to current iteration
489*
490 cmax = zero
491 DO 340 i = ilo, ihi
492 cor = alpha*work( i+n )
493 IF( abs( cor ).GT.cmax )
494 $ cmax = abs( cor )
495 lscale( i ) = lscale( i ) + cor
496 cor = alpha*work( i )
497 IF( abs( cor ).GT.cmax )
498 $ cmax = abs( cor )
499 rscale( i ) = rscale( i ) + cor
500 340 CONTINUE
501 IF( cmax.LT.half )
502 $ GO TO 350
503*
504 CALL daxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ), 1 )
505 CALL daxpy( nr, -alpha, work( ilo+3*n ), 1, work( ilo+5*n ), 1 )
506*
507 pgamma = gamma
508 it = it + 1
509 IF( it.LE.nrp2 )
510 $ GO TO 250
511*
512* End generalized conjugate gradient iteration
513*
514 350 CONTINUE
515 sfmin = dlamch( 'S' )
516 sfmax = one / sfmin
517 lsfmin = int( log10( sfmin ) / basl+one )
518 lsfmax = int( log10( sfmax ) / basl )
519 DO 360 i = ilo, ihi
520 irab = idamax( n-ilo+1, a( i, ilo ), lda )
521 rab = abs( a( i, irab+ilo-1 ) )
522 irab = idamax( n-ilo+1, b( i, ilo ), ldb )
523 rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
524 lrab = int( log10( rab+sfmin ) / basl+one )
525 ir = int(lscale( i ) + sign( half, lscale( i ) ))
526 ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
527 lscale( i ) = sclfac**ir
528 icab = idamax( ihi, a( 1, i ), 1 )
529 cab = abs( a( icab, i ) )
530 icab = idamax( ihi, b( 1, i ), 1 )
531 cab = max( cab, abs( b( icab, i ) ) )
532 lcab = int( log10( cab+sfmin ) / basl+one )
533 jc = int(rscale( i ) + sign( half, rscale( i ) ))
534 jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
535 rscale( i ) = sclfac**jc
536 360 CONTINUE
537*
538* Row scaling of matrices A and B
539*
540 DO 370 i = ilo, ihi
541 CALL dscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
542 CALL dscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
543 370 CONTINUE
544*
545* Column scaling of matrices A and B
546*
547 DO 380 j = ilo, ihi
548 CALL dscal( ihi, rscale( j ), a( 1, j ), 1 )
549 CALL dscal( ihi, rscale( j ), b( 1, j ), 1 )
550 380 CONTINUE
551*
552 RETURN
553*
554* End of DGGBAL
555*
alpha
#define alpha
Definition eval.h:35
jc
subroutine jc(p, t, a, b, cm, cn, tref, tm, epsm, sigmam, jc_yield, tan_jc)
Definition sigeps106.F:339

◆ dla_gbamv()

subroutine dla_gbamv ( integer trans,
integer m,
integer n,
integer kl,
integer ku,
double precision alpha,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( * ) x,
integer incx,
double precision beta,
double precision, dimension( * ) y,
integer incy )

DLA_GBAMV performs a matrix-vector operation to calculate error bounds.

Download DLA_GBAMV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DLA_GBAMV  performs one of the matrix-vector operations
!>
!>         y := alpha*abs(A)*abs(x) + beta*abs(y),
!>    or   y := alpha*abs(A)**T*abs(x) + beta*abs(y),
!>
!> where alpha and beta are scalars, x and y are vectors and A is an
!> m by n matrix.
!>
!> This function is primarily used in calculating error bounds.
!> To protect against underflow during evaluation, components in
!> the resulting vector are perturbed away from zero by (N+1)
!> times the underflow threshold.  To prevent unnecessarily large
!> errors for block-structure embedded in general matrices,
!>  zero components are not perturbed.  A zero
!> entry is considered  if all multiplications involved
!> in computing that entry have at least one zero multiplicand.
!>
Parameters
[in]TRANS
!>          TRANS is INTEGER
!>           On entry, TRANS specifies the operation to be performed as
!>           follows:
!>
!>             BLAS_NO_TRANS      y := alpha*abs(A)*abs(x) + beta*abs(y)
!>             BLAS_TRANS         y := alpha*abs(A**T)*abs(x) + beta*abs(y)
!>             BLAS_CONJ_TRANS    y := alpha*abs(A**T)*abs(x) + beta*abs(y)
!>
!>           Unchanged on exit.
!>
[in]M
!>          M is INTEGER
!>           On entry, M specifies the number of rows of the matrix A.
!>           M must be at least zero.
!>           Unchanged on exit.
!>
[in]N
!>          N is INTEGER
!>           On entry, N specifies the number of columns of the matrix A.
!>           N must be at least zero.
!>           Unchanged on exit.
!>
[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]ALPHA
!>          ALPHA is DOUBLE PRECISION
!>           On entry, ALPHA specifies the scalar alpha.
!>           Unchanged on exit.
!>
[in]AB
!>          AB is DOUBLE PRECISION array, dimension ( LDAB, n )
!>           Before entry, the leading m by n part of the array AB must
!>           contain the matrix of coefficients.
!>           Unchanged on exit.
!>
[in]LDAB
!>          LDAB is INTEGER
!>           On entry, LDA specifies the first dimension of AB as declared
!>           in the calling (sub) program. LDAB must be at least
!>           max( 1, m ).
!>           Unchanged on exit.
!>
[in]X
!>          X is DOUBLE PRECISION array, dimension
!>           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
!>           and at least
!>           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
!>           Before entry, the incremented array X must contain the
!>           vector x.
!>           Unchanged on exit.
!>
[in]INCX
!>          INCX is INTEGER
!>           On entry, INCX specifies the increment for the elements of
!>           X. INCX must not be zero.
!>           Unchanged on exit.
!>
[in]BETA
!>          BETA is DOUBLE PRECISION
!>           On entry, BETA specifies the scalar beta. When BETA is
!>           supplied as zero then Y need not be set on input.
!>           Unchanged on exit.
!>
[in,out]Y
!>          Y is DOUBLE PRECISION array, dimension
!>           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
!>           and at least
!>           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
!>           Before entry with BETA non-zero, the incremented array Y
!>           must contain the vector y. On exit, Y is overwritten by the
!>           updated vector y.
!>
[in]INCY
!>          INCY is INTEGER
!>           On entry, INCY specifies the increment for the elements of
!>           Y. INCY must not be zero.
!>           Unchanged on exit.
!>
!>  Level 2 Blas routine.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 183 of file dla_gbamv.f.

185*
186* -- LAPACK computational routine --
187* -- LAPACK is a software package provided by Univ. of Tennessee, --
188* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189*
190* .. Scalar Arguments ..
191 DOUBLE PRECISION ALPHA, BETA
192 INTEGER INCX, INCY, LDAB, M, N, KL, KU, TRANS
193* ..
194* .. Array Arguments ..
195 DOUBLE PRECISION AB( LDAB, * ), X( * ), Y( * )
196* ..
197*
198* =====================================================================
199*
200* .. Parameters ..
201 DOUBLE PRECISION ONE, ZERO
202 parameter( one = 1.0d+0, zero = 0.0d+0 )
203* ..
204* .. Local Scalars ..
205 LOGICAL SYMB_ZERO
206 DOUBLE PRECISION TEMP, SAFE1
207 INTEGER I, INFO, IY, J, JX, KX, KY, LENX, LENY, KD, KE
208* ..
209* .. External Subroutines ..
210 EXTERNAL xerbla, dlamch
211 DOUBLE PRECISION DLAMCH
212* ..
213* .. External Functions ..
214 EXTERNAL ilatrans
215 INTEGER ILATRANS
216* ..
217* .. Intrinsic Functions ..
218 INTRINSIC max, abs, sign
219* ..
220* .. Executable Statements ..
221*
222* Test the input parameters.
223*
224 info = 0
225 IF ( .NOT.( ( trans.EQ.ilatrans( 'N' ) )
226 $ .OR. ( trans.EQ.ilatrans( 'T' ) )
227 $ .OR. ( trans.EQ.ilatrans( 'C' ) ) ) ) THEN
228 info = 1
229 ELSE IF( m.LT.0 )THEN
230 info = 2
231 ELSE IF( n.LT.0 )THEN
232 info = 3
233 ELSE IF( kl.LT.0 .OR. kl.GT.m-1 ) THEN
234 info = 4
235 ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
236 info = 5
237 ELSE IF( ldab.LT.kl+ku+1 )THEN
238 info = 6
239 ELSE IF( incx.EQ.0 )THEN
240 info = 8
241 ELSE IF( incy.EQ.0 )THEN
242 info = 11
243 END IF
244 IF( info.NE.0 )THEN
245 CALL xerbla( 'DLA_GBAMV ', info )
246 RETURN
247 END IF
248*
249* Quick return if possible.
250*
251 IF( ( m.EQ.0 ).OR.( n.EQ.0 ).OR.
252 $ ( ( alpha.EQ.zero ).AND.( beta.EQ.one ) ) )
253 $ RETURN
254*
255* Set LENX and LENY, the lengths of the vectors x and y, and set
256* up the start points in X and Y.
257*
258 IF( trans.EQ.ilatrans( 'N' ) )THEN
259 lenx = n
260 leny = m
261 ELSE
262 lenx = m
263 leny = n
264 END IF
265 IF( incx.GT.0 )THEN
266 kx = 1
267 ELSE
268 kx = 1 - ( lenx - 1 )*incx
269 END IF
270 IF( incy.GT.0 )THEN
271 ky = 1
272 ELSE
273 ky = 1 - ( leny - 1 )*incy
274 END IF
275*
276* Set SAFE1 essentially to be the underflow threshold times the
277* number of additions in each row.
278*
279 safe1 = dlamch( 'Safe minimum' )
280 safe1 = (n+1)*safe1
281*
282* Form y := alpha*abs(A)*abs(x) + beta*abs(y).
283*
284* The O(M*N) SYMB_ZERO tests could be replaced by O(N) queries to
285* the inexact flag. Still doesn't help change the iteration order
286* to per-column.
287*
288 kd = ku + 1
289 ke = kl + 1
290 iy = ky
291 IF ( incx.EQ.1 ) THEN
292 IF( trans.EQ.ilatrans( 'N' ) )THEN
293 DO i = 1, leny
294 IF ( beta .EQ. zero ) THEN
295 symb_zero = .true.
296 y( iy ) = 0.0d+0
297 ELSE IF ( y( iy ) .EQ. zero ) THEN
298 symb_zero = .true.
299 ELSE
300 symb_zero = .false.
301 y( iy ) = beta * abs( y( iy ) )
302 END IF
303 IF ( alpha .NE. zero ) THEN
304 DO j = max( i-kl, 1 ), min( i+ku, lenx )
305 temp = abs( ab( kd+i-j, j ) )
306 symb_zero = symb_zero .AND.
307 $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
308
309 y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
310 END DO
311 END IF
312
313 IF ( .NOT.symb_zero )
314 $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
315 iy = iy + incy
316 END DO
317 ELSE
318 DO i = 1, leny
319 IF ( beta .EQ. zero ) THEN
320 symb_zero = .true.
321 y( iy ) = 0.0d+0
322 ELSE IF ( y( iy ) .EQ. zero ) THEN
323 symb_zero = .true.
324 ELSE
325 symb_zero = .false.
326 y( iy ) = beta * abs( y( iy ) )
327 END IF
328 IF ( alpha .NE. zero ) THEN
329 DO j = max( i-kl, 1 ), min( i+ku, lenx )
330 temp = abs( ab( ke-i+j, i ) )
331 symb_zero = symb_zero .AND.
332 $ ( x( j ) .EQ. zero .OR. temp .EQ. zero )
333
334 y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
335 END DO
336 END IF
337
338 IF ( .NOT.symb_zero )
339 $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
340 iy = iy + incy
341 END DO
342 END IF
343 ELSE
344 IF( trans.EQ.ilatrans( 'N' ) )THEN
345 DO i = 1, leny
346 IF ( beta .EQ. zero ) THEN
347 symb_zero = .true.
348 y( iy ) = 0.0d+0
349 ELSE IF ( y( iy ) .EQ. zero ) THEN
350 symb_zero = .true.
351 ELSE
352 symb_zero = .false.
353 y( iy ) = beta * abs( y( iy ) )
354 END IF
355 IF ( alpha .NE. zero ) THEN
356 jx = kx
357 DO j = max( i-kl, 1 ), min( i+ku, lenx )
358 temp = abs( ab( kd+i-j, j ) )
359 symb_zero = symb_zero .AND.
360 $ ( x( jx ) .EQ. zero .OR. temp .EQ. zero )
361
362 y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
363 jx = jx + incx
364 END DO
365 END IF
366
367 IF ( .NOT.symb_zero )
368 $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
369
370 iy = iy + incy
371 END DO
372 ELSE
373 DO i = 1, leny
374 IF ( beta .EQ. zero ) THEN
375 symb_zero = .true.
376 y( iy ) = 0.0d+0
377 ELSE IF ( y( iy ) .EQ. zero ) THEN
378 symb_zero = .true.
379 ELSE
380 symb_zero = .false.
381 y( iy ) = beta * abs( y( iy ) )
382 END IF
383 IF ( alpha .NE. zero ) THEN
384 jx = kx
385 DO j = max( i-kl, 1 ), min( i+ku, lenx )
386 temp = abs( ab( ke-i+j, i ) )
387 symb_zero = symb_zero .AND.
388 $ ( x( jx ) .EQ. zero .OR. temp .EQ. zero )
389
390 y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
391 jx = jx + incx
392 END DO
393 END IF
394
395 IF ( .NOT.symb_zero )
396 $ y( iy ) = y( iy ) + sign( safe1, y( iy ) )
397
398 iy = iy + incy
399 END DO
400 END IF
401
402 END IF
403*
404 RETURN
405*
406* End of DLA_GBAMV
407*

◆ dla_gbrcond()

double precision function dla_gbrcond ( character trans,
integer n,
integer kl,
integer ku,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( ldafb, * ) afb,
integer ldafb,
integer, dimension( * ) ipiv,
integer cmode,
double precision, dimension( * ) c,
integer info,
double precision, dimension( * ) work,
integer, dimension( * ) iwork )

DLA_GBRCOND estimates the Skeel condition number for a general banded matrix.

Download DLA_GBRCOND + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!>    DLA_GBRCOND Estimates the Skeel condition number of  op(A) * op2(C)
!>    where op2 is determined by CMODE as follows
!>    CMODE =  1    op2(C) = C
!>    CMODE =  0    op2(C) = I
!>    CMODE = -1    op2(C) = inv(C)
!>    The Skeel condition number  cond(A) = norminf( |inv(A)||A| )
!>    is computed by computing scaling factors R such that
!>    diag(R)*A*op2(C) is row equilibrated and computing the standard
!>    infinity-norm condition number.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>     Specifies the form of the system of equations:
!>       = 'N':  A * X = B     (No transpose)
!>       = 'T':  A**T * X = B  (Transpose)
!>       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose)
!>
[in]N
!>          N is INTEGER
!>     The number of linear equations, i.e., the order 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]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>     On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
!>     The j-th column of A is stored in the j-th column of the
!>     array AB as follows:
!>     AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
!>
[in]LDAB
!>          LDAB is INTEGER
!>     The leading dimension of the array AB.  LDAB >= KL+KU+1.
!>
[in]AFB
!>          AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
!>     Details of the LU factorization of the band matrix A, 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.
!>
[in]LDAFB
!>          LDAFB is INTEGER
!>     The leading dimension of the array AFB.  LDAFB >= 2*KL+KU+1.
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>     The pivot indices from the factorization A = P*L*U
!>     as computed by DGBTRF; row i of the matrix was interchanged
!>     with row IPIV(i).
!>
[in]CMODE
!>          CMODE is INTEGER
!>     Determines op2(C) in the formula op(A) * op2(C) as follows:
!>     CMODE =  1    op2(C) = C
!>     CMODE =  0    op2(C) = I
!>     CMODE = -1    op2(C) = inv(C)
!>
[in]C
!>          C is DOUBLE PRECISION array, dimension (N)
!>     The vector C in the formula op(A) * op2(C).
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit.
!>     i > 0:  The ith argument is invalid.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (5*N).
!>     Workspace.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N).
!>     Workspace.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 167 of file dla_gbrcond.f.

170*
171* -- LAPACK computational routine --
172* -- LAPACK is a software package provided by Univ. of Tennessee, --
173* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
174*
175* .. Scalar Arguments ..
176 CHARACTER TRANS
177 INTEGER N, LDAB, LDAFB, INFO, KL, KU, CMODE
178* ..
179* .. Array Arguments ..
180 INTEGER IWORK( * ), IPIV( * )
181 DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ),
182 $ C( * )
183* ..
184*
185* =====================================================================
186*
187* .. Local Scalars ..
188 LOGICAL NOTRANS
189 INTEGER KASE, I, J, KD, KE
190 DOUBLE PRECISION AINVNM, TMP
191* ..
192* .. Local Arrays ..
193 INTEGER ISAVE( 3 )
194* ..
195* .. External Functions ..
196 LOGICAL LSAME
197 EXTERNAL lsame
198* ..
199* .. External Subroutines ..
200 EXTERNAL dlacn2, dgbtrs, xerbla
201* ..
202* .. Intrinsic Functions ..
203 INTRINSIC abs, max
204* ..
205* .. Executable Statements ..
206*
207 dla_gbrcond = 0.0d+0
208*
209 info = 0
210 notrans = lsame( trans, 'N' )
211 IF ( .NOT. notrans .AND. .NOT. lsame(trans, 'T')
212 $ .AND. .NOT. lsame(trans, 'C') ) THEN
213 info = -1
214 ELSE IF( n.LT.0 ) THEN
215 info = -2
216 ELSE IF( kl.LT.0 .OR. kl.GT.n-1 ) THEN
217 info = -3
218 ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
219 info = -4
220 ELSE IF( ldab.LT.kl+ku+1 ) THEN
221 info = -6
222 ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
223 info = -8
224 END IF
225 IF( info.NE.0 ) THEN
226 CALL xerbla( 'DLA_GBRCOND', -info )
227 RETURN
228 END IF
229 IF( n.EQ.0 ) THEN
230 dla_gbrcond = 1.0d+0
231 RETURN
232 END IF
233*
234* Compute the equilibration matrix R such that
235* inv(R)*A*C has unit 1-norm.
236*
237 kd = ku + 1
238 ke = kl + 1
239 IF ( notrans ) THEN
240 DO i = 1, n
241 tmp = 0.0d+0
242 IF ( cmode .EQ. 1 ) THEN
243 DO j = max( i-kl, 1 ), min( i+ku, n )
244 tmp = tmp + abs( ab( kd+i-j, j ) * c( j ) )
245 END DO
246 ELSE IF ( cmode .EQ. 0 ) THEN
247 DO j = max( i-kl, 1 ), min( i+ku, n )
248 tmp = tmp + abs( ab( kd+i-j, j ) )
249 END DO
250 ELSE
251 DO j = max( i-kl, 1 ), min( i+ku, n )
252 tmp = tmp + abs( ab( kd+i-j, j ) / c( j ) )
253 END DO
254 END IF
255 work( 2*n+i ) = tmp
256 END DO
257 ELSE
258 DO i = 1, n
259 tmp = 0.0d+0
260 IF ( cmode .EQ. 1 ) THEN
261 DO j = max( i-kl, 1 ), min( i+ku, n )
262 tmp = tmp + abs( ab( ke-i+j, i ) * c( j ) )
263 END DO
264 ELSE IF ( cmode .EQ. 0 ) THEN
265 DO j = max( i-kl, 1 ), min( i+ku, n )
266 tmp = tmp + abs( ab( ke-i+j, i ) )
267 END DO
268 ELSE
269 DO j = max( i-kl, 1 ), min( i+ku, n )
270 tmp = tmp + abs( ab( ke-i+j, i ) / c( j ) )
271 END DO
272 END IF
273 work( 2*n+i ) = tmp
274 END DO
275 END IF
276*
277* Estimate the norm of inv(op(A)).
278*
279 ainvnm = 0.0d+0
280
281 kase = 0
282 10 CONTINUE
283 CALL dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
284 IF( kase.NE.0 ) THEN
285 IF( kase.EQ.2 ) THEN
286*
287* Multiply by R.
288*
289 DO i = 1, n
290 work( i ) = work( i ) * work( 2*n+i )
291 END DO
292
293 IF ( notrans ) THEN
294 CALL dgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
295 $ ipiv, work, n, info )
296 ELSE
297 CALL dgbtrs( 'Transpose', n, kl, ku, 1, afb, ldafb, ipiv,
298 $ work, n, info )
299 END IF
300*
301* Multiply by inv(C).
302*
303 IF ( cmode .EQ. 1 ) THEN
304 DO i = 1, n
305 work( i ) = work( i ) / c( i )
306 END DO
307 ELSE IF ( cmode .EQ. -1 ) THEN
308 DO i = 1, n
309 work( i ) = work( i ) * c( i )
310 END DO
311 END IF
312 ELSE
313*
314* Multiply by inv(C**T).
315*
316 IF ( cmode .EQ. 1 ) THEN
317 DO i = 1, n
318 work( i ) = work( i ) / c( i )
319 END DO
320 ELSE IF ( cmode .EQ. -1 ) THEN
321 DO i = 1, n
322 work( i ) = work( i ) * c( i )
323 END DO
324 END IF
325
326 IF ( notrans ) THEN
327 CALL dgbtrs( 'Transpose', n, kl, ku, 1, afb, ldafb, ipiv,
328 $ work, n, info )
329 ELSE
330 CALL dgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
331 $ ipiv, work, n, info )
332 END IF
333*
334* Multiply by R.
335*
336 DO i = 1, n
337 work( i ) = work( i ) * work( 2*n+i )
338 END DO
339 END IF
340 GO TO 10
341 END IF
342*
343* Compute the estimate of the reciprocal condition number.
344*
345 IF( ainvnm .NE. 0.0d+0 )
346 $ dla_gbrcond = ( 1.0d+0 / ainvnm )
347*
348 RETURN
349*
350* End of DLA_GBRCOND
351*

◆ dla_gbrfsx_extended()

subroutine dla_gbrfsx_extended ( integer prec_type,
integer trans_type,
integer n,
integer kl,
integer ku,
integer nrhs,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( ldafb, * ) afb,
integer ldafb,
integer, dimension( * ) ipiv,
logical colequ,
double precision, dimension( * ) c,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( ldy, * ) y,
integer ldy,
double precision, dimension(*) berr_out,
integer n_norms,
double precision, dimension( nrhs, * ) err_bnds_norm,
double precision, dimension( nrhs, * ) err_bnds_comp,
double precision, dimension(*) res,
double precision, dimension(*) ayb,
double precision, dimension(*) dy,
double precision, dimension(*) y_tail,
double precision rcond,
integer ithresh,
double precision rthresh,
double precision dz_ub,
logical ignore_cwise,
integer info )

DLA_GBRFSX_EXTENDED improves the computed solution to a system of linear equations for general banded matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.

Download DLA_GBRFSX_EXTENDED + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!>
!> DLA_GBRFSX_EXTENDED improves the computed solution to a system of
!> linear equations by performing extra-precise iterative refinement
!> and provides error bounds and backward error estimates for the solution.
!> This subroutine is called by DGBRFSX to perform iterative refinement.
!> In addition to normwise error bound, the code provides maximum
!> componentwise error bound if possible. See comments for ERR_BNDS_NORM
!> and ERR_BNDS_COMP for details of the error bounds. Note that this
!> subroutine is only responsible for setting the second fields of
!> ERR_BNDS_NORM and ERR_BNDS_COMP.
!>
Parameters
[in]PREC_TYPE
!>          PREC_TYPE is INTEGER
!>     Specifies the intermediate precision to be used in refinement.
!>     The value is defined by ILAPREC(P) where P is a CHARACTER and P
!>          = 'S':  Single
!>          = 'D':  Double
!>          = 'I':  Indigenous
!>          = 'X' or 'E':  Extra
!>
[in]TRANS_TYPE
!>          TRANS_TYPE is INTEGER
!>     Specifies the transposition operation on A.
!>     The value is defined by ILATRANS(T) where T is a CHARACTER and T
!>          = 'N':  No transpose
!>          = 'T':  Transpose
!>          = 'C':  Conjugate transpose
!>
[in]N
!>          N is INTEGER
!>     The number of linear equations, i.e., the order 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]NRHS
!>          NRHS is INTEGER
!>     The number of right-hand-sides, i.e., the number of columns of the
!>     matrix B.
!>
[in]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>          On entry, the N-by-N matrix AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDBA >= max(1,N).
!>
[in]AFB
!>          AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
!>     The factors L and U from the factorization
!>     A = P*L*U as computed by DGBTRF.
!>
[in]LDAFB
!>          LDAFB is INTEGER
!>     The leading dimension of the array AF.  LDAFB >= max(1,N).
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>     The pivot indices from the factorization A = P*L*U
!>     as computed by DGBTRF; row i of the matrix was interchanged
!>     with row IPIV(i).
!>
[in]COLEQU
!>          COLEQU is LOGICAL
!>     If .TRUE. then column equilibration was done to A before calling
!>     this routine. This is needed to compute the solution and error
!>     bounds correctly.
!>
[in]C
!>          C is DOUBLE PRECISION array, dimension (N)
!>     The column scale factors for A. If COLEQU = .FALSE., C
!>     is not accessed. If C is input, each element of C should be a power
!>     of the radix to ensure a reliable solution and error estimates.
!>     Scaling by powers of the radix does not cause rounding errors unless
!>     the result underflows or overflows. Rounding errors during scaling
!>     lead to refining with a matrix that is not equivalent to the
!>     input matrix, producing error estimates that may not be
!>     reliable.
!>
[in]B
!>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
!>     The right-hand-side matrix B.
!>
[in]LDB
!>          LDB is INTEGER
!>     The leading dimension of the array B.  LDB >= max(1,N).
!>
[in,out]Y
!>          Y is DOUBLE PRECISION array, dimension (LDY,NRHS)
!>     On entry, the solution matrix X, as computed by DGBTRS.
!>     On exit, the improved solution matrix Y.
!>
[in]LDY
!>          LDY is INTEGER
!>     The leading dimension of the array Y.  LDY >= max(1,N).
!>
[out]BERR_OUT
!>          BERR_OUT is DOUBLE PRECISION array, dimension (NRHS)
!>     On exit, BERR_OUT(j) contains the componentwise relative backward
!>     error for right-hand-side j from the formula
!>         max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
!>     where abs(Z) is the componentwise absolute value of the matrix
!>     or vector Z. This is computed by DLA_LIN_BERR.
!>
[in]N_NORMS
!>          N_NORMS is INTEGER
!>     Determines which error bounds to return (see ERR_BNDS_NORM
!>     and ERR_BNDS_COMP).
!>     If N_NORMS >= 1 return normwise error bounds.
!>     If N_NORMS >= 2 return componentwise error bounds.
!>
[in,out]ERR_BNDS_NORM
!>          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     normwise relative error, which is defined as follows:
!>
!>     Normwise relative error in the ith solution vector:
!>             max_j (abs(XTRUE(j,i) - X(j,i)))
!>            ------------------------------
!>                  max_j abs(X(j,i))
!>
!>     The array is indexed by the type of error information as described
!>     below. There currently are up to three pieces of information
!>     returned.
!>
!>     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_NORM(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated normwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*A, where S scales each row by a power of the
!>              radix so all absolute row sums of Z are approximately 1.
!>
!>     This subroutine is only responsible for setting the second field
!>     above.
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in,out]ERR_BNDS_COMP
!>          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
!>     For each right-hand side, this array contains information about
!>     various error bounds and condition numbers corresponding to the
!>     componentwise relative error, which is defined as follows:
!>
!>     Componentwise relative error in the ith solution vector:
!>                    abs(XTRUE(j,i) - X(j,i))
!>             max_j ----------------------
!>                         abs(X(j,i))
!>
!>     The array is indexed by the right-hand side i (on which the
!>     componentwise relative error depends), and the type of error
!>     information as described below. There currently are up to three
!>     pieces of information returned for each right-hand side. If
!>     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
!>     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS < 3, then at most
!>     the first (:,N_ERR_BNDS) entries are returned.
!>
!>     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
!>     right-hand side.
!>
!>     The second index in ERR_BNDS_COMP(:,err) contains the following
!>     three fields:
!>     err = 1  boolean. Trust the answer if the
!>              reciprocal condition number is less than the threshold
!>              sqrt(n) * slamch('Epsilon').
!>
!>     err = 2  error bound: The estimated forward error,
!>              almost certainly within a factor of 10 of the true error
!>              so long as the next entry is greater than the threshold
!>              sqrt(n) * slamch('Epsilon'). This error bound should only
!>              be trusted if the previous boolean is true.
!>
!>     err = 3  Reciprocal condition number: Estimated componentwise
!>              reciprocal condition number.  Compared with the threshold
!>              sqrt(n) * slamch('Epsilon') to determine if the error
!>              estimate is . These reciprocal condition
!>              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
!>              appropriately scaled matrix Z.
!>              Let Z = S*(A*diag(x)), where x is the solution for the
!>              current right-hand side and S scales each row of
!>              A*diag(x) by a power of the radix so all absolute row
!>              sums of Z are approximately 1.
!>
!>     This subroutine is only responsible for setting the second field
!>     above.
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[in]RES
!>          RES is DOUBLE PRECISION array, dimension (N)
!>     Workspace to hold the intermediate residual.
!>
[in]AYB
!>          AYB is DOUBLE PRECISION array, dimension (N)
!>     Workspace. This can be the same workspace passed for Y_TAIL.
!>
[in]DY
!>          DY is DOUBLE PRECISION array, dimension (N)
!>     Workspace to hold the intermediate solution.
!>
[in]Y_TAIL
!>          Y_TAIL is DOUBLE PRECISION array, dimension (N)
!>     Workspace to hold the trailing bits of the intermediate solution.
!>
[in]RCOND
!>          RCOND is DOUBLE PRECISION
!>     Reciprocal scaled condition number.  This is an estimate of the
!>     reciprocal Skeel condition number of the matrix A after
!>     equilibration (if done).  If this is less than the machine
!>     precision (in particular, if it is zero), the matrix is singular
!>     to working precision.  Note that the error may still be small even
!>     if this number is very small and the matrix appears ill-
!>     conditioned.
!>
[in]ITHRESH
!>          ITHRESH is INTEGER
!>     The maximum number of residual computations allowed for
!>     refinement. The default is 10. For 'aggressive' set to 100 to
!>     permit convergence using approximate factorizations or
!>     factorizations other than LU. If the factorization uses a
!>     technique other than Gaussian elimination, the guarantees in
!>     ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy.
!>
[in]RTHRESH
!>          RTHRESH is DOUBLE PRECISION
!>     Determines when to stop refinement if the error estimate stops
!>     decreasing. Refinement will stop when the next solution no longer
!>     satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is
!>     the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The
!>     default value is 0.5. For 'aggressive' set to 0.9 to permit
!>     convergence on extremely ill-conditioned matrices. See LAWN 165
!>     for more details.
!>
[in]DZ_UB
!>          DZ_UB is DOUBLE PRECISION
!>     Determines when to start considering componentwise convergence.
!>     Componentwise convergence is only considered after each component
!>     of the solution Y is stable, which we define as the relative
!>     change in each component being less than DZ_UB. The default value
!>     is 0.25, requiring the first bit to be stable. See LAWN 165 for
!>     more details.
!>
[in]IGNORE_CWISE
!>          IGNORE_CWISE is LOGICAL
!>     If .TRUE. then ignore componentwise convergence. Default value
!>     is .FALSE..
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit.
!>       < 0:  if INFO = -i, the ith argument to DGBTRS had an illegal
!>             value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 404 of file dla_gbrfsx_extended.f.

411*
412* -- LAPACK computational routine --
413* -- LAPACK is a software package provided by Univ. of Tennessee, --
414* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
415*
416* .. Scalar Arguments ..
417 INTEGER INFO, LDAB, LDAFB, LDB, LDY, N, KL, KU, NRHS,
418 $ PREC_TYPE, TRANS_TYPE, N_NORMS, ITHRESH
419 LOGICAL COLEQU, IGNORE_CWISE
420 DOUBLE PRECISION RTHRESH, DZ_UB
421* ..
422* .. Array Arguments ..
423 INTEGER IPIV( * )
424 DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
425 $ Y( LDY, * ), RES(*), DY(*), Y_TAIL(*)
426 DOUBLE PRECISION C( * ), AYB(*), RCOND, BERR_OUT(*),
427 $ ERR_BNDS_NORM( NRHS, * ),
428 $ ERR_BNDS_COMP( NRHS, * )
429* ..
430*
431* =====================================================================
432*
433* .. Local Scalars ..
434 CHARACTER TRANS
435 INTEGER CNT, I, J, M, X_STATE, Z_STATE, Y_PREC_STATE
436 DOUBLE PRECISION YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
437 $ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX,
438 $ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z,
439 $ EPS, HUGEVAL, INCR_THRESH
440 LOGICAL INCR_PREC
441* ..
442* .. Parameters ..
443 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
444 $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL,
445 $ EXTRA_Y
446 parameter( unstable_state = 0, working_state = 1,
447 $ conv_state = 2, noprog_state = 3 )
448 parameter( base_residual = 0, extra_residual = 1,
449 $ extra_y = 2 )
450 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
451 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
452 INTEGER CMP_ERR_I, PIV_GROWTH_I
453 parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
454 $ berr_i = 3 )
455 parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
456 parameter( cmp_rcond_i = 7, cmp_err_i = 8,
457 $ piv_growth_i = 9 )
458 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
459 $ LA_LINRX_CWISE_I
460 parameter( la_linrx_itref_i = 1,
461 $ la_linrx_ithresh_i = 2 )
462 parameter( la_linrx_cwise_i = 3 )
463 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
464 $ LA_LINRX_RCOND_I
465 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
466 parameter( la_linrx_rcond_i = 3 )
467* ..
468* .. External Subroutines ..
469 EXTERNAL daxpy, dcopy, dgbtrs, dgbmv, blas_dgbmv_x,
470 $ blas_dgbmv2_x, dla_gbamv, dla_wwaddw, dlamch,
471 $ chla_transtype, dla_lin_berr
472 DOUBLE PRECISION DLAMCH
473 CHARACTER CHLA_TRANSTYPE
474* ..
475* .. Intrinsic Functions ..
476 INTRINSIC abs, max, min
477* ..
478* .. Executable Statements ..
479*
480 IF (info.NE.0) RETURN
481 trans = chla_transtype(trans_type)
482 eps = dlamch( 'Epsilon' )
483 hugeval = dlamch( 'Overflow' )
484* Force HUGEVAL to Inf
485 hugeval = hugeval * hugeval
486* Using HUGEVAL may lead to spurious underflows.
487 incr_thresh = dble( n ) * eps
488 m = kl+ku+1
489
490 DO j = 1, nrhs
491 y_prec_state = extra_residual
492 IF ( y_prec_state .EQ. extra_y ) THEN
493 DO i = 1, n
494 y_tail( i ) = 0.0d+0
495 END DO
496 END IF
497
498 dxrat = 0.0d+0
499 dxratmax = 0.0d+0
500 dzrat = 0.0d+0
501 dzratmax = 0.0d+0
502 final_dx_x = hugeval
503 final_dz_z = hugeval
504 prevnormdx = hugeval
505 prev_dz_z = hugeval
506 dz_z = hugeval
507 dx_x = hugeval
508
509 x_state = working_state
510 z_state = unstable_state
511 incr_prec = .false.
512
513 DO cnt = 1, ithresh
514*
515* Compute residual RES = B_s - op(A_s) * Y,
516* op(A) = A, A**T, or A**H depending on TRANS (and type).
517*
518 CALL dcopy( n, b( 1, j ), 1, res, 1 )
519 IF ( y_prec_state .EQ. base_residual ) THEN
520 CALL dgbmv( trans, m, n, kl, ku, -1.0d+0, ab, ldab,
521 $ y( 1, j ), 1, 1.0d+0, res, 1 )
522 ELSE IF ( y_prec_state .EQ. extra_residual ) THEN
523 CALL blas_dgbmv_x( trans_type, n, n, kl, ku,
524 $ -1.0d+0, ab, ldab, y( 1, j ), 1, 1.0d+0, res, 1,
525 $ prec_type )
526 ELSE
527 CALL blas_dgbmv2_x( trans_type, n, n, kl, ku, -1.0d+0,
528 $ ab, ldab, y( 1, j ), y_tail, 1, 1.0d+0, res, 1,
529 $ prec_type )
530 END IF
531
532! XXX: RES is no longer needed.
533 CALL dcopy( n, res, 1, dy, 1 )
534 CALL dgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv, dy, n,
535 $ info )
536*
537* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
538*
539 normx = 0.0d+0
540 normy = 0.0d+0
541 normdx = 0.0d+0
542 dz_z = 0.0d+0
543 ymin = hugeval
544
545 DO i = 1, n
546 yk = abs( y( i, j ) )
547 dyk = abs( dy( i ) )
548
549 IF ( yk .NE. 0.0d+0 ) THEN
550 dz_z = max( dz_z, dyk / yk )
551 ELSE IF ( dyk .NE. 0.0d+0 ) THEN
552 dz_z = hugeval
553 END IF
554
555 ymin = min( ymin, yk )
556
557 normy = max( normy, yk )
558
559 IF ( colequ ) THEN
560 normx = max( normx, yk * c( i ) )
561 normdx = max( normdx, dyk * c( i ) )
562 ELSE
563 normx = normy
564 normdx = max( normdx, dyk )
565 END IF
566 END DO
567
568 IF ( normx .NE. 0.0d+0 ) THEN
569 dx_x = normdx / normx
570 ELSE IF ( normdx .EQ. 0.0d+0 ) THEN
571 dx_x = 0.0d+0
572 ELSE
573 dx_x = hugeval
574 END IF
575
576 dxrat = normdx / prevnormdx
577 dzrat = dz_z / prev_dz_z
578*
579* Check termination criteria.
580*
581 IF ( .NOT.ignore_cwise
582 $ .AND. ymin*rcond .LT. incr_thresh*normy
583 $ .AND. y_prec_state .LT. extra_y )
584 $ incr_prec = .true.
585
586 IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
587 $ x_state = working_state
588 IF ( x_state .EQ. working_state ) THEN
589 IF ( dx_x .LE. eps ) THEN
590 x_state = conv_state
591 ELSE IF ( dxrat .GT. rthresh ) THEN
592 IF ( y_prec_state .NE. extra_y ) THEN
593 incr_prec = .true.
594 ELSE
595 x_state = noprog_state
596 END IF
597 ELSE
598 IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
599 END IF
600 IF ( x_state .GT. working_state ) final_dx_x = dx_x
601 END IF
602
603 IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
604 $ z_state = working_state
605 IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
606 $ z_state = working_state
607 IF ( z_state .EQ. working_state ) THEN
608 IF ( dz_z .LE. eps ) THEN
609 z_state = conv_state
610 ELSE IF ( dz_z .GT. dz_ub ) THEN
611 z_state = unstable_state
612 dzratmax = 0.0d+0
613 final_dz_z = hugeval
614 ELSE IF ( dzrat .GT. rthresh ) THEN
615 IF ( y_prec_state .NE. extra_y ) THEN
616 incr_prec = .true.
617 ELSE
618 z_state = noprog_state
619 END IF
620 ELSE
621 IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
622 END IF
623 IF ( z_state .GT. working_state ) final_dz_z = dz_z
624 END IF
625*
626* Exit if both normwise and componentwise stopped working,
627* but if componentwise is unstable, let it go at least two
628* iterations.
629*
630 IF ( x_state.NE.working_state ) THEN
631 IF ( ignore_cwise ) GOTO 666
632 IF ( z_state.EQ.noprog_state .OR. z_state.EQ.conv_state )
633 $ GOTO 666
634 IF ( z_state.EQ.unstable_state .AND. cnt.GT.1 ) GOTO 666
635 END IF
636
637 IF ( incr_prec ) THEN
638 incr_prec = .false.
639 y_prec_state = y_prec_state + 1
640 DO i = 1, n
641 y_tail( i ) = 0.0d+0
642 END DO
643 END IF
644
645 prevnormdx = normdx
646 prev_dz_z = dz_z
647*
648* Update soluton.
649*
650 IF (y_prec_state .LT. extra_y) THEN
651 CALL daxpy( n, 1.0d+0, dy, 1, y(1,j), 1 )
652 ELSE
653 CALL dla_wwaddw( n, y(1,j), y_tail, dy )
654 END IF
655
656 END DO
657* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
658 666 CONTINUE
659*
660* Set final_* when cnt hits ithresh.
661*
662 IF ( x_state .EQ. working_state ) final_dx_x = dx_x
663 IF ( z_state .EQ. working_state ) final_dz_z = dz_z
664*
665* Compute error bounds.
666*
667 IF ( n_norms .GE. 1 ) THEN
668 err_bnds_norm( j, la_linrx_err_i ) =
669 $ final_dx_x / (1 - dxratmax)
670 END IF
671 IF (n_norms .GE. 2) THEN
672 err_bnds_comp( j, la_linrx_err_i ) =
673 $ final_dz_z / (1 - dzratmax)
674 END IF
675*
676* Compute componentwise relative backward error from formula
677* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
678* where abs(Z) is the componentwise absolute value of the matrix
679* or vector Z.
680*
681* Compute residual RES = B_s - op(A_s) * Y,
682* op(A) = A, A**T, or A**H depending on TRANS (and type).
683*
684 CALL dcopy( n, b( 1, j ), 1, res, 1 )
685 CALL dgbmv(trans, n, n, kl, ku, -1.0d+0, ab, ldab, y(1,j),
686 $ 1, 1.0d+0, res, 1 )
687
688 DO i = 1, n
689 ayb( i ) = abs( b( i, j ) )
690 END DO
691*
692* Compute abs(op(A_s))*abs(Y) + abs(B_s).
693*
694 CALL dla_gbamv( trans_type, n, n, kl, ku, 1.0d+0,
695 $ ab, ldab, y(1, j), 1, 1.0d+0, ayb, 1 )
696
697 CALL dla_lin_berr( n, n, 1, res, ayb, berr_out( j ) )
698*
699* End of loop for each RHS
700*
701 END DO
702*
703 RETURN
704*
705* End of DLA_GBRFSX_EXTENDED
706*
chla_transtype
character *1 function chla_transtype(trans)
CHLA_TRANSTYPE
Definition chla_transtype.f:58
dla_gbamv
subroutine dla_gbamv(trans, m, n, kl, ku, alpha, ab, ldab, x, incx, beta, y, incy)
DLA_GBAMV performs a matrix-vector operation to calculate error bounds.
Definition dla_gbamv.f:185
dla_lin_berr
subroutine dla_lin_berr(n, nz, nrhs, res, ayb, berr)
DLA_LIN_BERR computes a component-wise relative backward error.
Definition dla_lin_berr.f:101
dla_wwaddw
subroutine dla_wwaddw(n, x, y, w)
DLA_WWADDW adds a vector into a doubled-single vector.
Definition dla_wwaddw.f:81

◆ dla_gbrpvgrw()

double precision function dla_gbrpvgrw ( integer n,
integer kl,
integer ku,
integer ncols,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( ldafb, * ) afb,
integer ldafb )

DLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix.

Download DLA_GBRPVGRW + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DLA_GBRPVGRW computes the reciprocal pivot growth factor
!> norm(A)/norm(U). The  norm is used. If this is
!> much less than 1, the stability of the LU factorization of the
!> (equilibrated) matrix A could be poor. This also means that the
!> solution X, estimated condition numbers, and error bounds could be
!> unreliable.
!>
Parameters
[in]N
!>          N is INTEGER
!>     The number of linear equations, i.e., the order 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]NCOLS
!>          NCOLS is INTEGER
!>     The number of columns of the matrix A.  NCOLS >= 0.
!>
[in]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>     On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
!>     The j-th column of A is stored in the j-th column of the
!>     array AB as follows:
!>     AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
!>
[in]LDAB
!>          LDAB is INTEGER
!>     The leading dimension of the array AB.  LDAB >= KL+KU+1.
!>
[in]AFB
!>          AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
!>     Details of the LU factorization of the band matrix A, 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.
!>
[in]LDAFB
!>          LDAFB is INTEGER
!>     The leading dimension of the array AFB.  LDAFB >= 2*KL+KU+1.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 115 of file dla_gbrpvgrw.f.

117*
118* -- LAPACK computational routine --
119* -- LAPACK is a software package provided by Univ. of Tennessee, --
120* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
121*
122* .. Scalar Arguments ..
123 INTEGER N, KL, KU, NCOLS, LDAB, LDAFB
124* ..
125* .. Array Arguments ..
126 DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * )
127* ..
128*
129* =====================================================================
130*
131* .. Local Scalars ..
132 INTEGER I, J, KD
133 DOUBLE PRECISION AMAX, UMAX, RPVGRW
134* ..
135* .. Intrinsic Functions ..
136 INTRINSIC abs, max, min
137* ..
138* .. Executable Statements ..
139*
140 rpvgrw = 1.0d+0
141
142 kd = ku + 1
143 DO j = 1, ncols
144 amax = 0.0d+0
145 umax = 0.0d+0
146 DO i = max( j-ku, 1 ), min( j+kl, n )
147 amax = max( abs( ab( kd+i-j, j)), amax )
148 END DO
149 DO i = max( j-ku, 1 ), j
150 umax = max( abs( afb( kd+i-j, j ) ), umax )
151 END DO
152 IF ( umax /= 0.0d+0 ) THEN
153 rpvgrw = min( amax / umax, rpvgrw )
154 END IF
155 END DO
156 dla_gbrpvgrw = rpvgrw
157*
158* End of DLA_GBRPVGRW
159*
dla_gbrpvgrw
double precision function dla_gbrpvgrw(n, kl, ku, ncols, ab, ldab, afb, ldafb)
DLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix.
Definition dla_gbrpvgrw.f:117

◆ dorgbr()

subroutine dorgbr ( character vect,
integer m,
integer n,
integer k,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) tau,
double precision, dimension( * ) work,
integer lwork,
integer info )

DORGBR

Download DORGBR + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DORGBR generates one of the real orthogonal matrices Q or P**T
!> determined by DGEBRD when reducing a real matrix A to bidiagonal
!> form: A = Q * B * P**T.  Q and P**T are defined as products of
!> elementary reflectors H(i) or G(i) respectively.
!>
!> If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q
!> is of order M:
!> if m >= k, Q = H(1) H(2) . . . H(k) and DORGBR returns the first n
!> columns of Q, where m >= n >= k;
!> if m < k, Q = H(1) H(2) . . . H(m-1) and DORGBR returns Q as an
!> M-by-M matrix.
!>
!> If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T
!> is of order N:
!> if k < n, P**T = G(k) . . . G(2) G(1) and DORGBR returns the first m
!> rows of P**T, where n >= m >= k;
!> if k >= n, P**T = G(n-1) . . . G(2) G(1) and DORGBR returns P**T as
!> an N-by-N matrix.
!>
Parameters
[in]VECT
!>          VECT is CHARACTER*1
!>          Specifies whether the matrix Q or the matrix P**T is
!>          required, as defined in the transformation applied by DGEBRD:
!>          = 'Q':  generate Q;
!>          = 'P':  generate P**T.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix Q or P**T to be returned.
!>          M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix Q or P**T to be returned.
!>          N >= 0.
!>          If VECT = 'Q', M >= N >= min(M,K);
!>          if VECT = 'P', N >= M >= min(N,K).
!>
[in]K
!>          K is INTEGER
!>          If VECT = 'Q', the number of columns in the original M-by-K
!>          matrix reduced by DGEBRD.
!>          If VECT = 'P', the number of rows in the original K-by-N
!>          matrix reduced by DGEBRD.
!>          K >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>          On entry, the vectors which define the elementary reflectors,
!>          as returned by DGEBRD.
!>          On exit, the M-by-N matrix Q or P**T.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,M).
!>
[in]TAU
!>          TAU is DOUBLE PRECISION array, dimension
!>                                (min(M,K)) if VECT = 'Q'
!>                                (min(N,K)) if VECT = 'P'
!>          TAU(i) must contain the scalar factor of the elementary
!>          reflector H(i) or G(i), which determines Q or P**T, as
!>          returned by DGEBRD in its array argument TAUQ or TAUP.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK. LWORK >= max(1,min(M,N)).
!>          For optimum performance LWORK >= min(M,N)*NB, where NB
!>          is the optimal blocksize.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 156 of file dorgbr.f.

157*
158* -- LAPACK computational routine --
159* -- LAPACK is a software package provided by Univ. of Tennessee, --
160* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
161*
162* .. Scalar Arguments ..
163 CHARACTER VECT
164 INTEGER INFO, K, LDA, LWORK, M, N
165* ..
166* .. Array Arguments ..
167 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
168* ..
169*
170* =====================================================================
171*
172* .. Parameters ..
173 DOUBLE PRECISION ZERO, ONE
174 parameter( zero = 0.0d+0, one = 1.0d+0 )
175* ..
176* .. Local Scalars ..
177 LOGICAL LQUERY, WANTQ
178 INTEGER I, IINFO, J, LWKOPT, MN
179* ..
180* .. External Functions ..
181 LOGICAL LSAME
182 EXTERNAL lsame
183* ..
184* .. External Subroutines ..
185 EXTERNAL dorglq, dorgqr, xerbla
186* ..
187* .. Intrinsic Functions ..
188 INTRINSIC max, min
189* ..
190* .. Executable Statements ..
191*
192* Test the input arguments
193*
194 info = 0
195 wantq = lsame( vect, 'Q' )
196 mn = min( m, n )
197 lquery = ( lwork.EQ.-1 )
198 IF( .NOT.wantq .AND. .NOT.lsame( vect, 'P' ) ) THEN
199 info = -1
200 ELSE IF( m.LT.0 ) THEN
201 info = -2
202 ELSE IF( n.LT.0 .OR. ( wantq .AND. ( n.GT.m .OR. n.LT.min( m,
203 $ k ) ) ) .OR. ( .NOT.wantq .AND. ( m.GT.n .OR. m.LT.
204 $ min( n, k ) ) ) ) THEN
205 info = -3
206 ELSE IF( k.LT.0 ) THEN
207 info = -4
208 ELSE IF( lda.LT.max( 1, m ) ) THEN
209 info = -6
210 ELSE IF( lwork.LT.max( 1, mn ) .AND. .NOT.lquery ) THEN
211 info = -9
212 END IF
213*
214 IF( info.EQ.0 ) THEN
215 work( 1 ) = 1
216 IF( wantq ) THEN
217 IF( m.GE.k ) THEN
218 CALL dorgqr( m, n, k, a, lda, tau, work, -1, iinfo )
219 ELSE
220 IF( m.GT.1 ) THEN
221 CALL dorgqr( m-1, m-1, m-1, a, lda, tau, work, -1,
222 $ iinfo )
223 END IF
224 END IF
225 ELSE
226 IF( k.LT.n ) THEN
227 CALL dorglq( m, n, k, a, lda, tau, work, -1, iinfo )
228 ELSE
229 IF( n.GT.1 ) THEN
230 CALL dorglq( n-1, n-1, n-1, a, lda, tau, work, -1,
231 $ iinfo )
232 END IF
233 END IF
234 END IF
235 lwkopt = work( 1 )
236 lwkopt = max(lwkopt, mn)
237 END IF
238*
239 IF( info.NE.0 ) THEN
240 CALL xerbla( 'DORGBR', -info )
241 RETURN
242 ELSE IF( lquery ) THEN
243 work( 1 ) = lwkopt
244 RETURN
245 END IF
246*
247* Quick return if possible
248*
249 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
250 work( 1 ) = 1
251 RETURN
252 END IF
253*
254 IF( wantq ) THEN
255*
256* Form Q, determined by a call to DGEBRD to reduce an m-by-k
257* matrix
258*
259 IF( m.GE.k ) THEN
260*
261* If m >= k, assume m >= n >= k
262*
263 CALL dorgqr( m, n, k, a, lda, tau, work, lwork, iinfo )
264*
265 ELSE
266*
267* If m < k, assume m = n
268*
269* Shift the vectors which define the elementary reflectors one
270* column to the right, and set the first row and column of Q
271* to those of the unit matrix
272*
273 DO 20 j = m, 2, -1
274 a( 1, j ) = zero
275 DO 10 i = j + 1, m
276 a( i, j ) = a( i, j-1 )
277 10 CONTINUE
278 20 CONTINUE
279 a( 1, 1 ) = one
280 DO 30 i = 2, m
281 a( i, 1 ) = zero
282 30 CONTINUE
283 IF( m.GT.1 ) THEN
284*
285* Form Q(2:m,2:m)
286*
287 CALL dorgqr( m-1, m-1, m-1, a( 2, 2 ), lda, tau, work,
288 $ lwork, iinfo )
289 END IF
290 END IF
291 ELSE
292*
293* Form P**T, determined by a call to DGEBRD to reduce a k-by-n
294* matrix
295*
296 IF( k.LT.n ) THEN
297*
298* If k < n, assume k <= m <= n
299*
300 CALL dorglq( m, n, k, a, lda, tau, work, lwork, iinfo )
301*
302 ELSE
303*
304* If k >= n, assume m = n
305*
306* Shift the vectors which define the elementary reflectors one
307* row downward, and set the first row and column of P**T to
308* those of the unit matrix
309*
310 a( 1, 1 ) = one
311 DO 40 i = 2, n
312 a( i, 1 ) = zero
313 40 CONTINUE
314 DO 60 j = 2, n
315 DO 50 i = j - 1, 2, -1
316 a( i, j ) = a( i-1, j )
317 50 CONTINUE
318 a( 1, j ) = zero
319 60 CONTINUE
320 IF( n.GT.1 ) THEN
321*
322* Form P**T(2:n,2:n)
323*
324 CALL dorglq( n-1, n-1, n-1, a( 2, 2 ), lda, tau, work,
325 $ lwork, iinfo )
326 END IF
327 END IF
328 END IF
329 work( 1 ) = lwkopt
330 RETURN
331*
332* End of DORGBR
333*
dorgqr
subroutine dorgqr(m, n, k, a, lda, tau, work, lwork, info)
DORGQR
Definition dorgqr.f:128
dorglq
subroutine dorglq(m, n, k, a, lda, tau, work, lwork, info)
DORGLQ
Definition dorglq.f:127