Return to dpbtrf.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack vers la version 3.3.0.

    1:       SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
    2: *
    3: *  -- LAPACK routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       CHARACTER          UPLO
   10:       INTEGER            INFO, KD, LDAB, N
   11: *     ..
   12: *     .. Array Arguments ..
   13:       DOUBLE PRECISION   AB( LDAB, * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  DPBTRF computes the Cholesky factorization of a real symmetric
   20: *  positive definite band matrix A.
   21: *
   22: *  The factorization has the form
   23: *     A = U**T * U,  if UPLO = 'U', or
   24: *     A = L  * L**T,  if UPLO = 'L',
   25: *  where U is an upper triangular matrix and L is lower triangular.
   26: *
   27: *  Arguments
   28: *  =========
   29: *
   30: *  UPLO    (input) CHARACTER*1
   31: *          = 'U':  Upper triangle of A is stored;
   32: *          = 'L':  Lower triangle of A is stored.
   33: *
   34: *  N       (input) INTEGER
   35: *          The order of the matrix A.  N >= 0.
   36: *
   37: *  KD      (input) INTEGER
   38: *          The number of superdiagonals of the matrix A if UPLO = 'U',
   39: *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
   40: *
   41: *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
   42: *          On entry, the upper or lower triangle of the symmetric band
   43: *          matrix A, stored in the first KD+1 rows of the array.  The
   44: *          j-th column of A is stored in the j-th column of the array AB
   45: *          as follows:
   46: *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
   47: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
   48: *
   49: *          On exit, if INFO = 0, the triangular factor U or L from the
   50: *          Cholesky factorization A = U**T*U or A = L*L**T of the band
   51: *          matrix A, in the same storage format as A.
   52: *
   53: *  LDAB    (input) INTEGER
   54: *          The leading dimension of the array AB.  LDAB >= KD+1.
   55: *
   56: *  INFO    (output) INTEGER
   57: *          = 0:  successful exit
   58: *          < 0:  if INFO = -i, the i-th argument had an illegal value
   59: *          > 0:  if INFO = i, the leading minor of order i is not
   60: *                positive definite, and the factorization could not be
   61: *                completed.
   62: *
   63: *  Further Details
   64: *  ===============
   65: *
   66: *  The band storage scheme is illustrated by the following example, when
   67: *  N = 6, KD = 2, and UPLO = 'U':
   68: *
   69: *  On entry:                       On exit:
   70: *
   71: *      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
   72: *      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
   73: *     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
   74: *
   75: *  Similarly, if UPLO = 'L' the format of A is as follows:
   76: *
   77: *  On entry:                       On exit:
   78: *
   79: *     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
   80: *     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
   81: *     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
   82: *
   83: *  Array elements marked * are not used by the routine.
   84: *
   85: *  Contributed by
   86: *  Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989
   87: *
   88: *  =====================================================================
   89: *
   90: *     .. Parameters ..
   91:       DOUBLE PRECISION   ONE, ZERO
   92:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
   93:       INTEGER            NBMAX, LDWORK
   94:       PARAMETER          ( NBMAX = 32, LDWORK = NBMAX+1 )
   95: *     ..
   96: *     .. Local Scalars ..
   97:       INTEGER            I, I2, I3, IB, II, J, JJ, NB
   98: *     ..
   99: *     .. Local Arrays ..
  100:       DOUBLE PRECISION   WORK( LDWORK, NBMAX )
  101: *     ..
  102: *     .. External Functions ..
  103:       LOGICAL            LSAME
  104:       INTEGER            ILAENV
  105:       EXTERNAL           LSAME, ILAENV
  106: *     ..
  107: *     .. External Subroutines ..
  108:       EXTERNAL           DGEMM, DPBTF2, DPOTF2, DSYRK, DTRSM, XERBLA
  109: *     ..
  110: *     .. Intrinsic Functions ..
  111:       INTRINSIC          MIN
  112: *     ..
  113: *     .. Executable Statements ..
  114: *
  115: *     Test the input parameters.
  116: *
  117:       INFO = 0
  118:       IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND.
  119:      $    ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
  120:          INFO = -1
  121:       ELSE IF( N.LT.0 ) THEN
  122:          INFO = -2
  123:       ELSE IF( KD.LT.0 ) THEN
  124:          INFO = -3
  125:       ELSE IF( LDAB.LT.KD+1 ) THEN
  126:          INFO = -5
  127:       END IF
  128:       IF( INFO.NE.0 ) THEN
  129:          CALL XERBLA( 'DPBTRF', -INFO )
  130:          RETURN
  131:       END IF
  132: *
  133: *     Quick return if possible
  134: *
  135:       IF( N.EQ.0 )
  136:      $   RETURN
  137: *
  138: *     Determine the block size for this environment
  139: *
  140:       NB = ILAENV( 1, 'DPBTRF', UPLO, N, KD, -1, -1 )
  141: *
  142: *     The block size must not exceed the semi-bandwidth KD, and must not
  143: *     exceed the limit set by the size of the local array WORK.
  144: *
  145:       NB = MIN( NB, NBMAX )
  146: *
  147:       IF( NB.LE.1 .OR. NB.GT.KD ) THEN
  148: *
  149: *        Use unblocked code
  150: *
  151:          CALL DPBTF2( UPLO, N, KD, AB, LDAB, INFO )
  152:       ELSE
  153: *
  154: *        Use blocked code
  155: *
  156:          IF( LSAME( UPLO, 'U' ) ) THEN
  157: *
  158: *           Compute the Cholesky factorization of a symmetric band
  159: *           matrix, given the upper triangle of the matrix in band
  160: *           storage.
  161: *
  162: *           Zero the upper triangle of the work array.
  163: *
  164:             DO 20 J = 1, NB
  165:                DO 10 I = 1, J - 1
  166:                   WORK( I, J ) = ZERO
  167:    10          CONTINUE
  168:    20       CONTINUE
  169: *
  170: *           Process the band matrix one diagonal block at a time.
  171: *
  172:             DO 70 I = 1, N, NB
  173:                IB = MIN( NB, N-I+1 )
  174: *
  175: *              Factorize the diagonal block
  176: *
  177:                CALL DPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II )
  178:                IF( II.NE.0 ) THEN
  179:                   INFO = I + II - 1
  180:                   GO TO 150
  181:                END IF
  182:                IF( I+IB.LE.N ) THEN
  183: *
  184: *                 Update the relevant part of the trailing submatrix.
  185: *                 If A11 denotes the diagonal block which has just been
  186: *                 factorized, then we need to update the remaining
  187: *                 blocks in the diagram:
  188: *
  189: *                    A11   A12   A13
  190: *                          A22   A23
  191: *                                A33
  192: *
  193: *                 The numbers of rows and columns in the partitioning
  194: *                 are IB, I2, I3 respectively. The blocks A12, A22 and
  195: *                 A23 are empty if IB = KD. The upper triangle of A13
  196: *                 lies outside the band.
  197: *
  198:                   I2 = MIN( KD-IB, N-I-IB+1 )
  199:                   I3 = MIN( IB, N-I-KD+1 )
  200: *
  201:                   IF( I2.GT.0 ) THEN
  202: *
  203: *                    Update A12
  204: *
  205:                      CALL DTRSM( 'Left', 'Upper', 'Transpose',
  206:      $                           'Non-unit', IB, I2, ONE, AB( KD+1, I ),
  207:      $                           LDAB-1, AB( KD+1-IB, I+IB ), LDAB-1 )
  208: *
  209: *                    Update A22
  210: *
  211:                      CALL DSYRK( 'Upper', 'Transpose', I2, IB, -ONE,
  212:      $                           AB( KD+1-IB, I+IB ), LDAB-1, ONE,
  213:      $                           AB( KD+1, I+IB ), LDAB-1 )
  214:                   END IF
  215: *
  216:                   IF( I3.GT.0 ) THEN
  217: *
  218: *                    Copy the lower triangle of A13 into the work array.
  219: *
  220:                      DO 40 JJ = 1, I3
  221:                         DO 30 II = JJ, IB
  222:                            WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 )
  223:    30                   CONTINUE
  224:    40                CONTINUE
  225: *
  226: *                    Update A13 (in the work array).
  227: *
  228:                      CALL DTRSM( 'Left', 'Upper', 'Transpose',
  229:      $                           'Non-unit', IB, I3, ONE, AB( KD+1, I ),
  230:      $                           LDAB-1, WORK, LDWORK )
  231: *
  232: *                    Update A23
  233: *
  234:                      IF( I2.GT.0 )
  235:      $                  CALL DGEMM( 'Transpose', 'No Transpose', I2, I3,
  236:      $                              IB, -ONE, AB( KD+1-IB, I+IB ),
  237:      $                              LDAB-1, WORK, LDWORK, ONE,
  238:      $                              AB( 1+IB, I+KD ), LDAB-1 )
  239: *
  240: *                    Update A33
  241: *
  242:                      CALL DSYRK( 'Upper', 'Transpose', I3, IB, -ONE,
  243:      $                           WORK, LDWORK, ONE, AB( KD+1, I+KD ),
  244:      $                           LDAB-1 )
  245: *
  246: *                    Copy the lower triangle of A13 back into place.
  247: *
  248:                      DO 60 JJ = 1, I3
  249:                         DO 50 II = JJ, IB
  250:                            AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ )
  251:    50                   CONTINUE
  252:    60                CONTINUE
  253:                   END IF
  254:                END IF
  255:    70       CONTINUE
  256:          ELSE
  257: *
  258: *           Compute the Cholesky factorization of a symmetric band
  259: *           matrix, given the lower triangle of the matrix in band
  260: *           storage.
  261: *
  262: *           Zero the lower triangle of the work array.
  263: *
  264:             DO 90 J = 1, NB
  265:                DO 80 I = J + 1, NB
  266:                   WORK( I, J ) = ZERO
  267:    80          CONTINUE
  268:    90       CONTINUE
  269: *
  270: *           Process the band matrix one diagonal block at a time.
  271: *
  272:             DO 140 I = 1, N, NB
  273:                IB = MIN( NB, N-I+1 )
  274: *
  275: *              Factorize the diagonal block
  276: *
  277:                CALL DPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II )
  278:                IF( II.NE.0 ) THEN
  279:                   INFO = I + II - 1
  280:                   GO TO 150
  281:                END IF
  282:                IF( I+IB.LE.N ) THEN
  283: *
  284: *                 Update the relevant part of the trailing submatrix.
  285: *                 If A11 denotes the diagonal block which has just been
  286: *                 factorized, then we need to update the remaining
  287: *                 blocks in the diagram:
  288: *
  289: *                    A11
  290: *                    A21   A22
  291: *                    A31   A32   A33
  292: *
  293: *                 The numbers of rows and columns in the partitioning
  294: *                 are IB, I2, I3 respectively. The blocks A21, A22 and
  295: *                 A32 are empty if IB = KD. The lower triangle of A31
  296: *                 lies outside the band.
  297: *
  298:                   I2 = MIN( KD-IB, N-I-IB+1 )
  299:                   I3 = MIN( IB, N-I-KD+1 )
  300: *
  301:                   IF( I2.GT.0 ) THEN
  302: *
  303: *                    Update A21
  304: *
  305:                      CALL DTRSM( 'Right', 'Lower', 'Transpose',
  306:      $                           'Non-unit', I2, IB, ONE, AB( 1, I ),
  307:      $                           LDAB-1, AB( 1+IB, I ), LDAB-1 )
  308: *
  309: *                    Update A22
  310: *
  311:                      CALL DSYRK( 'Lower', 'No Transpose', I2, IB, -ONE,
  312:      $                           AB( 1+IB, I ), LDAB-1, ONE,
  313:      $                           AB( 1, I+IB ), LDAB-1 )
  314:                   END IF
  315: *
  316:                   IF( I3.GT.0 ) THEN
  317: *
  318: *                    Copy the upper triangle of A31 into the work array.
  319: *
  320:                      DO 110 JJ = 1, IB
  321:                         DO 100 II = 1, MIN( JJ, I3 )
  322:                            WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 )
  323:   100                   CONTINUE
  324:   110                CONTINUE
  325: *
  326: *                    Update A31 (in the work array).
  327: *
  328:                      CALL DTRSM( 'Right', 'Lower', 'Transpose',
  329:      $                           'Non-unit', I3, IB, ONE, AB( 1, I ),
  330:      $                           LDAB-1, WORK, LDWORK )
  331: *
  332: *                    Update A32
  333: *
  334:                      IF( I2.GT.0 )
  335:      $                  CALL DGEMM( 'No transpose', 'Transpose', I3, I2,
  336:      $                              IB, -ONE, WORK, LDWORK,
  337:      $                              AB( 1+IB, I ), LDAB-1, ONE,
  338:      $                              AB( 1+KD-IB, I+IB ), LDAB-1 )
  339: *
  340: *                    Update A33
  341: *
  342:                      CALL DSYRK( 'Lower', 'No Transpose', I3, IB, -ONE,
  343:      $                           WORK, LDWORK, ONE, AB( 1, I+KD ),
  344:      $                           LDAB-1 )
  345: *
  346: *                    Copy the upper triangle of A31 back into place.
  347: *
  348:                      DO 130 JJ = 1, IB
  349:                         DO 120 II = 1, MIN( JJ, I3 )
  350:                            AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ )
  351:   120                   CONTINUE
  352:   130                CONTINUE
  353:                   END IF
  354:                END IF
  355:   140       CONTINUE
  356:          END IF
  357:       END IF
  358:       RETURN
  359: *
  360:   150 CONTINUE
  361:       RETURN
  362: *
  363: *     End of DPBTRF
  364: *
  365:       END