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