Annotation of rpl/lapack/lapack/dlasyf.f, revision 1.11
1.9 bertrand 1: *> \brief \b DLASYF
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
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">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
22: *
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: * ..
31: *
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: *> If UPLO = 'U', only the last KB elements of IPIV are set;
113: *> if UPLO = 'L', only the first KB elements are set.
114: *>
115: *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
116: *> interchanged and D(k,k) is a 1-by-1 diagonal block.
117: *> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
118: *> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
119: *> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
120: *> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
121: *> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
122: *> \endverbatim
123: *>
124: *> \param[out] W
125: *> \verbatim
126: *> W is DOUBLE PRECISION array, dimension (LDW,NB)
127: *> \endverbatim
128: *>
129: *> \param[in] LDW
130: *> \verbatim
131: *> LDW is INTEGER
132: *> The leading dimension of the array W. LDW >= max(1,N).
133: *> \endverbatim
134: *>
135: *> \param[out] INFO
136: *> \verbatim
137: *> INFO is INTEGER
138: *> = 0: successful exit
139: *> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
140: *> has been completed, but the block diagonal matrix D is
141: *> exactly singular.
142: *> \endverbatim
143: *
144: * Authors:
145: * ========
146: *
147: *> \author Univ. of Tennessee
148: *> \author Univ. of California Berkeley
149: *> \author Univ. of Colorado Denver
150: *> \author NAG Ltd.
151: *
152: *> \date November 2011
153: *
154: *> \ingroup doubleSYcomputational
155: *
156: * =====================================================================
1.1 bertrand 157: SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
158: *
1.9 bertrand 159: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 160: * -- LAPACK is a software package provided by Univ. of Tennessee, --
161: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 bertrand 162: * November 2011
1.1 bertrand 163: *
164: * .. Scalar Arguments ..
165: CHARACTER UPLO
166: INTEGER INFO, KB, LDA, LDW, N, NB
167: * ..
168: * .. Array Arguments ..
169: INTEGER IPIV( * )
170: DOUBLE PRECISION A( LDA, * ), W( LDW, * )
171: * ..
172: *
173: * =====================================================================
174: *
175: * .. Parameters ..
176: DOUBLE PRECISION ZERO, ONE
177: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
178: DOUBLE PRECISION EIGHT, SEVTEN
179: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
180: * ..
181: * .. Local Scalars ..
182: INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
183: $ KSTEP, KW
184: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1,
185: $ ROWMAX, T
186: * ..
187: * .. External Functions ..
188: LOGICAL LSAME
189: INTEGER IDAMAX
190: EXTERNAL LSAME, IDAMAX
191: * ..
192: * .. External Subroutines ..
193: EXTERNAL DCOPY, DGEMM, DGEMV, DSCAL, DSWAP
194: * ..
195: * .. Intrinsic Functions ..
196: INTRINSIC ABS, MAX, MIN, SQRT
197: * ..
198: * .. Executable Statements ..
199: *
200: INFO = 0
201: *
202: * Initialize ALPHA for use in choosing pivot block size.
203: *
204: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
205: *
206: IF( LSAME( UPLO, 'U' ) ) THEN
207: *
208: * Factorize the trailing columns of A using the upper triangle
209: * of A and working backwards, and compute the matrix W = U12*D
210: * for use in updating A11
211: *
212: * K is the main loop index, decreasing from N in steps of 1 or 2
213: *
214: * KW is the column of W which corresponds to column K of A
215: *
216: K = N
217: 10 CONTINUE
218: KW = NB + K - N
219: *
220: * Exit from loop
221: *
222: IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
223: $ GO TO 30
224: *
225: * Copy column K of A to column KW of W and update it
226: *
227: CALL DCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 )
228: IF( K.LT.N )
229: $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ), LDA,
230: $ W( K, KW+1 ), LDW, ONE, W( 1, KW ), 1 )
231: *
232: KSTEP = 1
233: *
234: * Determine rows and columns to be interchanged and whether
235: * a 1-by-1 or 2-by-2 pivot block will be used
236: *
237: ABSAKK = ABS( W( K, KW ) )
238: *
239: * IMAX is the row-index of the largest off-diagonal element in
240: * column K, and COLMAX is its absolute value
241: *
242: IF( K.GT.1 ) THEN
243: IMAX = IDAMAX( K-1, W( 1, KW ), 1 )
244: COLMAX = ABS( W( IMAX, KW ) )
245: ELSE
246: COLMAX = ZERO
247: END IF
248: *
249: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
250: *
251: * Column K is zero: set INFO and continue
252: *
253: IF( INFO.EQ.0 )
254: $ INFO = K
255: KP = K
256: ELSE
257: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
258: *
259: * no interchange, use 1-by-1 pivot block
260: *
261: KP = K
262: ELSE
263: *
264: * Copy column IMAX to column KW-1 of W and update it
265: *
266: CALL DCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
267: CALL DCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
268: $ W( IMAX+1, KW-1 ), 1 )
269: IF( K.LT.N )
270: $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ),
271: $ LDA, W( IMAX, KW+1 ), LDW, ONE,
272: $ W( 1, KW-1 ), 1 )
273: *
274: * JMAX is the column-index of the largest off-diagonal
275: * element in row IMAX, and ROWMAX is its absolute value
276: *
277: JMAX = IMAX + IDAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
278: ROWMAX = ABS( W( JMAX, KW-1 ) )
279: IF( IMAX.GT.1 ) THEN
280: JMAX = IDAMAX( IMAX-1, W( 1, KW-1 ), 1 )
281: ROWMAX = MAX( ROWMAX, ABS( W( JMAX, KW-1 ) ) )
282: END IF
283: *
284: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
285: *
286: * no interchange, use 1-by-1 pivot block
287: *
288: KP = K
289: ELSE IF( ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
290: *
291: * interchange rows and columns K and IMAX, use 1-by-1
292: * pivot block
293: *
294: KP = IMAX
295: *
296: * copy column KW-1 of W to column KW
297: *
298: CALL DCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
299: ELSE
300: *
301: * interchange rows and columns K-1 and IMAX, use 2-by-2
302: * pivot block
303: *
304: KP = IMAX
305: KSTEP = 2
306: END IF
307: END IF
308: *
309: KK = K - KSTEP + 1
310: KKW = NB + KK - N
311: *
312: * Updated column KP is already stored in column KKW of W
313: *
314: IF( KP.NE.KK ) THEN
315: *
316: * Copy non-updated column KK to column KP
317: *
318: A( KP, K ) = A( KK, K )
319: CALL DCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
320: $ LDA )
321: CALL DCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
322: *
323: * Interchange rows KK and KP in last KK columns of A and W
324: *
325: CALL DSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
326: CALL DSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
327: $ LDW )
328: END IF
329: *
330: IF( KSTEP.EQ.1 ) THEN
331: *
332: * 1-by-1 pivot block D(k): column KW of W now holds
333: *
334: * W(k) = U(k)*D(k)
335: *
336: * where U(k) is the k-th column of U
337: *
338: * Store U(k) in column k of A
339: *
340: CALL DCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
341: R1 = ONE / A( K, K )
342: CALL DSCAL( K-1, R1, A( 1, K ), 1 )
343: ELSE
344: *
345: * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
346: * hold
347: *
348: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
349: *
350: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
351: * of U
352: *
353: IF( K.GT.2 ) THEN
354: *
355: * Store U(k) and U(k-1) in columns k and k-1 of A
356: *
357: D21 = W( K-1, KW )
358: D11 = W( K, KW ) / D21
359: D22 = W( K-1, KW-1 ) / D21
360: T = ONE / ( D11*D22-ONE )
361: D21 = T / D21
362: DO 20 J = 1, K - 2
363: A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
364: A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
365: 20 CONTINUE
366: END IF
367: *
368: * Copy D(k) to A
369: *
370: A( K-1, K-1 ) = W( K-1, KW-1 )
371: A( K-1, K ) = W( K-1, KW )
372: A( K, K ) = W( K, KW )
373: END IF
374: END IF
375: *
376: * Store details of the interchanges in IPIV
377: *
378: IF( KSTEP.EQ.1 ) THEN
379: IPIV( K ) = KP
380: ELSE
381: IPIV( K ) = -KP
382: IPIV( K-1 ) = -KP
383: END IF
384: *
385: * Decrease K and return to the start of the main loop
386: *
387: K = K - KSTEP
388: GO TO 10
389: *
390: 30 CONTINUE
391: *
392: * Update the upper triangle of A11 (= A(1:k,1:k)) as
393: *
1.8 bertrand 394: * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
1.1 bertrand 395: *
396: * computing blocks of NB columns at a time
397: *
398: DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
399: JB = MIN( NB, K-J+1 )
400: *
401: * Update the upper triangle of the diagonal block
402: *
403: DO 40 JJ = J, J + JB - 1
404: CALL DGEMV( 'No transpose', JJ-J+1, N-K, -ONE,
405: $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
406: $ A( J, JJ ), 1 )
407: 40 CONTINUE
408: *
409: * Update the rectangular superdiagonal block
410: *
411: CALL DGEMM( 'No transpose', 'Transpose', J-1, JB, N-K, -ONE,
412: $ A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, ONE,
413: $ A( 1, J ), LDA )
414: 50 CONTINUE
415: *
416: * Put U12 in standard form by partially undoing the interchanges
417: * in columns k+1:n
418: *
419: J = K + 1
420: 60 CONTINUE
421: JJ = J
422: JP = IPIV( J )
423: IF( JP.LT.0 ) THEN
424: JP = -JP
425: J = J + 1
426: END IF
427: J = J + 1
428: IF( JP.NE.JJ .AND. J.LE.N )
429: $ CALL DSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
430: IF( J.LE.N )
431: $ GO TO 60
432: *
433: * Set KB to the number of columns factorized
434: *
435: KB = N - K
436: *
437: ELSE
438: *
439: * Factorize the leading columns of A using the lower triangle
440: * of A and working forwards, and compute the matrix W = L21*D
441: * for use in updating A22
442: *
443: * K is the main loop index, increasing from 1 in steps of 1 or 2
444: *
445: K = 1
446: 70 CONTINUE
447: *
448: * Exit from loop
449: *
450: IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
451: $ GO TO 90
452: *
453: * Copy column K of A to column K of W and update it
454: *
455: CALL DCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
456: CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ), LDA,
457: $ W( K, 1 ), LDW, ONE, W( K, K ), 1 )
458: *
459: KSTEP = 1
460: *
461: * Determine rows and columns to be interchanged and whether
462: * a 1-by-1 or 2-by-2 pivot block will be used
463: *
464: ABSAKK = ABS( W( K, K ) )
465: *
466: * IMAX is the row-index of the largest off-diagonal element in
467: * column K, and COLMAX is its absolute value
468: *
469: IF( K.LT.N ) THEN
470: IMAX = K + IDAMAX( N-K, W( K+1, K ), 1 )
471: COLMAX = ABS( W( IMAX, K ) )
472: ELSE
473: COLMAX = ZERO
474: END IF
475: *
476: IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
477: *
478: * Column K is zero: set INFO and continue
479: *
480: IF( INFO.EQ.0 )
481: $ INFO = K
482: KP = K
483: ELSE
484: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
485: *
486: * no interchange, use 1-by-1 pivot block
487: *
488: KP = K
489: ELSE
490: *
491: * Copy column IMAX to column K+1 of W and update it
492: *
493: CALL DCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
494: CALL DCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
495: $ 1 )
496: CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ),
497: $ LDA, W( IMAX, 1 ), LDW, ONE, W( K, K+1 ), 1 )
498: *
499: * JMAX is the column-index of the largest off-diagonal
500: * element in row IMAX, and ROWMAX is its absolute value
501: *
502: JMAX = K - 1 + IDAMAX( IMAX-K, W( K, K+1 ), 1 )
503: ROWMAX = ABS( W( JMAX, K+1 ) )
504: IF( IMAX.LT.N ) THEN
505: JMAX = IMAX + IDAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
506: ROWMAX = MAX( ROWMAX, ABS( W( JMAX, K+1 ) ) )
507: END IF
508: *
509: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
510: *
511: * no interchange, use 1-by-1 pivot block
512: *
513: KP = K
514: ELSE IF( ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
515: *
516: * interchange rows and columns K and IMAX, use 1-by-1
517: * pivot block
518: *
519: KP = IMAX
520: *
521: * copy column K+1 of W to column K
522: *
523: CALL DCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
524: ELSE
525: *
526: * interchange rows and columns K+1 and IMAX, use 2-by-2
527: * pivot block
528: *
529: KP = IMAX
530: KSTEP = 2
531: END IF
532: END IF
533: *
534: KK = K + KSTEP - 1
535: *
536: * Updated column KP is already stored in column KK of W
537: *
538: IF( KP.NE.KK ) THEN
539: *
540: * Copy non-updated column KK to column KP
541: *
542: A( KP, K ) = A( KK, K )
543: CALL DCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
544: CALL DCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
545: *
546: * Interchange rows KK and KP in first KK columns of A and W
547: *
548: CALL DSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
549: CALL DSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
550: END IF
551: *
552: IF( KSTEP.EQ.1 ) THEN
553: *
554: * 1-by-1 pivot block D(k): column k of W now holds
555: *
556: * W(k) = L(k)*D(k)
557: *
558: * where L(k) is the k-th column of L
559: *
560: * Store L(k) in column k of A
561: *
562: CALL DCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
563: IF( K.LT.N ) THEN
564: R1 = ONE / A( K, K )
565: CALL DSCAL( N-K, R1, A( K+1, K ), 1 )
566: END IF
567: ELSE
568: *
569: * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
570: *
571: * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
572: *
573: * where L(k) and L(k+1) are the k-th and (k+1)-th columns
574: * of L
575: *
576: IF( K.LT.N-1 ) THEN
577: *
578: * Store L(k) and L(k+1) in columns k and k+1 of A
579: *
580: D21 = W( K+1, K )
581: D11 = W( K+1, K+1 ) / D21
582: D22 = W( K, K ) / D21
583: T = ONE / ( D11*D22-ONE )
584: D21 = T / D21
585: DO 80 J = K + 2, N
586: A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
587: A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
588: 80 CONTINUE
589: END IF
590: *
591: * Copy D(k) to A
592: *
593: A( K, K ) = W( K, K )
594: A( K+1, K ) = W( K+1, K )
595: A( K+1, K+1 ) = W( K+1, K+1 )
596: END IF
597: END IF
598: *
599: * Store details of the interchanges in IPIV
600: *
601: IF( KSTEP.EQ.1 ) THEN
602: IPIV( K ) = KP
603: ELSE
604: IPIV( K ) = -KP
605: IPIV( K+1 ) = -KP
606: END IF
607: *
608: * Increase K and return to the start of the main loop
609: *
610: K = K + KSTEP
611: GO TO 70
612: *
613: 90 CONTINUE
614: *
615: * Update the lower triangle of A22 (= A(k:n,k:n)) as
616: *
1.8 bertrand 617: * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
1.1 bertrand 618: *
619: * computing blocks of NB columns at a time
620: *
621: DO 110 J = K, N, NB
622: JB = MIN( NB, N-J+1 )
623: *
624: * Update the lower triangle of the diagonal block
625: *
626: DO 100 JJ = J, J + JB - 1
627: CALL DGEMV( 'No transpose', J+JB-JJ, K-1, -ONE,
628: $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
629: $ A( JJ, JJ ), 1 )
630: 100 CONTINUE
631: *
632: * Update the rectangular subdiagonal block
633: *
634: IF( J+JB.LE.N )
635: $ CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
636: $ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
637: $ ONE, A( J+JB, J ), LDA )
638: 110 CONTINUE
639: *
640: * Put L21 in standard form by partially undoing the interchanges
641: * in columns 1:k-1
642: *
643: J = K - 1
644: 120 CONTINUE
645: JJ = J
646: JP = IPIV( J )
647: IF( JP.LT.0 ) THEN
648: JP = -JP
649: J = J - 1
650: END IF
651: J = J - 1
652: IF( JP.NE.JJ .AND. J.GE.1 )
653: $ CALL DSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
654: IF( J.GE.1 )
655: $ GO TO 120
656: *
657: * Set KB to the number of columns factorized
658: *
659: KB = K - 1
660: *
661: END IF
662: RETURN
663: *
664: * End of DLASYF
665: *
666: 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