Annotation of rpl/lapack/lapack/zpbtf2.f, revision 1.19

1.12      bertrand    1: *> \brief \b ZPBTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite band matrix (unblocked algorithm).
1.9       bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.16      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.9       bertrand    7: *
                      8: *> \htmlonly
1.16      bertrand    9: *> Download ZPBTF2 + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpbtf2.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpbtf2.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpbtf2.f">
1.9       bertrand   15: *> [TXT]</a>
1.16      bertrand   16: *> \endhtmlonly
1.9       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE ZPBTF2( UPLO, N, KD, AB, LDAB, INFO )
1.16      bertrand   22: *
1.9       bertrand   23: *       .. Scalar Arguments ..
                     24: *       CHARACTER          UPLO
                     25: *       INTEGER            INFO, KD, LDAB, N
                     26: *       ..
                     27: *       .. Array Arguments ..
                     28: *       COMPLEX*16         AB( LDAB, * )
                     29: *       ..
1.16      bertrand   30: *
1.9       bertrand   31: *
                     32: *> \par Purpose:
                     33: *  =============
                     34: *>
                     35: *> \verbatim
                     36: *>
                     37: *> ZPBTF2 computes the Cholesky factorization of a complex Hermitian
                     38: *> positive definite band matrix A.
                     39: *>
                     40: *> The factorization has the form
                     41: *>    A = U**H * U ,  if UPLO = 'U', or
                     42: *>    A = L  * L**H,  if UPLO = 'L',
                     43: *> where U is an upper triangular matrix, U**H is the conjugate transpose
                     44: *> of U, and L is lower triangular.
                     45: *>
                     46: *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
                     47: *> \endverbatim
                     48: *
                     49: *  Arguments:
                     50: *  ==========
                     51: *
                     52: *> \param[in] UPLO
                     53: *> \verbatim
                     54: *>          UPLO is CHARACTER*1
                     55: *>          Specifies whether the upper or lower triangular part of the
                     56: *>          Hermitian matrix A is stored:
                     57: *>          = 'U':  Upper triangular
                     58: *>          = 'L':  Lower triangular
                     59: *> \endverbatim
                     60: *>
                     61: *> \param[in] N
                     62: *> \verbatim
                     63: *>          N is INTEGER
                     64: *>          The order of the matrix A.  N >= 0.
                     65: *> \endverbatim
                     66: *>
                     67: *> \param[in] KD
                     68: *> \verbatim
                     69: *>          KD is INTEGER
                     70: *>          The number of super-diagonals of the matrix A if UPLO = 'U',
                     71: *>          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
                     72: *> \endverbatim
                     73: *>
                     74: *> \param[in,out] AB
                     75: *> \verbatim
                     76: *>          AB is COMPLEX*16 array, dimension (LDAB,N)
                     77: *>          On entry, the upper or lower triangle of the Hermitian band
                     78: *>          matrix A, stored in the first KD+1 rows of the array.  The
                     79: *>          j-th column of A is stored in the j-th column of the array AB
                     80: *>          as follows:
                     81: *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
                     82: *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
                     83: *>
                     84: *>          On exit, if INFO = 0, the triangular factor U or L from the
                     85: *>          Cholesky factorization A = U**H *U or A = L*L**H of the band
                     86: *>          matrix A, in the same storage format as A.
                     87: *> \endverbatim
                     88: *>
                     89: *> \param[in] LDAB
                     90: *> \verbatim
                     91: *>          LDAB is INTEGER
                     92: *>          The leading dimension of the array AB.  LDAB >= KD+1.
                     93: *> \endverbatim
                     94: *>
                     95: *> \param[out] INFO
                     96: *> \verbatim
                     97: *>          INFO is INTEGER
                     98: *>          = 0: successful exit
                     99: *>          < 0: if INFO = -k, the k-th argument had an illegal value
                    100: *>          > 0: if INFO = k, the leading minor of order k is not
                    101: *>               positive definite, and the factorization could not be
                    102: *>               completed.
                    103: *> \endverbatim
                    104: *
                    105: *  Authors:
                    106: *  ========
                    107: *
1.16      bertrand  108: *> \author Univ. of Tennessee
                    109: *> \author Univ. of California Berkeley
                    110: *> \author Univ. of Colorado Denver
                    111: *> \author NAG Ltd.
1.9       bertrand  112: *
                    113: *> \ingroup complex16OTHERcomputational
                    114: *
                    115: *> \par Further Details:
                    116: *  =====================
                    117: *>
                    118: *> \verbatim
                    119: *>
                    120: *>  The band storage scheme is illustrated by the following example, when
                    121: *>  N = 6, KD = 2, and UPLO = 'U':
                    122: *>
                    123: *>  On entry:                       On exit:
                    124: *>
                    125: *>      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
                    126: *>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
                    127: *>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
                    128: *>
                    129: *>  Similarly, if UPLO = 'L' the format of A is as follows:
                    130: *>
                    131: *>  On entry:                       On exit:
                    132: *>
                    133: *>     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
                    134: *>     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
                    135: *>     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
                    136: *>
                    137: *>  Array elements marked * are not used by the routine.
                    138: *> \endverbatim
                    139: *>
                    140: *  =====================================================================
1.1       bertrand  141:       SUBROUTINE ZPBTF2( UPLO, N, KD, AB, LDAB, INFO )
                    142: *
