LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
Collaboration diagram for complex16:

Functions

subroutine zgbbrd (VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q, LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO)
 ZGBBRD More...
 
subroutine zgbcon (NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND, WORK, RWORK, INFO)
 ZGBCON More...
 
subroutine zgbequ (M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, AMAX, INFO)
 ZGBEQU More...
 
subroutine zgbequb (M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, AMAX, INFO)
 ZGBEQUB More...
 
subroutine zgbrfs (TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
 ZGBRFS More...
 
subroutine zgbrfsx (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)
 ZGBRFSX More...
 
subroutine zgbtf2 (M, N, KL, KU, AB, LDAB, IPIV, INFO)
 ZGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algorithm. More...
 
subroutine zgbtrf (M, N, KL, KU, AB, LDAB, IPIV, INFO)
 ZGBTRF More...
 
subroutine zgbtrs (TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
 ZGBTRS More...
 
subroutine zggbak (JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
 ZGGBAK More...
 
subroutine zggbal (JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
 ZGGBAL More...
 
subroutine zla_gbamv (TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X, INCX, BETA, Y, INCY)
 ZLA_GBAMV performs a matrix-vector operation to calculate error bounds. More...
 
double precision function zla_gbrcond_c (TRANS, N, KL, KU, AB, LDAB, AFB, LDAFB, IPIV, C, CAPPLY, INFO, WORK, RWORK)
 ZLA_GBRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general banded matrices. More...
 
double precision function zla_gbrcond_x (TRANS, N, KL, KU, AB, LDAB, AFB, LDAFB, IPIV, X, INFO, WORK, RWORK)
 ZLA_GBRCOND_X computes the infinity norm condition number of op(A)*diag(x) for general banded matrices. More...
 
subroutine zla_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)
 ZLA_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. More...
 
double precision function zla_gbrpvgrw (N, KL, KU, NCOLS, AB, LDAB, AFB, LDAFB)
 ZLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix. More...
 
subroutine zungbr (VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
 ZUNGBR More...
 

Detailed Description

This is the group of complex16 computational functions for GB matrices

Function Documentation

subroutine zgbbrd ( character  VECT,
integer  M,
integer  N,
integer  NCC,
integer  KL,
integer  KU,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
complex*16, dimension( ldq, * )  Q,
integer  LDQ,
complex*16, dimension( ldpt, * )  PT,
integer  LDPT,
complex*16, dimension( ldc, * )  C,
integer  LDC,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGBBRD

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

Purpose:
 ZGBBRD 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*16 array, dimension (LDAB,N)
          On entry, the m-by-n band matrix A, stored in rows 1 to
          KL+KU+1. The j-th column of A is stored in the j-th column of
          the array AB as follows:
          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
          On exit, A is overwritten by values generated during the
          reduction.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array A. LDAB >= KL+KU+1.
[out]D
          D is DOUBLE PRECISION array, dimension (min(M,N))
          The diagonal elements of the bidiagonal matrix B.
[out]E
          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
          The superdiagonal elements of the bidiagonal matrix B.
[out]Q
          Q is COMPLEX*16 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*16 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*16 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*16 array, dimension (max(M,N))
[out]RWORK
          RWORK is DOUBLE PRECISION 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.
Date
November 2011

Definition at line 195 of file zgbbrd.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgbcon ( character  NORM,
integer  N,
integer  KL,
integer  KU,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
integer, dimension( * )  IPIV,
double precision  ANORM,
double precision  RCOND,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGBCON

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

Purpose:
 ZGBCON 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 ZGBTRF.

 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*16 array, dimension (LDAB,N)
          Details of the LU factorization of the band matrix A, as
          computed by ZGBTRF.  U is stored as an upper triangular band
          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
          the multipliers used during the factorization are stored in
          rows KL+KU+2 to 2*KL+KU+1.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          The pivot indices; for 1 <= i <= N, row i of the matrix was
          interchanged with row IPIV(i).
[in]ANORM
          ANORM is DOUBLE PRECISION
          If NORM = '1' or 'O', the 1-norm of the original matrix A.
          If NORM = 'I', the infinity-norm of the original matrix A.
[out]RCOND
          RCOND is DOUBLE PRECISION
          The reciprocal of the condition number of the matrix A,
          computed as RCOND = 1/(norm(A) * norm(inv(A))).
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION 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.
Date
November 2011

Definition at line 149 of file zgbcon.f.

149 *
150 * -- LAPACK computational routine (version 3.4.0) --
151 * -- LAPACK is a software package provided by Univ. of Tennessee, --
152 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
153 * November 2011
154 *
155 * .. Scalar Arguments ..
156  CHARACTER norm
157  INTEGER info, kl, ku, ldab, n
158  DOUBLE PRECISION anorm, rcond
159 * ..
160 * .. Array Arguments ..
161  INTEGER ipiv( * )
162  DOUBLE PRECISION rwork( * )
163  COMPLEX*16 ab( ldab, * ), work( * )
164 * ..
165 *
166 * =====================================================================
167 *
168 * .. Parameters ..
169  DOUBLE PRECISION one, zero
170  parameter( one = 1.0d+0, zero = 0.0d+0 )
171 * ..
172 * .. Local Scalars ..
173  LOGICAL lnoti, onenrm
174  CHARACTER normin
175  INTEGER ix, j, jp, kase, kase1, kd, lm
176  DOUBLE PRECISION ainvnm, scale, smlnum
177  COMPLEX*16 t, zdum
178 * ..
179 * .. Local Arrays ..
180  INTEGER isave( 3 )
181 * ..
182 * .. External Functions ..
183  LOGICAL lsame
184  INTEGER izamax
185  DOUBLE PRECISION dlamch
186  COMPLEX*16 zdotc
187  EXTERNAL lsame, izamax, dlamch, zdotc
188 * ..
189 * .. External Subroutines ..
190  EXTERNAL xerbla, zaxpy, zdrscl, zlacn2, zlatbs
191 * ..
192 * .. Intrinsic Functions ..
193  INTRINSIC abs, dble, dimag, min
194 * ..
195 * .. Statement Functions ..
196  DOUBLE PRECISION cabs1
197 * ..
198 * .. Statement Function definitions ..
199  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
200 * ..
201 * .. Executable Statements ..
202 *
203 * Test the input parameters.
204 *
205  info = 0
206  onenrm = norm.EQ.'1' .OR. lsame( norm, 'O' )
207  IF( .NOT.onenrm .AND. .NOT.lsame( norm, 'I' ) ) THEN
208  info = -1
209  ELSE IF( n.LT.0 ) THEN
210  info = -2
211  ELSE IF( kl.LT.0 ) THEN
212  info = -3
213  ELSE IF( ku.LT.0 ) THEN
214  info = -4
215  ELSE IF( ldab.LT.2*kl+ku+1 ) THEN
216  info = -6
217  ELSE IF( anorm.LT.zero ) THEN
218  info = -8
219  END IF
220  IF( info.NE.0 ) THEN
221  CALL xerbla( 'ZGBCON', -info )
222  RETURN
223  END IF
224 *
225 * Quick return if possible
226 *
227  rcond = zero
228  IF( n.EQ.0 ) THEN
229  rcond = one
230  RETURN
231  ELSE IF( anorm.EQ.zero ) THEN
232  RETURN
233  END IF
234 *
235  smlnum = dlamch( 'Safe minimum' )
236 *
237 * Estimate the norm of inv(A).
238 *
239  ainvnm = zero
240  normin = 'N'
241  IF( onenrm ) THEN
242  kase1 = 1
243  ELSE
244  kase1 = 2
245  END IF
246  kd = kl + ku + 1
247  lnoti = kl.GT.0
248  kase = 0
249  10 CONTINUE
250  CALL zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
251  IF( kase.NE.0 ) THEN
252  IF( kase.EQ.kase1 ) THEN
253 *
254 * Multiply by inv(L).
255 *
256  IF( lnoti ) THEN
257  DO 20 j = 1, n - 1
258  lm = min( kl, n-j )
259  jp = ipiv( j )
260  t = work( jp )
261  IF( jp.NE.j ) THEN
262  work( jp ) = work( j )
263  work( j ) = t
264  END IF
265  CALL zaxpy( lm, -t, ab( kd+1, j ), 1, work( j+1 ), 1 )
266  20 CONTINUE
267  END IF
268 *
269 * Multiply by inv(U).
270 *
271  CALL zlatbs( 'Upper', 'No transpose', 'Non-unit', normin, n,
272  $ kl+ku, ab, ldab, work, scale, rwork, info )
273  ELSE
274 *
275 * Multiply by inv(U**H).
276 *
277  CALL zlatbs( 'Upper', 'Conjugate transpose', 'Non-unit',
278  $ normin, n, kl+ku, ab, ldab, work, scale, rwork,
279  $ info )
280 *
281 * Multiply by inv(L**H).
282 *
283  IF( lnoti ) THEN
284  DO 30 j = n - 1, 1, -1
285  lm = min( kl, n-j )
286  work( j ) = work( j ) - zdotc( lm, ab( kd+1, j ), 1,
287  $ work( j+1 ), 1 )
288  jp = ipiv( j )
289  IF( jp.NE.j ) THEN
290  t = work( jp )
291  work( jp ) = work( j )
292  work( j ) = t
293  END IF
294  30 CONTINUE
295  END IF
296  END IF
297 *
298 * Divide X by 1/SCALE if doing so will not cause overflow.
299 *
300  normin = 'Y'
301  IF( scale.NE.one ) THEN
302  ix = izamax( n, work, 1 )
303  IF( scale.LT.cabs1( work( ix ) )*smlnum .OR. scale.EQ.zero )
304  $ GO TO 40
305  CALL zdrscl( n, scale, work, 1 )
306  END IF
307  GO TO 10
308  END IF
309 *
310 * Compute the estimate of the reciprocal condition number.
311 *
312  IF( ainvnm.NE.zero )
313  $ rcond = ( one / ainvnm ) / anorm
314 *
315  40 CONTINUE
316  RETURN
317 *
318 * End of ZGBCON
319 *
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:135
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53
subroutine zlatbs(UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, SCALE, CNORM, INFO)
ZLATBS solves a triangular banded system of equations.
Definition: zlatbs.f:245
subroutine zdrscl(N, SA, SX, INCX)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: zdrscl.f:86
complex *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
Definition: zdotc.f:54

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgbequ ( integer  M,
integer  N,
integer  KL,
integer  KU,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
double precision, dimension( * )  R,
double precision, dimension( * )  C,
double precision  ROWCND,
double precision  COLCND,
double precision  AMAX,
integer  INFO 
)

ZGBEQU

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

Purpose:
 ZGBEQU 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*16 array, dimension (LDAB,N)
          The band matrix A, stored in rows 1 to KL+KU+1.  The j-th
          column of A is stored in the j-th column of the array AB as
          follows:
          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KL+KU+1.
[out]R
          R is DOUBLE PRECISION array, dimension (M)
          If INFO = 0, or INFO > M, R contains the row scale factors
          for A.
[out]C
          C is DOUBLE PRECISION array, dimension (N)
          If INFO = 0, C contains the column scale factors for A.
[out]ROWCND
          ROWCND is DOUBLE PRECISION
          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
          AMAX is neither too large nor too small, it is not worth
          scaling by R.
[out]COLCND
          COLCND is DOUBLE PRECISION
          If INFO = 0, COLCND contains the ratio of the smallest
          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
          worth scaling by C.
[out]AMAX
          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix element.  If AMAX is very
          close to overflow or very close to underflow, the matrix
          should be scaled.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, and i is
                <= M:  the i-th row of A is exactly zero
                >  M:  the (i-M)-th column of A is exactly zero
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 156 of file zgbequ.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgbequb ( integer  M,
integer  N,
integer  KL,
integer  KU,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
double precision, dimension( * )  R,
double precision, dimension( * )  C,
double precision  ROWCND,
double precision  COLCND,
double precision  AMAX,
integer  INFO 
)

ZGBEQUB

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

Purpose:
 ZGBEQUB 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 ZGEEQU by restricting the scaling factors
 to a power of the radix.  Baring over- and underflow, scaling by
 these factors introduces no additional rounding errors.  However, the
 scaled entries' magnitured are no longer approximately 1 but lie
 between sqrt(radix) and 1/sqrt(radix).
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in]KL
          KL is INTEGER
          The number of subdiagonals within the band of A.  KL >= 0.
[in]KU
          KU is INTEGER
          The number of superdiagonals within the band of A.  KU >= 0.
[in]AB
          AB is DOUBLE PRECISION array, dimension (LDAB,N)
          On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
          The j-th column of A is stored in the j-th column of the
          array AB as follows:
          AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array A.  LDAB >= max(1,M).
[out]R
          R is DOUBLE PRECISION array, dimension (M)
          If INFO = 0 or INFO > M, R contains the row scale factors
          for A.
[out]C
          C is DOUBLE PRECISION array, dimension (N)
          If INFO = 0,  C contains the column scale factors for A.
[out]ROWCND
          ROWCND is DOUBLE PRECISION
          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
          AMAX is neither too large nor too small, it is not worth
          scaling by R.
[out]COLCND
          COLCND is DOUBLE PRECISION
          If INFO = 0, COLCND contains the ratio of the smallest
          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
          worth scaling by C.
[out]AMAX
          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix element.  If AMAX is very
          close to overflow or very close to underflow, the matrix
          should be scaled.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i,  and i is
                <= M:  the i-th row of A is exactly zero
                >  M:  the (i-M)-th column of A is exactly zero
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 163 of file zgbequb.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgbrfs ( character  TRANS,
integer  N,
integer  KL,
integer  KU,
integer  NRHS,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
complex*16, dimension( ldafb, * )  AFB,
integer  LDAFB,
integer, dimension( * )  IPIV,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldx, * )  X,
integer  LDX,
double precision, dimension( * )  FERR,
double precision, dimension( * )  BERR,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGBRFS

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

Purpose:
 ZGBRFS 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*16 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*16 array, dimension (LDAFB,N)
          Details of the LU factorization of the band matrix A, as
          computed by ZGBTRF.  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 ZGBTRF; for 1<=i<=N, row i of the
          matrix was interchanged with row IPIV(i).
[in]B
          B is COMPLEX*16 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*16 array, dimension (LDX,NRHS)
          On entry, the solution matrix X, as computed by ZGBTRS.
          On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.  LDX >= max(1,N).
[out]FERR
          FERR is DOUBLE PRECISION array, dimension (NRHS)
          The estimated forward error bound for each solution vector
          X(j) (the j-th column of the solution matrix X).
          If XTRUE is the true solution corresponding to X(j), FERR(j)
          is an estimated upper bound for the magnitude of the largest
          element in (X(j) - XTRUE) divided by the magnitude of the
          largest element in X(j).  The estimate is as reliable as
          the estimate for RCOND, and is almost always a slight
          overestimate of the true error.
[out]BERR
          BERR is DOUBLE PRECISION array, dimension (NRHS)
          The componentwise relative backward error of each solution
          vector X(j) (i.e., the smallest relative change in
          any element of A or B that makes X(j) an exact solution).
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION 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.
Date
November 2011

Definition at line 208 of file zgbrfs.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgbrfsx ( character  TRANS,
character  EQUED,
integer  N,
integer  KL,
integer  KU,
integer  NRHS,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
complex*16, dimension( ldafb, * )  AFB,
integer  LDAFB,
integer, dimension( * )  IPIV,
double precision, dimension( * )  R,
double precision, dimension( * )  C,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldx , * )  X,
integer  LDX,
double precision  RCOND,
double precision, dimension( * )  BERR,
integer  N_ERR_BNDS,
double precision, dimension( nrhs, * )  ERR_BNDS_NORM,
double precision, dimension( nrhs, * )  ERR_BNDS_COMP,
integer  NPARAMS,
double precision, dimension( * )  PARAMS,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGBRFSX

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

Purpose:
    ZGBRFSX improves the computed solution to a system of linear
    equations and provides error bounds and backward error estimates
    for the solution.  In addition to normwise error bound, the code
    provides maximum componentwise error bound if possible.  See
    comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the
    error bounds.

    The original system of linear equations may have been equilibrated
    before calling this routine, as described by arguments EQUED, R
    and C below. In this case, the solution and error bounds returned
    are for the original unequilibrated system.
     Some optional parameters are bundled in the PARAMS array.  These
     settings determine how refinement is performed, but often the
     defaults are acceptable.  If the defaults are acceptable, users
     can pass NPARAMS = 0 which prevents the source code from accessing
     the PARAMS argument.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
     Specifies the form of the system of equations:
       = 'N':  A * X = B     (No transpose)
       = 'T':  A**T * X = B  (Transpose)
       = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
[in]EQUED
          EQUED is CHARACTER*1
     Specifies the form of equilibration that was done to A
     before calling this routine. This is needed to compute
     the solution and error bounds correctly.
       = 'N':  No equilibration
       = 'R':  Row equilibration, i.e., A has been premultiplied by
               diag(R).
       = 'C':  Column equilibration, i.e., A has been postmultiplied
               by diag(C).
       = 'B':  Both row and column equilibration, i.e., A has been
               replaced by diag(R) * A * diag(C).
               The right hand side B has been changed accordingly.
[in]N
          N is INTEGER
     The order of the matrix A.  N >= 0.
[in]KL
          KL is INTEGER
     The number of subdiagonals within the band of A.  KL >= 0.
[in]KU
          KU is INTEGER
     The number of superdiagonals within the band of A.  KU >= 0.
[in]NRHS
          NRHS is INTEGER
     The number of right hand sides, i.e., the number of columns
     of the matrices B and X.  NRHS >= 0.
[in]AB
          AB is COMPLEX*16 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*16 array, dimension (LDAFB,N)
     Details of the LU factorization of the band matrix A, as
     computed by DGBTRF.  U is stored as an upper triangular band
     matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
     the multipliers used during the factorization are stored in
     rows KL+KU+2 to 2*KL+KU+1.
[in]LDAFB
          LDAFB is INTEGER
     The leading dimension of the array AFB.  LDAFB >= 2*KL*KU+1.
[in]IPIV
          IPIV is INTEGER array, dimension (N)
     The pivot indices from DGETRF; for 1<=i<=N, row i of the
     matrix was interchanged with row IPIV(i).
[in,out]R
          R is DOUBLE PRECISION array, dimension (N)
     The row scale factors for A.  If EQUED = 'R' or 'B', A is
     multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
     is not accessed.  R is an input argument if FACT = 'F';
     otherwise, R is an output argument.  If FACT = 'F' and
     EQUED = 'R' or 'B', each element of R must be positive.
     If R is output, each element of R is a power of the radix.
     If R is input, each element of R should be a power of the radix
     to ensure a reliable solution and error estimates. Scaling by
     powers of the radix does not cause rounding errors unless the
     result underflows or overflows. Rounding errors during scaling
     lead to refining with a matrix that is not equivalent to the
     input matrix, producing error estimates that may not be
     reliable.
[in,out]C
          C is DOUBLE PRECISION array, dimension (N)
     The column scale factors for A.  If EQUED = 'C' or 'B', A is
     multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
     is not accessed.  C is an input argument if FACT = 'F';
     otherwise, C is an output argument.  If FACT = 'F' and
     EQUED = 'C' or 'B', each element of C must be positive.
     If C is output, each element of C is a power of the radix.
     If C is input, each element of C should be a power of the radix
     to ensure a reliable solution and error estimates. Scaling by
     powers of the radix does not cause rounding errors unless the
     result underflows or overflows. Rounding errors during scaling
     lead to refining with a matrix that is not equivalent to the
     input matrix, producing error estimates that may not be
     reliable.
[in]B
          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
     The right hand side matrix B.
[in]LDB
          LDB is INTEGER
     The leading dimension of the array B.  LDB >= max(1,N).
[in,out]X
          X is COMPLEX*16 array, dimension (LDX,NRHS)
     On entry, the solution matrix X, as computed by DGETRS.
     On exit, the improved solution matrix X.
[in]LDX
          LDX is INTEGER
     The leading dimension of the array X.  LDX >= max(1,N).
[out]RCOND
          RCOND is DOUBLE PRECISION
     Reciprocal scaled condition number.  This is an estimate of the
     reciprocal Skeel condition number of the matrix A after
     equilibration (if done).  If this is less than the machine
     precision (in particular, if it is zero), the matrix is singular
     to working precision.  Note that the error may still be small even
     if this number is very small and the matrix appears ill-
     conditioned.
[out]BERR
          BERR is COMPLEX*16 array, dimension (NRHS)
     Componentwise relative backward error.  This is the
     componentwise relative backward error of each solution vector X(j)
     (i.e., the smallest relative change in any element of A or B that
     makes X(j) an exact solution).
[in]N_ERR_BNDS
          N_ERR_BNDS is INTEGER
     Number of error bounds to return for each right hand side
     and each type (normwise or componentwise).  See ERR_BNDS_NORM and
     ERR_BNDS_COMP below.
[out]ERR_BNDS_NORM
          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     normwise relative error, which is defined as follows:

     Normwise relative error in the ith solution vector:
             max_j (abs(XTRUE(j,i) - X(j,i)))
            ------------------------------
                  max_j abs(X(j,i))

     The array is indexed by the type of error information as described
     below. There currently are up to three pieces of information
     returned.

     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_NORM(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * dlamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * dlamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated normwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * dlamch('Epsilon') to determine if the error
              estimate is "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*A, where S scales each row by a power of the
              radix so all absolute row sums of Z are approximately 1.

     See Lapack Working Note 165 for further details and extra
     cautions.
[out]ERR_BNDS_COMP
          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     componentwise relative error, which is defined as follows:

     Componentwise relative error in the ith solution vector:
                    abs(XTRUE(j,i) - X(j,i))
             max_j ----------------------
                         abs(X(j,i))

     The array is indexed by the right-hand side i (on which the
     componentwise relative error depends), and the type of error
     information as described below. There currently are up to three
     pieces of information returned for each right-hand side. If
     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 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 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * dlamch('Epsilon').

     err = 2 "Guaranteed" error bound: The estimated forward error,
              almost certainly within a factor of 10 of the true error
              so long as the next entry is greater than the threshold
              sqrt(n) * dlamch('Epsilon'). This error bound should only
              be trusted if the previous boolean is true.

     err = 3  Reciprocal condition number: Estimated componentwise
              reciprocal condition number.  Compared with the threshold
              sqrt(n) * dlamch('Epsilon') to determine if the error
              estimate is "guaranteed". 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 .LE. 0, the
     PARAMS array is never referenced and default values are used.
[in,out]PARAMS
          PARAMS is DOUBLE PRECISION array, dimension NPARAMS
     Specifies algorithm parameters.  If an entry is .LT. 0.0, then
     that entry will be filled with default value used for that
     parameter.  Only positions up to NPARAMS are accessed; defaults
     are used for higher-numbered parameters.

       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
            refinement or not.
         Default: 1.0D+0
            = 0.0 : No refinement is performed, and no error bounds are
                    computed.
            = 1.0 : Use the double-precision refinement algorithm,
                    possibly with doubled-single computations if the
                    compilation environment does not support DOUBLE
                    PRECISION.
              (other values are reserved for future use)

       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
            computations allowed for refinement.
         Default: 10
         Aggressive: Set to 100 to permit convergence using approximate
                     factorizations or factorizations other than LU. If
                     the factorization uses a technique other than
                     Gaussian elimination, the guarantees in
                     err_bnds_norm and err_bnds_comp may no longer be
                     trustworthy.

       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
            will attempt to find a solution with small componentwise
            relative error in the double-precision algorithm.  Positive
            is true, 0.0 is false.
         Default: 1.0 (attempt componentwise convergence)
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]RWORK
          RWORK is DOUBLE PRECISION 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.
Date
April 2012

Definition at line 442 of file zgbrfsx.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgbtf2 ( integer  M,
integer  N,
integer  KL,
integer  KU,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
integer, dimension( * )  IPIV,
integer  INFO 
)

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

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

Purpose:
 ZGBTF2 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*16 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.
Date
September 2012
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 147 of file zgbtf2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgbtrf ( integer  M,
integer  N,
integer  KL,
integer  KU,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
integer, dimension( * )  IPIV,
integer  INFO 
)

ZGBTRF

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

Purpose:
 ZGBTRF 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*16 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.
Date
November 2011
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 146 of file zgbtrf.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zgbtrs ( character  TRANS,
integer  N,
integer  KL,
integer  KU,
integer  NRHS,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
integer, dimension( * )  IPIV,
complex*16, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)

ZGBTRS

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

Purpose:
 ZGBTRS 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 ZGBTRF.
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*16 array, dimension (LDAB,N)
          Details of the LU factorization of the band matrix A, as
          computed by ZGBTRF.  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*16 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.
Date
November 2011

Definition at line 140 of file zgbtrs.f.

140 *
141 * -- LAPACK computational routine (version 3.4.0) --
142 * -- LAPACK is a software package provided by Univ. of Tennessee, --
143 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144 * November 2011
145 *
146 * .. Scalar Arguments ..
147  CHARACTER trans
148  INTEGER info, kl, ku, ldab, ldb, n, nrhs
149 * ..
150 * .. Array Arguments ..
151  INTEGER ipiv( * )
152  COMPLEX*16 ab( ldab, * ), b( ldb, * )
153 * ..
154 *
155 * =====================================================================
156 *
157 * .. Parameters ..
158  COMPLEX*16 one
159  parameter( one = ( 1.0d+0, 0.0d+0 ) )
160 * ..
161 * .. Local Scalars ..
162  LOGICAL lnoti, notran
163  INTEGER i, j, kd, l, lm
164 * ..
165 * .. External Functions ..
166  LOGICAL lsame
167  EXTERNAL lsame
168 * ..
169 * .. External Subroutines ..
170  EXTERNAL xerbla, zgemv, zgeru, zlacgv, zswap, ztbsv
171 * ..
172 * .. Intrinsic Functions ..
173  INTRINSIC max, min
174 * ..
175 * .. Executable Statements ..
176 *
177 * Test the input parameters.
178 *
179  info = 0
180  notran = lsame( trans, 'N' )
181  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
182  $ lsame( trans, 'C' ) ) THEN
183  info = -1
184  ELSE IF( n.LT.0 ) THEN
185  info = -2
186  ELSE IF( kl.LT.0 ) THEN
187  info = -3
188  ELSE IF( ku.LT.0 ) THEN
189  info = -4
190  ELSE IF( nrhs.LT.0 ) THEN
191  info = -5
192  ELSE IF( ldab.LT.( 2*kl+ku+1 ) ) THEN
193  info = -7
194  ELSE IF( ldb.LT.max( 1, n ) ) THEN
195  info = -10
196  END IF
197  IF( info.NE.0 ) THEN
198  CALL xerbla( 'ZGBTRS', -info )
199  RETURN
200  END IF
201 *
202 * Quick return if possible
203 *
204  IF( n.EQ.0 .OR. nrhs.EQ.0 )
205  $ RETURN
206 *
207  kd = ku + kl + 1
208  lnoti = kl.GT.0
209 *
210  IF( notran ) THEN
211 *
212 * Solve A*X = B.
213 *
214 * Solve L*X = B, overwriting B with X.
215 *
216 * L is represented as a product of permutations and unit lower
217 * triangular matrices L = P(1) * L(1) * ... * P(n-1) * L(n-1),
218 * where each transformation L(i) is a rank-one modification of
219 * the identity matrix.
220 *
221  IF( lnoti ) THEN
222  DO 10 j = 1, n - 1
223  lm = min( kl, n-j )
224  l = ipiv( j )
225  IF( l.NE.j )
226  $ CALL zswap( nrhs, b( l, 1 ), ldb, b( j, 1 ), ldb )
227  CALL zgeru( lm, nrhs, -one, ab( kd+1, j ), 1, b( j, 1 ),
228  $ ldb, b( j+1, 1 ), ldb )
229  10 CONTINUE
230  END IF
231 *
232  DO 20 i = 1, nrhs
233 *
234 * Solve U*X = B, overwriting B with X.
235 *
236  CALL ztbsv( 'Upper', 'No transpose', 'Non-unit', n, kl+ku,
237  $ ab, ldab, b( 1, i ), 1 )
238  20 CONTINUE
239 *
240  ELSE IF( lsame( trans, 'T' ) ) THEN
241 *
242 * Solve A**T * X = B.
243 *
244  DO 30 i = 1, nrhs
245 *
246 * Solve U**T * X = B, overwriting B with X.
247 *
248  CALL ztbsv( 'Upper', 'Transpose', 'Non-unit', n, kl+ku, ab,
249  $ ldab, b( 1, i ), 1 )
250  30 CONTINUE
251 *
252 * Solve L**T * X = B, overwriting B with X.
253 *
254  IF( lnoti ) THEN
255  DO 40 j = n - 1, 1, -1
256  lm = min( kl, n-j )
257  CALL zgemv( 'Transpose', lm, nrhs, -one, b( j+1, 1 ),
258  $ ldb, ab( kd+1, j ), 1, one, b( j, 1 ), ldb )
259  l = ipiv( j )
260  IF( l.NE.j )
261  $ CALL zswap( nrhs, b( l, 1 ), ldb, b( j, 1 ), ldb )
262  40 CONTINUE
263  END IF
264 *
265  ELSE
266 *
267 * Solve A**H * X = B.
268 *
269  DO 50 i = 1, nrhs
270 *
271 * Solve U**H * X = B, overwriting B with X.
272 *
273  CALL ztbsv( 'Upper', 'Conjugate transpose', 'Non-unit', n,
274  $ kl+ku, ab, ldab, b( 1, i ), 1 )
275  50 CONTINUE
276 *
277 * Solve L**H * X = B, overwriting B with X.
278 *
279  IF( lnoti ) THEN
280  DO 60 j = n - 1, 1, -1
281  lm = min( kl, n-j )
282  CALL zlacgv( nrhs, b( j, 1 ), ldb )
283  CALL zgemv( 'Conjugate transpose', lm, nrhs, -one,
284  $ b( j+1, 1 ), ldb, ab( kd+1, j ), 1, one,
285  $ b( j, 1 ), ldb )
286  CALL zlacgv( nrhs, b( j, 1 ), ldb )
287  l = ipiv( j )
288  IF( l.NE.j )
289  $ CALL zswap( nrhs, b( l, 1 ), ldb, b( j, 1 ), ldb )
290  60 CONTINUE
291  END IF
292  END IF
293  RETURN
294 *
295 * End of ZGBTRS
296 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:76
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERU
Definition: zgeru.f:132
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine ztbsv(UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX)
ZTBSV
Definition: ztbsv.f:191

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zggbak ( character  JOB,
character  SIDE,
integer  N,
integer  ILO,
integer  IHI,
double precision, dimension( * )  LSCALE,
double precision, dimension( * )  RSCALE,
integer  M,
complex*16, dimension( ldv, * )  V,
integer  LDV,
integer  INFO 
)

ZGGBAK

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

Purpose:
 ZGGBAK 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
 ZGGBAL.
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 ZGGBAL.
[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 ZGGBAL.
          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
[in]LSCALE
          LSCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and/or scaling factors applied
          to the left side of A and B, as returned by ZGGBAL.
[in]RSCALE
          RSCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and/or scaling factors applied
          to the right side of A and B, as returned by ZGGBAL.
[in]M
          M is INTEGER
          The number of columns of the matrix V.  M >= 0.
[in,out]V
          V is COMPLEX*16 array, dimension (LDV,M)
          On entry, the matrix of right or left eigenvectors to be
          transformed, as returned by ZTGEVC.
          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.
Date
November 2015
Further Details:
  See R.C. Ward, Balancing the generalized eigenvalue problem,
                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.

Definition at line 150 of file zggbak.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zggbal ( character  JOB,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
integer  ILO,
integer  IHI,
double precision, dimension( * )  LSCALE,
double precision, dimension( * )  RSCALE,
double precision, dimension( * )  WORK,
integer  INFO 
)

ZGGBAL

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

Purpose:
 ZGGBAL 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*16 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*16 array, dimension (LDB,N)
          On entry, the input matrix B.
          On exit, B is overwritten by the balanced matrix.
          If JOB = 'N', B is not referenced.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,N).
[out]ILO
          ILO is INTEGER
[out]IHI
          IHI is INTEGER
          ILO and IHI are set to integers such that on exit
          A(i,j) = 0 and B(i,j) = 0 if i > j and
          j = 1,...,ILO-1 or i = IHI+1,...,N.
          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
[out]LSCALE
          LSCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and scaling factors applied
          to the left side of A and B.  If P(j) is the index of the
          row interchanged with row j, and D(j) is the scaling factor
          applied to row j, then
            LSCALE(j) = P(j)    for J = 1,...,ILO-1
                      = D(j)    for J = ILO,...,IHI
                      = P(j)    for J = IHI+1,...,N.
          The order in which the interchanges are made is N to IHI+1,
          then 1 to ILO-1.
[out]RSCALE
          RSCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and scaling factors applied
          to the right side of A and B.  If P(j) is the index of the
          column interchanged with column j, and D(j) is the scaling
          factor applied to column j, then
            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.
Date
November 2015
Further Details:
  See R.C. WARD, Balancing the generalized eigenvalue problem,
                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.

Definition at line 179 of file zggbal.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zla_gbamv ( integer  TRANS,
integer  M,
integer  N,
integer  KL,
integer  KU,
double precision  ALPHA,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
complex*16, dimension( * )  X,
integer  INCX,
double precision  BETA,
double precision, dimension( * )  Y,
integer  INCY 
)

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

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

Purpose:
 ZLA_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,
 "symbolically" zero components are not perturbed.  A zero
 entry is considered "symbolic" if all multiplications involved
 in computing that entry have at least one zero multiplicand.
Parameters
[in]TRANS
          TRANS is INTEGER
           On entry, TRANS specifies the operation to be performed as
           follows:

             BLAS_NO_TRANS      y := alpha*abs(A)*abs(x) + beta*abs(y)
             BLAS_TRANS         y := alpha*abs(A**T)*abs(x) + beta*abs(y)
             BLAS_CONJ_TRANS    y := alpha*abs(A**T)*abs(x) + beta*abs(y)

           Unchanged on exit.
[in]M
          M is INTEGER
           On entry, M specifies the number of rows of the matrix A.
           M must be at least zero.
           Unchanged on exit.
[in]N
          N is INTEGER
           On entry, N specifies the number of columns of the matrix A.
           N must be at least zero.
           Unchanged on exit.
[in]KL
          KL is INTEGER
           The number of subdiagonals within the band of A.  KL >= 0.
[in]KU
          KU is INTEGER
           The number of superdiagonals within the band of A.  KU >= 0.
[in]ALPHA
          ALPHA is DOUBLE PRECISION
           On entry, ALPHA specifies the scalar alpha.
           Unchanged on exit.
[in]AB
          AB is COMPLEX*16 array of 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*16 array, dimension
           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
           and at least
           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
           Before entry, the incremented array X must contain the
           vector x.
           Unchanged on exit.
[in]INCX
          INCX is INTEGER
           On entry, INCX specifies the increment for the elements of
           X. INCX must not be zero.
           Unchanged on exit.
[in]BETA
          BETA is DOUBLE PRECISION
           On entry, BETA specifies the scalar beta. When BETA is
           supplied as zero then Y need not be set on input.
           Unchanged on exit.
[in,out]Y
          Y is DOUBLE PRECISION array, dimension
           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
           and at least
           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
           Before entry with BETA non-zero, the incremented array Y
           must contain the vector y. On exit, Y is overwritten by the
           updated vector y.
[in]INCY
          INCY is INTEGER
           On entry, INCY specifies the increment for the elements of
           Y. INCY must not be zero.
           Unchanged on exit.

  Level 2 Blas routine.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 188 of file zla_gbamv.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

double precision function zla_gbrcond_c ( character  TRANS,
integer  N,
integer  KL,
integer  KU,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
complex*16, dimension( ldafb, * )  AFB,
integer  LDAFB,
integer, dimension( * )  IPIV,
double precision, dimension( * )  C,
logical  CAPPLY,
integer  INFO,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK 
)

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

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

Purpose:
    ZLA_GBRCOND_C Computes the infinity norm condition number of
    op(A) * inv(diag(C)) where C is a DOUBLE PRECISION 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*16 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*16 array, dimension (LDAFB,N)
     Details of the LU factorization of the band matrix A, as
     computed by ZGBTRF.  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 ZGBTRF; row i of the matrix was interchanged
     with row IPIV(i).
[in]C
          C is DOUBLE PRECISION 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.
[in]WORK
          WORK is COMPLEX*16 array, dimension (2*N).
     Workspace.
[in]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N).
     Workspace.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 165 of file zla_gbrcond_c.f.

165 *
166 * -- LAPACK computational routine (version 3.4.2) --
167 * -- LAPACK is a software package provided by Univ. of Tennessee, --
168 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169 * September 2012
170 *
171 * .. Scalar Arguments ..
172  CHARACTER trans
173  LOGICAL capply
174  INTEGER n, kl, ku, kd, ke, ldab, ldafb, info
175 * ..
176 * .. Array Arguments ..
177  INTEGER ipiv( * )
178  COMPLEX*16 ab( ldab, * ), afb( ldafb, * ), work( * )
179  DOUBLE PRECISION c( * ), rwork( * )
180 *
181 *
182 * =====================================================================
183 *
184 * .. Local Scalars ..
185  LOGICAL notrans
186  INTEGER kase, i, j
187  DOUBLE PRECISION ainvnm, anorm, tmp
188  COMPLEX*16 zdum
189 * ..
190 * .. Local Arrays ..
191  INTEGER isave( 3 )
192 * ..
193 * .. External Functions ..
194  LOGICAL lsame
195  EXTERNAL lsame
196 * ..
197 * .. External Subroutines ..
198  EXTERNAL zlacn2, zgbtrs, xerbla
199 * ..
200 * .. Intrinsic Functions ..
201  INTRINSIC abs, max
202 * ..
203 * .. Statement Functions ..
204  DOUBLE PRECISION cabs1
205 * ..
206 * .. Statement Function Definitions ..
207  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
208 * ..
209 * .. Executable Statements ..
210  zla_gbrcond_c = 0.0d+0
211 *
212  info = 0
213  notrans = lsame( trans, 'N' )
214  IF ( .NOT. notrans .AND. .NOT. lsame( trans, 'T' ) .AND. .NOT.
215  $ lsame( trans, 'C' ) ) THEN
216  info = -1
217  ELSE IF( n.LT.0 ) THEN
218  info = -2
219  ELSE IF( kl.LT.0 .OR. kl.GT.n-1 ) THEN
220  info = -3
221  ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
222  info = -4
223  ELSE IF( ldab.LT.kl+ku+1 ) THEN
224  info = -6
225  ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
226  info = -8
227  END IF
228  IF( info.NE.0 ) THEN
229  CALL xerbla( 'ZLA_GBRCOND_C', -info )
230  RETURN
231  END IF
232 *
233 * Compute norm of op(A)*op2(C).
234 *
235  anorm = 0.0d+0
236  kd = ku + 1
237  ke = kl + 1
238  IF ( notrans ) THEN
239  DO i = 1, n
240  tmp = 0.0d+0
241  IF ( capply ) THEN
242  DO j = max( i-kl, 1 ), min( i+ku, n )
243  tmp = tmp + cabs1( ab( kd+i-j, j ) ) / c( j )
244  END DO
245  ELSE
246  DO j = max( i-kl, 1 ), min( i+ku, n )
247  tmp = tmp + cabs1( ab( kd+i-j, j ) )
248  END DO
249  END IF
250  rwork( i ) = tmp
251  anorm = max( anorm, tmp )
252  END DO
253  ELSE
254  DO i = 1, n
255  tmp = 0.0d+0
256  IF ( capply ) THEN
257  DO j = max( i-kl, 1 ), min( i+ku, n )
258  tmp = tmp + cabs1( ab( ke-i+j, i ) ) / c( j )
259  END DO
260  ELSE
261  DO j = max( i-kl, 1 ), min( i+ku, n )
262  tmp = tmp + cabs1( ab( ke-i+j, i ) )
263  END DO
264  END IF
265  rwork( i ) = tmp
266  anorm = max( anorm, tmp )
267  END DO
268  END IF
269 *
270 * Quick return if possible.
271 *
272  IF( n.EQ.0 ) THEN
273  zla_gbrcond_c = 1.0d+0
274  RETURN
275  ELSE IF( anorm .EQ. 0.0d+0 ) THEN
276  RETURN
277  END IF
278 *
279 * Estimate the norm of inv(op(A)).
280 *
281  ainvnm = 0.0d+0
282 *
283  kase = 0
284  10 CONTINUE
285  CALL zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
286  IF( kase.NE.0 ) THEN
287  IF( kase.EQ.2 ) THEN
288 *
289 * Multiply by R.
290 *
291  DO i = 1, n
292  work( i ) = work( i ) * rwork( i )
293  END DO
294 *
295  IF ( notrans ) THEN
296  CALL zgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
297  $ ipiv, work, n, info )
298  ELSE
299  CALL zgbtrs( 'Conjugate transpose', n, kl, ku, 1, afb,
300  $ ldafb, ipiv, work, n, info )
301  ENDIF
302 *
303 * Multiply by inv(C).
304 *
305  IF ( capply ) THEN
306  DO i = 1, n
307  work( i ) = work( i ) * c( i )
308  END DO
309  END IF
310  ELSE
311 *
312 * Multiply by inv(C**H).
313 *
314  IF ( capply ) THEN
315  DO i = 1, n
316  work( i ) = work( i ) * c( i )
317  END DO
318  END IF
319 *
320  IF ( notrans ) THEN
321  CALL zgbtrs( 'Conjugate transpose', n, kl, ku, 1, afb,
322  $ ldafb, ipiv, work, n, info )
323  ELSE
324  CALL zgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
325  $ ipiv, work, n, info )
326  END IF
327 *
328 * Multiply by R.
329 *
330  DO i = 1, n
331  work( i ) = work( i ) * rwork( i )
332  END DO
333  END IF
334  GO TO 10
335  END IF
336 *
337 * Compute the estimate of the reciprocal condition number.
338 *
339  IF( ainvnm .NE. 0.0d+0 )
340  $ zla_gbrcond_c = 1.0d+0 / ainvnm
341 *
342  RETURN
343 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
ZGBTRS
Definition: zgbtrs.f:140
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:135
double precision function zla_gbrcond_c(TRANS, N, KL, KU, AB, LDAB, AFB, LDAFB, IPIV, C, CAPPLY, INFO, WORK, RWORK)
ZLA_GBRCOND_C computes the infinity norm condition number of op(A)*inv(diag(c)) for general banded ma...
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function:

double precision function zla_gbrcond_x ( character  TRANS,
integer  N,
integer  KL,
integer  KU,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
complex*16, dimension( ldafb, * )  AFB,
integer  LDAFB,
integer, dimension( * )  IPIV,
complex*16, dimension( * )  X,
integer  INFO,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK 
)

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

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

Purpose:
    ZLA_GBRCOND_X Computes the infinity norm condition number of
    op(A) * diag(X) where X is a COMPLEX*16 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*16 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*16 array, dimension (LDAFB,N)
     Details of the LU factorization of the band matrix A, as
     computed by ZGBTRF.  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 ZGBTRF; row i of the matrix was interchanged
     with row IPIV(i).
[in]X
          X is COMPLEX*16 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.
[in]WORK
          WORK is COMPLEX*16 array, dimension (2*N).
     Workspace.
[in]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N).
     Workspace.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 157 of file zla_gbrcond_x.f.

