Annotation of rpl/lapack/lapack/dpbstf.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DPBSTF( 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: * DPBSTF computes a split Cholesky factorization of a real
! 20: * symmetric positive definite band matrix A.
! 21: *
! 22: * This routine is designed to be used in conjunction with DSBGST.
! 23: *
! 24: * The factorization has the form A = S**T*S where S is a band matrix
! 25: * of the same bandwidth as A and the following structure:
! 26: *
! 27: * S = ( U )
! 28: * ( M L )
! 29: *
! 30: * where U is upper triangular of order m = (n+kd)/2, and L is lower
! 31: * triangular of order n-m.
! 32: *
! 33: * Arguments
! 34: * =========
! 35: *
! 36: * UPLO (input) CHARACTER*1
! 37: * = 'U': Upper triangle of A is stored;
! 38: * = 'L': Lower triangle of A is stored.
! 39: *
! 40: * N (input) INTEGER
! 41: * The order of the matrix A. N >= 0.
! 42: *
! 43: * KD (input) INTEGER
! 44: * The number of superdiagonals of the matrix A if UPLO = 'U',
! 45: * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
! 46: *
! 47: * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
! 48: * On entry, the upper or lower triangle of the symmetric band
! 49: * matrix A, stored in the first kd+1 rows of the array. The
! 50: * j-th column of A is stored in the j-th column of the array AB
! 51: * as follows:
! 52: * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
! 53: * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
! 54: *
! 55: * On exit, if INFO = 0, the factor S from the split Cholesky
! 56: * factorization A = S**T*S. See Further Details.
! 57: *
! 58: * LDAB (input) INTEGER
! 59: * The leading dimension of the array AB. LDAB >= KD+1.
! 60: *
! 61: * INFO (output) INTEGER
! 62: * = 0: successful exit
! 63: * < 0: if INFO = -i, the i-th argument had an illegal value
! 64: * > 0: if INFO = i, the factorization could not be completed,
! 65: * because the updated element a(i,i) was negative; the
! 66: * matrix A is not positive definite.
! 67: *
! 68: * Further Details
! 69: * ===============
! 70: *
! 71: * The band storage scheme is illustrated by the following example, when
! 72: * N = 7, KD = 2:
! 73: *
! 74: * S = ( s11 s12 s13 )
! 75: * ( s22 s23 s24 )
! 76: * ( s33 s34 )
! 77: * ( s44 )
! 78: * ( s53 s54 s55 )
! 79: * ( s64 s65 s66 )
! 80: * ( s75 s76 s77 )
! 81: *
! 82: * If UPLO = 'U', the array AB holds:
! 83: *
! 84: * on entry: on exit:
! 85: *
! 86: * * * a13 a24 a35 a46 a57 * * s13 s24 s53 s64 s75
! 87: * * a12 a23 a34 a45 a56 a67 * s12 s23 s34 s54 s65 s76
! 88: * a11 a22 a33 a44 a55 a66 a77 s11 s22 s33 s44 s55 s66 s77
! 89: *
! 90: * If UPLO = 'L', the array AB holds:
! 91: *
! 92: * on entry: on exit:
! 93: *
! 94: * a11 a22 a33 a44 a55 a66 a77 s11 s22 s33 s44 s55 s66 s77
! 95: * a21 a32 a43 a54 a65 a76 * s12 s23 s34 s54 s65 s76 *
! 96: * a31 a42 a53 a64 a64 * * s13 s24 s53 s64 s75 * *
! 97: *
! 98: * Array elements marked * are not used by the routine.
! 99: *
! 100: * =====================================================================
! 101: *
! 102: * .. Parameters ..
! 103: DOUBLE PRECISION ONE, ZERO
! 104: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! 105: * ..
! 106: * .. Local Scalars ..
! 107: LOGICAL UPPER
! 108: INTEGER J, KLD, KM, M
! 109: DOUBLE PRECISION AJJ
! 110: * ..
! 111: * .. External Functions ..
! 112: LOGICAL LSAME
! 113: EXTERNAL LSAME
! 114: * ..
! 115: * .. External Subroutines ..
! 116: EXTERNAL DSCAL, DSYR, XERBLA
! 117: * ..
! 118: * .. Intrinsic Functions ..
! 119: INTRINSIC MAX, MIN, SQRT
! 120: * ..
! 121: * .. Executable Statements ..
! 122: *
! 123: * Test the input parameters.
! 124: *
! 125: INFO = 0
! 126: UPPER = LSAME( UPLO, 'U' )
! 127: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 128: INFO = -1
! 129: ELSE IF( N.LT.0 ) THEN
! 130: INFO = -2
! 131: ELSE IF( KD.LT.0 ) THEN
! 132: INFO = -3
! 133: ELSE IF( LDAB.LT.KD+1 ) THEN
! 134: INFO = -5
! 135: END IF
! 136: IF( INFO.NE.0 ) THEN
! 137: CALL XERBLA( 'DPBSTF', -INFO )
! 138: RETURN
! 139: END IF
! 140: *
! 141: * Quick return if possible
! 142: *
! 143: IF( N.EQ.0 )
! 144: $ RETURN
! 145: *
! 146: KLD = MAX( 1, LDAB-1 )
! 147: *
! 148: * Set the splitting point m.
! 149: *
! 150: M = ( N+KD ) / 2
! 151: *
! 152: IF( UPPER ) THEN
! 153: *
! 154: * Factorize A(m+1:n,m+1:n) as L**T*L, and update A(1:m,1:m).
! 155: *
! 156: DO 10 J = N, M + 1, -1
! 157: *
! 158: * Compute s(j,j) and test for non-positive-definiteness.
! 159: *
! 160: AJJ = AB( KD+1, J )
! 161: IF( AJJ.LE.ZERO )
! 162: $ GO TO 50
! 163: AJJ = SQRT( AJJ )
! 164: AB( KD+1, J ) = AJJ
! 165: KM = MIN( J-1, KD )
! 166: *
! 167: * Compute elements j-km:j-1 of the j-th column and update the
! 168: * the leading submatrix within the band.
! 169: *
! 170: CALL DSCAL( KM, ONE / AJJ, AB( KD+1-KM, J ), 1 )
! 171: CALL DSYR( 'Upper', KM, -ONE, AB( KD+1-KM, J ), 1,
! 172: $ AB( KD+1, J-KM ), KLD )
! 173: 10 CONTINUE
! 174: *
! 175: * Factorize the updated submatrix A(1:m,1:m) as U**T*U.
! 176: *
! 177: DO 20 J = 1, M
! 178: *
! 179: * Compute s(j,j) and test for non-positive-definiteness.
! 180: *
! 181: AJJ = AB( KD+1, J )
! 182: IF( AJJ.LE.ZERO )
! 183: $ GO TO 50
! 184: AJJ = SQRT( AJJ )
! 185: AB( KD+1, J ) = AJJ
! 186: KM = MIN( KD, M-J )
! 187: *
! 188: * Compute elements j+1:j+km of the j-th row and update the
! 189: * trailing submatrix within the band.
! 190: *
! 191: IF( KM.GT.0 ) THEN
! 192: CALL DSCAL( KM, ONE / AJJ, AB( KD, J+1 ), KLD )
! 193: CALL DSYR( 'Upper', KM, -ONE, AB( KD, J+1 ), KLD,
! 194: $ AB( KD+1, J+1 ), KLD )
! 195: END IF
! 196: 20 CONTINUE
! 197: ELSE
! 198: *
! 199: * Factorize A(m+1:n,m+1:n) as L**T*L, and update A(1:m,1:m).
! 200: *
! 201: DO 30 J = N, M + 1, -1
! 202: *
! 203: * Compute s(j,j) and test for non-positive-definiteness.
! 204: *
! 205: AJJ = AB( 1, J )
! 206: IF( AJJ.LE.ZERO )
! 207: $ GO TO 50
! 208: AJJ = SQRT( AJJ )
! 209: AB( 1, J ) = AJJ
! 210: KM = MIN( J-1, KD )
! 211: *
! 212: * Compute elements j-km:j-1 of the j-th row and update the
! 213: * trailing submatrix within the band.
! 214: *
! 215: CALL DSCAL( KM, ONE / AJJ, AB( KM+1, J-KM ), KLD )
! 216: CALL DSYR( 'Lower', KM, -ONE, AB( KM+1, J-KM ), KLD,
! 217: $ AB( 1, J-KM ), KLD )
! 218: 30 CONTINUE
! 219: *
! 220: * Factorize the updated submatrix A(1:m,1:m) as U**T*U.
! 221: *
! 222: DO 40 J = 1, M
! 223: *
! 224: * Compute s(j,j) and test for non-positive-definiteness.
! 225: *
! 226: AJJ = AB( 1, J )
! 227: IF( AJJ.LE.ZERO )
! 228: $ GO TO 50
! 229: AJJ = SQRT( AJJ )
! 230: AB( 1, J ) = AJJ
! 231: KM = MIN( KD, M-J )
! 232: *
! 233: * Compute elements j+1:j+km of the j-th column and update the
! 234: * trailing submatrix within the band.
! 235: *
! 236: IF( KM.GT.0 ) THEN
! 237: CALL DSCAL( KM, ONE / AJJ, AB( 2, J ), 1 )
! 238: CALL DSYR( 'Lower', KM, -ONE, AB( 2, J ), 1,
! 239: $ AB( 1, J+1 ), KLD )
! 240: END IF
! 241: 40 CONTINUE
! 242: END IF
! 243: RETURN
! 244: *
! 245: 50 CONTINUE
! 246: INFO = J
! 247: RETURN
! 248: *
! 249: * End of DPBSTF
! 250: *
! 251: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dpbstf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack