Annotation of rpl/lapack/lapack/dlasyf.f, revision 1.19
1.14 bertrand 1: *> \brief \b DLASYF computes a partial factorization of a real symmetric matrix using the Bunch-Kaufman diagonal pivoting method.
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.14 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.14 bertrand 9: *> Download DLASYF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasyf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasyf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasyf.f">
1.9 bertrand 15: *> [TXT]</a>
1.14 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
1.14 bertrand 22: *
1.9 bertrand 23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER INFO, KB, LDA, LDW, N, NB
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * )
29: * DOUBLE PRECISION A( LDA, * ), W( LDW, * )
30: * ..
1.14 bertrand 31: *
1.9 bertrand 32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DLASYF computes a partial factorization of a real symmetric matrix A
39: *> using the Bunch-Kaufman diagonal pivoting method. The partial
40: *> factorization has the form:
41: *>
42: *> A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or:
43: *> ( 0 U22 ) ( 0 D ) ( U12**T U22**T )
44: *>
45: *> A = ( L11 0 ) ( D 0 ) ( L11**T L21**T ) if UPLO = 'L'
46: *> ( L21 I ) ( 0 A22 ) ( 0 I )
47: *>
48: *> where the order of D is at most NB. The actual order is returned in
49: *> the argument KB, and is either NB or NB-1, or N if N <= NB.
50: *>
51: *> DLASYF is an auxiliary routine called by DSYTRF. It uses blocked code
52: *> (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
53: *> A22 (if UPLO = 'L').
54: *> \endverbatim
55: *
56: * Arguments:
57: * ==========
58: *
59: *> \param[in] UPLO
60: *> \verbatim
61: *> UPLO is CHARACTER*1
62: *> Specifies whether the upper or lower triangular part of the
63: *> symmetric matrix A is stored:
64: *> = 'U': Upper triangular
65: *> = 'L': Lower triangular
66: *> \endverbatim
67: *>
68: *> \param[in] N
69: *> \verbatim
70: *> N is INTEGER
71: *> The order of the matrix A. N >= 0.
72: *> \endverbatim
73: *>
74: *> \param[in] NB
75: *> \verbatim
76: *> NB is INTEGER
77: *> The maximum number of columns of the matrix A that should be
78: *> factored. NB should be at least 2 to allow for 2-by-2 pivot
79: *> blocks.
80: *> \endverbatim
81: *>
82: *> \param[out] KB
83: *> \verbatim
84: *> KB is INTEGER
85: *> The number of columns of A that were actually factored.
86: *> KB is either NB-1 or NB, or N if N <= NB.
87: *> \endverbatim
88: *>
89: *> \param[in,out] A
90: *> \verbatim
91: *> A is DOUBLE PRECISION array, dimension (LDA,N)
92: *> On entry, the symmetric matrix A. If UPLO = 'U', the leading
93: *> n-by-n upper triangular part of A contains the upper
94: *> triangular part of the matrix A, and the strictly lower
95: *> triangular part of A is not referenced. If UPLO = 'L', the
96: *> leading n-by-n lower triangular part of A contains the lower
97: *> triangular part of the matrix A, and the strictly upper
98: *> triangular part of A is not referenced.
99: *> On exit, A contains details of the partial factorization.
100: *> \endverbatim
101: *>
102: *> \param[in] LDA
103: *> \verbatim
104: *> LDA is INTEGER
105: *> The leading dimension of the array A. LDA >= max(1,N).
106: *> \endverbatim
107: *>
108: *> \param[out] IPIV
109: *> \verbatim
110: *> IPIV is INTEGER array, dimension (N)
111: *> Details of the interchanges and the block structure of D.
112: *>
1.14 bertrand 113: *> If UPLO = 'U':
114: *> Only the last KB elements of IPIV are set.
115: *>
116: *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
117: *> interchanged and D(k,k) is a 1-by-1 diagonal block.
118: *>
119: *> If IPIV(k) = IPIV(k-1) < 0, then rows and columns
120: *> k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
121: *> is a 2-by-2 diagonal block.
122: *>
123: *> If UPLO = 'L':
124: *> Only the first KB elements of IPIV are set.
125: *>
126: *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
127: *> interchanged and D(k,k) is a 1-by-1 diagonal block.
128: *>
129: *> If IPIV(k) = IPIV(k+1) < 0, then rows and columns
130: *> k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
131: *> is a 2-by-2 diagonal block.
1.9 bertrand 132: *> \endverbatim
133: *>
134: *> \param[out] W
135: *> \verbatim
136: *> W is DOUBLE PRECISION array, dimension (LDW,NB)
137: *> \endverbatim
138: *>
139: *> \param[in] LDW
140: *> \verbatim
141: *> LDW is INTEGER
142: *> The leading dimension of the array W. LDW >= max(1,N).
143: *> \endverbatim
144: *>
145: *> \param[out] INFO
146: *> \verbatim
147: *> INFO is INTEGER
148: *> = 0: successful exit
149: *> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
150: *> has been completed, but the block diagonal matrix D is
151: *> exactly singular.
152: *> \endverbatim
153: *
154: * Authors:
155: * ========
156: *
1.14 bertrand 157: *> \author Univ. of Tennessee
158: *> \author Univ. of California Berkeley
159: *> \author Univ. of Colorado Denver
160: *> \author NAG Ltd.
1.9 bertrand 161: *
162: *> \ingroup doubleSYcomputational
163: *
1.14 bertrand 164: *> \par Contributors:
165: * ==================
166: *>
167: *> \verbatim
168: *>
169: *> November 2013, Igor Kozachenko,
170: *> Computer Science Division,
171: *> University of California, Berkeley
172: *> \endverbatim
173: *
1.9 bertrand 174: * =====================================================================
1.1 bertrand 175: SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
176: *
1.19 ! bertrand 177: * -- LAPACK computational routine --
1.1 bertrand 178: * -- LAPACK is a software package provided by Univ. of Tennessee, --
179: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180: *
181: * .. Scalar Arguments ..
182: CHARACTER UPLO
183: INTEGER INFO, KB, LDA, LDW, N, NB
184: * ..
185: * .. Array Arguments ..
186: INTEGER IPIV( * )
187: DOUBLE PRECISION A( LDA, * ), W( LDW, * )
188: * ..
189: *
190: * =====================================================================
191: *
192: * .. Parameters ..
193: DOUBLE PRECISION ZERO, ONE
194: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
195: DOUBLE PRECISION EIGHT, SEVTEN
196: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
197: * ..
198: * .. Local Scalars ..
199: INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
200: $ KSTEP, KW
201: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1,
202: $ ROWMAX, T
203: * ..
204: * .. External Functions ..
205: LOGICAL LSAME
206: INTEGER IDAMAX
207: EXTERNAL LSAME, IDAMAX
208: * ..
209: * .. External Subroutines ..
210: EXTERNAL DCOPY, DGEMM, DGEMV, DSCAL, DSWAP
211: * ..
212: * .. Intrinsic Functions ..
213: INTRINSIC ABS, MAX, MIN, SQRT
214: * ..
215: * .. Executable Statements ..
216: *
217: INFO = 0
218: *
219: * Initialize ALPHA for use in choosing pivot block size.
220: *
221: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
222: *
223: IF( LSAME( UPLO, 'U' ) ) THEN
224: *
225: * Factorize the trailing columns of A using the upper triangle
226: * of A and working backwards, and compute the matrix W = U12*D
227: * for use in updating A11
228: *
229: * K is the main loop index, decreasing from N in steps of 1 or 2
230: *
231: * KW is the column of W which corresponds to column K of A
232: *
233: K = N
234: 10 CONTINUE
235: KW = NB + K - N
236: *
237: * Exit from loop
238: *
239: IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
240: $ GO TO 30
241: *
242: * Copy column K of A to column KW of W and update it
243: *
244: CALL DCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 )
245: IF( K.LT.N )
246: $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ), LDA,
247: $ W( K, KW+1 ), LDW, ONE, W( 1, KW ), 1 )
248: *
249: KSTEP = 1
250: *
251: * Determine rows and columns to be interchanged and whether
252: * a 1-by-1 or 2-by-2 pivot block will be used
253: *
254: ABSAKK = ABS( W( K, KW ) )
255: *
256: * IMAX is the row-index of the largest off-diagonal element in
1.14 bertrand 257: * column K, and COLMAX is its absolute value.
258: * Determine both COLMAX and IMAX.
1.1 bertrand 259: *
260: IF( K.GT.1 ) THEN
261: IMAX = IDAMAX( K-1, W( 1, KW ), 1 )
262: COLMAX = ABS( W( IMAX, KW ) )
263: ELSE
264: COLMAX = ZERO
265: END IF
266: *
267: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
268: *
1.14 bertrand 269: * Column K is zero or underflow: set INFO and continue
1.1 bertrand 270: *
271: IF( INFO.EQ.0 )
272: $ INFO = K
273: KP = K
274: ELSE
275: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
276: *
277: * no interchange, use 1-by-1 pivot block
278: *
279: KP = K
280: ELSE
281: *
282: * Copy column IMAX to column KW-1 of W and update it
283: *
284: CALL DCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
285: CALL DCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
286: $ W( IMAX+1, KW-1 ), 1 )
287: IF( K.LT.N )
288: $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ),
289: $ LDA, W( IMAX, KW+1 ), LDW, ONE,
290: $ W( 1, KW-1 ), 1 )
291: *
292: * JMAX is the column-index of the largest off-diagonal
293: * element in row IMAX, and ROWMAX is its absolute value
294: *
295: JMAX = IMAX + IDAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
296: ROWMAX = ABS( W( JMAX, KW-1 ) )
297: IF( IMAX.GT.1 ) THEN
298: JMAX = IDAMAX( IMAX-1, W( 1, KW-1 ), 1 )
299: ROWMAX = MAX( ROWMAX, ABS( W( JMAX, KW-1 ) ) )
300: END IF
301: *
302: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
303: *
304: * no interchange, use 1-by-1 pivot block
305: *
306: KP = K
307: ELSE IF( ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
308: *
309: * interchange rows and columns K and IMAX, use 1-by-1
310: * pivot block
311: *
312: KP = IMAX
313: *
1.14 bertrand 314: * copy column KW-1 of W to column KW of W
1.1 bertrand 315: *
316: CALL DCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
317: ELSE
318: *
319: * interchange rows and columns K-1 and IMAX, use 2-by-2
320: * pivot block
321: *
322: KP = IMAX
323: KSTEP = 2
324: END IF
325: END IF
326: *
1.14 bertrand 327: * ============================================================
328: *
329: * KK is the column of A where pivoting step stopped
330: *
1.1 bertrand 331: KK = K - KSTEP + 1
1.14 bertrand 332: *
333: * KKW is the column of W which corresponds to column KK of A
334: *
1.1 bertrand 335: KKW = NB + KK - N
336: *
1.14 bertrand 337: * Interchange rows and columns KP and KK.
338: * Updated column KP is already stored in column KKW of W.
1.1 bertrand 339: *
340: IF( KP.NE.KK ) THEN
341: *
1.14 bertrand 342: * Copy non-updated column KK to column KP of submatrix A
343: * at step K. No need to copy element into column K
344: * (or K and K-1 for 2-by-2 pivot) of A, since these columns
345: * will be later overwritten.
1.1 bertrand 346: *
1.14 bertrand 347: A( KP, KP ) = A( KK, KK )
348: CALL DCOPY( KK-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
1.1 bertrand 349: $ LDA )
1.14 bertrand 350: IF( KP.GT.1 )
351: $ CALL DCOPY( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
1.1 bertrand 352: *
1.14 bertrand 353: * Interchange rows KK and KP in last K+1 to N columns of A
354: * (columns K (or K and K-1 for 2-by-2 pivot) of A will be
355: * later overwritten). Interchange rows KK and KP
356: * in last KKW to NB columns of W.
1.1 bertrand 357: *
1.14 bertrand 358: IF( K.LT.N )
359: $ CALL DSWAP( N-K, A( KK, K+1 ), LDA, A( KP, K+1 ),
360: $ LDA )
1.1 bertrand 361: CALL DSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
362: $ LDW )
363: END IF
364: *
365: IF( KSTEP.EQ.1 ) THEN
366: *
1.14 bertrand 367: * 1-by-1 pivot block D(k): column kw of W now holds
1.1 bertrand 368: *
1.14 bertrand 369: * W(kw) = U(k)*D(k),
1.1 bertrand 370: *
371: * where U(k) is the k-th column of U
372: *
1.14 bertrand 373: * Store subdiag. elements of column U(k)
374: * and 1-by-1 block D(k) in column k of A.
375: * NOTE: Diagonal element U(k,k) is a UNIT element
376: * and not stored.
377: * A(k,k) := D(k,k) = W(k,kw)
378: * A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k)
1.1 bertrand 379: *
380: CALL DCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
381: R1 = ONE / A( K, K )
382: CALL DSCAL( K-1, R1, A( 1, K ), 1 )
1.14 bertrand 383: *
1.1 bertrand 384: ELSE
385: *
1.14 bertrand 386: * 2-by-2 pivot block D(k): columns kw and kw-1 of W now hold
1.1 bertrand 387: *
1.14 bertrand 388: * ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k)
1.1 bertrand 389: *
390: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
391: * of U
392: *
1.14 bertrand 393: * Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2
394: * block D(k-1:k,k-1:k) in columns k-1 and k of A.
395: * NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT
396: * block and not stored.
397: * A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw)
398: * A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) =
399: * = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) )
400: *
1.1 bertrand 401: IF( K.GT.2 ) THEN
402: *
1.14 bertrand 403: * Compose the columns of the inverse of 2-by-2 pivot
404: * block D in the following way to reduce the number
405: * of FLOPS when we myltiply panel ( W(kw-1) W(kw) ) by
406: * this inverse
407: *
408: * D**(-1) = ( d11 d21 )**(-1) =
409: * ( d21 d22 )
410: *
411: * = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
412: * ( (-d21 ) ( d11 ) )
413: *
414: * = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
415: *
416: * * ( ( d22/d21 ) ( -1 ) ) =
417: * ( ( -1 ) ( d11/d21 ) )
418: *
419: * = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) =
420: * ( ( -1 ) ( D22 ) )
421: *
422: * = 1/d21 * T * ( ( D11 ) ( -1 ) )
423: * ( ( -1 ) ( D22 ) )
424: *
425: * = D21 * ( ( D11 ) ( -1 ) )
426: * ( ( -1 ) ( D22 ) )
1.1 bertrand 427: *
428: D21 = W( K-1, KW )
429: D11 = W( K, KW ) / D21
430: D22 = W( K-1, KW-1 ) / D21
431: T = ONE / ( D11*D22-ONE )
432: D21 = T / D21
1.14 bertrand 433: *
434: * Update elements in columns A(k-1) and A(k) as
435: * dot products of rows of ( W(kw-1) W(kw) ) and columns
436: * of D**(-1)
437: *
1.1 bertrand 438: DO 20 J = 1, K - 2
439: A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
440: A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
441: 20 CONTINUE
442: END IF
443: *
444: * Copy D(k) to A
445: *
446: A( K-1, K-1 ) = W( K-1, KW-1 )
447: A( K-1, K ) = W( K-1, KW )
448: A( K, K ) = W( K, KW )
1.14 bertrand 449: *
1.1 bertrand 450: END IF
1.14 bertrand 451: *
1.1 bertrand 452: END IF
453: *
454: * Store details of the interchanges in IPIV
455: *
456: IF( KSTEP.EQ.1 ) THEN
457: IPIV( K ) = KP
458: ELSE
459: IPIV( K ) = -KP
460: IPIV( K-1 ) = -KP
461: END IF
462: *
463: * Decrease K and return to the start of the main loop
464: *
465: K = K - KSTEP
466: GO TO 10
467: *
468: 30 CONTINUE
469: *
470: * Update the upper triangle of A11 (= A(1:k,1:k)) as
471: *
1.8 bertrand 472: * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
1.1 bertrand 473: *
474: * computing blocks of NB columns at a time
475: *
476: DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
477: JB = MIN( NB, K-J+1 )
478: *
479: * Update the upper triangle of the diagonal block
480: *
481: DO 40 JJ = J, J + JB - 1
482: CALL DGEMV( 'No transpose', JJ-J+1, N-K, -ONE,
483: $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
484: $ A( J, JJ ), 1 )
485: 40 CONTINUE
486: *
487: * Update the rectangular superdiagonal block
488: *
489: CALL DGEMM( 'No transpose', 'Transpose', J-1, JB, N-K, -ONE,
490: $ A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, ONE,
491: $ A( 1, J ), LDA )
492: 50 CONTINUE
493: *
494: * Put U12 in standard form by partially undoing the interchanges
1.14 bertrand 495: * in columns k+1:n looping backwards from k+1 to n
1.1 bertrand 496: *
497: J = K + 1
498: 60 CONTINUE
1.14 bertrand 499: *
500: * Undo the interchanges (if any) of rows JJ and JP at each
501: * step J
502: *
503: * (Here, J is a diagonal index)
504: JJ = J
505: JP = IPIV( J )
506: IF( JP.LT.0 ) THEN
507: JP = -JP
508: * (Here, J is a diagonal index)
509: J = J + 1
510: END IF
511: * (NOTE: Here, J is used to determine row length. Length N-J+1
512: * of the rows to swap back doesn't include diagonal element)
1.1 bertrand 513: J = J + 1
1.14 bertrand 514: IF( JP.NE.JJ .AND. J.LE.N )
515: $ CALL DSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
516: IF( J.LT.N )
1.1 bertrand 517: $ GO TO 60
518: *
519: * Set KB to the number of columns factorized
520: *
521: KB = N - K
522: *
523: ELSE
524: *
525: * Factorize the leading columns of A using the lower triangle
526: * of A and working forwards, and compute the matrix W = L21*D
527: * for use in updating A22
528: *
529: * K is the main loop index, increasing from 1 in steps of 1 or 2
530: *
531: K = 1
532: 70 CONTINUE
533: *
534: * Exit from loop
535: *
536: IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
537: $ GO TO 90
538: *
539: * Copy column K of A to column K of W and update it
540: *
541: CALL DCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
542: CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ), LDA,
543: $ W( K, 1 ), LDW, ONE, W( K, K ), 1 )
544: *
545: KSTEP = 1
546: *
547: * Determine rows and columns to be interchanged and whether
548: * a 1-by-1 or 2-by-2 pivot block will be used
549: *
550: ABSAKK = ABS( W( K, K ) )
551: *
552: * IMAX is the row-index of the largest off-diagonal element in
1.14 bertrand 553: * column K, and COLMAX is its absolute value.
554: * Determine both COLMAX and IMAX.
1.1 bertrand 555: *
556: IF( K.LT.N ) THEN
557: IMAX = K + IDAMAX( N-K, W( K+1, K ), 1 )
558: COLMAX = ABS( W( IMAX, K ) )
559: ELSE
560: COLMAX = ZERO
561: END IF
562: *
563: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
564: *
1.14 bertrand 565: * Column K is zero or underflow: set INFO and continue
1.1 bertrand 566: *
567: IF( INFO.EQ.0 )
568: $ INFO = K
569: KP = K
570: ELSE
571: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
572: *
573: * no interchange, use 1-by-1 pivot block
574: *
575: KP = K
576: ELSE
577: *
578: * Copy column IMAX to column K+1 of W and update it
579: *
580: CALL DCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
581: CALL DCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
582: $ 1 )
583: CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ),
584: $ LDA, W( IMAX, 1 ), LDW, ONE, W( K, K+1 ), 1 )
585: *
586: * JMAX is the column-index of the largest off-diagonal
587: * element in row IMAX, and ROWMAX is its absolute value
588: *
589: JMAX = K - 1 + IDAMAX( IMAX-K, W( K, K+1 ), 1 )
590: ROWMAX = ABS( W( JMAX, K+1 ) )
591: IF( IMAX.LT.N ) THEN
592: JMAX = IMAX + IDAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
593: ROWMAX = MAX( ROWMAX, ABS( W( JMAX, K+1 ) ) )
594: END IF
595: *
596: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
597: *
598: * no interchange, use 1-by-1 pivot block
599: *
600: KP = K
601: ELSE IF( ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
602: *
603: * interchange rows and columns K and IMAX, use 1-by-1
604: * pivot block
605: *
606: KP = IMAX
607: *
1.14 bertrand 608: * copy column K+1 of W to column K of W
1.1 bertrand 609: *
610: CALL DCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
611: ELSE
612: *
613: * interchange rows and columns K+1 and IMAX, use 2-by-2
614: * pivot block
615: *
616: KP = IMAX
617: KSTEP = 2
618: END IF
619: END IF
620: *
1.14 bertrand 621: * ============================================================
622: *
623: * KK is the column of A where pivoting step stopped
624: *
1.1 bertrand 625: KK = K + KSTEP - 1
626: *
1.14 bertrand 627: * Interchange rows and columns KP and KK.
628: * Updated column KP is already stored in column KK of W.
1.1 bertrand 629: *
630: IF( KP.NE.KK ) THEN
631: *
1.14 bertrand 632: * Copy non-updated column KK to column KP of submatrix A
633: * at step K. No need to copy element into column K
634: * (or K and K+1 for 2-by-2 pivot) of A, since these columns
635: * will be later overwritten.
1.1 bertrand 636: *
1.14 bertrand 637: A( KP, KP ) = A( KK, KK )
638: CALL DCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
639: $ LDA )
640: IF( KP.LT.N )
641: $ CALL DCOPY( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
1.1 bertrand 642: *
1.14 bertrand 643: * Interchange rows KK and KP in first K-1 columns of A
644: * (columns K (or K and K+1 for 2-by-2 pivot) of A will be
645: * later overwritten). Interchange rows KK and KP
646: * in first KK columns of W.
1.1 bertrand 647: *
1.14 bertrand 648: IF( K.GT.1 )
649: $ CALL DSWAP( K-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
1.1 bertrand 650: CALL DSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
651: END IF
652: *
653: IF( KSTEP.EQ.1 ) THEN
654: *
655: * 1-by-1 pivot block D(k): column k of W now holds
656: *
1.14 bertrand 657: * W(k) = L(k)*D(k),
1.1 bertrand 658: *
659: * where L(k) is the k-th column of L
660: *
1.14 bertrand 661: * Store subdiag. elements of column L(k)
662: * and 1-by-1 block D(k) in column k of A.
663: * (NOTE: Diagonal element L(k,k) is a UNIT element
664: * and not stored)
665: * A(k,k) := D(k,k) = W(k,k)
666: * A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k)
1.1 bertrand 667: *
668: CALL DCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
669: IF( K.LT.N ) THEN
670: R1 = ONE / A( K, K )
671: CALL DSCAL( N-K, R1, A( K+1, K ), 1 )
672: END IF
1.14 bertrand 673: *
1.1 bertrand 674: ELSE
675: *
676: * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
677: *
678: * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
679: *
680: * where L(k) and L(k+1) are the k-th and (k+1)-th columns
681: * of L
682: *
1.14 bertrand 683: * Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2
684: * block D(k:k+1,k:k+1) in columns k and k+1 of A.
685: * (NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT
686: * block and not stored)
687: * A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1)
688: * A(k+2:N,k:k+1) := L(k+2:N,k:k+1) =
689: * = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) )
690: *
1.1 bertrand 691: IF( K.LT.N-1 ) THEN
692: *
1.14 bertrand 693: * Compose the columns of the inverse of 2-by-2 pivot
694: * block D in the following way to reduce the number
695: * of FLOPS when we myltiply panel ( W(k) W(k+1) ) by
696: * this inverse
697: *
698: * D**(-1) = ( d11 d21 )**(-1) =
699: * ( d21 d22 )
700: *
701: * = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
702: * ( (-d21 ) ( d11 ) )
703: *
704: * = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
705: *
706: * * ( ( d22/d21 ) ( -1 ) ) =
707: * ( ( -1 ) ( d11/d21 ) )
708: *
709: * = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) =
710: * ( ( -1 ) ( D22 ) )
711: *
712: * = 1/d21 * T * ( ( D11 ) ( -1 ) )
713: * ( ( -1 ) ( D22 ) )
714: *
715: * = D21 * ( ( D11 ) ( -1 ) )
716: * ( ( -1 ) ( D22 ) )
1.1 bertrand 717: *
718: D21 = W( K+1, K )
719: D11 = W( K+1, K+1 ) / D21
720: D22 = W( K, K ) / D21
721: T = ONE / ( D11*D22-ONE )
722: D21 = T / D21
1.14 bertrand 723: *
724: * Update elements in columns A(k) and A(k+1) as
725: * dot products of rows of ( W(k) W(k+1) ) and columns
726: * of D**(-1)
727: *
1.1 bertrand 728: DO 80 J = K + 2, N
729: A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
730: A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
731: 80 CONTINUE
732: END IF
733: *
734: * Copy D(k) to A
735: *
736: A( K, K ) = W( K, K )
737: A( K+1, K ) = W( K+1, K )
738: A( K+1, K+1 ) = W( K+1, K+1 )
1.14 bertrand 739: *
1.1 bertrand 740: END IF
1.14 bertrand 741: *
1.1 bertrand 742: END IF
743: *
744: * Store details of the interchanges in IPIV
745: *
746: IF( KSTEP.EQ.1 ) THEN
747: IPIV( K ) = KP
748: ELSE
749: IPIV( K ) = -KP
750: IPIV( K+1 ) = -KP
751: END IF
752: *
753: * Increase K and return to the start of the main loop
754: *
755: K = K + KSTEP
756: GO TO 70
757: *
758: 90 CONTINUE
759: *
760: * Update the lower triangle of A22 (= A(k:n,k:n)) as
761: *
1.8 bertrand 762: * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
1.1 bertrand 763: *
764: * computing blocks of NB columns at a time
765: *
766: DO 110 J = K, N, NB
767: JB = MIN( NB, N-J+1 )
768: *
769: * Update the lower triangle of the diagonal block
770: *
771: DO 100 JJ = J, J + JB - 1
772: CALL DGEMV( 'No transpose', J+JB-JJ, K-1, -ONE,
773: $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
774: $ A( JJ, JJ ), 1 )
775: 100 CONTINUE
776: *
777: * Update the rectangular subdiagonal block
778: *
779: IF( J+JB.LE.N )
780: $ CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
781: $ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
782: $ ONE, A( J+JB, J ), LDA )
783: 110 CONTINUE
784: *
785: * Put L21 in standard form by partially undoing the interchanges
1.14 bertrand 786: * of rows in columns 1:k-1 looping backwards from k-1 to 1
1.1 bertrand 787: *
788: J = K - 1
789: 120 CONTINUE
1.14 bertrand 790: *
791: * Undo the interchanges (if any) of rows JJ and JP at each
792: * step J
793: *
794: * (Here, J is a diagonal index)
795: JJ = J
796: JP = IPIV( J )
797: IF( JP.LT.0 ) THEN
798: JP = -JP
799: * (Here, J is a diagonal index)
800: J = J - 1
801: END IF
802: * (NOTE: Here, J is used to determine row length. Length J
803: * of the rows to swap back doesn't include diagonal element)
1.1 bertrand 804: J = J - 1
1.14 bertrand 805: IF( JP.NE.JJ .AND. J.GE.1 )
806: $ CALL DSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
807: IF( J.GT.1 )
1.1 bertrand 808: $ GO TO 120
809: *
810: * Set KB to the number of columns factorized
811: *
812: KB = K - 1
813: *
814: END IF
815: RETURN
816: *
817: * End of DLASYF
818: *
819: 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