157 *
158 * -- LAPACK computational routine (version 3.4.2) --
159 * -- LAPACK is a software package provided by Univ. of Tennessee, --
160 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
161 * September 2012
162 *
163 * .. Scalar Arguments ..
164  CHARACTER trans
165  INTEGER n, kl, ku, kd, ke, ldab, ldafb, info
166 * ..
167 * .. Array Arguments ..
168  INTEGER ipiv( * )
169  COMPLEX*16 ab( ldab, * ), afb( ldafb, * ), work( * ),
170  $ x( * )
171  DOUBLE PRECISION rwork( * )
172 *
173 *
174 * =====================================================================
175 *
176 * .. Local Scalars ..
177  LOGICAL notrans
178  INTEGER kase, i, j
179  DOUBLE PRECISION ainvnm, anorm, tmp
180  COMPLEX*16 zdum
181 * ..
182 * .. Local Arrays ..
183  INTEGER isave( 3 )
184 * ..
185 * .. External Functions ..
186  LOGICAL lsame
187  EXTERNAL lsame
188 * ..
189 * .. External Subroutines ..
190  EXTERNAL zlacn2, zgbtrs, xerbla
191 * ..
192 * .. Intrinsic Functions ..
193  INTRINSIC abs, max
194 * ..
195 * .. Statement Functions ..
196  DOUBLE PRECISION cabs1
197 * ..
198 * .. Statement Function Definitions ..
199  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
200 * ..
201 * .. Executable Statements ..
202 *
203  zla_gbrcond_x = 0.0d+0
204 *
205  info = 0
206  notrans = lsame( trans, 'N' )
207  IF ( .NOT. notrans .AND. .NOT. lsame(trans, 'T') .AND. .NOT.
208  $ lsame( trans, 'C' ) ) THEN
209  info = -1
210  ELSE IF( n.LT.0 ) THEN
211  info = -2
212  ELSE IF( kl.LT.0 .OR. kl.GT.n-1 ) THEN
213  info = -3
214  ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
215  info = -4
216  ELSE IF( ldab.LT.kl+ku+1 ) THEN
217  info = -6
218  ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
219  info = -8
220  END IF
221  IF( info.NE.0 ) THEN
222  CALL xerbla( 'ZLA_GBRCOND_X', -info )
223  RETURN
224  END IF
225 *
226 * Compute norm of op(A)*op2(C).
227 *
228  kd = ku + 1
229  ke = kl + 1
230  anorm = 0.0d+0
231  IF ( notrans ) THEN
232  DO i = 1, n
233  tmp = 0.0d+0
234  DO j = max( i-kl, 1 ), min( i+ku, n )
235  tmp = tmp + cabs1( ab( kd+i-j, j) * x( j ) )
236  END DO
237  rwork( i ) = tmp
238  anorm = max( anorm, tmp )
239  END DO
240  ELSE
241  DO i = 1, n
242  tmp = 0.0d+0
243  DO j = max( i-kl, 1 ), min( i+ku, n )
244  tmp = tmp + cabs1( ab( ke-i+j, i ) * x( j ) )
245  END DO
246  rwork( i ) = tmp
247  anorm = max( anorm, tmp )
248  END DO
249  END IF
250 *
251 * Quick return if possible.
252 *
253  IF( n.EQ.0 ) THEN
254  zla_gbrcond_x = 1.0d+0
255  RETURN
256  ELSE IF( anorm .EQ. 0.0d+0 ) THEN
257  RETURN
258  END IF
259 *
260 * Estimate the norm of inv(op(A)).
261 *
262  ainvnm = 0.0d+0
263 *
264  kase = 0
265  10 CONTINUE
266  CALL zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
267  IF( kase.NE.0 ) THEN
268  IF( kase.EQ.2 ) THEN
269 *
270 * Multiply by R.
271 *
272  DO i = 1, n
273  work( i ) = work( i ) * rwork( i )
274  END DO
275 *
276  IF ( notrans ) THEN
277  CALL zgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
278  $ ipiv, work, n, info )
279  ELSE
280  CALL zgbtrs( 'Conjugate transpose', n, kl, ku, 1, afb,
281  $ ldafb, ipiv, work, n, info )
282  ENDIF
283 *
284 * Multiply by inv(X).
285 *
286  DO i = 1, n
287  work( i ) = work( i ) / x( i )
288  END DO
289  ELSE
290 *
291 * Multiply by inv(X**H).
292 *
293  DO i = 1, n
294  work( i ) = work( i ) / x( i )
295  END DO
296 *
297  IF ( notrans ) THEN
298  CALL zgbtrs( 'Conjugate transpose', n, kl, ku, 1, afb,
299  $ ldafb, ipiv, work, n, info )
300  ELSE
301  CALL zgbtrs( 'No transpose', n, kl, ku, 1, afb, ldafb,
302  $ ipiv, work, n, info )
303  END IF
304 *
305 * Multiply by R.
306 *
307  DO i = 1, n
308  work( i ) = work( i ) * rwork( i )
309  END DO
310  END IF
311  GO TO 10
312  END IF
313 *
314 * Compute the estimate of the reciprocal condition number.
315 *
316  IF( ainvnm .NE. 0.0d+0 )
317  $ zla_gbrcond_x = 1.0d+0 / ainvnm
318 *
319  RETURN
320 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
ZGBTRS
Definition: zgbtrs.f:140
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:135
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
double precision function zla_gbrcond_x(TRANS, N, KL, KU, AB, LDAB, AFB, LDAFB, IPIV, X, INFO, WORK, RWORK)
ZLA_GBRCOND_X computes the infinity norm condition number of op(A)*diag(x) for general banded matrice...

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine zla_gbrfsx_extended ( integer  PREC_TYPE,
integer  TRANS_TYPE,
integer  N,
integer  KL,
integer  KU,
integer  NRHS,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
complex*16, dimension( ldafb, * )  AFB,
integer  LDAFB,
integer, dimension( * )  IPIV,
logical  COLEQU,
double precision, dimension( * )  C,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldy, * )  Y,
integer  LDY,
double precision, dimension( * )  BERR_OUT,
integer  N_NORMS,
double precision, dimension( nrhs, * )  ERR_BNDS_NORM,
double precision, dimension( nrhs, * )  ERR_BNDS_COMP,
complex*16, dimension( * )  RES,
double precision, dimension(*)  AYB,
complex*16, dimension( * )  DY,
complex*16, dimension( * )  Y_TAIL,
double precision  RCOND,
integer  ITHRESH,
double precision  RTHRESH,
double precision  DZ_UB,
logical  IGNORE_CWISE,
integer  INFO 
)

ZLA_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 ZLA_GBRFSX_EXTENDED + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZLA_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 ZGBRFSX 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 resonsible 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', '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*16 array, dimension (LDAB,N)
     On entry, the N-by-N matrix A.
[in]LDAB
          LDAB is INTEGER
     The leading dimension of the array A.  LDAB >= max(1,N).
[in]AFB
          AFB is COMPLEX*16 array, dimension (LDAF,N)
     The factors L and U from the factorization
     A = P*L*U as computed by ZGBTRF.
[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 ZGBTRF; row i of the matrix was interchanged
     with row IPIV(i).
[in]COLEQU
          COLEQU is LOGICAL
     If .TRUE. then column equilibration was done to A before calling
     this routine. This is needed to compute the solution and error
     bounds correctly.
[in]C
          C is DOUBLE PRECISION array, dimension (N)
     The column scale factors for A. If COLEQU = .FALSE., C
     is not accessed. If C is input, each element of C should be a power
     of the radix to ensure a reliable solution and error estimates.
     Scaling by powers of the radix does not cause rounding errors unless
     the result underflows or overflows. Rounding errors during scaling
     lead to refining with a matrix that is not equivalent to the
     input matrix, producing error estimates that may not be
     reliable.
[in]B
          B is COMPLEX*16 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*16 array, dimension (LDY,NRHS)
     On entry, the solution matrix X, as computed by ZGBTRS.
     On exit, the improved solution matrix Y.
[in]LDY
          LDY is INTEGER
     The leading dimension of the array Y.  LDY >= max(1,N).
[out]BERR_OUT
          BERR_OUT is DOUBLE PRECISION array, dimension (NRHS)
     On exit, BERR_OUT(j) contains the componentwise relative backward
     error for right-hand-side j from the formula
         max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
     where abs(Z) is the componentwise absolute value of the matrix
     or vector Z. This is computed by ZLA_LIN_BERR.
[in]N_NORMS
          N_NORMS is INTEGER
     Determines which error bounds to return (see ERR_BNDS_NORM
     and ERR_BNDS_COMP).
     If N_NORMS >= 1 return normwise error bounds.
     If N_NORMS >= 2 return componentwise error bounds.
[in,out]ERR_BNDS_NORM
          ERR_BNDS_NORM is DOUBLE PRECISION array, dimension
                    (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     normwise relative error, which is defined as follows:

     Normwise relative error in the ith solution vector:
             max_j (abs(XTRUE(j,i) - X(j,i)))
            ------------------------------
                  max_j abs(X(j,i))

     The array is indexed by the type of error information as described
     below. There currently are up to three pieces of information
     returned.

     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
     right-hand side.

     The second index in ERR_BNDS_NORM(:,err) contains the following
     three fields:
     err = 1 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" 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 "guaranteed". These reciprocal condition
              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
              appropriately scaled matrix Z.
              Let Z = S*A, where S scales each row by a power of the
              radix so all absolute row sums of Z are approximately 1.

     This subroutine is only responsible for setting the second field
     above.
     See Lapack Working Note 165 for further details and extra
     cautions.
[in,out]ERR_BNDS_COMP
          ERR_BNDS_COMP is DOUBLE PRECISION array, dimension
                    (NRHS, N_ERR_BNDS)
     For each right-hand side, this array contains information about
     various error bounds and condition numbers corresponding to the
     componentwise relative error, which is defined as follows:

     Componentwise relative error in the ith solution vector:
                    abs(XTRUE(j,i) - X(j,i))
             max_j ----------------------
                         abs(X(j,i))

     The array is indexed by the right-hand side i (on which the
     componentwise relative error depends), and the type of error
     information as described below. There currently are up to three
     pieces of information returned for each right-hand side. If
     componentwise accuracy is not requested (PARAMS(3) = 0.0), then
     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 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 "Trust/don't trust" boolean. Trust the answer if the
              reciprocal condition number is less than the threshold
              sqrt(n) * slamch('Epsilon').

     err = 2 "Guaranteed" 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 "guaranteed". 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*16 array, dimension (N)
     Workspace to hold the intermediate residual.
[in]AYB
          AYB is DOUBLE PRECISION array, dimension (N)
     Workspace.
[in]DY
          DY is COMPLEX*16 array, dimension (N)
     Workspace to hold the intermediate solution.
[in]Y_TAIL
          Y_TAIL is COMPLEX*16 array, dimension (N)
     Workspace to hold the trailing bits of the intermediate solution.
[in]RCOND
          RCOND is DOUBLE PRECISION
     Reciprocal scaled condition number.  This is an estimate of the
     reciprocal Skeel condition number of the matrix A after
     equilibration (if done).  If this is less than the machine
     precision (in particular, if it is zero), the matrix is singular
     to working precision.  Note that the error may still be small even
     if this number is very small and the matrix appears ill-
     conditioned.
[in]ITHRESH
          ITHRESH is INTEGER
     The maximum number of residual computations allowed for
     refinement. The default is 10. For 'aggressive' set to 100 to
     permit convergence using approximate factorizations or
     factorizations other than LU. If the factorization uses a
     technique other than Gaussian elimination, the guarantees in
     ERR_BNDS_NORM and ERR_BNDS_COMP may no longer be trustworthy.
[in]RTHRESH
          RTHRESH is DOUBLE PRECISION
     Determines when to stop refinement if the error estimate stops
     decreasing. Refinement will stop when the next solution no longer
     satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is
     the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The
     default value is 0.5. For 'aggressive' set to 0.9 to permit
     convergence on extremely ill-conditioned matrices. See LAWN 165
     for more details.
[in]DZ_UB
          DZ_UB is DOUBLE PRECISION
     Determines when to start considering componentwise convergence.
     Componentwise convergence is only considered after each component
     of the solution Y is stable, which we definte 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 ZGBTRS had an illegal
             value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 414 of file zla_gbrfsx_extended.f.

414 *
415 * -- LAPACK computational routine (version 3.4.2) --
416 * -- LAPACK is a software package provided by Univ. of Tennessee, --
417 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
418 * September 2012
419 *
420 * .. Scalar Arguments ..
421  INTEGER info, ldab, ldafb, ldb, ldy, n, kl, ku, nrhs,
422  $ prec_type, trans_type, n_norms, ithresh
423  LOGICAL colequ, ignore_cwise
424  DOUBLE PRECISION rthresh, dz_ub
425 * ..
426 * .. Array Arguments ..
427  INTEGER ipiv( * )
428  COMPLEX*16 ab( ldab, * ), afb( ldafb, * ), b( ldb, * ),
429  $ y( ldy, * ), res( * ), dy( * ), y_tail( * )
430  DOUBLE PRECISION c( * ), ayb(*), rcond, berr_out( * ),
431  $ err_bnds_norm( nrhs, * ),
432  $ err_bnds_comp( nrhs, * )
433 * ..
434 *
435 * =====================================================================
436 *
437 * .. Local Scalars ..
438  CHARACTER trans
439  INTEGER cnt, i, j, m, x_state, z_state, y_prec_state
440  DOUBLE PRECISION yk, dyk, ymin, normy, normx, normdx, dxrat,
441  $ dzrat, prevnormdx, prev_dz_z, dxratmax,
442  $ dzratmax, dx_x, dz_z, final_dx_x, final_dz_z,
443  $ eps, hugeval, incr_thresh
444  LOGICAL incr_prec
445  COMPLEX*16 zdum
446 * ..
447 * .. Parameters ..
448  INTEGER unstable_state, working_state, conv_state,
449  $ noprog_state, base_residual, extra_residual,
450  $ extra_y
451  parameter( unstable_state = 0, working_state = 1,
452  $ conv_state = 2, noprog_state = 3 )
453  parameter( base_residual = 0, extra_residual = 1,
454  $ extra_y = 2 )
455  INTEGER final_nrm_err_i, final_cmp_err_i, berr_i
456  INTEGER rcond_i, nrm_rcond_i, nrm_err_i, cmp_rcond_i
457  INTEGER cmp_err_i, piv_growth_i
458  parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
459  $ berr_i = 3 )
460  parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
461  parameter( cmp_rcond_i = 7, cmp_err_i = 8,
462  $ piv_growth_i = 9 )
463  INTEGER la_linrx_itref_i, la_linrx_ithresh_i,
464  $ la_linrx_cwise_i
465  parameter( la_linrx_itref_i = 1,
466  $ la_linrx_ithresh_i = 2 )
467  parameter( la_linrx_cwise_i = 3 )
468  INTEGER la_linrx_trust_i, la_linrx_err_i,
469  $ la_linrx_rcond_i
470  parameter( la_linrx_trust_i = 1, la_linrx_err_i = 2 )
471  parameter( la_linrx_rcond_i = 3 )
472 * ..
473 * .. External Subroutines ..
474  EXTERNAL zaxpy, zcopy, zgbtrs, zgbmv, blas_zgbmv_x,
475  $ blas_zgbmv2_x, zla_gbamv, zla_wwaddw, dlamch,
477  DOUBLE PRECISION dlamch
478  CHARACTER chla_transtype
479 * ..
480 * .. Intrinsic Functions..
481  INTRINSIC abs, max, min
482 * ..
483 * .. Statement Functions ..
484  DOUBLE PRECISION cabs1
485 * ..
486 * .. Statement Function Definitions ..
487  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
488 * ..
489 * .. Executable Statements ..
490 *
491  IF (info.NE.0) RETURN
492  trans = chla_transtype(trans_type)
493  eps = dlamch( 'Epsilon' )
494  hugeval = dlamch( 'Overflow' )
495 * Force HUGEVAL to Inf
496  hugeval = hugeval * hugeval
497 * Using HUGEVAL may lead to spurious underflows.
498  incr_thresh = dble( n ) * eps
499  m = kl+ku+1
500 
501  DO j = 1, nrhs
502  y_prec_state = extra_residual
503  IF ( y_prec_state .EQ. extra_y ) THEN
504  DO i = 1, n
505  y_tail( i ) = 0.0d+0
506  END DO
507  END IF
508 
509  dxrat = 0.0d+0
510  dxratmax = 0.0d+0
511  dzrat = 0.0d+0
512  dzratmax = 0.0d+0
513  final_dx_x = hugeval
514  final_dz_z = hugeval
515  prevnormdx = hugeval
516  prev_dz_z = hugeval
517  dz_z = hugeval
518  dx_x = hugeval
519 
520  x_state = working_state
521  z_state = unstable_state
522  incr_prec = .false.
523 
524  DO cnt = 1, ithresh
525 *
526 * Compute residual RES = B_s - op(A_s) * Y,
527 * op(A) = A, A**T, or A**H depending on TRANS (and type).
528 *
529  CALL zcopy( n, b( 1, j ), 1, res, 1 )
530  IF ( y_prec_state .EQ. base_residual ) THEN
531  CALL zgbmv( trans, m, n, kl, ku, (-1.0d+0,0.0d+0), ab,
532  $ ldab, y( 1, j ), 1, (1.0d+0,0.0d+0), res, 1 )
533  ELSE IF ( y_prec_state .EQ. extra_residual ) THEN
534  CALL blas_zgbmv_x( trans_type, n, n, kl, ku,
535  $ (-1.0d+0,0.0d+0), ab, ldab, y( 1, j ), 1,
536  $ (1.0d+0,0.0d+0), res, 1, prec_type )
537  ELSE
538  CALL blas_zgbmv2_x( trans_type, n, n, kl, ku,
539  $ (-1.0d+0,0.0d+0), ab, ldab, y( 1, j ), y_tail, 1,
540  $ (1.0d+0,0.0d+0), res, 1, prec_type )
541  END IF
542 
543 ! XXX: RES is no longer needed.
544  CALL zcopy( n, res, 1, dy, 1 )
545  CALL zgbtrs( trans, n, kl, ku, 1, afb, ldafb, ipiv, dy, n,
546  $ info )
547 *
548 * Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
549 *
550  normx = 0.0d+0
551  normy = 0.0d+0
552  normdx = 0.0d+0
553  dz_z = 0.0d+0
554  ymin = hugeval
555 
556  DO i = 1, n
557  yk = cabs1( y( i, j ) )
558  dyk = cabs1( dy( i ) )
559 
560  IF (yk .NE. 0.0d+0) THEN
561  dz_z = max( dz_z, dyk / yk )
562  ELSE IF ( dyk .NE. 0.0d+0 ) THEN
563  dz_z = hugeval
564  END IF
565 
566  ymin = min( ymin, yk )
567 
568  normy = max( normy, yk )
569 
570  IF ( colequ ) THEN
571  normx = max( normx, yk * c( i ) )
572  normdx = max(normdx, dyk * c(i))
573  ELSE
574  normx = normy
575  normdx = max( normdx, dyk )
576  END IF
577  END DO
578 
579  IF ( normx .NE. 0.0d+0 ) THEN
580  dx_x = normdx / normx
581  ELSE IF ( normdx .EQ. 0.0d+0 ) THEN
582  dx_x = 0.0d+0
583  ELSE
584  dx_x = hugeval
585  END IF
586 
587  dxrat = normdx / prevnormdx
588  dzrat = dz_z / prev_dz_z
589 *
590 * Check termination criteria.
591 *
592  IF (.NOT.ignore_cwise
593  $ .AND. ymin*rcond .LT. incr_thresh*normy
594  $ .AND. y_prec_state .LT. extra_y )
595  $ incr_prec = .true.
596 
597  IF ( x_state .EQ. noprog_state .AND. dxrat .LE. rthresh )
598  $ x_state = working_state
599  IF ( x_state .EQ. working_state ) THEN
600  IF ( dx_x .LE. eps ) THEN
601  x_state = conv_state
602  ELSE IF ( dxrat .GT. rthresh ) THEN
603  IF ( y_prec_state .NE. extra_y ) THEN
604  incr_prec = .true.
605  ELSE
606  x_state = noprog_state
607  END IF
608  ELSE
609  IF ( dxrat .GT. dxratmax ) dxratmax = dxrat
610  END IF
611  IF ( x_state .GT. working_state ) final_dx_x = dx_x
612  END IF
613 
614  IF ( z_state .EQ. unstable_state .AND. dz_z .LE. dz_ub )
615  $ z_state = working_state
616  IF ( z_state .EQ. noprog_state .AND. dzrat .LE. rthresh )
617  $ z_state = working_state
618  IF ( z_state .EQ. working_state ) THEN
619  IF ( dz_z .LE. eps ) THEN
620  z_state = conv_state
621  ELSE IF ( dz_z .GT. dz_ub ) THEN
622  z_state = unstable_state
623  dzratmax = 0.0d+0
624  final_dz_z = hugeval
625  ELSE IF ( dzrat .GT. rthresh ) THEN
626  IF ( y_prec_state .NE. extra_y ) THEN
627  incr_prec = .true.
628  ELSE
629  z_state = noprog_state
630  END IF
631  ELSE
632  IF ( dzrat .GT. dzratmax ) dzratmax = dzrat
633  END IF
634  IF ( z_state .GT. working_state ) final_dz_z = dz_z
635  END IF
636 *
637 * Exit if both normwise and componentwise stopped working,
638 * but if componentwise is unstable, let it go at least two
639 * iterations.
640 *
641  IF ( x_state.NE.working_state ) THEN
642  IF ( ignore_cwise ) GOTO 666
643  IF ( z_state.EQ.noprog_state .OR. z_state.EQ.conv_state )
644  $ GOTO 666
645  IF ( z_state.EQ.unstable_state .AND. cnt.GT.1 ) GOTO 666
646  END IF
647 
648  IF ( incr_prec ) THEN
649  incr_prec = .false.
650  y_prec_state = y_prec_state + 1
651  DO i = 1, n
652  y_tail( i ) = 0.0d+0
653  END DO
654  END IF
655 
656  prevnormdx = normdx
657  prev_dz_z = dz_z
658 *
659 * Update soluton.
660 *
661  IF ( y_prec_state .LT. extra_y ) THEN
662  CALL zaxpy( n, (1.0d+0,0.0d+0), dy, 1, y(1,j), 1 )
663  ELSE
664  CALL zla_wwaddw( n, y(1,j), y_tail, dy )
665  END IF
666 
667  END DO
668 * Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
669  666 CONTINUE
670 *
671 * Set final_* when cnt hits ithresh.
672 *
673  IF ( x_state .EQ. working_state ) final_dx_x = dx_x
674  IF ( z_state .EQ. working_state ) final_dz_z = dz_z
675 *
676 * Compute error bounds.
677 *
678  IF ( n_norms .GE. 1 ) THEN
679  err_bnds_norm( j, la_linrx_err_i ) =
680  $ final_dx_x / (1 - dxratmax)
681  END IF
682  IF ( n_norms .GE. 2 ) THEN
683  err_bnds_comp( j, la_linrx_err_i ) =
684  $ final_dz_z / (1 - dzratmax)
685  END IF
686 *
687 * Compute componentwise relative backward error from formula
688 * max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
689 * where abs(Z) is the componentwise absolute value of the matrix
690 * or vector Z.
691 *
692 * Compute residual RES = B_s - op(A_s) * Y,
693 * op(A) = A, A**T, or A**H depending on TRANS (and type).
694 *
695  CALL zcopy( n, b( 1, j ), 1, res, 1 )
696  CALL zgbmv( trans, n, n, kl, ku, (-1.0d+0,0.0d+0), ab, ldab,
697  $ y(1,j), 1, (1.0d+0,0.0d+0), res, 1 )
698 
699  DO i = 1, n
700  ayb( i ) = cabs1( b( i, j ) )
701  END DO
702 *
703 * Compute abs(op(A_s))*abs(Y) + abs(B_s).
704 *
705  CALL zla_gbamv( trans_type, n, n, kl, ku, 1.0d+0,
706  $ ab, ldab, y(1, j), 1, 1.0d+0, ayb, 1 )
707 
708  CALL zla_lin_berr( n, n, 1, res, ayb, berr_out( j ) )
709 *
710 * End of loop for each RHS.
711 *
712  END DO
713 *
714  RETURN
subroutine zgbmv(TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGBMV
Definition: zgbmv.f:189
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
ZGBTRS
Definition: zgbtrs.f:140
subroutine zla_lin_berr(N, NZ, NRHS, RES, AYB, BERR)
ZLA_LIN_BERR computes a component-wise relative backward error.
Definition: zla_lin_berr.f:103
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zla_wwaddw(N, X, Y, W)
ZLA_WWADDW adds a vector into a doubled-single vector.
Definition: zla_wwaddw.f:83
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53
subroutine zla_gbamv(TRANS, M, N, KL, KU, ALPHA, AB, LDAB, X, INCX, BETA, Y, INCY)
ZLA_GBAMV performs a matrix-vector operation to calculate error bounds.
Definition: zla_gbamv.f:188
character *1 function chla_transtype(TRANS)
CHLA_TRANSTYPE

Here is the call graph for this function:

Here is the caller graph for this function:

double precision function zla_gbrpvgrw ( integer  N,
integer  KL,
integer  KU,
integer  NCOLS,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
complex*16, dimension( ldafb, * )  AFB,
integer  LDAFB 
)

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

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

Purpose:
 ZLA_GBRPVGRW computes the reciprocal pivot growth factor
 norm(A)/norm(U). The "max absolute element" 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*16 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*16 array, dimension (LDAFB,N)
     Details of the LU factorization of the band matrix A, as
     computed by ZGBTRF.  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.
Date
September 2012

Definition at line 119 of file zla_gbrpvgrw.f.

119 *
120 * -- LAPACK computational routine (version 3.4.2) --
121 * -- LAPACK is a software package provided by Univ. of Tennessee, --
122 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
123 * September 2012
124 *
125 * .. Scalar Arguments ..
126  INTEGER n, kl, ku, ncols, ldab, ldafb
127 * ..
128 * .. Array Arguments ..
129  COMPLEX*16 ab( ldab, * ), afb( ldafb, * )
130 * ..
131 *
132 * =====================================================================
133 *
134 * .. Local Scalars ..
135  INTEGER i, j, kd
136  DOUBLE PRECISION amax, umax, rpvgrw
137  COMPLEX*16 zdum
138 * ..
139 * .. Intrinsic Functions ..
140  INTRINSIC abs, max, min, REAL, dimag
141 * ..
142 * .. Statement Functions ..
143  DOUBLE PRECISION cabs1
144 * ..
145 * .. Statement Function Definitions ..
146  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
147 * ..
148 * .. Executable Statements ..
149 *
150  rpvgrw = 1.0d+0
151 
152  kd = ku + 1
153  DO j = 1, ncols
154  amax = 0.0d+0
155  umax = 0.0d+0
156  DO i = max( j-ku, 1 ), min( j+kl, n )
157  amax = max( cabs1( ab( kd+i-j, j ) ), amax )
158  END DO
159  DO i = max( j-ku, 1 ), j
160  umax = max( cabs1( afb( kd+i-j, j ) ), umax )
161  END DO
162  IF ( umax /= 0.0d+0 ) THEN
163  rpvgrw = min( amax / umax, rpvgrw )
164  END IF
165  END DO
166  zla_gbrpvgrw = rpvgrw
double precision function zla_gbrpvgrw(N, KL, KU, NCOLS, AB, LDAB, AFB, LDAFB)
ZLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix...
Definition: zla_gbrpvgrw.f:119

Here is the caller graph for this function:

subroutine zungbr ( character  VECT,
integer  M,
integer  N,
integer  K,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  TAU,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZUNGBR

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

Purpose:
 ZUNGBR generates one of the complex unitary matrices Q or P**H
 determined by ZGEBRD 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 ZUNGBR returns the first n
 columns of Q, where m >= n >= k;
 if m < k, Q = H(1) H(2) . . . H(m-1) and ZUNGBR 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 ZUNGBR 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 ZUNGBR 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 ZGEBRD:
          = '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 ZGEBRD.
          If VECT = 'P', the number of rows in the original K-by-N
          matrix reduced by ZGEBRD.
          K >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the vectors which define the elementary reflectors,
          as returned by ZGEBRD.
          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*16 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 ZGEBRD in its array argument TAUQ or TAUP.
[out]WORK
          WORK is COMPLEX*16 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.
Date
April 2012

Definition at line 159 of file zungbr.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: