1: *> \brief \b DSPTRF
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DSPTRF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsptrf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsptrf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsptrf.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO )
22: *
23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER INFO, N
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * )
29: * DOUBLE PRECISION AP( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DSPTRF computes the factorization of a real symmetric matrix A stored
39: *> in packed format using the Bunch-Kaufman diagonal pivoting method:
40: *>
41: *> A = U*D*U**T or A = L*D*L**T
42: *>
43: *> where U (or L) is a product of permutation and unit upper (lower)
44: *> triangular matrices, and D is symmetric and block diagonal with
45: *> 1-by-1 and 2-by-2 diagonal blocks.
46: *> \endverbatim
47: *
48: * Arguments:
49: * ==========
50: *
51: *> \param[in] UPLO
52: *> \verbatim
53: *> UPLO is CHARACTER*1
54: *> = 'U': Upper triangle of A is stored;
55: *> = 'L': Lower triangle of A is stored.
56: *> \endverbatim
57: *>
58: *> \param[in] N
59: *> \verbatim
60: *> N is INTEGER
61: *> The order of the matrix A. N >= 0.
62: *> \endverbatim
63: *>
64: *> \param[in,out] AP
65: *> \verbatim
66: *> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
67: *> On entry, the upper or lower triangle of the symmetric matrix
68: *> A, packed columnwise in a linear array. The j-th column of A
69: *> is stored in the array AP as follows:
70: *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
71: *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
72: *>
73: *> On exit, the block diagonal matrix D and the multipliers used
74: *> to obtain the factor U or L, stored as a packed triangular
75: *> matrix overwriting A (see below for further details).
76: *> \endverbatim
77: *>
78: *> \param[out] IPIV
79: *> \verbatim
80: *> IPIV is INTEGER array, dimension (N)
81: *> Details of the interchanges and the block structure of D.
82: *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
83: *> interchanged and D(k,k) is a 1-by-1 diagonal block.
84: *> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
85: *> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
86: *> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
87: *> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
88: *> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
89: *> \endverbatim
90: *>
91: *> \param[out] INFO
92: *> \verbatim
93: *> INFO is INTEGER
94: *> = 0: successful exit
95: *> < 0: if INFO = -i, the i-th argument had an illegal value
96: *> > 0: if INFO = i, D(i,i) is exactly zero. The factorization
97: *> has been completed, but the block diagonal matrix D is
98: *> exactly singular, and division by zero will occur if it
99: *> is used to solve a system of equations.
100: *> \endverbatim
101: *
102: * Authors:
103: * ========
104: *
105: *> \author Univ. of Tennessee
106: *> \author Univ. of California Berkeley
107: *> \author Univ. of Colorado Denver
108: *> \author NAG Ltd.
109: *
110: *> \ingroup doubleOTHERcomputational
111: *
112: *> \par Further Details:
113: * =====================
114: *>
115: *> \verbatim
116: *>
117: *> If UPLO = 'U', then A = U*D*U**T, where
118: *> U = P(n)*U(n)* ... *P(k)U(k)* ...,
119: *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
120: *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
121: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
122: *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
123: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
124: *>
125: *> ( I v 0 ) k-s
126: *> U(k) = ( 0 I 0 ) s
127: *> ( 0 0 I ) n-k
128: *> k-s s n-k
129: *>
130: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
131: *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
132: *> and A(k,k), and v overwrites A(1:k-2,k-1:k).
133: *>
134: *> If UPLO = 'L', then A = L*D*L**T, where
135: *> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
136: *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
137: *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
138: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
139: *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
140: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
141: *>
142: *> ( I 0 0 ) k-1
143: *> L(k) = ( 0 I 0 ) s
144: *> ( 0 v I ) n-k-s+1
145: *> k-1 s n-k-s+1
146: *>
147: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
148: *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
149: *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
150: *> \endverbatim
151: *
152: *> \par Contributors:
153: * ==================
154: *>
155: *> J. Lewis, Boeing Computer Services Company
156: *>
157: * =====================================================================
158: SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO )
159: *
160: * -- LAPACK computational routine --
161: * -- LAPACK is a software package provided by Univ. of Tennessee, --
162: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
163: *
164: * .. Scalar Arguments ..
165: CHARACTER UPLO
166: INTEGER INFO, N
167: * ..
168: * .. Array Arguments ..
169: INTEGER IPIV( * )
170: DOUBLE PRECISION AP( * )
171: * ..
172: *
173: * =====================================================================
174: *
175: * .. Parameters ..
176: DOUBLE PRECISION ZERO, ONE
177: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
178: DOUBLE PRECISION EIGHT, SEVTEN
179: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
180: * ..
181: * .. Local Scalars ..
182: LOGICAL UPPER
183: INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
184: $ KSTEP, KX, NPP
185: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
186: $ ROWMAX, T, WK, WKM1, WKP1
187: * ..
188: * .. External Functions ..
189: LOGICAL LSAME
190: INTEGER IDAMAX
191: EXTERNAL LSAME, IDAMAX
192: * ..
193: * .. External Subroutines ..
194: EXTERNAL DSCAL, DSPR, DSWAP, XERBLA
195: * ..
196: * .. Intrinsic Functions ..
197: INTRINSIC ABS, MAX, SQRT
198: * ..
199: * .. Executable Statements ..
200: *
201: * Test the input parameters.
202: *
203: INFO = 0
204: UPPER = LSAME( UPLO, 'U' )
205: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
206: INFO = -1
207: ELSE IF( N.LT.0 ) THEN
208: INFO = -2
209: END IF
210: IF( INFO.NE.0 ) THEN
211: CALL XERBLA( 'DSPTRF', -INFO )
212: RETURN
213: END IF
214: *
215: * Initialize ALPHA for use in choosing pivot block size.
216: *
217: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
218: *
219: IF( UPPER ) THEN
220: *
221: * Factorize A as U*D*U**T using the upper triangle of A
222: *
223: * K is the main loop index, decreasing from N to 1 in steps of
224: * 1 or 2
225: *
226: K = N
227: KC = ( N-1 )*N / 2 + 1
228: 10 CONTINUE
229: KNC = KC
230: *
231: * If K < 1, exit from loop
232: *
233: IF( K.LT.1 )
234: $ GO TO 110
235: KSTEP = 1
236: *
237: * Determine rows and columns to be interchanged and whether
238: * a 1-by-1 or 2-by-2 pivot block will be used
239: *
240: ABSAKK = ABS( AP( KC+K-1 ) )
241: *
242: * IMAX is the row-index of the largest off-diagonal element in
243: * column K, and COLMAX is its absolute value
244: *
245: IF( K.GT.1 ) THEN
246: IMAX = IDAMAX( K-1, AP( KC ), 1 )
247: COLMAX = ABS( AP( KC+IMAX-1 ) )
248: ELSE
249: COLMAX = ZERO
250: END IF
251: *
252: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
253: *
254: * Column K is zero: set INFO and continue
255: *
256: IF( INFO.EQ.0 )
257: $ INFO = K
258: KP = K
259: ELSE
260: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
261: *
262: * no interchange, use 1-by-1 pivot block
263: *
264: KP = K
265: ELSE
266: *
267: ROWMAX = ZERO
268: JMAX = IMAX
269: KX = IMAX*( IMAX+1 ) / 2 + IMAX
270: DO 20 J = IMAX + 1, K
271: IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
272: ROWMAX = ABS( AP( KX ) )
273: JMAX = J
274: END IF
275: KX = KX + J
276: 20 CONTINUE
277: KPC = ( IMAX-1 )*IMAX / 2 + 1
278: IF( IMAX.GT.1 ) THEN
279: JMAX = IDAMAX( IMAX-1, AP( KPC ), 1 )
280: ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-1 ) ) )
281: END IF
282: *
283: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
284: *
285: * no interchange, use 1-by-1 pivot block
286: *
287: KP = K
288: ELSE IF( ABS( AP( KPC+IMAX-1 ) ).GE.ALPHA*ROWMAX ) THEN
289: *
290: * interchange rows and columns K and IMAX, use 1-by-1
291: * pivot block
292: *
293: KP = IMAX
294: ELSE
295: *
296: * interchange rows and columns K-1 and IMAX, use 2-by-2
297: * pivot block
298: *
299: KP = IMAX
300: KSTEP = 2
301: END IF
302: END IF
303: *
304: KK = K - KSTEP + 1
305: IF( KSTEP.EQ.2 )
306: $ KNC = KNC - K + 1
307: IF( KP.NE.KK ) THEN
308: *
309: * Interchange rows and columns KK and KP in the leading
310: * submatrix A(1:k,1:k)
311: *
312: CALL DSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 )
313: KX = KPC + KP - 1
314: DO 30 J = KP + 1, KK - 1
315: KX = KX + J - 1
316: T = AP( KNC+J-1 )
317: AP( KNC+J-1 ) = AP( KX )
318: AP( KX ) = T
319: 30 CONTINUE
320: T = AP( KNC+KK-1 )
321: AP( KNC+KK-1 ) = AP( KPC+KP-1 )
322: AP( KPC+KP-1 ) = T
323: IF( KSTEP.EQ.2 ) THEN
324: T = AP( KC+K-2 )
325: AP( KC+K-2 ) = AP( KC+KP-1 )
326: AP( KC+KP-1 ) = T
327: END IF
328: END IF
329: *
330: * Update the leading submatrix
331: *
332: IF( KSTEP.EQ.1 ) THEN
333: *
334: * 1-by-1 pivot block D(k): column k now holds
335: *
336: * W(k) = U(k)*D(k)
337: *
338: * where U(k) is the k-th column of U
339: *
340: * Perform a rank-1 update of A(1:k-1,1:k-1) as
341: *
342: * A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
343: *
344: R1 = ONE / AP( KC+K-1 )
345: CALL DSPR( UPLO, K-1, -R1, AP( KC ), 1, AP )
346: *
347: * Store U(k) in column k
348: *
349: CALL DSCAL( K-1, R1, AP( KC ), 1 )
350: ELSE
351: *
352: * 2-by-2 pivot block D(k): columns k and k-1 now hold
353: *
354: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
355: *
356: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
357: * of U
358: *
359: * Perform a rank-2 update of A(1:k-2,1:k-2) as
360: *
361: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
362: * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
363: *
364: IF( K.GT.2 ) THEN
365: *
366: D12 = AP( K-1+( K-1 )*K / 2 )
367: D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12
368: D11 = AP( K+( K-1 )*K / 2 ) / D12
369: T = ONE / ( D11*D22-ONE )
370: D12 = T / D12
371: *
372: DO 50 J = K - 2, 1, -1
373: WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
374: $ AP( J+( K-1 )*K / 2 ) )
375: WK = D12*( D22*AP( J+( K-1 )*K / 2 )-
376: $ AP( J+( K-2 )*( K-1 ) / 2 ) )
377: DO 40 I = J, 1, -1
378: AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) -
379: $ AP( I+( K-1 )*K / 2 )*WK -
380: $ AP( I+( K-2 )*( K-1 ) / 2 )*WKM1
381: 40 CONTINUE
382: AP( J+( K-1 )*K / 2 ) = WK
383: AP( J+( K-2 )*( K-1 ) / 2 ) = WKM1
384: 50 CONTINUE
385: *
386: END IF
387: *
388: END IF
389: END IF
390: *
391: * Store details of the interchanges in IPIV
392: *
393: IF( KSTEP.EQ.1 ) THEN
394: IPIV( K ) = KP
395: ELSE
396: IPIV( K ) = -KP
397: IPIV( K-1 ) = -KP
398: END IF
399: *
400: * Decrease K and return to the start of the main loop
401: *
402: K = K - KSTEP
403: KC = KNC - K
404: GO TO 10
405: *
406: ELSE
407: *
408: * Factorize A as L*D*L**T using the lower triangle of A
409: *
410: * K is the main loop index, increasing from 1 to N in steps of
411: * 1 or 2
412: *
413: K = 1
414: KC = 1
415: NPP = N*( N+1 ) / 2
416: 60 CONTINUE
417: KNC = KC
418: *
419: * If K > N, exit from loop
420: *
421: IF( K.GT.N )
422: $ GO TO 110
423: KSTEP = 1
424: *
425: * Determine rows and columns to be interchanged and whether
426: * a 1-by-1 or 2-by-2 pivot block will be used
427: *
428: ABSAKK = ABS( AP( KC ) )
429: *
430: * IMAX is the row-index of the largest off-diagonal element in
431: * column K, and COLMAX is its absolute value
432: *
433: IF( K.LT.N ) THEN
434: IMAX = K + IDAMAX( N-K, AP( KC+1 ), 1 )
435: COLMAX = ABS( AP( KC+IMAX-K ) )
436: ELSE
437: COLMAX = ZERO
438: END IF
439: *
440: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
441: *
442: * Column K is zero: set INFO and continue
443: *
444: IF( INFO.EQ.0 )
445: $ INFO = K
446: KP = K
447: ELSE
448: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
449: *
450: * no interchange, use 1-by-1 pivot block
451: *
452: KP = K
453: ELSE
454: *
455: * JMAX is the column-index of the largest off-diagonal
456: * element in row IMAX, and ROWMAX is its absolute value
457: *
458: ROWMAX = ZERO
459: KX = KC + IMAX - K
460: DO 70 J = K, IMAX - 1
461: IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
462: ROWMAX = ABS( AP( KX ) )
463: JMAX = J
464: END IF
465: KX = KX + N - J
466: 70 CONTINUE
467: KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1
468: IF( IMAX.LT.N ) THEN
469: JMAX = IMAX + IDAMAX( N-IMAX, AP( KPC+1 ), 1 )
470: ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-IMAX ) ) )
471: END IF
472: *
473: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
474: *
475: * no interchange, use 1-by-1 pivot block
476: *
477: KP = K
478: ELSE IF( ABS( AP( KPC ) ).GE.ALPHA*ROWMAX ) THEN
479: *
480: * interchange rows and columns K and IMAX, use 1-by-1
481: * pivot block
482: *
483: KP = IMAX
484: ELSE
485: *
486: * interchange rows and columns K+1 and IMAX, use 2-by-2
487: * pivot block
488: *
489: KP = IMAX
490: KSTEP = 2
491: END IF
492: END IF
493: *
494: KK = K + KSTEP - 1
495: IF( KSTEP.EQ.2 )
496: $ KNC = KNC + N - K + 1
497: IF( KP.NE.KK ) THEN
498: *
499: * Interchange rows and columns KK and KP in the trailing
500: * submatrix A(k:n,k:n)
501: *
502: IF( KP.LT.N )
503: $ CALL DSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ),
504: $ 1 )
505: KX = KNC + KP - KK
506: DO 80 J = KK + 1, KP - 1
507: KX = KX + N - J + 1
508: T = AP( KNC+J-KK )
509: AP( KNC+J-KK ) = AP( KX )
510: AP( KX ) = T
511: 80 CONTINUE
512: T = AP( KNC )
513: AP( KNC ) = AP( KPC )
514: AP( KPC ) = T
515: IF( KSTEP.EQ.2 ) THEN
516: T = AP( KC+1 )
517: AP( KC+1 ) = AP( KC+KP-K )
518: AP( KC+KP-K ) = T
519: END IF
520: END IF
521: *
522: * Update the trailing submatrix
523: *
524: IF( KSTEP.EQ.1 ) THEN
525: *
526: * 1-by-1 pivot block D(k): column k now holds
527: *
528: * W(k) = L(k)*D(k)
529: *
530: * where L(k) is the k-th column of L
531: *
532: IF( K.LT.N ) THEN
533: *
534: * Perform a rank-1 update of A(k+1:n,k+1:n) as
535: *
536: * A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
537: *
538: R1 = ONE / AP( KC )
539: CALL DSPR( UPLO, N-K, -R1, AP( KC+1 ), 1,
540: $ AP( KC+N-K+1 ) )
541: *
542: * Store L(k) in column K
543: *
544: CALL DSCAL( N-K, R1, AP( KC+1 ), 1 )
545: END IF
546: ELSE
547: *
548: * 2-by-2 pivot block D(k): columns K and K+1 now hold
549: *
550: * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
551: *
552: * where L(k) and L(k+1) are the k-th and (k+1)-th columns
553: * of L
554: *
555: IF( K.LT.N-1 ) THEN
556: *
557: * Perform a rank-2 update of A(k+2:n,k+2:n) as
558: *
559: * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T
560: * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T
561: *
562: * where L(k) and L(k+1) are the k-th and (k+1)-th
563: * columns of L
564: *
565: D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 )
566: D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21
567: D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21
568: T = ONE / ( D11*D22-ONE )
569: D21 = T / D21
570: *
571: DO 100 J = K + 2, N
572: WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-
573: $ AP( J+K*( 2*N-K-1 ) / 2 ) )
574: WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
575: $ AP( J+( K-1 )*( 2*N-K ) / 2 ) )
576: *
577: DO 90 I = J, N
578: AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )*
579: $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) /
580: $ 2 )*WK - AP( I+K*( 2*N-K-1 ) / 2 )*WKP1
581: 90 CONTINUE
582: *
583: AP( J+( K-1 )*( 2*N-K ) / 2 ) = WK
584: AP( J+K*( 2*N-K-1 ) / 2 ) = WKP1
585: *
586: 100 CONTINUE
587: END IF
588: END IF
589: END IF
590: *
591: * Store details of the interchanges in IPIV
592: *
593: IF( KSTEP.EQ.1 ) THEN
594: IPIV( K ) = KP
595: ELSE
596: IPIV( K ) = -KP
597: IPIV( K+1 ) = -KP
598: END IF
599: *
600: * Increase K and return to the start of the main loop
601: *
602: K = K + KSTEP
603: KC = KNC + N - K + 2
604: GO TO 60
605: *
606: END IF
607: *
608: 110 CONTINUE
609: RETURN
610: *
611: * End of DSPTRF
612: *
613: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsptrf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack