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

Functions

subroutine cgbbrd (vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, c, ldc, work, rwork, info)
 CGBBRD
subroutine cgbcon (norm, n, kl, ku, ab, ldab, ipiv, anorm, rcond, work, rwork, info)
 CGBCON
subroutine cgbequ (m, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd, amax, info)
 CGBEQU
subroutine cgbequb (m, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd, amax, info)
 CGBEQUB
subroutine cgbrfs (trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, b, ldb, x, ldx, ferr, berr, work, rwork, info)
 CGBRFS
subroutine cgbrfsx (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, rwork, info)
 CGBRFSX
subroutine cgbtf2 (m, n, kl, ku, ab, ldab, ipiv, info)
 CGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algorithm.
subroutine cgbtrf (m, n, kl, ku, ab, ldab, ipiv, info)
 CGBTRF
subroutine cgbtrs (trans, n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
 CGBTRS
subroutine cggbak (job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
 CGGBAK
subroutine cggbal (job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
 CGGBAL
subroutine cla_gbamv (trans, m, n, kl, ku, alpha, ab, ldab, x, incx, beta, y, incy)
 CLA_GBAMV performs a matrix-vector operation to calculate error bounds.
real function cla_gbrcond_c (trans, n, kl, ku, ab, ldab, afb, ldafb, ipiv, c, capply, info, work, rwork)
 CLA_GBRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general banded matrices.
real function cla_gbrcond_x (trans, n, kl, ku, ab, ldab, afb, ldafb, ipiv, x, info, work, rwork)
 CLA_GBRCOND_X computes the infinity norm condition number of op(A)*diag(x) for general banded matrices.
subroutine cla_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)
 CLA_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 cla_gbrpvgrw (n, kl, ku, ncols, ab, ldab, afb, ldafb)
 CLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix.
subroutine cungbr (vect, m, n, k, a, lda, tau, work, lwork, info)
 CUNGBR

Detailed Description

This is the group of complex computational functions for GB matrices

Function Documentation

◆ cgbbrd()

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

CGBBRD

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

Purpose:
!>
!> CGBBRD reduces a complex general m-by-n band matrix A to real upper
!> bidiagonal form B by a unitary transformation: Q**H * A * P = B.
!>
!> The routine computes B, and optionally forms Q or P**H, or computes
!> Q**H*C for a given matrix C.
!>
Parameters
[in]VECT
!>          VECT is CHARACTER*1
!>          Specifies whether or not the matrices Q and P**H are to be
!>          formed.
!>          = 'N': do not form Q or P**H;
!>          = 'Q': form Q only;
!>          = 'P': form P**H 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 COMPLEX 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 COMPLEX array, dimension (LDQ,M)
!>          If VECT = 'Q' or 'B', the m-by-m unitary 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 COMPLEX array, dimension (LDPT,N)
!>          If VECT = 'P' or 'B', the n-by-n unitary 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 COMPLEX array, dimension (LDC,NCC)
!>          On entry, an m-by-ncc matrix C.
!>          On exit, C is overwritten by Q**H*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 COMPLEX array, dimension (max(M,N))
!>
[out]RWORK
!>          RWORK is REAL array, dimension (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 191 of file cgbbrd.f.

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

◆ cgbcon()

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

CGBCON

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

Purpose:
!>
!> CGBCON estimates the reciprocal of the condition number of a complex
!> general band matrix A, in either the 1-norm or the infinity-norm,
!> using the LU factorization computed by CGBTRF.
!>
!> 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 COMPLEX array, dimension (LDAB,N)
!>          Details of the LU factorization of the band matrix A, as
!>          computed by CGBTRF.  U is stored as an upper triangular band
!>          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
!>          the multipliers used during the factorization are stored in
!>          rows KL+KU+2 to 2*KL+KU+1.
!>
[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 COMPLEX array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is REAL 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 145 of file cgbcon.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 NORM
154 INTEGER INFO, KL, KU, LDAB, N
155 REAL ANORM, RCOND
156* ..
157* .. Array Arguments ..
158 INTEGER IPIV( * )
159 REAL RWORK( * )
160 COMPLEX AB( LDAB, * ), WORK( * )
161* ..
162*
163* =====================================================================
164*
165* .. Parameters ..
166 REAL ONE, ZERO
167 parameter( one = 1.0e+0, zero = 0.0e+0 )
168* ..
169* .. Local Scalars ..
170 LOGICAL LNOTI, ONENRM
171 CHARACTER NORMIN
172 INTEGER IX, J, JP, KASE, KASE1, KD, LM
173 REAL AINVNM, SCALE, SMLNUM
174 COMPLEX T, ZDUM
175* ..
176* .. Local Arrays ..
177 INTEGER ISAVE( 3 )
178* ..
179* .. External Functions ..
180 LOGICAL LSAME
181 INTEGER ICAMAX
182 REAL SLAMCH
183 COMPLEX CDOTC
184 EXTERNAL lsame, icamax, slamch, cdotc
185* ..
186* .. External Subroutines ..
187 EXTERNAL caxpy, clacn2, clatbs, csrscl, xerbla
188* ..
189* .. Intrinsic Functions ..
190 INTRINSIC abs, aimag, min, real
191* ..
192* .. Statement Functions ..
193 REAL CABS1
194* ..
195* .. Statement Function definitions ..
196 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
197* ..
198* .. Executable Statements ..
199*
200* Test the input parameters.
201*
202 info = 0
203 onenrm = norm.EQ.'1' .OR. lsame( norm, 'O' )
204 IF( .NOT.onenrm .AND. .NOT.lsame( norm, 'I' ) ) THEN
205 info = -1
206 ELSE IF( n.LT.0 ) THEN
207 info = -2
208 ELSE IF( kl.LT.0 ) THEN
209 info = -3
210 ELSE IF( ku.LT.0 ) THEN
211 info = -4
212 ELSE IF( ldab.LT.2*kl+ku+1 ) THEN
213 info = -6
214 ELSE IF( anorm.LT.zero ) THEN
215 info = -8
216 END IF
217 IF( info.NE.0 ) THEN
218 CALL xerbla( 'CGBCON', -info )
219 RETURN
220 END IF
221*
222* Quick return if possible
223*
224 rcond = zero
225 IF( n.EQ.0 ) THEN
226 rcond = one
227 RETURN
228 ELSE IF( anorm.EQ.zero ) THEN
229 RETURN
230 END IF
231*
232 smlnum = slamch( 'Safe minimum' )
233*
234* Estimate the norm of inv(A).
235*
236 ainvnm = zero
237 normin = 'N'
238 IF( onenrm ) THEN
239 kase1 = 1
240 ELSE
241 kase1 = 2
242 END IF
243 kd = kl + ku + 1
244 lnoti = kl.GT.0
245 kase = 0
246 10 CONTINUE
247 CALL clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
248 IF( kase.NE.0 ) THEN
249 IF( kase.EQ.kase1 ) THEN
250*
251* Multiply by inv(L).
252*
253 IF( lnoti ) THEN
254 DO 20 j = 1, n - 1
255 lm = min( kl, n-j )
256 jp = ipiv( j )
257 t = work( jp )
258 IF( jp.NE.j ) THEN
259 work( jp ) = work( j )
260 work( j ) = t
261 END IF
262 CALL caxpy( lm, -t, ab( kd+1, j ), 1, work( j+1 ), 1 )
263 20 CONTINUE
264 END IF
265*
266* Multiply by inv(U).
267*
268 CALL clatbs( 'Upper', 'No transpose', 'Non-unit', normin, n,
269 $ kl+ku, ab, ldab, work, scale, rwork, info )
270 ELSE
271*
272* Multiply by inv(U**H).
273*
274 CALL clatbs( 'Upper', 'Conjugate transpose', 'Non-unit',
275 $ normin, n, kl+ku, ab, ldab, work, scale, rwork,
276 $ info )
277*
278* Multiply by inv(L**H).
279*
280 IF( lnoti ) THEN
281 DO 30 j = n - 1, 1, -1
282 lm = min( kl, n-j )
283 work( j ) = work( j ) - cdotc( lm, ab( kd+1, j ), 1,
284 $ work( j+1 ), 1 )
285 jp = ipiv( j )
286 IF( jp.NE.j ) THEN
287 t = work( jp )
288 work( jp ) = work( j )
289 work( j ) = t
290 END IF
291 30 CONTINUE
292 END IF
293 END IF
294*
295* Divide X by 1/SCALE if doing so will not cause overflow.
296*
297 normin = 'Y'
298 IF( scale.NE.one ) THEN
299 ix = icamax( n, work, 1 )
300 IF( scale.LT.cabs1( work( ix ) )*smlnum .OR. scale.EQ.zero )
301 $ GO TO 40
302 CALL csrscl( n, scale, work, 1 )
303 END IF
304 GO TO 10
305 END IF
306*
307* Compute the estimate of the reciprocal condition number.
308*
309 IF( ainvnm.NE.zero )
310 $ rcond = ( one / ainvnm ) / anorm
311*
312 40 CONTINUE
313 RETURN
314*
315* End of CGBCON
316*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
icamax
integer function icamax(n, cx, incx)
ICAMAX
Definition icamax.f:71
clatbs
subroutine clatbs(uplo, trans, diag, normin, n, kd, ab, ldab, x, scale, cnorm, info)
CLATBS solves a triangular banded system of equations.
Definition clatbs.f:243
csrscl
subroutine csrscl(n, sa, sx, incx)
CSRSCL multiplies a vector by the reciprocal of a real scalar.
Definition csrscl.f:84
clacn2
subroutine clacn2(n, v, x, est, kase, isave)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition clacn2.f:133
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
cdotc
complex function cdotc(n, cx, incx, cy, incy)
CDOTC
Definition cdotc.f:83
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68

◆ cgbequ()

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

CGBEQU

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

Purpose:
!>
!> CGBEQU 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 COMPLEX 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 152 of file cgbequ.f.

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

◆ cgbequb()

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

CGBEQUB

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

Purpose:
!>
!> CGBEQUB 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 CGEEQU 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 COMPLEX 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 159 of file cgbequb.f.

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

◆ cgbrfs()

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

CGBRFS

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

Purpose:
!>
!> CGBRFS 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)
!>
[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 COMPLEX 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 COMPLEX array, dimension (LDAFB,N)
!>          Details of the LU factorization of the band matrix A, as
!>          computed by CGBTRF.  U is stored as an upper triangular band
!>          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
!>          the multipliers used during the factorization are stored in
!>          rows KL+KU+2 to 2*KL+KU+1.
!>
[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 CGBTRF; for 1<=i<=N, row i of the
!>          matrix was interchanged with row IPIV(i).
!>
[in]B
!>          B is COMPLEX 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 COMPLEX array, dimension (LDX,NRHS)
!>          On entry, the solution matrix X, as computed by CGBTRS.
!>          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 COMPLEX array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is REAL 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 203 of file cgbrfs.f.

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

◆ cgbrfsx()

subroutine cgbrfsx ( character trans,
character equed,
integer n,
integer kl,
integer ku,
integer nrhs,
complex, dimension( ldab, * ) ab,
integer ldab,
complex, dimension( ldafb, * ) afb,
integer ldafb,
integer, dimension( * ) ipiv,
real, dimension( * ) r,
real, dimension( * ) c,
complex, dimension( ldb, * ) b,
integer ldb,
complex, 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,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer info )

CGBRFSX

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

Purpose:
!>
!>    CGBRFSX 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)
!>
[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 COMPLEX 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 COMPLEX array, dimension (LDAFB,N)
!>     Details of the LU factorization of the band matrix A, as
!>     computed by CGBTRF.  U is stored as an upper triangular band
!>     matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
!>     the multipliers used during the factorization are stored in
!>     rows KL+KU+2 to 2*KL+KU+1.
!>
[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 CGETRF; 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 COMPLEX 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 COMPLEX array, dimension (LDX,NRHS)
!>     On entry, the solution matrix X, as computed by CGETRS.
!>     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 COMPLEX array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is REAL array, dimension (2*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 cgbrfsx.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( * )
453 COMPLEX 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, * ), RWORK( * )
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, IGNORE_CWISE
486 INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE, N_NORMS,
487 $ ITHRESH
488 REAL ANORM, RCOND_TMP, ILLRCOND_THRESH, ERR_LBND,
489 $ CWISE_WRONG, RTHRESH, UNSTABLE_THRESH
490* ..
491* .. External Subroutines ..
492 EXTERNAL xerbla, cgbcon, cla_gbrfsx_extended
493* ..
494* .. Intrinsic Functions ..
495 INTRINSIC max, sqrt, transfer
496* ..
497* .. External Functions ..
498 EXTERNAL lsame, ilatrans, ilaprec
499 EXTERNAL slamch, clangb, cla_gbrcond_x, cla_gbrcond_c
500 REAL SLAMCH, CLANGB, CLA_GBRCOND_X, CLA_GBRCOND_C
501 LOGICAL LSAME
502 INTEGER ILATRANS, ILAPREC
503* ..
504* .. Executable Statements ..
505*
506* Check the input parameters.
507*
508 info = 0
509 trans_type = ilatrans( trans )
510 ref_type = int( itref_default )
511 IF ( nparams .GE. la_linrx_itref_i ) THEN
512 IF ( params( la_linrx_itref_i ) .LT. 0.0 ) THEN
513 params( la_linrx_itref_i ) = itref_default
514 ELSE
515 ref_type = params( la_linrx_itref_i )
516 END IF
517 END IF
518*
519* Set default parameters.
520*
521 illrcond_thresh = real( n ) * slamch( 'Epsilon' )
522 ithresh = int( ithresh_default )
523 rthresh = rthresh_default
524 unstable_thresh = dzthresh_default
525 ignore_cwise = componentwise_default .EQ. 0.0
526*
527 IF ( nparams.GE.la_linrx_ithresh_i ) THEN
528 IF ( params( la_linrx_ithresh_i ).LT.0.0 ) THEN
529 params( la_linrx_ithresh_i ) = ithresh
530 ELSE
531 ithresh = int( params( la_linrx_ithresh_i ) )
532 END IF
533 END IF
534 IF ( nparams.GE.la_linrx_cwise_i ) THEN
535 IF ( params( la_linrx_cwise_i ).LT.0.0 ) THEN
536 IF ( ignore_cwise ) THEN
537 params( la_linrx_cwise_i ) = 0.0
538 ELSE
539 params( la_linrx_cwise_i ) = 1.0
540 END IF
541 ELSE
542 ignore_cwise = params( la_linrx_cwise_i ) .EQ. 0.0
543 END IF
544 END IF
545 IF ( ref_type .EQ. 0 .OR. n_err_bnds .EQ. 0 ) THEN
546 n_norms = 0
547 ELSE IF ( ignore_cwise ) THEN
548 n_norms = 1
549 ELSE
550 n_norms = 2
551 END IF
552*
553 notran = lsame( trans, 'N' )
554 rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
555 colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
556*
557* Test input parameters.
558*
559 IF( trans_type.EQ.-1 ) THEN
560 info = -1
561 ELSE IF( .NOT.rowequ .AND. .NOT.colequ .AND.
562 $ .NOT.lsame( equed, 'N' ) ) THEN
563 info = -2
564 ELSE IF( n.LT.0 ) THEN
565 info = -3
566 ELSE IF( kl.LT.0 ) THEN
567 info = -4
568 ELSE IF( ku.LT.0 ) THEN
569 info = -5
570 ELSE IF( nrhs.LT.0 ) THEN
571 info = -6
572 ELSE IF( ldab.LT.kl+ku+1 ) THEN
573 info = -8
574 ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
575 info = -10
576 ELSE IF( ldb.LT.max( 1, n ) ) THEN
577 info = -13
578 ELSE IF( ldx.LT.max( 1, n ) ) THEN
579 info = -15
580 END IF
581 IF( info.NE.0 ) THEN
582 CALL xerbla( 'CGBRFSX', -info )
583 RETURN
584 END IF
585*
586* Quick return if possible.
587*
588 IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
589 rcond = 1.0
590 DO j = 1, nrhs
591 berr( j ) = 0.0
592 IF ( n_err_bnds .GE. 1 ) THEN
593 err_bnds_norm( j, la_linrx_trust_i ) = 1.0
594 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
595 END IF
596 IF ( n_err_bnds .GE. 2 ) THEN
597 err_bnds_norm( j, la_linrx_err_i ) = 0.0
598 err_bnds_comp( j, la_linrx_err_i ) = 0.0
599 END IF
600 IF ( n_err_bnds .GE. 3 ) THEN
601 err_bnds_norm( j, la_linrx_rcond_i ) = 1.0
602 err_bnds_comp( j, la_linrx_rcond_i ) = 1.0
603 END IF
604 END DO
605 RETURN
606 END IF
607*
608* Default to failure.
609*
610 rcond = 0.0
611 DO j = 1, nrhs
612 berr( j ) = 1.0
613 IF ( n_err_bnds .GE. 1 ) THEN
614 err_bnds_norm( j, la_linrx_trust_i ) = 1.0
615 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
616 END IF
617 IF ( n_err_bnds .GE. 2 ) THEN
618 err_bnds_norm( j, la_linrx_err_i ) = 1.0
619 err_bnds_comp( j, la_linrx_err_i ) = 1.0
620 END IF
621 IF ( n_err_bnds .GE. 3 ) THEN
622 err_bnds_norm( j, la_linrx_rcond_i ) = 0.0
623 err_bnds_comp( j, la_linrx_rcond_i ) = 0.0
624 END IF
625 END DO
626*
627* Compute the norm of A and the reciprocal of the condition
628* number of A.
629*
630 IF( notran ) THEN
631 norm = 'I'
632 ELSE
633 norm = '1'
634 END IF
635 anorm = clangb( norm, n, kl, ku, ab, ldab, rwork )
636 CALL cgbcon( norm, n, kl, ku, afb, ldafb, ipiv, anorm, rcond,
637 $ work, rwork, info )
638*
639* Perform refinement on each right-hand side
640*
641 IF ( ref_type .NE. 0 .AND. info .EQ. 0 ) THEN
642
643 prec_type = ilaprec( 'D' )
644
645 IF ( notran ) THEN
646 CALL cla_gbrfsx_extended( prec_type, trans_type, n, kl, ku,
647 $ nrhs, ab, ldab, afb, ldafb, ipiv, colequ, c, b,
648 $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
649 $ err_bnds_comp, work, rwork, work(n+1),
650 $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
651 $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
652 $ info )
653 ELSE
654 CALL cla_gbrfsx_extended( prec_type, trans_type, n, kl, ku,
655 $ nrhs, ab, ldab, afb, ldafb, ipiv, rowequ, r, b,
656 $ ldb, x, ldx, berr, n_norms, err_bnds_norm,
657 $ err_bnds_comp, work, rwork, work(n+1),
658 $ transfer(rwork(1:2*n), (/ (zero, zero) /), n),
659 $ rcond, ithresh, rthresh, unstable_thresh, ignore_cwise,
660 $ info )
661 END IF
662 END IF
663
664 err_lbnd = max( 10.0, sqrt( real( n ) ) ) * slamch( 'Epsilon' )
665 IF (n_err_bnds .GE. 1 .AND. n_norms .GE. 1) THEN
666*
667* Compute scaled normwise condition number cond(A*C).
668*
669 IF ( colequ .AND. notran ) THEN
670 rcond_tmp = cla_gbrcond_c( trans, n, kl, ku, ab, ldab, afb,
671 $ ldafb, ipiv, c, .true., info, work, rwork )
672 ELSE IF ( rowequ .AND. .NOT. notran ) THEN
673 rcond_tmp = cla_gbrcond_c( trans, n, kl, ku, ab, ldab, afb,
674 $ ldafb, ipiv, r, .true., info, work, rwork )
675 ELSE
676 rcond_tmp = cla_gbrcond_c( trans, n, kl, ku, ab, ldab, afb,
677 $ ldafb, ipiv, c, .false., info, work, rwork )
678 END IF
679 DO j = 1, nrhs
680*
681* Cap the error at 1.0.
682*
683 IF ( n_err_bnds .GE. la_linrx_err_i
684 $ .AND. err_bnds_norm( j, la_linrx_err_i ) .GT. 1.0)
685 $ err_bnds_norm( j, la_linrx_err_i ) = 1.0
686*
687* Threshold the error (see LAWN).
688*
689 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
690 err_bnds_norm( j, la_linrx_err_i ) = 1.0
691 err_bnds_norm( j, la_linrx_trust_i ) = 0.0
692 IF ( info .LE. n ) info = n + j
693 ELSE IF ( err_bnds_norm( j, la_linrx_err_i ) .LT. err_lbnd )
694 $ THEN
695 err_bnds_norm( j, la_linrx_err_i ) = err_lbnd
696 err_bnds_norm( j, la_linrx_trust_i ) = 1.0
697 END IF
698*
699* Save the condition number.
700*
701 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
702 err_bnds_norm( j, la_linrx_rcond_i ) = rcond_tmp
703 END IF
704
705 END DO
706 END IF
707
708 IF (n_err_bnds .GE. 1 .AND. n_norms .GE. 2) THEN
709*
710* Compute componentwise condition number cond(A*diag(Y(:,J))) for
711* each right-hand side using the current solution as an estimate of
712* the true solution. If the componentwise error estimate is too
713* large, then the solution is a lousy estimate of truth and the
714* estimated RCOND may be too optimistic. To avoid misleading users,
715* the inverse condition number is set to 0.0 when the estimated
716* cwise error is at least CWISE_WRONG.
717*
718 cwise_wrong = sqrt( slamch( 'Epsilon' ) )
719 DO j = 1, nrhs
720 IF (err_bnds_comp( j, la_linrx_err_i ) .LT. cwise_wrong )
721 $ THEN
722 rcond_tmp = cla_gbrcond_x( trans, n, kl, ku, ab, ldab,
723 $ afb, ldafb, ipiv, x( 1, j ), info, work, rwork )
724 ELSE
725 rcond_tmp = 0.0
726 END IF
727*
728* Cap the error at 1.0.
729*
730 IF ( n_err_bnds .GE. la_linrx_err_i
731 $ .AND. err_bnds_comp( j, la_linrx_err_i ) .GT. 1.0 )
732 $ err_bnds_comp( j, la_linrx_err_i ) = 1.0
733*
734* Threshold the error (see LAWN).
735*
736 IF ( rcond_tmp .LT. illrcond_thresh ) THEN
737 err_bnds_comp( j, la_linrx_err_i ) = 1.0
738 err_bnds_comp( j, la_linrx_trust_i ) = 0.0
739 IF ( params( la_linrx_cwise_i ) .EQ. 1.0
740 $ .AND. info.LT.n + j ) info = n + j
741 ELSE IF ( err_bnds_comp( j, la_linrx_err_i )
742 $ .LT. err_lbnd ) THEN
743 err_bnds_comp( j, la_linrx_err_i ) = err_lbnd
744 err_bnds_comp( j, la_linrx_trust_i ) = 1.0
745 END IF
746*
747* Save the condition number.
748*
749 IF ( n_err_bnds .GE. la_linrx_rcond_i ) THEN
750 err_bnds_comp( j, la_linrx_rcond_i ) = rcond_tmp
751 END IF
752
753 END DO
754 END IF
755*
756 RETURN
757*
758* End of CGBRFSX
759*
ilaprec
integer function ilaprec(prec)
ILAPREC
Definition ilaprec.f:58
ilatrans
integer function ilatrans(trans)
ILATRANS
Definition ilatrans.f:58
clangb
real function clangb(norm, n, kl, ku, ab, ldab, work)
CLANGB returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clangb.f:125
cla_gbrfsx_extended
subroutine cla_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)
CLA_GBRFSX_EXTENDED improves the computed solution to a system of linear equations for general banded...
Definition cla_gbrfsx_extended.f:410
cla_gbrcond_c
real function cla_gbrcond_c(trans, n, kl, ku, ab, ldab, afb, ldafb, ipiv, c, capply, info, work, rwork)
CLA_GBRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general banded ma...
Definition cla_gbrcond_c.f:161
cgbcon
subroutine cgbcon(norm, n, kl, ku, ab, ldab, ipiv, anorm, rcond, work, rwork, info)
CGBCON
Definition cgbcon.f:147
cla_gbrcond_x
real function cla_gbrcond_x(trans, n, kl, ku, ab, ldab, afb, ldafb, ipiv, x, info, work, rwork)
CLA_GBRCOND_X computes the infinity norm condition number of op(A)*diag(x) for general banded matrice...
Definition cla_gbrcond_x.f:153

◆ cgbtf2()

subroutine cgbtf2 ( integer m,
integer n,
integer kl,
integer ku,
complex, dimension( ldab, * ) ab,
integer ldab,
integer, dimension( * ) ipiv,
integer info )

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

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

Purpose:
!>
!> CGBTF2 computes an LU factorization of a complex 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 COMPLEX 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 cgbtf2.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 COMPLEX AB( LDAB, * )
156* ..
157*
158* =====================================================================
159*
160* .. Parameters ..
161 COMPLEX ONE, ZERO
162 parameter( one = ( 1.0e+0, 0.0e+0 ),
163 $ zero = ( 0.0e+0, 0.0e+0 ) )
164* ..
165* .. Local Scalars ..
166 INTEGER I, J, JP, JU, KM, KV
167* ..
168* .. External Functions ..
169 INTEGER ICAMAX
170 EXTERNAL icamax
171* ..
172* .. External Subroutines ..
173 EXTERNAL cgeru, cscal, cswap, xerbla
174* ..
175* .. Intrinsic Functions ..
176 INTRINSIC max, min
177* ..
178* .. Executable Statements ..
179*
180* KV is the number of superdiagonals in the factor U, allowing for
181* fill-in.
182*
183 kv = ku + kl
184*
185* Test the input parameters.
186*
187 info = 0
188 IF( m.LT.0 ) THEN
189 info = -1
190 ELSE IF( n.LT.0 ) THEN
191 info = -2
192 ELSE IF( kl.LT.0 ) THEN
193 info = -3
194 ELSE IF( ku.LT.0 ) THEN
195 info = -4
196 ELSE IF( ldab.LT.kl+kv+1 ) THEN
197 info = -6
198 END IF
199 IF( info.NE.0 ) THEN
200 CALL xerbla( 'CGBTF2', -info )
201 RETURN
202 END IF
203*
204* Quick return if possible
205*
206 IF( m.EQ.0 .OR. n.EQ.0 )
207 $ RETURN
208*
209* Gaussian elimination with partial pivoting
210*
211* Set fill-in elements in columns KU+2 to KV to zero.
212*
213 DO 20 j = ku + 2, min( kv, n )
214 DO 10 i = kv - j + 2, kl
215 ab( i, j ) = zero
216 10 CONTINUE
217 20 CONTINUE
218*
219* JU is the index of the last column affected by the current stage
220* of the factorization.
221*
222 ju = 1
223*
224 DO 40 j = 1, min( m, n )
225*
226* Set fill-in elements in column J+KV to zero.
227*
228 IF( j+kv.LE.n ) THEN
229 DO 30 i = 1, kl
230 ab( i, j+kv ) = zero
231 30 CONTINUE
232 END IF
233*
234* Find pivot and test for singularity. KM is the number of
235* subdiagonal elements in the current column.
236*
237 km = min( kl, m-j )
238 jp = icamax( km+1, ab( kv+1, j ), 1 )
239 ipiv( j ) = jp + j - 1
240 IF( ab( kv+jp, j ).NE.zero ) THEN
241 ju = max( ju, min( j+ku+jp-1, n ) )
242*
243* Apply interchange to columns J to JU.
244*
245 IF( jp.NE.1 )
246 $ CALL cswap( ju-j+1, ab( kv+jp, j ), ldab-1,
247 $ ab( kv+1, j ), ldab-1 )
248 IF( km.GT.0 ) THEN
249*
250* Compute multipliers.
251*
252 CALL cscal( 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 cgeru( 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 CGBTF2
273*
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
cgeru
subroutine cgeru(m, n, alpha, x, incx, y, incy, a, lda)
CGERU
Definition cgeru.f:130

◆ cgbtrf()

subroutine cgbtrf ( integer m,
integer n,
integer kl,
integer ku,
complex, dimension( ldab, * ) ab,
integer ldab,
integer, dimension( * ) ipiv,
integer info )

CGBTRF

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

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

◆ cgbtrs()

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

CGBTRS

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

Purpose:
!>
!> CGBTRS solves a system of linear equations
!>    A * X = B,  A**T * X = B,  or  A**H * X = B
!> with a general band matrix A using the LU factorization computed
!> by CGBTRF.
!>
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)
!>
[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 COMPLEX array, dimension (LDAB,N)
!>          Details of the LU factorization of the band matrix A, as
!>          computed by CGBTRF.  U is stored as an upper triangular band
!>          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
!>          the multipliers used during the factorization are stored in
!>          rows KL+KU+2 to 2*KL+KU+1.
!>
[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 COMPLEX 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 cgbtrs.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 COMPLEX AB( LDAB, * ), B( LDB, * )
150* ..
151*
152* =====================================================================
153*
154* .. Parameters ..
155 COMPLEX ONE
156 parameter( one = ( 1.0e+0, 0.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 cgemv, cgeru, clacgv, cswap, ctbsv, 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( 'CGBTRS', -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 cswap( nrhs, b( l, 1 ), ldb, b( j, 1 ), ldb )
224 CALL cgeru( 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 ctbsv( 'Upper', 'No transpose', 'Non-unit', n, kl+ku,
234 $ ab, ldab, b( 1, i ), 1 )
235 20 CONTINUE
236*
237 ELSE IF( lsame( trans, 'T' ) ) THEN
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 ctbsv( '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 cgemv( '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 cswap( nrhs, b( l, 1 ), ldb, b( j, 1 ), ldb )
259 40 CONTINUE
260 END IF
261*
262 ELSE
263*
264* Solve A**H * X = B.
265*
266 DO 50 i = 1, nrhs
267*
268* Solve U**H * X = B, overwriting B with X.
269*
270 CALL ctbsv( 'Upper', 'Conjugate transpose', 'Non-unit', n,
271 $ kl+ku, ab, ldab, b( 1, i ), 1 )
272 50 CONTINUE
273*
274* Solve L**H * X = B, overwriting B with X.
275*
276 IF( lnoti ) THEN
277 DO 60 j = n - 1, 1, -1
278 lm = min( kl, n-j )
279 CALL clacgv( nrhs, b( j, 1 ), ldb )
280 CALL cgemv( 'Conjugate transpose', lm, nrhs, -one,
281 $ b( j+1, 1 ), ldb, ab( kd+1, j ), 1, one,
282 $ b( j, 1 ), ldb )
283 CALL clacgv( nrhs, b( j, 1 ), ldb )
284 l = ipiv( j )
285 IF( l.NE.j )
286 $ CALL cswap( nrhs, b( l, 1 ), ldb, b( j, 1 ), ldb )
287 60 CONTINUE
288 END IF
289 END IF
290 RETURN
291*
292* End of CGBTRS
293*
clacgv
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:74
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:158
ctbsv
subroutine ctbsv(uplo, trans, diag, n, k, a, lda, x, incx)
CTBSV
Definition ctbsv.f:189

◆ cggbak()

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

CGGBAK

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

Purpose:
!>
!> CGGBAK forms the right or left eigenvectors of a complex generalized
!> eigenvalue problem A*x = lambda*B*x, by backward transformation on
!> the computed eigenvectors of the balanced pair of matrices output by
!> CGGBAL.
!>
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 CGGBAL.
!>
[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 CGGBAL.
!>          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 CGGBAL.
!>
[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 CGGBAL.
!>
[in]M
!>          M is INTEGER
!>          The number of columns of the matrix V.  M >= 0.
!>
[in,out]V
!>          V is COMPLEX array, dimension (LDV,M)
!>          On entry, the matrix of right or left eigenvectors to be
!>          transformed, as returned by CTGEVC.
!>          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 146 of file cggbak.f.

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

◆ cggbal()

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

CGGBAL

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

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

◆ cla_gbamv()

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

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

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

Purpose:
!>
!> CLA_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 COMPLEX 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, LDAB 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 COMPLEX 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 184 of file cla_gbamv.f.

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

◆ cla_gbrcond_c()

real function cla_gbrcond_c ( character trans,
integer n,
integer kl,
integer ku,
complex, dimension( ldab, * ) ab,
integer ldab,
complex, dimension( ldafb, * ) afb,
integer ldafb,
integer, dimension( * ) ipiv,
real, dimension( * ) c,
logical capply,
integer info,
complex, dimension( * ) work,
real, dimension( * ) rwork )

CLA_GBRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general banded matrices.

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

Purpose:
!>
!>    CLA_GBRCOND_C Computes the infinity norm condition number of
!>    op(A) * inv(diag(C)) where C is a REAL vector.
!>
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 COMPLEX 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 COMPLEX array, dimension (LDAFB,N)
!>     Details of the LU factorization of the band matrix A, as
!>     computed by CGBTRF.  U is stored as an upper triangular
!>     band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
!>     and the multipliers used during the factorization are stored
!>     in rows KL+KU+2 to 2*KL+KU+1.
!>
[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 CGBTRF; row i of the matrix was interchanged
!>     with row IPIV(i).
!>
[in]C
!>          C is REAL array, dimension (N)
!>     The vector C in the formula op(A) * inv(diag(C)).
!>
[in]CAPPLY
!>          CAPPLY is LOGICAL
!>     If .TRUE. then access the vector C in the formula above.
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit.
!>     i > 0:  The ith argument is invalid.
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (2*N).
!>     Workspace.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (N).
!>     Workspace.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 158 of file cla_gbrcond_c.f.

161*
162* -- LAPACK computational routine --
163* -- LAPACK is a software package provided by Univ. of Tennessee, --
164* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165*
166* .. Scalar Arguments ..
167 CHARACTER TRANS
168 LOGICAL CAPPLY
169 INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO
170* ..
171* .. Array Arguments ..
172 INTEGER IPIV( * )
173 COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), WORK( * )
174 REAL C( * ), RWORK( * )
175* ..
176*
177* =====================================================================
178*
179* .. Local Scalars ..
180 LOGICAL NOTRANS
181 INTEGER KASE, I, J
182 REAL AINVNM, ANORM, TMP
183 COMPLEX ZDUM
184* ..
185* .. Local Arrays ..
186 INTEGER ISAVE( 3 )
187* ..
188* .. External Functions ..
189 LOGICAL LSAME
190 EXTERNAL lsame
191* ..
192* .. External Subroutines ..
193 EXTERNAL clacn2, cgbtrs, xerbla
194* ..
195* .. Intrinsic Functions ..
196 INTRINSIC abs, max
197* ..
198* .. Statement Functions ..
199 REAL CABS1
200* ..
201* .. Statement Function Definitions ..
202 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
203* ..
204* .. Executable Statements ..
205 cla_gbrcond_c = 0.0e+0
206*
207 info = 0
208 notrans = lsame( trans, 'N' )
209 IF ( .NOT. notrans .AND. .NOT. lsame( trans, 'T' ) .AND. .NOT.
210 $ 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( 'CLA_GBRCOND_C', -info )
225 RETURN
226 END IF
227*
228* Compute norm of op(A)*op2(C).
229*
230 anorm = 0.0e+0
231 kd = ku + 1
232 ke = kl + 1
233 IF ( notrans ) THEN
234 DO i = 1, n
235 tmp = 0.0e+0
236 IF ( capply ) THEN
237 DO j = max( i-kl, 1 ), min( i+ku, n )
238 tmp = tmp + cabs1( ab( kd+i-j, j ) ) / c( j )
239 END DO
240 ELSE
241 DO j = max( i-kl, 1 ), min( i+ku, n )
242 tmp = tmp + cabs1( ab( kd+i-j, j ) )
243 END DO
244 END IF
245 rwork( i ) = tmp
246 anorm = max( anorm, tmp )
247 END DO
248 ELSE
249 DO i = 1, n
250 tmp = 0.0e+0
251 IF ( capply ) THEN
252 DO j = max( i-kl, 1 ), min( i+ku, n )
253 tmp = tmp + cabs1( ab( ke-i+j, i ) ) / c( j )
254 END DO
255 ELSE
256 DO j = max( i-kl, 1 ), min( i+ku, n )
257 tmp = tmp + cabs1( ab( ke-i+j, i ) )
258 END DO
259 END IF
260 rwork( i ) = tmp
261 anorm = max( anorm, tmp )
262 END DO
263 END IF
264*
265* Quick return if possible.
266*
267 IF( n.EQ.0 ) THEN
268 cla_gbrcond_c = 1.0e+0
269 RETURN
270 ELSE IF( anorm .EQ. 0.0e+0 ) THEN
271 RETURN
272 END IF
273*
274* Estimate the norm of inv(op(A)).
275*
276 ainvnm = 0.0e+0
277*
278 kase = 0
279 10 CONTINUE
280 CALL clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
281 IF( kase.NE.0 ) THEN
282 IF( kase.EQ.2 ) THEN
283*
284* Multiply by R.
285*
286 DO i = 1, n
287 work( i ) = work( i ) * rwork( i )
288 END DO
289*
290 IF ( notrans ) THEN
291 CALL cgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
292 $ ipiv, work, n, info )
293 ELSE
294 CALL cgbtrs( 'Conjugate transpose', n, kl, ku, 1, afb,
295 $ ldafb, ipiv, work, n, info )
296 ENDIF
297*
298* Multiply by inv(C).
299*
300 IF ( capply ) THEN
301 DO i = 1, n
302 work( i ) = work( i ) * c( i )
303 END DO
304 END IF
305 ELSE
306*
307* Multiply by inv(C**H).
308*
309 IF ( capply ) THEN
310 DO i = 1, n
311 work( i ) = work( i ) * c( i )
312 END DO
313 END IF
314*
315 IF ( notrans ) THEN
316 CALL cgbtrs( 'Conjugate transpose', n, kl, ku, 1, afb,
317 $ ldafb, ipiv, work, n, info )
318 ELSE
319 CALL cgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
320 $ ipiv, work, n, info )
321 END IF
322*
323* Multiply by R.
324*
325 DO i = 1, n
326 work( i ) = work( i ) * rwork( i )
327 END DO
328 END IF
329 GO TO 10
330 END IF
331*
332* Compute the estimate of the reciprocal condition number.
333*
334 IF( ainvnm .NE. 0.0e+0 )
335 $ cla_gbrcond_c = 1.0e+0 / ainvnm
336*
337 RETURN
338*
339* End of CLA_GBRCOND_C
340*

◆ cla_gbrcond_x()

real function cla_gbrcond_x ( character trans,
integer n,
integer kl,
integer ku,
complex, dimension( ldab, * ) ab,
integer ldab,
complex, dimension( ldafb, * ) afb,
integer ldafb,
integer, dimension( * ) ipiv,
complex, dimension( * ) x,
integer info,
complex, dimension( * ) work,
real, dimension( * ) rwork )

CLA_GBRCOND_X computes the infinity norm condition number of op(A)*diag(x) for general banded matrices.

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

Purpose:
!>
!>    CLA_GBRCOND_X Computes the infinity norm condition number of
!>    op(A) * diag(X) where X is a COMPLEX vector.
!>
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 COMPLEX 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 COMPLEX array, dimension (LDAFB,N)
!>     Details of the LU factorization of the band matrix A, as
!>     computed by CGBTRF.  U is stored as an upper triangular
!>     band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
!>     and the multipliers used during the factorization are stored
!>     in rows KL+KU+2 to 2*KL+KU+1.
!>
[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 CGBTRF; row i of the matrix was interchanged
!>     with row IPIV(i).
!>
[in]X
!>          X is COMPLEX array, dimension (N)
!>     The vector X in the formula op(A) * diag(X).
!>
[out]INFO
!>          INFO is INTEGER
!>       = 0:  Successful exit.
!>     i > 0:  The ith argument is invalid.
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (2*N).
!>     Workspace.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (N).
!>     Workspace.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 151 of file cla_gbrcond_x.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 CHARACTER TRANS
160 INTEGER N, KL, KU, KD, KE, LDAB, LDAFB, INFO
161* ..
162* .. Array Arguments ..
163 INTEGER IPIV( * )
164 COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), WORK( * ),
165 $ X( * )
166 REAL RWORK( * )
167* ..
168*
169* =====================================================================
170*
171* .. Local Scalars ..
172 LOGICAL NOTRANS
173 INTEGER KASE, I, J
174 REAL AINVNM, ANORM, TMP
175 COMPLEX ZDUM
176* ..
177* .. Local Arrays ..
178 INTEGER ISAVE( 3 )
179* ..
180* .. External Functions ..
181 LOGICAL LSAME
182 EXTERNAL lsame
183* ..
184* .. External Subroutines ..
185 EXTERNAL clacn2, cgbtrs, xerbla
186* ..
187* .. Intrinsic Functions ..
188 INTRINSIC abs, max
189* ..
190* .. Statement Functions ..
191 REAL CABS1
192* ..
193* .. Statement Function Definitions ..
194 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
195* ..
196* .. Executable Statements ..
197*
198 cla_gbrcond_x = 0.0e+0
199*
200 info = 0
201 notrans = lsame( trans, 'N' )
202 IF ( .NOT. notrans .AND. .NOT. lsame(trans, 'T') .AND. .NOT.
203 $ lsame( trans, 'C' ) ) THEN
204 info = -1
205 ELSE IF( n.LT.0 ) THEN
206 info = -2
207 ELSE IF( kl.LT.0 .OR. kl.GT.n-1 ) THEN
208 info = -3
209 ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
210 info = -4
211 ELSE IF( ldab.LT.kl+ku+1 ) THEN
212 info = -6
213 ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
214 info = -8
215 END IF
216 IF( info.NE.0 ) THEN
217 CALL xerbla( 'CLA_GBRCOND_X', -info )
218 RETURN
219 END IF
220*
221* Compute norm of op(A)*op2(C).
222*
223 kd = ku + 1
224 ke = kl + 1
225 anorm = 0.0
226 IF ( notrans ) THEN
227 DO i = 1, n
228 tmp = 0.0e+0
229 DO j = max( i-kl, 1 ), min( i+ku, n )
230 tmp = tmp + cabs1( ab( kd+i-j, j) * x( j ) )
231 END DO
232 rwork( i ) = tmp
233 anorm = max( anorm, tmp )
234 END DO
235 ELSE
236 DO i = 1, n
237 tmp = 0.0e+0
238 DO j = max( i-kl, 1 ), min( i+ku, n )
239 tmp = tmp + cabs1( ab( ke-i+j, i ) * x( j ) )
240 END DO
241 rwork( i ) = tmp
242 anorm = max( anorm, tmp )
243 END DO
244 END IF
245*
246* Quick return if possible.
247*
248 IF( n.EQ.0 ) THEN
249 cla_gbrcond_x = 1.0e+0
250 RETURN
251 ELSE IF( anorm .EQ. 0.0e+0 ) THEN
252 RETURN
253 END IF
254*
255* Estimate the norm of inv(op(A)).
256*
257 ainvnm = 0.0e+0
258*
259 kase = 0
260 10 CONTINUE
261 CALL clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
262 IF( kase.NE.0 ) THEN
263 IF( kase.EQ.2 ) THEN
264*
265* Multiply by R.
266*
267 DO i = 1, n
268 work( i ) = work( i ) * rwork( i )
269 END DO
270*
271 IF ( notrans ) THEN
272 CALL cgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
273 $ ipiv, work, n, info )
274 ELSE
275 CALL cgbtrs( 'Conjugate transpose', n, kl, ku, 1, afb,
276 $ ldafb, ipiv, work, n, info )
277 ENDIF
278*
279* Multiply by inv(X).
280*
281 DO i = 1, n
282 work( i ) = work( i ) / x( i )
283 END DO
284 ELSE
285*
286* Multiply by inv(X**H).
287*
288 DO i = 1, n
289 work( i ) = work( i ) / x( i )
290 END DO
291*
292 IF ( notrans ) THEN
293 CALL cgbtrs( 'Conjugate transpose', n, kl, ku, 1, afb,
294 $ ldafb, ipiv, work, n, info )
295 ELSE
296 CALL cgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
297 $ ipiv, work, n, info )
298 END IF
299*
300* Multiply by R.
301*
302 DO i = 1, n
303 work( i ) = work( i ) * rwork( i )
304 END DO
305 END IF
306 GO TO 10
307 END IF
308*
309* Compute the estimate of the reciprocal condition number.
310*
311 IF( ainvnm .NE. 0.0e+0 )
312 $ cla_gbrcond_x = 1.0e+0 / ainvnm
313*
314 RETURN
315*
316* End of CLA_GBRCOND_X
317*

◆ cla_gbrfsx_extended()

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

CLA_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 CLA_GBRFSX_EXTENDED + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> CLA_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 CGBRFSX 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 COMPLEX 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 COMPLEX array, dimension (LDAF,N)
!>     The factors L and U from the factorization
!>     A = P*L*U as computed by CGBTRF.
!>
[in]LDAFB
!>          LDAFB is INTEGER
!>     The leading dimension of the array AF.  LDAF >= max(1,N).
!>
[in]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>     The pivot indices from the factorization A = P*L*U
!>     as computed by CGBTRF; 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 COMPLEX 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 COMPLEX array, dimension (LDY,NRHS)
!>     On entry, the solution matrix X, as computed by CGBTRS.
!>     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 CLA_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 COMPLEX array, dimension (N)
!>     Workspace to hold the intermediate residual.
!>
[in]AYB
!>          AYB is REAL array, dimension (N)
!>     Workspace.
!>
[in]DY
!>          DY is COMPLEX array, dimension (N)
!>     Workspace to hold the intermediate solution.
!>
[in]Y_TAIL
!>          Y_TAIL is COMPLEX 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 CGBTRS had an illegal
!>             value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 403 of file cla_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 COMPLEX 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 COMPLEX ZDUM
441* ..
442* .. Parameters ..
443 INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
444 $ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL,
445 $ EXTRA_Y
446 parameter( unstable_state = 0, working_state = 1,
447 $ conv_state = 2, noprog_state = 3 )
448 parameter( base_residual = 0, extra_residual = 1,
449 $ extra_y = 2 )
450 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
451 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
452 INTEGER CMP_ERR_I, PIV_GROWTH_I
453 parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
454 $ berr_i = 3 )
455 parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
456 parameter( cmp_rcond_i = 7, cmp_err_i = 8,
457 $ piv_growth_i = 9 )
458 INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
459 $ LA_LINRX_CWISE_I
460 parameter( la_linrx_itref_i = 1,
461 $ la_linrx_ithresh_i = 2 )
462 parameter( la_linrx_cwise_i = 3 )
463 INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
464 $ LA_LINRX_RCOND_I
465 parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
466 parameter( la_linrx_rcond_i = 3 )
467* ..
468* .. External Subroutines ..
469 EXTERNAL caxpy, ccopy, cgbtrs, cgbmv, blas_cgbmv_x,
470 $ blas_cgbmv2_x, cla_gbamv, cla_wwaddw, slamch,
471 $ chla_transtype, cla_lin_berr
472 REAL SLAMCH
473 CHARACTER CHLA_TRANSTYPE
474* ..
475* .. Intrinsic Functions..
476 INTRINSIC abs, max, min
477* ..
478* .. Statement Functions ..
479 REAL CABS1
480* ..
481* .. Statement Function Definitions ..
482 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
483* ..
484* .. Executable Statements ..
485*
486 IF (info.NE.0) RETURN
487 trans = chla_transtype(trans_type)
488 eps = slamch( 'Epsilon' )
489 hugeval = slamch( 'Overflow' )
490* Force HUGEVAL to Inf
491 hugeval = hugeval * hugeval
492* Using HUGEVAL may lead to spurious underflows.
493 incr_thresh = real( n ) * eps
494 m = kl+ku+1
495
496 DO j = 1, nrhs
497 y_prec_state = extra_residual
498 IF ( y_prec_state .EQ. extra_y ) then
499 DO i = 1, n
500 y_tail( i ) = 0.0
501 END DO
502 END IF
503
504 dxrat = 0.0e+0
505 dxratmax = 0.0e+0
506 dzrat = 0.0e+0
507 dzratmax = 0.0e+0
508 final_dx_x = hugeval
509 final_dz_z = hugeval
510 prevnormdx = hugeval
511 prev_dz_z = hugeval
512 dz_z = hugeval
513 dx_x = hugeval
514
515 x_state = working_state
516 z_state = unstable_state
517 incr_prec = .false.
518
519 DO cnt = 1, ithresh
520*
521* Compute residual RES = B_s - op(A_s) * Y,
522* op(A) = A, A**T, or A**H depending on TRANS (and type).
523*
524 CALL ccopy( n, b( 1, j ), 1, res, 1 )
525 IF ( y_prec_state .EQ. base_residual ) THEN
526 CALL cgbmv( trans, m, n, kl, ku, (-1.0e+0,0.0e+0), ab,
527 $ ldab, y( 1, j ), 1, (1.0e+0,0.0e+0), res, 1 )
528 ELSE IF ( y_prec_state .EQ. extra_residual ) THEN
529 CALL blas_cgbmv_x( trans_type, n, n, kl, ku,
530 $ (-1.0e+0,0.0e+0), ab, ldab, y( 1, j ), 1,
531 $ (1.0e+0,0.0e+0), res, 1, prec_type )
532 ELSE
533 CALL blas_cgbmv2_x( trans_type, n, n, kl, ku,
534 $ (-1.0e+0,0.0e+0), ab, ldab, y( 1, j ), y_tail, 1,
535 $ (1.0e+0,0.0e+0), res, 1, prec_type )
536 END IF
537
538! XXX: RES is no longer needed.
539 CALL ccopy( n, res, 1, dy, 1 )
540 CALL cgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv, dy, n,
541 $ info )
542*
543* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
544*
545 normx = 0.0e+0
546 normy = 0.0e+0
547 normdx = 0.0e+0
548 dz_z = 0.0e+0
549 ymin = hugeval
550
551 DO i = 1, n
552 yk = cabs1( y( i, j ) )
553 dyk = cabs1( dy( i ) )
554
555 IF (yk .NE. 0.0) THEN
556 dz_z = max( dz_z, dyk / yk )
557 ELSE IF ( dyk .NE. 0.0 ) THEN
558 dz_z = hugeval
559 END IF
560
561 ymin = min( ymin, yk )
562
563 normy = max( normy, yk )
564
565 IF ( colequ ) THEN
566 normx = max( normx, yk * c( i ) )
567 normdx = max(normdx, dyk * c(i))
568 ELSE
569 normx = normy
570 normdx = max( normdx, dyk )
571 END IF
572 END DO
573
574 IF ( normx .NE. 0.0 ) THEN
575 dx_x = normdx / normx
576 ELSE IF ( normdx .EQ. 0.0 ) THEN
577 dx_x = 0.0
578 ELSE
579 dx_x = hugeval
580 END IF
581
582 dxrat = normdx / prevnormdx
583 dzrat = dz_z / prev_dz_z
584*
585* Check termination criteria.
586*
587 IF (.NOT.ignore_cwise
588 $ .AND. ymin*rcond .LT. incr_thresh*normy
589 $ .AND. y_prec_state .LT. extra_y )
590 $ incr_prec = .true.
591
592 IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
593 $ x_state = working_state
594 IF ( x_state .EQ. working_state ) THEN
595 IF ( dx_x .LE. eps ) THEN
596 x_state = conv_state
597 ELSE IF ( dxrat .GT. rthresh ) THEN
598 IF ( y_prec_state .NE. extra_y ) THEN
599 incr_prec = .true.
600 ELSE
601 x_state = noprog_state
602 END IF
603 ELSE
604 IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
605 END IF
606 IF ( x_state .GT. working_state ) final_dx_x = dx_x
607 END IF
608
609 IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
610 $ z_state = working_state
611 IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
612 $ z_state = working_state
613 IF ( z_state .EQ. working_state ) THEN
614 IF ( dz_z .LE. eps ) THEN
615 z_state = conv_state
616 ELSE IF ( dz_z .GT. dz_ub ) THEN
617 z_state = unstable_state
618 dzratmax = 0.0
619 final_dz_z = hugeval
620 ELSE IF ( dzrat .GT. rthresh ) THEN
621 IF ( y_prec_state .NE. extra_y ) THEN
622 incr_prec = .true.
623 ELSE
624 z_state = noprog_state
625 END IF
626 ELSE
627 IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
628 END IF
629 IF ( z_state .GT. working_state ) final_dz_z = dz_z
630 END IF
631*
632* Exit if both normwise and componentwise stopped working,
633* but if componentwise is unstable, let it go at least two
634* iterations.
635*
636 IF ( x_state.NE.working_state ) THEN
637 IF ( ignore_cwise ) GOTO 666
638 IF ( z_state.EQ.noprog_state .OR. z_state.EQ.conv_state )
639 $ GOTO 666
640 IF ( z_state.EQ.unstable_state .AND. cnt.GT.1 ) GOTO 666
641 END IF
642
643 IF ( incr_prec ) THEN
644 incr_prec = .false.
645 y_prec_state = y_prec_state + 1
646 DO i = 1, n
647 y_tail( i ) = 0.0
648 END DO
649 END IF
650
651 prevnormdx = normdx
652 prev_dz_z = dz_z
653*
654* Update soluton.
655*
656 IF ( y_prec_state .LT. extra_y ) THEN
657 CALL caxpy( n, (1.0e+0,0.0e+0), dy, 1, y(1,j), 1 )
658 ELSE
659 CALL cla_wwaddw( n, y(1,j), y_tail, dy )
660 END IF
661
662 END DO
663* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
664 666 CONTINUE
665*
666* Set final_* when cnt hits ithresh.
667*
668 IF ( x_state .EQ. working_state ) final_dx_x = dx_x
669 IF ( z_state .EQ. working_state ) final_dz_z = dz_z
670*
671* Compute error bounds.
672*
673 IF ( n_norms .GE. 1 ) THEN
674 err_bnds_norm( j, la_linrx_err_i ) =
675 $ final_dx_x / (1 - dxratmax)
676 END IF
677 IF ( n_norms .GE. 2 ) THEN
678 err_bnds_comp( j, la_linrx_err_i ) =
679 $ final_dz_z / (1 - dzratmax)
680 END IF
681*
682* Compute componentwise relative backward error from formula
683* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
684* where abs(Z) is the componentwise absolute value of the matrix
685* or vector Z.
686*
687* Compute residual RES = B_s - op(A_s) * Y,
688* op(A) = A, A**T, or A**H depending on TRANS (and type).
689*
690 CALL ccopy( n, b( 1, j ), 1, res, 1 )
691 CALL cgbmv( trans, n, n, kl, ku, (-1.0e+0,0.0e+0), ab, ldab,
692 $ y(1,j), 1, (1.0e+0,0.0e+0), res, 1 )
693
694 DO i = 1, n
695 ayb( i ) = cabs1( b( i, j ) )
696 END DO
697*
698* Compute abs(op(A_s))*abs(Y) + abs(B_s).
699*
700 CALL cla_gbamv( trans_type, n, n, kl, ku, 1.0e+0,
701 $ ab, ldab, y(1, j), 1, 1.0e+0, ayb, 1 )
702
703 CALL cla_lin_berr( n, n, 1, res, ayb, berr_out( j ) )
704*
705* End of loop for each RHS.
706*
707 END DO
708*
709 RETURN
710*
711* End of CLA_GBRFSX_EXTENDED
712*
chla_transtype
character *1 function chla_transtype(trans)
CHLA_TRANSTYPE
Definition chla_transtype.f:58
cla_gbamv
subroutine cla_gbamv(trans, m, n, kl, ku, alpha, ab, ldab, x, incx, beta, y, incy)
CLA_GBAMV performs a matrix-vector operation to calculate error bounds.
Definition cla_gbamv.f:186
cla_lin_berr
subroutine cla_lin_berr(n, nz, nrhs, res, ayb, berr)
CLA_LIN_BERR computes a component-wise relative backward error.
Definition cla_lin_berr.f:101
cla_wwaddw
subroutine cla_wwaddw(n, x, y, w)
CLA_WWADDW adds a vector into a doubled-single vector.
Definition cla_wwaddw.f:81

◆ cla_gbrpvgrw()

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

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

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

Purpose:
!>
!> CLA_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 COMPLEX 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 COMPLEX array, dimension (LDAFB,N)
!>     Details of the LU factorization of the band matrix A, as
!>     computed by CGBTRF.  U is stored as an upper triangular
!>     band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
!>     and the multipliers used during the factorization are stored
!>     in rows KL+KU+2 to 2*KL+KU+1.
!>
[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 cla_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 COMPLEX AB( LDAB, * ), AFB( LDAFB, * )
127* ..
128*
129* =====================================================================
130*
131* .. Local Scalars ..
132 INTEGER I, J, KD
133 REAL AMAX, UMAX, RPVGRW
134 COMPLEX ZDUM
135* ..
136* .. Intrinsic Functions ..
137 INTRINSIC abs, max, min, real, aimag
138* ..
139* .. Statement Functions ..
140 REAL CABS1
141* ..
142* .. Statement Function Definitions ..
143 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
144* ..
145* .. Executable Statements ..
146*
147 rpvgrw = 1.0
148
149 kd = ku + 1
150 DO j = 1, ncols
151 amax = 0.0
152 umax = 0.0
153 DO i = max( j-ku, 1 ), min( j+kl, n )
154 amax = max( cabs1( ab( kd+i-j, j ) ), amax )
155 END DO
156 DO i = max( j-ku, 1 ), j
157 umax = max( cabs1( afb( kd+i-j, j ) ), umax )
158 END DO
159 IF ( umax /= 0.0 ) THEN
160 rpvgrw = min( amax / umax, rpvgrw )
161 END IF
162 END DO
163 cla_gbrpvgrw = rpvgrw
164*
165* End of CLA_GBRPVGRW
166*
cla_gbrpvgrw
real function cla_gbrpvgrw(n, kl, ku, ncols, ab, ldab, afb, ldafb)
CLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix.
Definition cla_gbrpvgrw.f:117

◆ cungbr()

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

CUNGBR

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

Purpose:
!>
!> CUNGBR generates one of the complex unitary matrices Q or P**H
!> determined by CGEBRD when reducing a complex matrix A to bidiagonal
!> form: A = Q * B * P**H.  Q and P**H 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 CUNGBR returns the first n
!> columns of Q, where m >= n >= k;
!> if m < k, Q = H(1) H(2) . . . H(m-1) and CUNGBR returns Q as an
!> M-by-M matrix.
!>
!> If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**H
!> is of order N:
!> if k < n, P**H = G(k) . . . G(2) G(1) and CUNGBR returns the first m
!> rows of P**H, where n >= m >= k;
!> if k >= n, P**H = G(n-1) . . . G(2) G(1) and CUNGBR returns P**H as
!> an N-by-N matrix.
!>
Parameters
[in]VECT
!>          VECT is CHARACTER*1
!>          Specifies whether the matrix Q or the matrix P**H is
!>          required, as defined in the transformation applied by CGEBRD:
!>          = 'Q':  generate Q;
!>          = 'P':  generate P**H.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix Q or P**H to be returned.
!>          M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix Q or P**H 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 CGEBRD.
!>          If VECT = 'P', the number of rows in the original K-by-N
!>          matrix reduced by CGEBRD.
!>          K >= 0.
!>
[in,out]A
!>          A is COMPLEX array, dimension (LDA,N)
!>          On entry, the vectors which define the elementary reflectors,
!>          as returned by CGEBRD.
!>          On exit, the M-by-N matrix Q or P**H.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= M.
!>
[in]TAU
!>          TAU is COMPLEX 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**H, as
!>          returned by CGEBRD in its array argument TAUQ or TAUP.
!>
[out]WORK
!>          WORK is COMPLEX 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 cungbr.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 COMPLEX A( LDA, * ), TAU( * ), WORK( * )
168* ..
169*
170* =====================================================================
171*
172* .. Parameters ..
173 COMPLEX ZERO, ONE
174 parameter( zero = ( 0.0e+0, 0.0e+0 ),
175 $ one = ( 1.0e+0, 0.0e+0 ) )
176* ..
177* .. Local Scalars ..
178 LOGICAL LQUERY, WANTQ
179 INTEGER I, IINFO, J, LWKOPT, MN
180* ..
181* .. External Functions ..
182 LOGICAL LSAME
183 EXTERNAL lsame
184* ..
185* .. External Subroutines ..
186 EXTERNAL cunglq, cungqr, xerbla
187* ..
188* .. Intrinsic Functions ..
189 INTRINSIC max, min
190* ..
191* .. Executable Statements ..
192*
193* Test the input arguments
194*
195 info = 0
196 wantq = lsame( vect, 'Q' )
197 mn = min( m, n )
198 lquery = ( lwork.EQ.-1 )
199 IF( .NOT.wantq .AND. .NOT.lsame( vect, 'P' ) ) THEN
200 info = -1
201 ELSE IF( m.LT.0 ) THEN
202 info = -2
203 ELSE IF( n.LT.0 .OR. ( wantq .AND. ( n.GT.m .OR. n.LT.min( m,
204 $ k ) ) ) .OR. ( .NOT.wantq .AND. ( m.GT.n .OR. m.LT.
205 $ min( n, k ) ) ) ) THEN
206 info = -3
207 ELSE IF( k.LT.0 ) THEN
208 info = -4
209 ELSE IF( lda.LT.max( 1, m ) ) THEN
210 info = -6
211 ELSE IF( lwork.LT.max( 1, mn ) .AND. .NOT.lquery ) THEN
212 info = -9
213 END IF
214*
215 IF( info.EQ.0 ) THEN
216 work( 1 ) = 1
217 IF( wantq ) THEN
218 IF( m.GE.k ) THEN
219 CALL cungqr( m, n, k, a, lda, tau, work, -1, iinfo )
220 ELSE
221 IF( m.GT.1 ) THEN
222 CALL cungqr( m-1, m-1, m-1, a, lda, tau, work, -1,
223 $ iinfo )
224 END IF
225 END IF
226 ELSE
227 IF( k.LT.n ) THEN
228 CALL cunglq( m, n, k, a, lda, tau, work, -1, iinfo )
229 ELSE
230 IF( n.GT.1 ) THEN
231 CALL cunglq( n-1, n-1, n-1, a, lda, tau, work, -1,
232 $ iinfo )
233 END IF
234 END IF
235 END IF
236 lwkopt = real( work( 1 ) )
237 lwkopt = max(lwkopt, mn)
238 END IF
239*
240 IF( info.NE.0 ) THEN
241 CALL xerbla( 'CUNGBR', -info )
242 RETURN
243 ELSE IF( lquery ) THEN
244 work( 1 ) = lwkopt
245 RETURN
246 END IF
247*
248* Quick return if possible
249*
250 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
251 work( 1 ) = 1
252 RETURN
253 END IF
254*
255 IF( wantq ) THEN
256*
257* Form Q, determined by a call to CGEBRD to reduce an m-by-k
258* matrix
259*
260 IF( m.GE.k ) THEN
261*
262* If m >= k, assume m >= n >= k
263*
264 CALL cungqr( m, n, k, a, lda, tau, work, lwork, iinfo )
265*
266 ELSE
267*
268* If m < k, assume m = n
269*
270* Shift the vectors which define the elementary reflectors one
271* column to the right, and set the first row and column of Q
272* to those of the unit matrix
273*
274 DO 20 j = m, 2, -1
275 a( 1, j ) = zero
276 DO 10 i = j + 1, m
277 a( i, j ) = a( i, j-1 )
278 10 CONTINUE
279 20 CONTINUE
280 a( 1, 1 ) = one
281 DO 30 i = 2, m
282 a( i, 1 ) = zero
283 30 CONTINUE
284 IF( m.GT.1 ) THEN
285*
286* Form Q(2:m,2:m)
287*
288 CALL cungqr( m-1, m-1, m-1, a( 2, 2 ), lda, tau, work,
289 $ lwork, iinfo )
290 END IF
291 END IF
292 ELSE
293*
294* Form P**H, determined by a call to CGEBRD to reduce a k-by-n
295* matrix
296*
297 IF( k.LT.n ) THEN
298*
299* If k < n, assume k <= m <= n
300*
301 CALL cunglq( m, n, k, a, lda, tau, work, lwork, iinfo )
302*
303 ELSE
304*
305* If k >= n, assume m = n
306*
307* Shift the vectors which define the elementary reflectors one
308* row downward, and set the first row and column of P**H to
309* those of the unit matrix
310*
311 a( 1, 1 ) = one
312 DO 40 i = 2, n
313 a( i, 1 ) = zero
314 40 CONTINUE
315 DO 60 j = 2, n
316 DO 50 i = j - 1, 2, -1
317 a( i, j ) = a( i-1, j )
318 50 CONTINUE
319 a( 1, j ) = zero
320 60 CONTINUE
321 IF( n.GT.1 ) THEN
322*
323* Form P**H(2:n,2:n)
324*
325 CALL cunglq( n-1, n-1, n-1, a( 2, 2 ), lda, tau, work,
326 $ lwork, iinfo )
327 END IF
328 END IF
329 END IF
330 work( 1 ) = lwkopt
331 RETURN
332*
333* End of CUNGBR
334*
cunglq
subroutine cunglq(m, n, k, a, lda, tau, work, lwork, info)
CUNGLQ
Definition cunglq.f:127
cungqr
subroutine cungqr(m, n, k, a, lda, tau, work, lwork, info)
CUNGQR
Definition cungqr.f:128