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: *> \date November 2011
111: *
112: *> \ingroup doubleOTHERcomputational
113: *
114: *> \par Further Details:
115: * =====================
116: *>
117: *> \verbatim
118: *>
119: *> If UPLO = 'U', then A = U*D*U**T, where
120: *> U = P(n)*U(n)* ... *P(k)U(k)* ...,
121: *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
122: *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
123: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
124: *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
125: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
126: *>
127: *> ( I v 0 ) k-s
128: *> U(k) = ( 0 I 0 ) s
129: *> ( 0 0 I ) n-k
130: *> k-s s n-k
131: *>
132: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
133: *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
134: *> and A(k,k), and v overwrites A(1:k-2,k-1:k).
135: *>
136: *> If UPLO = 'L', then A = L*D*L**T, where
137: *> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
138: *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
139: *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
140: *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
141: *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
142: *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
143: *>
144: *> ( I 0 0 ) k-1
145: *> L(k) = ( 0 I 0 ) s
146: *> ( 0 v I ) n-k-s+1
147: *> k-1 s n-k-s+1
148: *>
149: *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
150: *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
151: *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
152: *> \endverbatim
153: *
154: *> \par Contributors:
155: * ==================
156: *>
157: *> J. Lewis, Boeing Computer Services Company
158: *>
159: * =====================================================================
160: SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO )
161: *
162: * -- LAPACK computational routine (version 3.4.0) --
163: * -- LAPACK is a software package provided by Univ. of Tennessee, --
164: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165: * November 2011
166: *
167: * .. Scalar Arguments ..
168: CHARACTER UPLO
169: INTEGER INFO, N
170: * ..
171: * .. Array Arguments ..
172: INTEGER IPIV( * )
173: DOUBLE PRECISION AP( * )
174: * ..
175: *
176: * =====================================================================
177: *
178: * .. Parameters ..
179: DOUBLE PRECISION ZERO, ONE
180: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
181: DOUBLE PRECISION EIGHT, SEVTEN
182: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
183: * ..
184: * .. Local Scalars ..
185: LOGICAL UPPER
186: INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
187: $ KSTEP, KX, NPP
188: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
189: $ ROWMAX, T, WK, WKM1, WKP1
190: * ..
191: * .. External Functions ..
192: LOGICAL LSAME
193: INTEGER IDAMAX
194: EXTERNAL LSAME, IDAMAX
195: * ..
196: * .. External Subroutines ..
197: EXTERNAL DSCAL, DSPR, DSWAP, XERBLA
198: * ..
199: * .. Intrinsic Functions ..
200: INTRINSIC ABS, MAX, SQRT
201: * ..
202: * .. Executable Statements ..
203: *
204: * Test the input parameters.
205: *
206: INFO = 0
207: UPPER = LSAME( UPLO, 'U' )
208: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
209: INFO = -1
210: ELSE IF( N.LT.0 ) THEN
211: INFO = -2
212: END IF
213: IF( INFO.NE.0 ) THEN
214: CALL XERBLA( 'DSPTRF', -INFO )
215: RETURN
216: END IF
217: *
218: * Initialize ALPHA for use in choosing pivot block size.
219: *
220: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
221: *
222: IF( UPPER ) THEN
223: *
224: * Factorize A as U*D*U**T using the upper triangle of A
225: *
226: * K is the main loop index, decreasing from N to 1 in steps of
227: * 1 or 2
228: *
229: K = N
230: KC = ( N-1 )*N / 2 + 1
231: 10 CONTINUE
232: KNC = KC
233: *
234: * If K < 1, exit from loop
235: *
236: IF( K.LT.1 )
237: $ GO TO 110
238: KSTEP = 1
239: *
240: * Determine rows and columns to be interchanged and whether
241: * a 1-by-1 or 2-by-2 pivot block will be used
242: *
243: ABSAKK = ABS( AP( KC+K-1 ) )
244: *
245: * IMAX is the row-index of the largest off-diagonal element in
246: * column K, and COLMAX is its absolute value
247: *
248: IF( K.GT.1 ) THEN
249: IMAX = IDAMAX( K-1, AP( KC ), 1 )
250: COLMAX = ABS( AP( KC+IMAX-1 ) )
251: ELSE
252: COLMAX = ZERO
253: END IF
254: *
255: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
256: *
257: * Column K is zero: set INFO and continue
258: *
259: IF( INFO.EQ.0 )
260: $ INFO = K
261: KP = K
262: ELSE
263: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
264: *
265: * no interchange, use 1-by-1 pivot block
266: *
267: KP = K
268: ELSE
269: *
270: ROWMAX = ZERO
271: JMAX = IMAX
272: KX = IMAX*( IMAX+1 ) / 2 + IMAX
273: DO 20 J = IMAX + 1, K
274: IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
275: ROWMAX = ABS( AP( KX ) )
276: JMAX = J
277: END IF
278: KX = KX + J
279: 20 CONTINUE
280: KPC = ( IMAX-1 )*IMAX / 2 + 1
281: IF( IMAX.GT.1 ) THEN
282: JMAX = IDAMAX( IMAX-1, AP( KPC ), 1 )
283: ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-1 ) ) )
284: END IF
285: *
286: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
287: *
288: * no interchange, use 1-by-1 pivot block
289: *
290: KP = K
291: ELSE IF( ABS( AP( KPC+IMAX-1 ) ).GE.ALPHA*ROWMAX ) THEN
292: *
293: * interchange rows and columns K and IMAX, use 1-by-1
294: * pivot block
295: *
296: KP = IMAX
297: ELSE
298: *
299: * interchange rows and columns K-1 and IMAX, use 2-by-2
300: * pivot block
301: *
302: KP = IMAX
303: KSTEP = 2
304: END IF
305: END IF
306: *
307: KK = K - KSTEP + 1
308: IF( KSTEP.EQ.2 )
309: $ KNC = KNC - K + 1
310: IF( KP.NE.KK ) THEN
311: *
312: * Interchange rows and columns KK and KP in the leading
313: * submatrix A(1:k,1:k)
314: *
315: CALL DSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 )
316: KX = KPC + KP - 1
317: DO 30 J = KP + 1, KK - 1
318: KX = KX + J - 1
319: T = AP( KNC+J-1 )
320: AP( KNC+J-1 ) = AP( KX )
321: AP( KX ) = T
322: 30 CONTINUE
323: T = AP( KNC+KK-1 )
324: AP( KNC+KK-1 ) = AP( KPC+KP-1 )
325: AP( KPC+KP-1 ) = T
326: IF( KSTEP.EQ.2 ) THEN
327: T = AP( KC+K-2 )
328: AP( KC+K-2 ) = AP( KC+KP-1 )
329: AP( KC+KP-1 ) = T
330: END IF
331: END IF
332: *
333: * Update the leading submatrix
334: *
335: IF( KSTEP.EQ.1 ) THEN
336: *
337: * 1-by-1 pivot block D(k): column k now holds
338: *
339: * W(k) = U(k)*D(k)
340: *
341: * where U(k) is the k-th column of U
342: *
343: * Perform a rank-1 update of A(1:k-1,1:k-1) as
344: *
345: * A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T
346: *
347: R1 = ONE / AP( KC+K-1 )
348: CALL DSPR( UPLO, K-1, -R1, AP( KC ), 1, AP )
349: *
350: * Store U(k) in column k
351: *
352: CALL DSCAL( K-1, R1, AP( KC ), 1 )
353: ELSE
354: *
355: * 2-by-2 pivot block D(k): columns k and k-1 now hold
356: *
357: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
358: *
359: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
360: * of U
361: *
362: * Perform a rank-2 update of A(1:k-2,1:k-2) as
363: *
364: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
365: * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T
366: *
367: IF( K.GT.2 ) THEN
368: *
369: D12 = AP( K-1+( K-1 )*K / 2 )
370: D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12
371: D11 = AP( K+( K-1 )*K / 2 ) / D12
372: T = ONE / ( D11*D22-ONE )
373: D12 = T / D12
374: *
375: DO 50 J = K - 2, 1, -1
376: WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
377: $ AP( J+( K-1 )*K / 2 ) )
378: WK = D12*( D22*AP( J+( K-1 )*K / 2 )-
379: $ AP( J+( K-2 )*( K-1 ) / 2 ) )
380: DO 40 I = J, 1, -1
381: AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) -
382: $ AP( I+( K-1 )*K / 2 )*WK -
383: $ AP( I+( K-2 )*( K-1 ) / 2 )*WKM1
384: 40 CONTINUE
385: AP( J+( K-1 )*K / 2 ) = WK
386: AP( J+( K-2 )*( K-1 ) / 2 ) = WKM1
387: 50 CONTINUE
388: *
389: END IF
390: *
391: END IF
392: END IF
393: *
394: * Store details of the interchanges in IPIV
395: *
396: IF( KSTEP.EQ.1 ) THEN
397: IPIV( K ) = KP
398: ELSE
399: IPIV( K ) = -KP
400: IPIV( K-1 ) = -KP
401: END IF
402: *
403: * Decrease K and return to the start of the main loop
404: *
405: K = K - KSTEP
406: KC = KNC - K
407: GO TO 10
408: *
409: ELSE
410: *
411: * Factorize A as L*D*L**T using the lower triangle of A
412: *
413: * K is the main loop index, increasing from 1 to N in steps of
414: * 1 or 2
415: *
416: K = 1
417: KC = 1
418: NPP = N*( N+1 ) / 2
419: 60 CONTINUE
420: KNC = KC
421: *
422: * If K > N, exit from loop
423: *
424: IF( K.GT.N )
425: $ GO TO 110
426: KSTEP = 1
427: *
428: * Determine rows and columns to be interchanged and whether
429: * a 1-by-1 or 2-by-2 pivot block will be used
430: *
431: ABSAKK = ABS( AP( KC ) )
432: *
433: * IMAX is the row-index of the largest off-diagonal element in
434: * column K, and COLMAX is its absolute value
435: *
436: IF( K.LT.N ) THEN
437: IMAX = K + IDAMAX( N-K, AP( KC+1 ), 1 )
438: COLMAX = ABS( AP( KC+IMAX-K ) )
439: ELSE
440: COLMAX = ZERO
441: END IF
442: *
443: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
444: *
445: * Column K is zero: set INFO and continue
446: *
447: IF( INFO.EQ.0 )
448: $ INFO = K
449: KP = K
450: ELSE
451: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
452: *
453: * no interchange, use 1-by-1 pivot block
454: *
455: KP = K
456: ELSE
457: *
458: * JMAX is the column-index of the largest off-diagonal
459: * element in row IMAX, and ROWMAX is its absolute value
460: *
461: ROWMAX = ZERO
462: KX = KC + IMAX - K
463: DO 70 J = K, IMAX - 1
464: IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
465: ROWMAX = ABS( AP( KX ) )
466: JMAX = J
467: END IF
468: KX = KX + N - J
469: 70 CONTINUE
470: KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1
471: IF( IMAX.LT.N ) THEN
472: JMAX = IMAX + IDAMAX( N-IMAX, AP( KPC+1 ), 1 )
473: ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-IMAX ) ) )
474: END IF
475: *
476: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
477: *
478: * no interchange, use 1-by-1 pivot block
479: *
480: KP = K
481: ELSE IF( ABS( AP( KPC ) ).GE.ALPHA*ROWMAX ) THEN
482: *
483: * interchange rows and columns K and IMAX, use 1-by-1
484: * pivot block
485: *
486: KP = IMAX
487: ELSE
488: *
489: * interchange rows and columns K+1 and IMAX, use 2-by-2
490: * pivot block
491: *
492: KP = IMAX
493: KSTEP = 2
494: END IF
495: END IF
496: *
497: KK = K + KSTEP - 1
498: IF( KSTEP.EQ.2 )
499: $ KNC = KNC + N - K + 1
500: IF( KP.NE.KK ) THEN
501: *
502: * Interchange rows and columns KK and KP in the trailing
503: * submatrix A(k:n,k:n)
504: *
505: IF( KP.LT.N )
506: $ CALL DSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ),
507: $ 1 )
508: KX = KNC + KP - KK
509: DO 80 J = KK + 1, KP - 1
510: KX = KX + N - J + 1
511: T = AP( KNC+J-KK )
512: AP( KNC+J-KK ) = AP( KX )
513: AP( KX ) = T
514: 80 CONTINUE
515: T = AP( KNC )
516: AP( KNC ) = AP( KPC )
517: AP( KPC ) = T
518: IF( KSTEP.EQ.2 ) THEN
519: T = AP( KC+1 )
520: AP( KC+1 ) = AP( KC+KP-K )
521: AP( KC+KP-K ) = T
522: END IF
523: END IF
524: *
525: * Update the trailing submatrix
526: *
527: IF( KSTEP.EQ.1 ) THEN
528: *
529: * 1-by-1 pivot block D(k): column k now holds
530: *
531: * W(k) = L(k)*D(k)
532: *
533: * where L(k) is the k-th column of L
534: *
535: IF( K.LT.N ) THEN
536: *
537: * Perform a rank-1 update of A(k+1:n,k+1:n) as
538: *
539: * A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T
540: *
541: R1 = ONE / AP( KC )
542: CALL DSPR( UPLO, N-K, -R1, AP( KC+1 ), 1,
543: $ AP( KC+N-K+1 ) )
544: *
545: * Store L(k) in column K
546: *
547: CALL DSCAL( N-K, R1, AP( KC+1 ), 1 )
548: END IF
549: ELSE
550: *
551: * 2-by-2 pivot block D(k): columns K and K+1 now hold
552: *
553: * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
554: *
555: * where L(k) and L(k+1) are the k-th and (k+1)-th columns
556: * of L
557: *
558: IF( K.LT.N-1 ) THEN
559: *
560: * Perform a rank-2 update of A(k+2:n,k+2:n) as
561: *
562: * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**T
563: * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**T
564: *
565: * where L(k) and L(k+1) are the k-th and (k+1)-th
566: * columns of L
567: *
568: D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 )
569: D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21
570: D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21
571: T = ONE / ( D11*D22-ONE )
572: D21 = T / D21
573: *
574: DO 100 J = K + 2, N
575: WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-
576: $ AP( J+K*( 2*N-K-1 ) / 2 ) )
577: WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
578: $ AP( J+( K-1 )*( 2*N-K ) / 2 ) )
579: *
580: DO 90 I = J, N
581: AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )*
582: $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) /
583: $ 2 )*WK - AP( I+K*( 2*N-K-1 ) / 2 )*WKP1
584: 90 CONTINUE
585: *
586: AP( J+( K-1 )*( 2*N-K ) / 2 ) = WK
587: AP( J+K*( 2*N-K-1 ) / 2 ) = WKP1
588: *
589: 100 CONTINUE
590: END IF
591: END IF
592: END IF
593: *
594: * Store details of the interchanges in IPIV
595: *
596: IF( KSTEP.EQ.1 ) THEN
597: IPIV( K ) = KP
598: ELSE
599: IPIV( K ) = -KP
600: IPIV( K+1 ) = -KP
601: END IF
602: *
603: * Increase K and return to the start of the main loop
604: *
605: K = K + KSTEP
606: KC = KNC + N - K + 2
607: GO TO 60
608: *
609: END IF
610: *
611: 110 CONTINUE
612: RETURN
613: *
614: * End of DSPTRF
615: *
616: 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