1.19    ! bertrand  143: *  -- LAPACK computational routine --
1.1       bertrand  144: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    145: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    146: *
                    147: *     .. Scalar Arguments ..
                    148:       CHARACTER          UPLO
                    149:       INTEGER            INFO, KD, LDAB, N
                    150: *     ..
                    151: *     .. Array Arguments ..
                    152:       COMPLEX*16         AB( LDAB, * )
                    153: *     ..
                    154: *
                    155: *  =====================================================================
                    156: *
                    157: *     .. Parameters ..
                    158:       DOUBLE PRECISION   ONE, ZERO
                    159:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
                    160: *     ..
                    161: *     .. Local Scalars ..
                    162:       LOGICAL            UPPER
                    163:       INTEGER            J, KLD, KN
                    164:       DOUBLE PRECISION   AJJ
                    165: *     ..
                    166: *     .. External Functions ..
                    167:       LOGICAL            LSAME
                    168:       EXTERNAL           LSAME
                    169: *     ..
                    170: *     .. External Subroutines ..
                    171:       EXTERNAL           XERBLA, ZDSCAL, ZHER, ZLACGV
                    172: *     ..
                    173: *     .. Intrinsic Functions ..
                    174:       INTRINSIC          DBLE, MAX, MIN, SQRT
                    175: *     ..
                    176: *     .. Executable Statements ..
                    177: *
                    178: *     Test the input parameters.
                    179: *
                    180:       INFO = 0
                    181:       UPPER = LSAME( UPLO, 'U' )
                    182:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
                    183:          INFO = -1
                    184:       ELSE IF( N.LT.0 ) THEN
                    185:          INFO = -2
                    186:       ELSE IF( KD.LT.0 ) THEN
                    187:          INFO = -3
                    188:       ELSE IF( LDAB.LT.KD+1 ) THEN
                    189:          INFO = -5
                    190:       END IF
                    191:       IF( INFO.NE.0 ) THEN
                    192:          CALL XERBLA( 'ZPBTF2', -INFO )
                    193:          RETURN
                    194:       END IF
                    195: *
                    196: *     Quick return if possible
                    197: *
                    198:       IF( N.EQ.0 )
                    199:      $   RETURN
                    200: *
                    201:       KLD = MAX( 1, LDAB-1 )
                    202: *
                    203:       IF( UPPER ) THEN
                    204: *
1.8       bertrand  205: *        Compute the Cholesky factorization A = U**H * U.
1.1       bertrand  206: *
                    207:          DO 10 J = 1, N
                    208: *
                    209: *           Compute U(J,J) and test for non-positive-definiteness.
                    210: *
                    211:             AJJ = DBLE( AB( KD+1, J ) )
                    212:             IF( AJJ.LE.ZERO ) THEN
                    213:                AB( KD+1, J ) = AJJ
                    214:                GO TO 30
                    215:             END IF
                    216:             AJJ = SQRT( AJJ )
                    217:             AB( KD+1, J ) = AJJ
                    218: *
                    219: *           Compute elements J+1:J+KN of row J and update the
                    220: *           trailing submatrix within the band.
                    221: *
                    222:             KN = MIN( KD, N-J )
                    223:             IF( KN.GT.0 ) THEN
                    224:                CALL ZDSCAL( KN, ONE / AJJ, AB( KD, J+1 ), KLD )
                    225:                CALL ZLACGV( KN, AB( KD, J+1 ), KLD )
                    226:                CALL ZHER( 'Upper', KN, -ONE, AB( KD, J+1 ), KLD,
                    227:      $                    AB( KD+1, J+1 ), KLD )
                    228:                CALL ZLACGV( KN, AB( KD, J+1 ), KLD )
                    229:             END IF
                    230:    10    CONTINUE
                    231:       ELSE
                    232: *
1.8       bertrand  233: *        Compute the Cholesky factorization A = L*L**H.
1.1       bertrand  234: *
                    235:          DO 20 J = 1, N
                    236: *
                    237: *           Compute L(J,J) and test for non-positive-definiteness.
                    238: *
                    239:             AJJ = DBLE( AB( 1, J ) )
                    240:             IF( AJJ.LE.ZERO ) THEN
                    241:                AB( 1, J ) = AJJ
                    242:                GO TO 30
                    243:             END IF
                    244:             AJJ = SQRT( AJJ )
                    245:             AB( 1, J ) = AJJ
                    246: *
                    247: *           Compute elements J+1:J+KN of column J and update the
                    248: *           trailing submatrix within the band.
                    249: *
                    250:             KN = MIN( KD, N-J )
                    251:             IF( KN.GT.0 ) THEN
                    252:                CALL ZDSCAL( KN, ONE / AJJ, AB( 2, J ), 1 )
                    253:                CALL ZHER( 'Lower', KN, -ONE, AB( 2, J ), 1,
                    254:      $                    AB( 1, J+1 ), KLD )
                    255:             END IF
                    256:    20    CONTINUE
                    257:       END IF
                    258:       RETURN
                    259: *
                    260:    30 CONTINUE
                    261:       INFO = J
                    262:       RETURN
                    263: *
                    264: *     End of ZPBTF2
                    265: *
                    266:       END

CVSweb interface <joel.bertrand@systella.fr>