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

Functions

subroutine sgbbrd (vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, c, ldc, work, info)
 SGBBRD
subroutine sgbcon (norm, n, kl, ku, ab, ldab, ipiv, anorm, rcond, work, iwork, info)
 SGBCON
subroutine sgbequ (m, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd, amax, info)
 SGBEQU
subroutine sgbequb (m, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd, amax, info)
 SGBEQUB
subroutine sgbrfs (trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, b, ldb, x, ldx, ferr, berr, work, iwork, info)
 SGBRFS
subroutine sgbrfsx (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)
 SGBRFSX
subroutine sgbtf2 (m, n, kl, ku, ab, ldab, ipiv, info)
 SGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algorithm.
subroutine sgbtrf (m, n, kl, ku, ab, ldab, ipiv, info)
 SGBTRF
subroutine sgbtrs (trans, n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
 SGBTRS
subroutine sggbak (job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
 SGGBAK
subroutine sggbal (job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
 SGGBAL
subroutine sla_gbamv (trans, m, n, kl, ku, alpha, ab, ldab, x, incx, beta, y, incy)
 SLA_GBAMV performs a matrix-vector operation to calculate error bounds.
real function sla_gbrcond (trans, n, kl, ku, ab, ldab, afb, ldafb, ipiv, cmode, c, info, work, iwork)
 SLA_GBRCOND estimates the Skeel condition number for a general banded matrix.
subroutine sla_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)
 SLA_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.
real function sla_gbrpvgrw (n, kl, ku, ncols, ab, ldab, afb, ldafb)
 SLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix.
subroutine sorgbr (vect, m, n, k, a, lda, tau, work, lwork, info)
 SORGBR

Detailed Description

This is the group of real computational functions for GB matrices

Function Documentation

◆ sgbbrd()

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

SGBBRD

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

Purpose:
!>
!> SGBBRD 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 REAL 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 REAL array, dimension (min(M,N))
!>          The diagonal elements of the bidiagonal matrix B.
!>
[out]E
!>          E is REAL array, dimension (min(M,N)-1)
!>          The superdiagonal elements of the bidiagonal matrix B.
!>
[out]Q
!>          Q is REAL 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 REAL 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 REAL 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 REAL 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 sgbbrd.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 REAL AB( LDAB, * ), C( LDC, * ), D( * ), E( * ),
198 $ PT( LDPT, * ), Q( LDQ, * ), WORK( * )
199* ..
200*
201* =====================================================================
202*
203* .. Parameters ..
204 REAL ZERO, ONE
205 parameter( zero = 0.0e+0, one = 1.0e+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 REAL RA, RB, RC, RS
212* ..
213* .. External Subroutines ..
214 EXTERNAL slargv, slartg, slartv, slaset, srot, 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( 'SGBBRD', -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 slaset( 'Full', m, m, zero, one, q, ldq )
264 IF( wantpt )
265 $ CALL slaset( '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 slargv( 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 slartv( 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 slartg( 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 srot( 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 srot( 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 srot( 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 slargv( 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 slartv( 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 slartg( 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 srot( 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 srot( 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 slartg( 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 srot( m, q( 1, i ), 1, q( 1, i+1 ), 1, rc, rs )
489 IF( wantc )
490 $ CALL srot( 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 slartg( 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 srot( 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 SGBBRD
543*
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:113
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
slargv
subroutine slargv(n, x, incx, y, incy, c, incc)
SLARGV generates a vector of plane rotations with real cosines and real sines.
Definition slargv.f:104
slartv
subroutine slartv(n, x, incx, y, incy, c, s, incc)
SLARTV applies a vector of plane rotations with real cosines and real sines to the elements of a pair...
Definition slartv.f:108
srot
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ sgbcon()

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

SGBCON

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

Purpose:
!>
!> SGBCON 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 SGBTRF.
!>
!> 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 REAL array, dimension (LDAB,N)
!>          Details of the LU factorization of the band matrix A, as
!>          computed by SGBTRF.  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 REAL
!>          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 REAL
!>          The reciprocal of the condition number of the matrix A,
!>          computed as RCOND = 1/(norm(A) * norm(inv(A))).
!>
[out]WORK
!>          WORK is REAL 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 sgbcon.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 REAL ANORM, RCOND
155* ..
156* .. Array Arguments ..
157 INTEGER IPIV( * ), IWORK( * )
158 REAL AB( LDAB, * ), WORK( * )
159* ..
160*
161* =====================================================================
162*
163* .. Parameters ..
164 REAL ONE, ZERO
165 parameter( one = 1.0e+0, zero = 0.0e+0 )
166* ..
167* .. Local Scalars ..
168 LOGICAL LNOTI, ONENRM
169 CHARACTER NORMIN
170 INTEGER IX, J, JP, KASE, KASE1, KD, LM
171 REAL AINVNM, SCALE, SMLNUM, T
172* ..
173* .. Local Arrays ..
174 INTEGER ISAVE( 3 )
175* ..
176* .. External Functions ..
177 LOGICAL LSAME
178 INTEGER ISAMAX
179 REAL SDOT, SLAMCH
180 EXTERNAL lsame, isamax, sdot, slamch
181* ..
182* .. External Subroutines ..
183 EXTERNAL saxpy, slacn2, slatbs, srscl, 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( 'SGBCON', -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 = slamch( '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 slacn2( 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 saxpy( 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 slatbs( '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 slatbs( '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 ) - sdot( 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 = isamax( n, work, 1 )
291 IF( scale.LT.abs( work( ix ) )*smlnum .OR. scale.EQ.zero )
292 $ GO TO 40
293 CALL srscl( 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 SGBCON
307*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
slatbs
subroutine slatbs(uplo, trans, diag, normin, n, kd, ab, ldab, x, scale, cnorm, info)
SLATBS solves a triangular banded system of equations.
Definition slatbs.f:242
srscl
subroutine srscl(n, sa, sx, incx)
SRSCL multiplies a vector by the reciprocal of a real scalar.
Definition srscl.f:84
slacn2
subroutine slacn2(n, v, x, isgn, est, kase, isave)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition slacn2.f:136
sdot
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68

◆ sgbequ()

subroutine sgbequ ( integer m,
integer n,
integer kl,
integer ku,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( * ) r,
real, dimension( * ) c,
real rowcnd,
real colcnd,
real amax,
integer info )

SGBEQU

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

Purpose:
!>
!> SGBEQU 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 REAL 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 REAL array, dimension (M)
!>          If INFO = 0, or INFO > M, R contains the row scale factors
!>          for A.
!>
[out]C
!>          C is REAL array, dimension (N)
!>          If INFO = 0, C contains the column scale factors for A.
!>
[out]ROWCND
!>          ROWCND is REAL
!>          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 REAL
!>          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 REAL
!>          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 sgbequ.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 REAL AMAX, COLCND, ROWCND
161* ..
162* .. Array Arguments ..
163 REAL AB( LDAB, * ), C( * ), R( * )
164* ..
165*
166* =====================================================================
167*
168* .. Parameters ..
169 REAL ONE, ZERO
170 parameter( one = 1.0e+0, zero = 0.0e+0 )
171* ..
172* .. Local Scalars ..
173 INTEGER I, J, KD
174 REAL BIGNUM, RCMAX, RCMIN, SMLNUM
175* ..
176* .. External Functions ..
177 REAL SLAMCH
178 EXTERNAL slamch
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( 'SGBEQU', -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 = slamch( '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 SGBEQU
320*

◆ sgbequb()

subroutine sgbequb ( integer m,
integer n,
integer kl,
integer ku,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( * ) r,
real, dimension( * ) c,
real rowcnd,
real colcnd,
real amax,
integer info )

SGBEQUB

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

Purpose:
!>
!> SGBEQUB 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 SGEEQU 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 REAL 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 REAL array, dimension (M)
!>          If INFO = 0 or INFO > M, R contains the row scale factors
!>          for A.
!>
[out]C
!>          C is REAL array, dimension (N)
!>          If INFO = 0,  C contains the column scale factors for A.
!>
[out]ROWCND
!>          ROWCND is REAL
!>          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 REAL
!>          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 REAL
!>          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 sgbequb.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 REAL AMAX, COLCND, ROWCND
168* ..
169* .. Array Arguments ..
170 REAL AB( LDAB, * ), C( * ), R( * )
171* ..
172*
173* =====================================================================
174*
175* .. Parameters ..
176 REAL ONE, ZERO
177 parameter( one = 1.0e+0, zero = 0.0e+0 )
178* ..
179* .. Local Scalars ..
180 INTEGER I, J, KD
181 REAL BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX
182* ..
183* .. External Functions ..
184 REAL SLAMCH
185 EXTERNAL slamch
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( 'SGBEQUB', -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 = slamch( 'S' )
226 bignum = one / smlnum
227 radix = slamch( '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 SGBEQUB
336*

◆ sgbrfs()

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

SGBRFS

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

Purpose:
!>
!> SGBRFS 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 REAL 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 REAL array, dimension (LDAFB,N)
!>          Details of the LU factorization of the band matrix A, as
!>          computed by SGBTRF.  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 SGBTRF; for 1<=i<=N, row i of the
!>          matrix was interchanged with row IPIV(i).
!>
[in]B
!>          B is REAL 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 REAL array, dimension (LDX,NRHS)
!>          On entry, the solution matrix X, as computed by SGBTRS.
!>          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 REAL 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 REAL 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 REAL 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 sgbrfs.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 REAL 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 REAL ZERO
226 parameter( zero = 0.0e+0 )
227 REAL ONE
228 parameter( one = 1.0e+0 )
229 REAL TWO
230 parameter( two = 2.0e+0 )
231 REAL THREE
232 parameter( three = 3.0e+0 )
233* ..
234* .. Local Scalars ..
235 LOGICAL NOTRAN
236 CHARACTER TRANST
237 INTEGER COUNT, I, J, K, KASE, KK, NZ
238 REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
239* ..
240* .. Local Arrays ..
241 INTEGER ISAVE( 3 )
242* ..
243* .. External Subroutines ..
244 EXTERNAL saxpy, scopy, sgbmv, sgbtrs, slacn2, xerbla
245* ..
246* .. Intrinsic Functions ..
247 INTRINSIC abs, max, min
248* ..
249* .. External Functions ..
250 LOGICAL LSAME
251 REAL SLAMCH
252 EXTERNAL lsame, slamch
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( 'SGBRFS', -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 = slamch( 'Epsilon' )
305 safmin = slamch( '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 scopy( n, b( 1, j ), 1, work( n+1 ), 1 )
323 CALL sgbmv( 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 sgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv,
382 $ work( n+1 ), n, info )
383 CALL saxpy( 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 SLACN2 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 slacn2( 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 sgbtrs( 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 sgbtrs( 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 SGBRFS
460*
sgbtrs
subroutine sgbtrs(trans, n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
SGBTRS
Definition sgbtrs.f:138
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sgbmv
subroutine sgbmv(trans, m, n, kl, ku, alpha, a, lda, x, incx, beta, y, incy)
SGBMV
Definition sgbmv.f:185

◆ sgbrfsx()

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

SGBRFSX

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

Purpose:
!>
!>    SGBRFSX 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 REAL 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 REAL array, dimension (LDAFB,N)
!>     Details of the LU factorization of the band matrix A, as
!>     computed by SGBTRF.  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 SGETRF; for 1<=i<=N, row i of the
!>     matrix was interchanged with row IPIV(i).
!>
[in,out]R
!>          R is REAL 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 REAL 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 REAL 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 REAL array, dimension (LDX,NRHS)
!>     On entry, the solution matrix X, as computed by SGETRS.
!>     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 REAL
!>     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 REAL 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 REAL 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.
!>
!>     See Lapack Working Note 165 for further details and extra
!>     cautions.
!>
[out]ERR_BNDS_COMP
!>          ERR_BNDS_COMP is REAL 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.
!>
!>     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 REAL 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.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 REAL 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 sgbrfsx.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 REAL RCOND
450* ..
451* .. Array Arguments ..
452 INTEGER IPIV( * ), IWORK( * )
453 REAL AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
454 $ X( LDX , * ),WORK( * )
455 REAL R( * ), C( * ), PARAMS( * ), BERR( * ),
456 $ ERR_BNDS_NORM( NRHS, * ),
457 $ ERR_BNDS_COMP( NRHS, * )
458* ..
459*
460* ==================================================================
461*
462* .. Parameters ..
463 REAL ZERO, ONE
464 parameter( zero = 0.0e+0, one = 1.0e+0 )
465 REAL ITREF_DEFAULT, ITHRESH_DEFAULT,
466 $ COMPONENTWISE_DEFAULT
467 REAL RTHRESH_DEFAULT, DZTHRESH_DEFAULT
468 parameter( itref_default = 1.0 )
469 parameter( ithresh_default = 10.0 )
470 parameter( componentwise_default = 1.0 )
471 parameter( rthresh_default = 0.5 )
472 parameter( dzthresh_default = 0.25 )
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 REAL ANORM, RCOND_TMP
489 REAL ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
490 LOGICAL IGNORE_CWISE
491 INTEGER ITHRESH
492 REAL RTHRESH, UNSTABLE_THRESH
493* ..
494* .. External Subroutines ..
495 EXTERNAL xerbla, sgbcon
496 EXTERNAL sla_gbrfsx_extended
497* ..
498* .. Intrinsic Functions ..
499 INTRINSIC max, sqrt
500* ..
501* .. External Functions ..
502 EXTERNAL lsame, ilatrans, ilaprec
503 EXTERNAL slamch, slangb, sla_gbrcond
504 REAL SLAMCH, SLANGB, SLA_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.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 = real( n ) * slamch( 'Epsilon' )
526 ithresh = int( ithresh_default )
527 rthresh = rthresh_default
528 unstable_thresh = dzthresh_default
529 ignore_cwise = componentwise_default .EQ. 0.0
530*
531 IF ( nparams.GE.la_linrx_ithresh_i ) THEN
532 IF ( params( la_linrx_ithresh_i ).LT.0.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.0 ) THEN
540 IF ( ignore_cwise ) THEN
541 params( la_linrx_cwise_i ) = 0.0
542 ELSE
543 params( la_linrx_cwise_i ) = 1.0
544 END IF
545 ELSE
546 ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.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( 'SGBRFSX', -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.0
594 DO j = 1, nrhs
595 berr( j ) = 0.0
596 IF ( n_err_bnds .GE. 1 ) THEN
597 err_bnds_norm( j, la_linrx_trust_i ) = 1.0
598 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
599 END IF
600 IF ( n_err_bnds .GE. 2 ) THEN
601 err_bnds_norm( j, la_linrx_err_i ) = 0.0
602 err_bnds_comp( j, la_linrx_err_i ) = 0.0
603 END IF
604 IF ( n_err_bnds .GE. 3 ) THEN
605 err_bnds_norm( j, la_linrx_rcond_i ) = 1.0
606 err_bnds_comp( j, la_linrx_rcond_i ) = 1.0
607 END IF
608 END DO
609 RETURN
610 END IF
611*
612* Default to failure.
613*
614 rcond = 0.0
615 DO j = 1, nrhs
616 berr( j ) = 1.0
617 IF ( n_err_bnds .GE. 1 ) THEN
618 err_bnds_norm( j, la_linrx_trust_i ) = 1.0
619 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
620 END IF
621 IF ( n_err_bnds .GE. 2 ) THEN
622 err_bnds_norm( j, la_linrx_err_i ) = 1.0
623 err_bnds_comp( j, la_linrx_err_i ) = 1.0
624 END IF
625 IF ( n_err_bnds .GE. 3 ) THEN
626 err_bnds_norm( j, la_linrx_rcond_i ) = 0.0
627 err_bnds_comp( j, la_linrx_rcond_i ) = 0.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 = slangb( norm, n, kl, ku, ab, ldab, work )
640 CALL sgbcon( 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( 'D' )
648
649 IF ( notran ) THEN
650 CALL sla_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 sla_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.0, sqrt( real( n ) ) ) * slamch( '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 = sla_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 = sla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
676 $ ldafb, ipiv, -1, r, info, work, iwork )
677 ELSE
678 rcond_tmp = sla_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.0 )
687 $ err_bnds_norm( j, la_linrx_err_i ) = 1.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.0
693 err_bnds_norm( j, la_linrx_trust_i ) = 0.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.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( slamch( 'Epsilon' ) )
721 DO j = 1, nrhs
722 IF ( err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
723 $ THEN
724 rcond_tmp = sla_gbrcond( trans, n, kl, ku, ab, ldab, afb,
725 $ ldafb, ipiv, 1, x( 1, j ), info, work, iwork )
726 ELSE
727 rcond_tmp = 0.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.0 )
734 $ err_bnds_comp( j, la_linrx_err_i ) = 1.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.0
740 err_bnds_comp( j, la_linrx_trust_i ) = 0.0
741 IF ( params( la_linrx_cwise_i ) .EQ. 1.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.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 SGBRFSX
761*
ilaprec
integer function ilaprec(prec)
ILAPREC
Definition ilaprec.f:58
ilatrans
integer function ilatrans(trans)
ILATRANS
Definition ilatrans.f:58
slangb
real function slangb(norm, n, kl, ku, ab, ldab, work)
SLANGB returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slangb.f:124
sla_gbrfsx_extended
subroutine sla_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)
SLA_GBRFSX_EXTENDED improves the computed solution to a system of linear equations for general banded...
Definition sla_gbrfsx_extended.f:410
sgbcon
subroutine sgbcon(norm, n, kl, ku, ab, ldab, ipiv, anorm, rcond, work, iwork, info)
SGBCON
Definition sgbcon.f:146
sla_gbrcond
real function sla_gbrcond(trans, n, kl, ku, ab, ldab, afb, ldafb, ipiv, cmode, c, info, work, iwork)
SLA_GBRCOND estimates the Skeel condition number for a general banded matrix.
Definition sla_gbrcond.f:168

◆ sgbtf2()

subroutine sgbtf2 ( integer m,
integer n,
integer kl,
integer ku,
real, dimension( ldab, * ) ab,
integer ldab,
integer, dimension( * ) ipiv,
integer info )

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

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

Purpose:
!>
!> SGBTF2 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 REAL 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 sgbtf2.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 REAL AB( LDAB, * )
156* ..
157*
158* =====================================================================
159*
160* .. Parameters ..
161 REAL ONE, ZERO
162 parameter( one = 1.0e+0, zero = 0.0e+0 )
163* ..
164* .. Local Scalars ..
165 INTEGER I, J, JP, JU, KM, KV
166* ..
167* .. External Functions ..
168 INTEGER ISAMAX
169 EXTERNAL isamax
170* ..
171* .. External Subroutines ..
172 EXTERNAL sger, sscal, sswap, 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( 'SGBTF2', -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 = isamax( 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 sswap( 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 sscal( 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 sger( 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 SGBTF2
273*
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
sger
subroutine sger(m, n, alpha, x, incx, y, incy, a, lda)
SGER
Definition sger.f:130

◆ sgbtrf()

subroutine sgbtrf ( integer m,
integer n,
integer kl,
integer ku,
real, dimension( ldab, * ) ab,
integer ldab,
integer, dimension( * ) ipiv,
integer info )

SGBTRF

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

Purpose:
!>
!> SGBTRF 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 REAL 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 sgbtrf.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 REAL AB( LDAB, * )
155* ..
156*
157* =====================================================================
158*
159* .. Parameters ..
160 REAL ONE, ZERO
161 parameter( one = 1.0e+0, zero = 0.0e+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 REAL TEMP
169* ..
170* .. Local Arrays ..
171 REAL WORK13( LDWORK, NBMAX ),
172 $ WORK31( LDWORK, NBMAX )
173* ..
174* .. External Functions ..
175 INTEGER ILAENV, ISAMAX
176 EXTERNAL ilaenv, isamax
177* ..
178* .. External Subroutines ..
179 EXTERNAL scopy, sgbtf2, sgemm, sger, slaswp, sscal,
180 $ sswap, strsm, 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( 'SGBTRF', -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, 'SGBTRF', ' ', 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 sgbtf2( 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 = isamax( 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 sswap( 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 sswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
319 $ work31( jp+jj-j-kl, 1 ), ldwork )
320 CALL sswap( 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 sscal( 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 sger( 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 scopy( 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 SLASWP to apply the row interchanges to A12, A22, and
363* A32.
364*
365 CALL slaswp( 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 strsm( '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 sgemm( '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 sgemm( '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 strsm( '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 sgemm( '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 sgemm( '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 sswap( 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 sswap( 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 scopy( 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 SGBTRF
512*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
sgbtf2
subroutine sgbtf2(m, n, kl, ku, ab, ldab, ipiv, info)
SGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algo...
Definition sgbtf2.f:145
slaswp
subroutine slaswp(n, a, lda, k1, k2, ipiv, incx)
SLASWP performs a series of row interchanges on a general rectangular matrix.
Definition slaswp.f:115
strsm
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187

◆ sgbtrs()

subroutine sgbtrs ( character trans,
integer n,
integer kl,
integer ku,
integer nrhs,
real, dimension( ldab, * ) ab,
integer ldab,
integer, dimension( * ) ipiv,
real, dimension( ldb, * ) b,
integer ldb,
integer info )

SGBTRS

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

Purpose:
!>
!> SGBTRS 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 SGBTRF.
!>
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 REAL array, dimension (LDAB,N)
!>          Details of the LU factorization of the band matrix A, as
!>          computed by SGBTRF.  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 REAL 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 sgbtrs.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 REAL AB( LDAB, * ), B( LDB, * )
150* ..
151*
152* =====================================================================
153*
154* .. Parameters ..
155 REAL ONE
156 parameter( one = 1.0e+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 sgemv, sger, sswap, stbsv, 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( 'SGBTRS', -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 sswap( nrhs, b( l, 1 ), ldb, b( j, 1 ), ldb )
224 CALL sger( 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 stbsv( '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 stbsv( '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 sgemv( '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 sswap( nrhs, b( l, 1 ), ldb, b( j, 1 ), ldb )
259 40 CONTINUE
260 END IF
261 END IF
262 RETURN
263*
264* End of SGBTRS
265*
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:156
stbsv
subroutine stbsv(uplo, trans, diag, n, k, a, lda, x, incx)
STBSV
Definition stbsv.f:189

◆ sggbak()

subroutine sggbak ( character job,
character side,
integer n,
integer ilo,
integer ihi,
real, dimension( * ) lscale,
real, dimension( * ) rscale,
integer m,
real, dimension( ldv, * ) v,
integer ldv,
integer info )

SGGBAK

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

Purpose:
!>
!> SGGBAK 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
!> SGGBAL.
!>
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 SGGBAL.
!>
[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 SGGBAL.
!>          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
!>
[in]LSCALE
!>          LSCALE is REAL array, dimension (N)
!>          Details of the permutations and/or scaling factors applied
!>          to the left side of A and B, as returned by SGGBAL.
!>
[in]RSCALE
!>          RSCALE is REAL array, dimension (N)
!>          Details of the permutations and/or scaling factors applied
!>          to the right side of A and B, as returned by SGGBAL.
!>
[in]M
!>          M is INTEGER
!>          The number of columns of the matrix V.  M >= 0.
!>
[in,out]V
!>          V is REAL array, dimension (LDV,M)
!>          On entry, the matrix of right or left eigenvectors to be
!>          transformed, as returned by STGEVC.
!>          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 sggbak.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 REAL 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 sscal, sswap, xerbla
172* ..
173* .. Intrinsic Functions ..
174 INTRINSIC max
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( 'SGGBAK', -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 sscal( 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 sscal( 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 = rscale( i )
256 IF( k.EQ.i )
257 $ GO TO 40
258 CALL sswap( 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 = rscale( i )
266 IF( k.EQ.i )
267 $ GO TO 60
268 CALL sswap( 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 = lscale( i )
280 IF( k.EQ.i )
281 $ GO TO 80
282 CALL sswap( 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 = lscale( i )
290 IF( k.EQ.i )
291 $ GO TO 100
292 CALL sswap( 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 SGGBAK
302*

◆ sggbal()

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

SGGBAL

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

Purpose:
!>
!> SGGBAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 sggbal.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 REAL A( LDA, * ), B( LDB, * ), LSCALE( * ),
188 $ RSCALE( * ), WORK( * )
189* ..
190*
191* =====================================================================
192*
193* .. Parameters ..
194 REAL ZERO, HALF, ONE
195 parameter( zero = 0.0e+0, half = 0.5e+0, one = 1.0e+0 )
196 REAL THREE, SCLFAC
197 parameter( three = 3.0e+0, sclfac = 1.0e+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 REAL 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 ISAMAX
210 REAL SDOT, SLAMCH
211 EXTERNAL lsame, isamax, sdot, slamch
212* ..
213* .. External Subroutines ..
214 EXTERNAL saxpy, sscal, sswap, xerbla
215* ..
216* .. Intrinsic Functions ..
217 INTRINSIC abs, int, log10, max, min, real, 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( 'SGGBAL', -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 sswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
344 CALL sswap( 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 sswap( l, a( 1, j ), 1, a( 1, m ), 1 )
353 CALL sswap( 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 / real( 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 = sdot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
420 $ sdot( 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 sscal( nr, beta, work( ilo ), 1 )
438 CALL sscal( nr, beta, work( ilo+n ), 1 )
439*
440 CALL saxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
441 CALL saxpy( 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 ) = real( 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 ) = real( kount )*work( j ) + sum
482 330 CONTINUE
483*
484 sum = sdot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
485 $ sdot( 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 saxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ), 1 )
505 CALL saxpy( 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 = slamch( '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 = isamax( n-ilo+1, a( i, ilo ), lda )
521 rab = abs( a( i, irab+ilo-1 ) )
522 irab = isamax( 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 = lscale( i ) + sign( half, lscale( i ) )
526 ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
527 lscale( i ) = sclfac**ir
528 icab = isamax( ihi, a( 1, i ), 1 )
529 cab = abs( a( icab, i ) )
530 icab = isamax( ihi, b( 1, i ), 1 )
531 cab = max( cab, abs( b( icab, i ) ) )
532 lcab = int( log10( cab+sfmin ) / basl+one )
533 jc = 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 sscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
542 CALL sscal( 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 sscal( ihi, rscale( j ), a( 1, j ), 1 )
549 CALL sscal( ihi, rscale( j ), b( 1, j ), 1 )
550 380 CONTINUE
551*
552 RETURN
553*
554* End of SGGBAL
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

◆ sla_gbamv()

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

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

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

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

◆ sla_gbrcond()

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

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

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

Purpose:
!>
!>    SLA_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 REAL 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 REAL array, dimension (LDAFB,N)
!>     Details of the LU factorization of the band matrix A, as
!>     computed by SGBTRF.  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 SGBTRF; 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 REAL 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 REAL 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 166 of file sla_gbrcond.f.

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

◆ sla_gbrfsx_extended()

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

SLA_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 SLA_GBRFSX_EXTENDED + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SLA_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 SGBRFSX 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 REAL array, dimension (LDAB,N)
!>     On entry, the N-by-N matrix AB.
!>
[in]LDAB
!>          LDAB is INTEGER
!>     The leading dimension of the array AB.  LDAB >= max(1,N).
!>
[in]AFB
!>          AFB is REAL array, dimension (LDAFB,N)
!>     The factors L and U from the factorization
!>     A = P*L*U as computed by SGBTRF.
!>
[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 SGBTRF; 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 REAL 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 REAL 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 REAL array, dimension (LDY,NRHS)
!>     On entry, the solution matrix X, as computed by SGBTRS.
!>     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 REAL 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 SLA_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 REAL 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 REAL 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 REAL array, dimension (N)
!>     Workspace to hold the intermediate residual.
!>
[in]AYB
!>          AYB is REAL array, dimension (N)
!>     Workspace. This can be the same workspace passed for Y_TAIL.
!>
[in]DY
!>          DY is REAL array, dimension (N)
!>     Workspace to hold the intermediate solution.
!>
[in]Y_TAIL
!>          Y_TAIL is REAL array, dimension (N)
!>     Workspace to hold the trailing bits of the intermediate solution.
!>
[in]RCOND
!>          RCOND is REAL
!>     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 REAL
!>     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 REAL
!>     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 SGBTRS had an illegal
!>             value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 403 of file sla_gbrfsx_extended.f.

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

◆ sla_gbrpvgrw()

real function sla_gbrpvgrw ( integer n,
integer kl,
integer ku,
integer ncols,
real, dimension( ldab, * ) ab,
integer ldab,
real, dimension( ldafb, * ) afb,
integer ldafb )

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

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

Purpose:
!>
!> SLA_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 REAL 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 REAL array, dimension (LDAFB,N)
!>     Details of the LU factorization of the band matrix A, as
!>     computed by SGBTRF.  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 sla_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 REAL AB( LDAB, * ), AFB( LDAFB, * )
127* ..
128*
129* =====================================================================
130*
131* .. Local Scalars ..
132 INTEGER I, J, KD
133 REAL AMAX, UMAX, RPVGRW
134* ..
135* .. Intrinsic Functions ..
136 INTRINSIC abs, max, min
137* ..
138* .. Executable Statements ..
139*
140 rpvgrw = 1.0
141
142 kd = ku + 1
143 DO j = 1, ncols
144 amax = 0.0
145 umax = 0.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.0 ) THEN
153 rpvgrw = min( amax / umax, rpvgrw )
154 END IF
155 END DO
156 sla_gbrpvgrw = rpvgrw
157*
158* End of SLA_GBRPVGRW
159*
sla_gbrpvgrw
real function sla_gbrpvgrw(n, kl, ku, ncols, ab, ldab, afb, ldafb)
SLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix.
Definition sla_gbrpvgrw.f:117

◆ sorgbr()

subroutine sorgbr ( character vect,
integer m,
integer n,
integer k,
real, dimension( lda, * ) a,
integer lda,
real, dimension( * ) tau,
real, dimension( * ) work,
integer lwork,
integer info )

SORGBR

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

Purpose:
!>
!> SORGBR generates one of the real orthogonal matrices Q or P**T
!> determined by SGEBRD 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 SORGBR returns the first n
!> columns of Q, where m >= n >= k;
!> if m < k, Q = H(1) H(2) . . . H(m-1) and SORGBR 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 SORGBR 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 SORGBR 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 SGEBRD:
!>          = '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 SGEBRD.
!>          If VECT = 'P', the number of rows in the original K-by-N
!>          matrix reduced by SGEBRD.
!>          K >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA,N)
!>          On entry, the vectors which define the elementary reflectors,
!>          as returned by SGEBRD.
!>          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 REAL 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 SGEBRD in its array argument TAUQ or TAUP.
!>
[out]WORK
!>          WORK is REAL 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 sorgbr.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 REAL A( LDA, * ), TAU( * ), WORK( * )
168* ..
169*
170* =====================================================================
171*
172* .. Parameters ..
173 REAL ZERO, ONE
174 parameter( zero = 0.0e+0, one = 1.0e+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 sorglq, sorgqr, 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 sorgqr( m, n, k, a, lda, tau, work, -1, iinfo )
219 ELSE
220 IF( m.GT.1 ) THEN
221 CALL sorgqr( 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 sorglq( m, n, k, a, lda, tau, work, -1, iinfo )
228 ELSE
229 IF( n.GT.1 ) THEN
230 CALL sorglq( 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( 'SORGBR', -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 SGEBRD 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 sorgqr( 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 sorgqr( 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 SGEBRD 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 sorglq( 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 sorglq( 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 SORGBR
333*
sorglq
subroutine sorglq(m, n, k, a, lda, tau, work, lwork, info)
SORGLQ
Definition sorglq.f:127
sorgqr
subroutine sorgqr(m, n, k, a, lda, tau, work, lwork, info)
SORGQR
Definition sorgqr.f:128