1 SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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**T U22**T )
26 *
27 * A = ( L11 0 ) ( D 0 ) ( L11**T L21**T ) 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**T = A11 - U12*W**T
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**T = A22 - L21*W**T
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
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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**T U22**T )
26 *
27 * A = ( L11 0 ) ( D 0 ) ( L11**T L21**T ) 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**T = A11 - U12*W**T
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**T = A22 - L21*W**T
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