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: *> \ingroup doubleOTHERcomputational
114: *
115: *> \par Further Details:
116: * =====================
117: *>
118: *> \verbatim
119: *>
120: *> The band storage scheme is illustrated by the following example, when
121: *> N = 7, KD = 2:
122: *>
123: *> S = ( s11 s12 s13 )
124: *> ( s22 s23 s24 )
125: *> ( s33 s34 )
126: *> ( s44 )
127: *> ( s53 s54 s55 )
128: *> ( s64 s65 s66 )
129: *> ( s75 s76 s77 )
130: *>
131: *> If UPLO = 'U', the array AB holds:
132: *>
133: *> on entry: on exit:
134: *>
135: *> * * a13 a24 a35 a46 a57 * * s13 s24 s53 s64 s75
136: *> * a12 a23 a34 a45 a56 a67 * s12 s23 s34 s54 s65 s76
137: *> a11 a22 a33 a44 a55 a66 a77 s11 s22 s33 s44 s55 s66 s77
138: *>
139: *> If UPLO = 'L', the array AB holds:
140: *>
141: *> on entry: on exit:
142: *>
143: *> a11 a22 a33 a44 a55 a66 a77 s11 s22 s33 s44 s55 s66 s77
144: *> a21 a32 a43 a54 a65 a76 * s12 s23 s34 s54 s65 s76 *
145: *> a31 a42 a53 a64 a64 * * s13 s24 s53 s64 s75 * *
146: *>
147: *> Array elements marked * are not used by the routine.
148: *> \endverbatim
149: *>
150: * =====================================================================
151: SUBROUTINE DPBSTF( UPLO, N, KD, AB, LDAB, INFO )
152: *
153: * -- LAPACK computational routine --
154: * -- LAPACK is a software package provided by Univ. of Tennessee, --
155: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
156: *
157: * .. Scalar Arguments ..
158: CHARACTER UPLO
159: INTEGER INFO, KD, LDAB, N
160: * ..
161: * .. Array Arguments ..
162: DOUBLE PRECISION AB( LDAB, * )
163: * ..
164: *
165: * =====================================================================
166: *
167: * .. Parameters ..
168: DOUBLE PRECISION ONE, ZERO
169: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
170: * ..
171: * .. Local Scalars ..
172: LOGICAL UPPER
173: INTEGER J, KLD, KM, M
174: DOUBLE PRECISION AJJ
175: * ..
176: * .. External Functions ..
177: LOGICAL LSAME
178: EXTERNAL LSAME
179: * ..
180: * .. External Subroutines ..
181: EXTERNAL DSCAL, DSYR, XERBLA
182: * ..
183: * .. Intrinsic Functions ..
184: INTRINSIC MAX, MIN, SQRT
185: * ..
186: * .. Executable Statements ..
187: *
188: * Test the input parameters.
189: *
190: INFO = 0
191: UPPER = LSAME( UPLO, 'U' )
192: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
193: INFO = -1
194: ELSE IF( N.LT.0 ) THEN
195: INFO = -2
196: ELSE IF( KD.LT.0 ) THEN
197: INFO = -3
198: ELSE IF( LDAB.LT.KD+1 ) THEN
199: INFO = -5
200: END IF
201: IF( INFO.NE.0 ) THEN
202: CALL XERBLA( 'DPBSTF', -INFO )
203: RETURN
204: END IF
205: *
206: * Quick return if possible
207: *
208: IF( N.EQ.0 )
209: $ RETURN
210: *
211: KLD = MAX( 1, LDAB-1 )
212: *
213: * Set the splitting point m.
214: *
215: M = ( N+KD ) / 2
216: *
217: IF( UPPER ) THEN
218: *
219: * Factorize A(m+1:n,m+1:n) as L**T*L, and update A(1:m,1:m).
220: *
221: DO 10 J = N, M + 1, -1
222: *
223: * Compute s(j,j) and test for non-positive-definiteness.
224: *
225: AJJ = AB( KD+1, J )
226: IF( AJJ.LE.ZERO )
227: $ GO TO 50
228: AJJ = SQRT( AJJ )
229: AB( KD+1, J ) = AJJ
230: KM = MIN( J-1, KD )
231: *
232: * Compute elements j-km:j-1 of the j-th column and update the
233: * the leading submatrix within the band.
234: *
235: CALL DSCAL( KM, ONE / AJJ, AB( KD+1-KM, J ), 1 )
236: CALL DSYR( 'Upper', KM, -ONE, AB( KD+1-KM, J ), 1,
237: $ AB( KD+1, J-KM ), KLD )
238: 10 CONTINUE
239: *
240: * Factorize the updated submatrix A(1:m,1:m) as U**T*U.
241: *
242: DO 20 J = 1, M
243: *
244: * Compute s(j,j) and test for non-positive-definiteness.
245: *
246: AJJ = AB( KD+1, J )
247: IF( AJJ.LE.ZERO )
248: $ GO TO 50
249: AJJ = SQRT( AJJ )
250: AB( KD+1, J ) = AJJ
251: KM = MIN( KD, M-J )
252: *
253: * Compute elements j+1:j+km of the j-th row and update the
254: * trailing submatrix within the band.
255: *
256: IF( KM.GT.0 ) THEN
257: CALL DSCAL( KM, ONE / AJJ, AB( KD, J+1 ), KLD )
258: CALL DSYR( 'Upper', KM, -ONE, AB( KD, J+1 ), KLD,
259: $ AB( KD+1, J+1 ), KLD )
260: END IF
261: 20 CONTINUE
262: ELSE
263: *
264: * Factorize A(m+1:n,m+1:n) as L**T*L, and update A(1:m,1:m).
265: *
266: DO 30 J = N, M + 1, -1
267: *
268: * Compute s(j,j) and test for non-positive-definiteness.
269: *
270: AJJ = AB( 1, J )
271: IF( AJJ.LE.ZERO )
272: $ GO TO 50
273: AJJ = SQRT( AJJ )
274: AB( 1, J ) = AJJ
275: KM = MIN( J-1, KD )
276: *
277: * Compute elements j-km:j-1 of the j-th row and update the
278: * trailing submatrix within the band.
279: *
280: CALL DSCAL( KM, ONE / AJJ, AB( KM+1, J-KM ), KLD )
281: CALL DSYR( 'Lower', KM, -ONE, AB( KM+1, J-KM ), KLD,
282: $ AB( 1, J-KM ), KLD )
283: 30 CONTINUE
284: *
285: * Factorize the updated submatrix A(1:m,1:m) as U**T*U.
286: *
287: DO 40 J = 1, M
288: *
289: * Compute s(j,j) and test for non-positive-definiteness.
290: *
291: AJJ = AB( 1, J )
292: IF( AJJ.LE.ZERO )
293: $ GO TO 50
294: AJJ = SQRT( AJJ )
295: AB( 1, J ) = AJJ
296: KM = MIN( KD, M-J )
297: *
298: * Compute elements j+1:j+km of the j-th column and update the
299: * trailing submatrix within the band.
300: *
301: IF( KM.GT.0 ) THEN
302: CALL DSCAL( KM, ONE / AJJ, AB( 2, J ), 1 )
303: CALL DSYR( 'Lower', KM, -ONE, AB( 2, J ), 1,
304: $ AB( 1, J+1 ), KLD )
305: END IF
306: 40 CONTINUE
307: END IF
308: RETURN
309: *
310: 50 CONTINUE
311: INFO = J
312: RETURN
313: *
314: * End of DPBSTF
315: *
316: 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