Annotation of rpl/lapack/lapack/dpbstf.f, revision 1.8

1.8     ! bertrand    1: *> \brief \b DPBSTF
        !             2: *
        !             3: *  =========== DOCUMENTATION ===========
        !             4: *
        !             5: * Online html documentation available at 
        !             6: *            http://www.netlib.org/lapack/explore-html/ 
        !             7: *
        !             8: *> \htmlonly
        !             9: *> Download DPBSTF + dependencies 
        !            10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpbstf.f"> 
        !            11: *> [TGZ]</a> 
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpbstf.f"> 
        !            13: *> [ZIP]</a> 
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpbstf.f"> 
        !            15: *> [TXT]</a>
        !            16: *> \endhtmlonly 
        !            17: *
        !            18: *  Definition:
        !            19: *  ===========
        !            20: *
        !            21: *       SUBROUTINE DPBSTF( UPLO, N, KD, AB, LDAB, INFO )
        !            22: * 
        !            23: *       .. Scalar Arguments ..
        !            24: *       CHARACTER          UPLO
        !            25: *       INTEGER            INFO, KD, LDAB, N
        !            26: *       ..
        !            27: *       .. Array Arguments ..
        !            28: *       DOUBLE PRECISION   AB( LDAB, * )
        !            29: *       ..
        !            30: *  
        !            31: *
        !            32: *> \par Purpose:
        !            33: *  =============
        !            34: *>
        !            35: *> \verbatim
        !            36: *>
        !            37: *> DPBSTF computes a split Cholesky factorization of a real
        !            38: *> symmetric positive definite band matrix A.
        !            39: *>
        !            40: *> This routine is designed to be used in conjunction with DSBGST.
        !            41: *>
        !            42: *> The factorization has the form  A = S**T*S  where S is a band matrix
        !            43: *> of the same bandwidth as A and the following structure:
        !            44: *>
        !            45: *>   S = ( U    )
        !            46: *>       ( M  L )
        !            47: *>
        !            48: *> where U is upper triangular of order m = (n+kd)/2, and L is lower
        !            49: *> triangular of order n-m.
        !            50: *> \endverbatim
        !            51: *
        !            52: *  Arguments:
        !            53: *  ==========
        !            54: *
        !            55: *> \param[in] UPLO
        !            56: *> \verbatim
        !            57: *>          UPLO is CHARACTER*1
        !            58: *>          = 'U':  Upper triangle of A is stored;
        !            59: *>          = 'L':  Lower triangle of A is stored.
        !            60: *> \endverbatim
        !            61: *>
        !            62: *> \param[in] N
        !            63: *> \verbatim
        !            64: *>          N is INTEGER
        !            65: *>          The order of the matrix A.  N >= 0.
        !            66: *> \endverbatim
        !            67: *>
        !            68: *> \param[in] KD
        !            69: *> \verbatim
        !            70: *>          KD is INTEGER
        !            71: *>          The number of superdiagonals of the matrix A if UPLO = 'U',
        !            72: *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
        !            73: *> \endverbatim
        !            74: *>
        !            75: *> \param[in,out] AB
        !            76: *> \verbatim
        !            77: *>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
        !            78: *>          On entry, the upper or lower triangle of the symmetric band
        !            79: *>          matrix A, stored in the first kd+1 rows of the array.  The
        !            80: *>          j-th column of A is stored in the j-th column of the array AB
        !            81: *>          as follows:
        !            82: *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
        !            83: *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
        !            84: *>
        !            85: *>          On exit, if INFO = 0, the factor S from the split Cholesky
        !            86: *>          factorization A = S**T*S. See Further Details.
        !            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 = -i, the i-th argument had an illegal value
        !           100: *>          > 0: if INFO = i, the factorization could not be completed,
        !           101: *>               because the updated element a(i,i) was negative; the
        !           102: *>               matrix A is not positive definite.
        !           103: *> \endverbatim
        !           104: *
        !           105: *  Authors:
        !           106: *  ========
        !           107: *
        !           108: *> \author Univ. of Tennessee 
        !           109: *> \author Univ. of California Berkeley 
        !           110: *> \author Univ. of Colorado Denver 
        !           111: *> \author NAG Ltd. 
        !           112: *
        !           113: *> \date November 2011
        !           114: *
        !           115: *> \ingroup doubleOTHERcomputational
        !           116: *
        !           117: *> \par Further Details:
        !           118: *  =====================
        !           119: *>
        !           120: *> \verbatim
        !           121: *>
        !           122: *>  The band storage scheme is illustrated by the following example, when
        !           123: *>  N = 7, KD = 2:
        !           124: *>
        !           125: *>  S = ( s11  s12  s13                     )
        !           126: *>      (      s22  s23  s24                )
        !           127: *>      (           s33  s34                )
        !           128: *>      (                s44                )
        !           129: *>      (           s53  s54  s55           )
        !           130: *>      (                s64  s65  s66      )
        !           131: *>      (                     s75  s76  s77 )
        !           132: *>
        !           133: *>  If UPLO = 'U', the array AB holds:
        !           134: *>
        !           135: *>  on entry:                          on exit:
        !           136: *>
        !           137: *>   *    *   a13  a24  a35  a46  a57   *    *   s13  s24  s53  s64  s75
        !           138: *>   *   a12  a23  a34  a45  a56  a67   *   s12  s23  s34  s54  s65  s76
        !           139: *>  a11  a22  a33  a44  a55  a66  a77  s11  s22  s33  s44  s55  s66  s77
        !           140: *>
        !           141: *>  If UPLO = 'L', the array AB holds:
        !           142: *>
        !           143: *>  on entry:                          on exit:
        !           144: *>
        !           145: *>  a11  a22  a33  a44  a55  a66  a77  s11  s22  s33  s44  s55  s66  s77
        !           146: *>  a21  a32  a43  a54  a65  a76   *   s12  s23  s34  s54  s65  s76   *
        !           147: *>  a31  a42  a53  a64  a64   *    *   s13  s24  s53  s64  s75   *    *
        !           148: *>
        !           149: *>  Array elements marked * are not used by the routine.
        !           150: *> \endverbatim
        !           151: *>
        !           152: *  =====================================================================
1.1       bertrand  153:       SUBROUTINE DPBSTF( UPLO, N, KD, AB, LDAB, INFO )
                    154: *
1.8     ! bertrand  155: *  -- LAPACK computational routine (version 3.4.0) --
1.1       bertrand  156: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    157: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8     ! bertrand  158: *     November 2011
1.1       bertrand  159: *
                    160: *     .. Scalar Arguments ..
                    161:       CHARACTER          UPLO
                    162:       INTEGER            INFO, KD, LDAB, N
                    163: *     ..
                    164: *     .. Array Arguments ..
                    165:       DOUBLE PRECISION   AB( LDAB, * )
                    166: *     ..
                    167: *
                    168: *  =====================================================================
                    169: *
                    170: *     .. Parameters ..
                    171:       DOUBLE PRECISION   ONE, ZERO
                    172:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
                    173: *     ..
                    174: *     .. Local Scalars ..
                    175:       LOGICAL            UPPER
                    176:       INTEGER            J, KLD, KM, M
                    177:       DOUBLE PRECISION   AJJ
                    178: *     ..
                    179: *     .. External Functions ..
                    180:       LOGICAL            LSAME
                    181:       EXTERNAL           LSAME
                    182: *     ..
                    183: *     .. External Subroutines ..
                    184:       EXTERNAL           DSCAL, DSYR, XERBLA
                    185: *     ..
                    186: *     .. Intrinsic Functions ..
                    187:       INTRINSIC          MAX, MIN, SQRT
                    188: *     ..
                    189: *     .. Executable Statements ..
                    190: *
                    191: *     Test the input parameters.
                    192: *
                    193:       INFO = 0
                    194:       UPPER = LSAME( UPLO, 'U' )
                    195:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
                    196:          INFO = -1
                    197:       ELSE IF( N.LT.0 ) THEN
                    198:          INFO = -2
                    199:       ELSE IF( KD.LT.0 ) THEN
                    200:          INFO = -3
                    201:       ELSE IF( LDAB.LT.KD+1 ) THEN
                    202:          INFO = -5
                    203:       END IF
                    204:       IF( INFO.NE.0 ) THEN
                    205:          CALL XERBLA( 'DPBSTF', -INFO )
                    206:          RETURN
                    207:       END IF
                    208: *
                    209: *     Quick return if possible
                    210: *
                    211:       IF( N.EQ.0 )
                    212:      $   RETURN
                    213: *
                    214:       KLD = MAX( 1, LDAB-1 )
                    215: *
                    216: *     Set the splitting point m.
                    217: *
                    218:       M = ( N+KD ) / 2
                    219: *
                    220:       IF( UPPER ) THEN
                    221: *
                    222: *        Factorize A(m+1:n,m+1:n) as L**T*L, and update A(1:m,1:m).
                    223: *
                    224:          DO 10 J = N, M + 1, -1
                    225: *
                    226: *           Compute s(j,j) and test for non-positive-definiteness.
                    227: *
                    228:             AJJ = AB( KD+1, J )
                    229:             IF( AJJ.LE.ZERO )
                    230:      $         GO TO 50
                    231:             AJJ = SQRT( AJJ )
                    232:             AB( KD+1, J ) = AJJ
                    233:             KM = MIN( J-1, KD )
                    234: *
                    235: *           Compute elements j-km:j-1 of the j-th column and update the
                    236: *           the leading submatrix within the band.
                    237: *
                    238:             CALL DSCAL( KM, ONE / AJJ, AB( KD+1-KM, J ), 1 )
                    239:             CALL DSYR( 'Upper', KM, -ONE, AB( KD+1-KM, J ), 1,
                    240:      $                 AB( KD+1, J-KM ), KLD )
                    241:    10    CONTINUE
                    242: *
                    243: *        Factorize the updated submatrix A(1:m,1:m) as U**T*U.
                    244: *
                    245:          DO 20 J = 1, M
                    246: *
                    247: *           Compute s(j,j) and test for non-positive-definiteness.
                    248: *
                    249:             AJJ = AB( KD+1, J )
                    250:             IF( AJJ.LE.ZERO )
                    251:      $         GO TO 50
                    252:             AJJ = SQRT( AJJ )
                    253:             AB( KD+1, J ) = AJJ
                    254:             KM = MIN( KD, M-J )
                    255: *
                    256: *           Compute elements j+1:j+km of the j-th row and update the
                    257: *           trailing submatrix within the band.
                    258: *
                    259:             IF( KM.GT.0 ) THEN
                    260:                CALL DSCAL( KM, ONE / AJJ, AB( KD, J+1 ), KLD )
                    261:                CALL DSYR( 'Upper', KM, -ONE, AB( KD, J+1 ), KLD,
                    262:      $                    AB( KD+1, J+1 ), KLD )
                    263:             END IF
                    264:    20    CONTINUE
                    265:       ELSE
                    266: *
                    267: *        Factorize A(m+1:n,m+1:n) as L**T*L, and update A(1:m,1:m).
                    268: *
                    269:          DO 30 J = N, M + 1, -1
                    270: *
                    271: *           Compute s(j,j) and test for non-positive-definiteness.
                    272: *
                    273:             AJJ = AB( 1, J )
                    274:             IF( AJJ.LE.ZERO )
                    275:      $         GO TO 50
                    276:             AJJ = SQRT( AJJ )
                    277:             AB( 1, J ) = AJJ
                    278:             KM = MIN( J-1, KD )
                    279: *
                    280: *           Compute elements j-km:j-1 of the j-th row and update the
                    281: *           trailing submatrix within the band.
                    282: *
                    283:             CALL DSCAL( KM, ONE / AJJ, AB( KM+1, J-KM ), KLD )
                    284:             CALL DSYR( 'Lower', KM, -ONE, AB( KM+1, J-KM ), KLD,
                    285:      $                 AB( 1, J-KM ), KLD )
                    286:    30    CONTINUE
                    287: *
                    288: *        Factorize the updated submatrix A(1:m,1:m) as U**T*U.
                    289: *
                    290:          DO 40 J = 1, M
                    291: *
                    292: *           Compute s(j,j) and test for non-positive-definiteness.
                    293: *
                    294:             AJJ = AB( 1, J )
                    295:             IF( AJJ.LE.ZERO )
                    296:      $         GO TO 50
                    297:             AJJ = SQRT( AJJ )
                    298:             AB( 1, J ) = AJJ
                    299:             KM = MIN( KD, M-J )
                    300: *
                    301: *           Compute elements j+1:j+km of the j-th column and update the
                    302: *           trailing submatrix within the band.
                    303: *
                    304:             IF( KM.GT.0 ) THEN
                    305:                CALL DSCAL( KM, ONE / AJJ, AB( 2, J ), 1 )
                    306:                CALL DSYR( 'Lower', KM, -ONE, AB( 2, J ), 1,
                    307:      $                    AB( 1, J+1 ), KLD )
                    308:             END IF
                    309:    40    CONTINUE
                    310:       END IF
                    311:       RETURN
                    312: *
                    313:    50 CONTINUE
                    314:       INFO = J
                    315:       RETURN
                    316: *
                    317: *     End of DPBSTF
                    318: *
                    319:       END

CVSweb interface <joel.bertrand@systella.fr>