1: SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: CHARACTER UPLO
10: INTEGER INFO, KB, LDA, LDW, N, NB
11: * ..
12: * .. Array Arguments ..
13: INTEGER IPIV( * )
14: DOUBLE PRECISION A( LDA, * ), W( LDW, * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DLASYF computes a partial factorization of a real symmetric matrix A
21: * using the Bunch-Kaufman diagonal pivoting method. The partial
22: * factorization has the form:
23: *
24: * A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or:
25: * ( 0 U22 ) ( 0 D ) ( U12' U22' )
26: *
27: * A = ( L11 0 ) ( D 0 ) ( L11' L21' ) if UPLO = 'L'
28: * ( L21 I ) ( 0 A22 ) ( 0 I )
29: *
30: * where the order of D is at most NB. The actual order is returned in
31: * the argument KB, and is either NB or NB-1, or N if N <= NB.
32: *
33: * DLASYF is an auxiliary routine called by DSYTRF. It uses blocked code
34: * (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
35: * A22 (if UPLO = 'L').
36: *
37: * Arguments
38: * =========
39: *
40: * UPLO (input) CHARACTER*1
41: * Specifies whether the upper or lower triangular part of the
42: * symmetric matrix A is stored:
43: * = 'U': Upper triangular
44: * = 'L': Lower triangular
45: *
46: * N (input) INTEGER
47: * The order of the matrix A. N >= 0.
48: *
49: * NB (input) INTEGER
50: * The maximum number of columns of the matrix A that should be
51: * factored. NB should be at least 2 to allow for 2-by-2 pivot
52: * blocks.
53: *
54: * KB (output) INTEGER
55: * The number of columns of A that were actually factored.
56: * KB is either NB-1 or NB, or N if N <= NB.
57: *
58: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
59: * On entry, the symmetric matrix A. If UPLO = 'U', the leading
60: * n-by-n upper triangular part of A contains the upper
61: * triangular part of the matrix A, and the strictly lower
62: * triangular part of A is not referenced. If UPLO = 'L', the
63: * leading n-by-n lower triangular part of A contains the lower
64: * triangular part of the matrix A, and the strictly upper
65: * triangular part of A is not referenced.
66: * On exit, A contains details of the partial factorization.
67: *
68: * LDA (input) INTEGER
69: * The leading dimension of the array A. LDA >= max(1,N).
70: *
71: * IPIV (output) INTEGER array, dimension (N)
72: * Details of the interchanges and the block structure of D.
73: * If UPLO = 'U', only the last KB elements of IPIV are set;
74: * if UPLO = 'L', only the first KB elements are set.
75: *
76: * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
77: * interchanged and D(k,k) is a 1-by-1 diagonal block.
78: * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
79: * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
80: * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
81: * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
82: * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
83: *
84: * W (workspace) DOUBLE PRECISION array, dimension (LDW,NB)
85: *
86: * LDW (input) INTEGER
87: * The leading dimension of the array W. LDW >= max(1,N).
88: *
89: * INFO (output) INTEGER
90: * = 0: successful exit
91: * > 0: if INFO = k, D(k,k) is exactly zero. The factorization
92: * has been completed, but the block diagonal matrix D is
93: * exactly singular.
94: *
95: * =====================================================================
96: *
97: * .. Parameters ..
98: DOUBLE PRECISION ZERO, ONE
99: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
100: DOUBLE PRECISION EIGHT, SEVTEN
101: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
102: * ..
103: * .. Local Scalars ..
104: INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
105: $ KSTEP, KW
106: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1,
107: $ ROWMAX, T
108: * ..
109: * .. External Functions ..
110: LOGICAL LSAME
111: INTEGER IDAMAX
112: EXTERNAL LSAME, IDAMAX
113: * ..
114: * .. External Subroutines ..
115: EXTERNAL DCOPY, DGEMM, DGEMV, DSCAL, DSWAP
116: * ..
117: * .. Intrinsic Functions ..
118: INTRINSIC ABS, MAX, MIN, SQRT
119: * ..
120: * .. Executable Statements ..
121: *
122: INFO = 0
123: *
124: * Initialize ALPHA for use in choosing pivot block size.
125: *
126: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
127: *
128: IF( LSAME( UPLO, 'U' ) ) THEN
129: *
130: * Factorize the trailing columns of A using the upper triangle
131: * of A and working backwards, and compute the matrix W = U12*D
132: * for use in updating A11
133: *
134: * K is the main loop index, decreasing from N in steps of 1 or 2
135: *
136: * KW is the column of W which corresponds to column K of A
137: *
138: K = N
139: 10 CONTINUE
140: KW = NB + K - N
141: *
142: * Exit from loop
143: *
144: IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
145: $ GO TO 30
146: *
147: * Copy column K of A to column KW of W and update it
148: *
149: CALL DCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 )
150: IF( K.LT.N )
151: $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ), LDA,
152: $ W( K, KW+1 ), LDW, ONE, W( 1, KW ), 1 )
153: *
154: KSTEP = 1
155: *
156: * Determine rows and columns to be interchanged and whether
157: * a 1-by-1 or 2-by-2 pivot block will be used
158: *
159: ABSAKK = ABS( W( K, KW ) )
160: *
161: * IMAX is the row-index of the largest off-diagonal element in
162: * column K, and COLMAX is its absolute value
163: *
164: IF( K.GT.1 ) THEN
165: IMAX = IDAMAX( K-1, W( 1, KW ), 1 )
166: COLMAX = ABS( W( IMAX, KW ) )
167: ELSE
168: COLMAX = ZERO
169: END IF
170: *
171: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
172: *
173: * Column K is zero: set INFO and continue
174: *
175: IF( INFO.EQ.0 )
176: $ INFO = K
177: KP = K
178: ELSE
179: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
180: *
181: * no interchange, use 1-by-1 pivot block
182: *
183: KP = K
184: ELSE
185: *
186: * Copy column IMAX to column KW-1 of W and update it
187: *
188: CALL DCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
189: CALL DCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
190: $ W( IMAX+1, KW-1 ), 1 )
191: IF( K.LT.N )
192: $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ),
193: $ LDA, W( IMAX, KW+1 ), LDW, ONE,
194: $ W( 1, KW-1 ), 1 )
195: *
196: * JMAX is the column-index of the largest off-diagonal
197: * element in row IMAX, and ROWMAX is its absolute value
198: *
199: JMAX = IMAX + IDAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
200: ROWMAX = ABS( W( JMAX, KW-1 ) )
201: IF( IMAX.GT.1 ) THEN
202: JMAX = IDAMAX( IMAX-1, W( 1, KW-1 ), 1 )
203: ROWMAX = MAX( ROWMAX, ABS( W( JMAX, KW-1 ) ) )
204: END IF
205: *
206: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
207: *
208: * no interchange, use 1-by-1 pivot block
209: *
210: KP = K
211: ELSE IF( ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
212: *
213: * interchange rows and columns K and IMAX, use 1-by-1
214: * pivot block
215: *
216: KP = IMAX
217: *
218: * copy column KW-1 of W to column KW
219: *
220: CALL DCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
221: ELSE
222: *
223: * interchange rows and columns K-1 and IMAX, use 2-by-2
224: * pivot block
225: *
226: KP = IMAX
227: KSTEP = 2
228: END IF
229: END IF
230: *
231: KK = K - KSTEP + 1
232: KKW = NB + KK - N
233: *
234: * Updated column KP is already stored in column KKW of W
235: *
236: IF( KP.NE.KK ) THEN
237: *
238: * Copy non-updated column KK to column KP
239: *
240: A( KP, K ) = A( KK, K )
241: CALL DCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
242: $ LDA )
243: CALL DCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
244: *
245: * Interchange rows KK and KP in last KK columns of A and W
246: *
247: CALL DSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
248: CALL DSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
249: $ LDW )
250: END IF
251: *
252: IF( KSTEP.EQ.1 ) THEN
253: *
254: * 1-by-1 pivot block D(k): column KW of W now holds
255: *
256: * W(k) = U(k)*D(k)
257: *
258: * where U(k) is the k-th column of U
259: *
260: * Store U(k) in column k of A
261: *
262: CALL DCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
263: R1 = ONE / A( K, K )
264: CALL DSCAL( K-1, R1, A( 1, K ), 1 )
265: ELSE
266: *
267: * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
268: * hold
269: *
270: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
271: *
272: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
273: * of U
274: *
275: IF( K.GT.2 ) THEN
276: *
277: * Store U(k) and U(k-1) in columns k and k-1 of A
278: *
279: D21 = W( K-1, KW )
280: D11 = W( K, KW ) / D21
281: D22 = W( K-1, KW-1 ) / D21
282: T = ONE / ( D11*D22-ONE )
283: D21 = T / D21
284: DO 20 J = 1, K - 2
285: A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
286: A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
287: 20 CONTINUE
288: END IF
289: *
290: * Copy D(k) to A
291: *
292: A( K-1, K-1 ) = W( K-1, KW-1 )
293: A( K-1, K ) = W( K-1, KW )
294: A( K, K ) = W( K, KW )
295: END IF
296: END IF
297: *
298: * Store details of the interchanges in IPIV
299: *
300: IF( KSTEP.EQ.1 ) THEN
301: IPIV( K ) = KP
302: ELSE
303: IPIV( K ) = -KP
304: IPIV( K-1 ) = -KP
305: END IF
306: *
307: * Decrease K and return to the start of the main loop
308: *
309: K = K - KSTEP
310: GO TO 10
311: *
312: 30 CONTINUE
313: *
314: * Update the upper triangle of A11 (= A(1:k,1:k)) as
315: *
316: * A11 := A11 - U12*D*U12' = A11 - U12*W'
317: *
318: * computing blocks of NB columns at a time
319: *
320: DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
321: JB = MIN( NB, K-J+1 )
322: *
323: * Update the upper triangle of the diagonal block
324: *
325: DO 40 JJ = J, J + JB - 1
326: CALL DGEMV( 'No transpose', JJ-J+1, N-K, -ONE,
327: $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
328: $ A( J, JJ ), 1 )
329: 40 CONTINUE
330: *
331: * Update the rectangular superdiagonal block
332: *
333: CALL DGEMM( 'No transpose', 'Transpose', J-1, JB, N-K, -ONE,
334: $ A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, ONE,
335: $ A( 1, J ), LDA )
336: 50 CONTINUE
337: *
338: * Put U12 in standard form by partially undoing the interchanges
339: * in columns k+1:n
340: *
341: J = K + 1
342: 60 CONTINUE
343: JJ = J
344: JP = IPIV( J )
345: IF( JP.LT.0 ) THEN
346: JP = -JP
347: J = J + 1
348: END IF
349: J = J + 1
350: IF( JP.NE.JJ .AND. J.LE.N )
351: $ CALL DSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
352: IF( J.LE.N )
353: $ GO TO 60
354: *
355: * Set KB to the number of columns factorized
356: *
357: KB = N - K
358: *
359: ELSE
360: *
361: * Factorize the leading columns of A using the lower triangle
362: * of A and working forwards, and compute the matrix W = L21*D
363: * for use in updating A22
364: *
365: * K is the main loop index, increasing from 1 in steps of 1 or 2
366: *
367: K = 1
368: 70 CONTINUE
369: *
370: * Exit from loop
371: *
372: IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
373: $ GO TO 90
374: *
375: * Copy column K of A to column K of W and update it
376: *
377: CALL DCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
378: CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ), LDA,
379: $ W( K, 1 ), LDW, ONE, W( K, K ), 1 )
380: *
381: KSTEP = 1
382: *
383: * Determine rows and columns to be interchanged and whether
384: * a 1-by-1 or 2-by-2 pivot block will be used
385: *
386: ABSAKK = ABS( W( K, K ) )
387: *
388: * IMAX is the row-index of the largest off-diagonal element in
389: * column K, and COLMAX is its absolute value
390: *
391: IF( K.LT.N ) THEN
392: IMAX = K + IDAMAX( N-K, W( K+1, K ), 1 )
393: COLMAX = ABS( W( IMAX, K ) )
394: ELSE
395: COLMAX = ZERO
396: END IF
397: *
398: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
399: *
400: * Column K is zero: set INFO and continue
401: *
402: IF( INFO.EQ.0 )
403: $ INFO = K
404: KP = K
405: ELSE
406: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
407: *
408: * no interchange, use 1-by-1 pivot block
409: *
410: KP = K
411: ELSE
412: *
413: * Copy column IMAX to column K+1 of W and update it
414: *
415: CALL DCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
416: CALL DCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
417: $ 1 )
418: CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ),
419: $ LDA, W( IMAX, 1 ), LDW, ONE, W( K, K+1 ), 1 )
420: *
421: * JMAX is the column-index of the largest off-diagonal
422: * element in row IMAX, and ROWMAX is its absolute value
423: *
424: JMAX = K - 1 + IDAMAX( IMAX-K, W( K, K+1 ), 1 )
425: ROWMAX = ABS( W( JMAX, K+1 ) )
426: IF( IMAX.LT.N ) THEN
427: JMAX = IMAX + IDAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
428: ROWMAX = MAX( ROWMAX, ABS( W( JMAX, K+1 ) ) )
429: END IF
430: *
431: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
432: *
433: * no interchange, use 1-by-1 pivot block
434: *
435: KP = K
436: ELSE IF( ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
437: *
438: * interchange rows and columns K and IMAX, use 1-by-1
439: * pivot block
440: *
441: KP = IMAX
442: *
443: * copy column K+1 of W to column K
444: *
445: CALL DCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
446: ELSE
447: *
448: * interchange rows and columns K+1 and IMAX, use 2-by-2
449: * pivot block
450: *
451: KP = IMAX
452: KSTEP = 2
453: END IF
454: END IF
455: *
456: KK = K + KSTEP - 1
457: *
458: * Updated column KP is already stored in column KK of W
459: *
460: IF( KP.NE.KK ) THEN
461: *
462: * Copy non-updated column KK to column KP
463: *
464: A( KP, K ) = A( KK, K )
465: CALL DCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
466: CALL DCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
467: *
468: * Interchange rows KK and KP in first KK columns of A and W
469: *
470: CALL DSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
471: CALL DSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
472: END IF
473: *
474: IF( KSTEP.EQ.1 ) THEN
475: *
476: * 1-by-1 pivot block D(k): column k of W now holds
477: *
478: * W(k) = L(k)*D(k)
479: *
480: * where L(k) is the k-th column of L
481: *
482: * Store L(k) in column k of A
483: *
484: CALL DCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
485: IF( K.LT.N ) THEN
486: R1 = ONE / A( K, K )
487: CALL DSCAL( N-K, R1, A( K+1, K ), 1 )
488: END IF
489: ELSE
490: *
491: * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
492: *
493: * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
494: *
495: * where L(k) and L(k+1) are the k-th and (k+1)-th columns
496: * of L
497: *
498: IF( K.LT.N-1 ) THEN
499: *
500: * Store L(k) and L(k+1) in columns k and k+1 of A
501: *
502: D21 = W( K+1, K )
503: D11 = W( K+1, K+1 ) / D21
504: D22 = W( K, K ) / D21
505: T = ONE / ( D11*D22-ONE )
506: D21 = T / D21
507: DO 80 J = K + 2, N
508: A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
509: A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
510: 80 CONTINUE
511: END IF
512: *
513: * Copy D(k) to A
514: *
515: A( K, K ) = W( K, K )
516: A( K+1, K ) = W( K+1, K )
517: A( K+1, K+1 ) = W( K+1, K+1 )
518: END IF
519: END IF
520: *
521: * Store details of the interchanges in IPIV
522: *
523: IF( KSTEP.EQ.1 ) THEN
524: IPIV( K ) = KP
525: ELSE
526: IPIV( K ) = -KP
527: IPIV( K+1 ) = -KP
528: END IF
529: *
530: * Increase K and return to the start of the main loop
531: *
532: K = K + KSTEP
533: GO TO 70
534: *
535: 90 CONTINUE
536: *
537: * Update the lower triangle of A22 (= A(k:n,k:n)) as
538: *
539: * A22 := A22 - L21*D*L21' = A22 - L21*W'
540: *
541: * computing blocks of NB columns at a time
542: *
543: DO 110 J = K, N, NB
544: JB = MIN( NB, N-J+1 )
545: *
546: * Update the lower triangle of the diagonal block
547: *
548: DO 100 JJ = J, J + JB - 1
549: CALL DGEMV( 'No transpose', J+JB-JJ, K-1, -ONE,
550: $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
551: $ A( JJ, JJ ), 1 )
552: 100 CONTINUE
553: *
554: * Update the rectangular subdiagonal block
555: *
556: IF( J+JB.LE.N )
557: $ CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
558: $ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
559: $ ONE, A( J+JB, J ), LDA )
560: 110 CONTINUE
561: *
562: * Put L21 in standard form by partially undoing the interchanges
563: * in columns 1:k-1
564: *
565: J = K - 1
566: 120 CONTINUE
567: JJ = J
568: JP = IPIV( J )
569: IF( JP.LT.0 ) THEN
570: JP = -JP
571: J = J - 1
572: END IF
573: J = J - 1
574: IF( JP.NE.JJ .AND. J.GE.1 )
575: $ CALL DSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
576: IF( J.GE.1 )
577: $ GO TO 120
578: *
579: * Set KB to the number of columns factorized
580: *
581: KB = K - 1
582: *
583: END IF
584: RETURN
585: *
586: * End of DLASYF
587: *
588: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasyf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack