|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 |
SUBROUTINE DPBTRS( UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO )
* * -- LAPACK routine (version 3.3.1) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- April 2011 -- * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, KD, LDAB, LDB, N, NRHS * .. * .. Array Arguments .. DOUBLE PRECISION AB( LDAB, * ), B( LDB, * ) * .. * * Purpose * ======= * * DPBTRS solves a system of linear equations A*X = B with a symmetric * positive definite band matrix A using the Cholesky factorization * A = U**T*U or A = L*L**T computed by DPBTRF. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * = 'U': Upper triangular factor stored in AB; * = 'L': Lower triangular factor stored in AB. * * N (input) INTEGER * The order of the matrix A. N >= 0. * * KD (input) INTEGER * The number of superdiagonals of the matrix A if UPLO = 'U', * or the number of subdiagonals if UPLO = 'L'. KD >= 0. * * NRHS (input) INTEGER * The number of right hand sides, i.e., the number of columns * of the matrix B. NRHS >= 0. * * AB (input) DOUBLE PRECISION array, dimension (LDAB,N) * The triangular factor U or L from the Cholesky factorization * A = U**T*U or A = L*L**T of the band matrix A, stored in the * first KD+1 rows of the array. The j-th column of U or L is * stored in the j-th column of the array AB as follows: * if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j; * if UPLO ='L', AB(1+i-j,j) = L(i,j) for j<=i<=min(n,j+kd). * * LDAB (input) INTEGER * The leading dimension of the array AB. LDAB >= KD+1. * * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) * On entry, the right hand side matrix B. * On exit, the solution matrix X. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * ===================================================================== * * .. Local Scalars .. LOGICAL UPPER INTEGER J * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DTBSV, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( KD.LT.0 ) THEN INFO = -3 ELSE IF( NRHS.LT.0 ) THEN INFO = -4 ELSE IF( LDAB.LT.KD+1 ) THEN INFO = -6 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -8 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DPBTRS', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 .OR. NRHS.EQ.0 ) $ RETURN * IF( UPPER ) THEN * * Solve A*X = B where A = U**T *U. * DO 10 J = 1, NRHS * * Solve U**T *X = B, overwriting B with X. * CALL DTBSV( 'Upper', 'Transpose', 'Non-unit', N, KD, AB, $ LDAB, B( 1, J ), 1 ) * * Solve U*X = B, overwriting B with X. * CALL DTBSV( 'Upper', 'No transpose', 'Non-unit', N, KD, AB, $ LDAB, B( 1, J ), 1 ) 10 CONTINUE ELSE * * Solve A*X = B where A = L*L**T. * DO 20 J = 1, NRHS * * Solve L*X = B, overwriting B with X. * CALL DTBSV( 'Lower', 'No transpose', 'Non-unit', N, KD, AB, $ LDAB, B( 1, J ), 1 ) * * Solve L**T *X = B, overwriting B with X. * CALL DTBSV( 'Lower', 'Transpose', 'Non-unit', N, KD, AB, $ LDAB, B( 1, J ), 1 ) 20 CONTINUE END IF * RETURN * * End of DPBTRS * END |