Annotation of rpl/lapack/lapack/dlasyf_aa.f, revision 1.5
1.1 bertrand 1: *> \brief \b DLASYF_AA
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLASYF_AA + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasyf_aa.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasyf_aa.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasyf_aa.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
1.3 bertrand 22: * H, LDH, WORK )
1.1 bertrand 23: *
24: * .. Scalar Arguments ..
25: * CHARACTER UPLO
1.3 bertrand 26: * INTEGER J1, M, NB, LDA, LDH
1.1 bertrand 27: * ..
28: * .. Array Arguments ..
29: * INTEGER IPIV( * )
30: * DOUBLE PRECISION A( LDA, * ), H( LDH, * ), WORK( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DLATRF_AA factorizes a panel of a real symmetric matrix A using
40: *> the Aasen's algorithm. The panel consists of a set of NB rows of A
41: *> when UPLO is U, or a set of NB columns when UPLO is L.
42: *>
43: *> In order to factorize the panel, the Aasen's algorithm requires the
44: *> last row, or column, of the previous panel. The first row, or column,
45: *> of A is set to be the first row, or column, of an identity matrix,
46: *> which is used to factorize the first panel.
47: *>
48: *> The resulting J-th row of U, or J-th column of L, is stored in the
49: *> (J-1)-th row, or column, of A (without the unit diagonals), while
50: *> the diagonal and subdiagonal of A are overwritten by those of T.
51: *>
52: *> \endverbatim
53: *
54: * Arguments:
55: * ==========
56: *
57: *> \param[in] UPLO
58: *> \verbatim
59: *> UPLO is CHARACTER*1
60: *> = 'U': Upper triangle of A is stored;
61: *> = 'L': Lower triangle of A is stored.
62: *> \endverbatim
63: *>
64: *> \param[in] J1
65: *> \verbatim
66: *> J1 is INTEGER
67: *> The location of the first row, or column, of the panel
68: *> within the submatrix of A, passed to this routine, e.g.,
69: *> when called by DSYTRF_AA, for the first panel, J1 is 1,
70: *> while for the remaining panels, J1 is 2.
71: *> \endverbatim
72: *>
73: *> \param[in] M
74: *> \verbatim
75: *> M is INTEGER
76: *> The dimension of the submatrix. M >= 0.
77: *> \endverbatim
78: *>
79: *> \param[in] NB
80: *> \verbatim
81: *> NB is INTEGER
82: *> The dimension of the panel to be facotorized.
83: *> \endverbatim
84: *>
85: *> \param[in,out] A
86: *> \verbatim
87: *> A is DOUBLE PRECISION array, dimension (LDA,M) for
88: *> the first panel, while dimension (LDA,M+1) for the
89: *> remaining panels.
90: *>
91: *> On entry, A contains the last row, or column, of
92: *> the previous panel, and the trailing submatrix of A
93: *> to be factorized, except for the first panel, only
94: *> the panel is passed.
95: *>
96: *> On exit, the leading panel is factorized.
97: *> \endverbatim
98: *>
99: *> \param[in] LDA
100: *> \verbatim
101: *> LDA is INTEGER
1.3 bertrand 102: *> The leading dimension of the array A. LDA >= max(1,M).
1.1 bertrand 103: *> \endverbatim
104: *>
105: *> \param[out] IPIV
106: *> \verbatim
1.3 bertrand 107: *> IPIV is INTEGER array, dimension (M)
1.1 bertrand 108: *> Details of the row and column interchanges,
109: *> the row and column k were interchanged with the row and
110: *> column IPIV(k).
111: *> \endverbatim
112: *>
113: *> \param[in,out] H
114: *> \verbatim
115: *> H is DOUBLE PRECISION workspace, dimension (LDH,NB).
116: *>
117: *> \endverbatim
118: *>
119: *> \param[in] LDH
120: *> \verbatim
121: *> LDH is INTEGER
122: *> The leading dimension of the workspace H. LDH >= max(1,M).
123: *> \endverbatim
124: *>
125: *> \param[out] WORK
126: *> \verbatim
127: *> WORK is DOUBLE PRECISION workspace, dimension (M).
128: *> \endverbatim
129: *>
130: *
131: * Authors:
132: * ========
133: *
134: *> \author Univ. of Tennessee
135: *> \author Univ. of California Berkeley
136: *> \author Univ. of Colorado Denver
137: *> \author NAG Ltd.
138: *
1.3 bertrand 139: *> \date November 2017
1.1 bertrand 140: *
141: *> \ingroup doubleSYcomputational
142: *
143: * =====================================================================
144: SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
1.3 bertrand 145: $ H, LDH, WORK )
1.1 bertrand 146: *
1.3 bertrand 147: * -- LAPACK computational routine (version 3.8.0) --
1.1 bertrand 148: * -- LAPACK is a software package provided by Univ. of Tennessee, --
149: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.3 bertrand 150: * November 2017
1.1 bertrand 151: *
152: IMPLICIT NONE
153: *
154: * .. Scalar Arguments ..
155: CHARACTER UPLO
1.3 bertrand 156: INTEGER M, NB, J1, LDA, LDH
1.1 bertrand 157: * ..
158: * .. Array Arguments ..
159: INTEGER IPIV( * )
160: DOUBLE PRECISION A( LDA, * ), H( LDH, * ), WORK( * )
161: * ..
162: *
163: * =====================================================================
164: * .. Parameters ..
165: DOUBLE PRECISION ZERO, ONE
166: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
167: *
168: * .. Local Scalars ..
1.3 bertrand 169: INTEGER J, K, K1, I1, I2, MJ
1.1 bertrand 170: DOUBLE PRECISION PIV, ALPHA
171: * ..
172: * .. External Functions ..
173: LOGICAL LSAME
174: INTEGER IDAMAX, ILAENV
175: EXTERNAL LSAME, ILAENV, IDAMAX
176: * ..
177: * .. External Subroutines ..
1.3 bertrand 178: EXTERNAL DGEMV, DAXPY, DCOPY, DSWAP, DSCAL, DLASET,
179: $ XERBLA
1.1 bertrand 180: * ..
181: * .. Intrinsic Functions ..
182: INTRINSIC MAX
183: * ..
184: * .. Executable Statements ..
185: *
186: J = 1
187: *
188: * K1 is the first column of the panel to be factorized
189: * i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks
190: *
191: K1 = (2-J1)+1
192: *
193: IF( LSAME( UPLO, 'U' ) ) THEN
194: *
195: * .....................................................
196: * Factorize A as U**T*D*U using the upper triangle of A
197: * .....................................................
198: *
199: 10 CONTINUE
200: IF ( J.GT.MIN(M, NB) )
201: $ GO TO 20
202: *
203: * K is the column to be factorized
204: * when being called from DSYTRF_AA,
205: * > for the first block column, J1 is 1, hence J1+J-1 is J,
206: * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
207: *
208: K = J1+J-1
1.3 bertrand 209: IF( J.EQ.M ) THEN
210: *
211: * Only need to compute T(J, J)
212: *
213: MJ = 1
214: ELSE
215: MJ = M-J+1
216: END IF
1.1 bertrand 217: *
1.3 bertrand 218: * H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J),
219: * where H(J:M, J) has been initialized to be A(J, J:M)
1.1 bertrand 220: *
221: IF( K.GT.2 ) THEN
222: *
223: * K is the column to be factorized
224: * > for the first block column, K is J, skipping the first two
225: * columns
226: * > for the rest of the columns, K is J+1, skipping only the
227: * first column
228: *
1.3 bertrand 229: CALL DGEMV( 'No transpose', MJ, J-K1,
1.1 bertrand 230: $ -ONE, H( J, K1 ), LDH,
231: $ A( 1, J ), 1,
232: $ ONE, H( J, J ), 1 )
233: END IF
234: *
1.3 bertrand 235: * Copy H(i:M, i) into WORK
1.1 bertrand 236: *
1.3 bertrand 237: CALL DCOPY( MJ, H( J, J ), 1, WORK( 1 ), 1 )
1.1 bertrand 238: *
239: IF( J.GT.K1 ) THEN
240: *
1.3 bertrand 241: * Compute WORK := WORK - L(J-1, J:M) * T(J-1,J),
242: * where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M)
1.1 bertrand 243: *
244: ALPHA = -A( K-1, J )
1.3 bertrand 245: CALL DAXPY( MJ, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )
1.1 bertrand 246: END IF
247: *
248: * Set A(J, J) = T(J, J)
249: *
250: A( K, J ) = WORK( 1 )
251: *
252: IF( J.LT.M ) THEN
253: *
1.3 bertrand 254: * Compute WORK(2:M) = T(J, J) L(J, (J+1):M)
255: * where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M)
1.1 bertrand 256: *
257: IF( K.GT.1 ) THEN
258: ALPHA = -A( K, J )
259: CALL DAXPY( M-J, ALPHA, A( K-1, J+1 ), LDA,
260: $ WORK( 2 ), 1 )
261: ENDIF
262: *
1.3 bertrand 263: * Find max(|WORK(2:M)|)
1.1 bertrand 264: *
265: I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1
266: PIV = WORK( I2 )
267: *
268: * Apply symmetric pivot
269: *
270: IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN
271: *
272: * Swap WORK(I1) and WORK(I2)
273: *
274: I1 = 2
275: WORK( I2 ) = WORK( I1 )
276: WORK( I1 ) = PIV
277: *
1.3 bertrand 278: * Swap A(I1, I1+1:M) with A(I1+1:M, I2)
1.1 bertrand 279: *
280: I1 = I1+J-1
281: I2 = I2+J-1
282: CALL DSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA,
283: $ A( J1+I1, I2 ), 1 )
284: *
1.3 bertrand 285: * Swap A(I1, I2+1:M) with A(I2, I2+1:M)
1.1 bertrand 286: *
1.5 ! bertrand 287: IF( I2.LT.M )
! 288: $ CALL DSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
! 289: $ A( J1+I2-1, I2+1 ), LDA )
1.1 bertrand 290: *
291: * Swap A(I1, I1) with A(I2,I2)
292: *
293: PIV = A( I1+J1-1, I1 )
294: A( J1+I1-1, I1 ) = A( J1+I2-1, I2 )
295: A( J1+I2-1, I2 ) = PIV
296: *
297: * Swap H(I1, 1:J1) with H(I2, 1:J1)
298: *
299: CALL DSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH )
300: IPIV( I1 ) = I2
301: *
302: IF( I1.GT.(K1-1) ) THEN
303: *
304: * Swap L(1:I1-1, I1) with L(1:I1-1, I2),
305: * skipping the first column
306: *
307: CALL DSWAP( I1-K1+1, A( 1, I1 ), 1,
308: $ A( 1, I2 ), 1 )
309: END IF
310: ELSE
311: IPIV( J+1 ) = J+1
312: ENDIF
313: *
314: * Set A(J, J+1) = T(J, J+1)
315: *
316: A( K, J+1 ) = WORK( 2 )
317: *
318: IF( J.LT.NB ) THEN
319: *
1.3 bertrand 320: * Copy A(J+1:M, J+1) into H(J:M, J),
1.1 bertrand 321: *
322: CALL DCOPY( M-J, A( K+1, J+1 ), LDA,
323: $ H( J+1, J+1 ), 1 )
324: END IF
325: *
1.3 bertrand 326: * Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
327: * where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
1.1 bertrand 328: *
1.5 ! bertrand 329: IF( J.LT.(M-1) ) THEN
! 330: IF( A( K, J+1 ).NE.ZERO ) THEN
! 331: ALPHA = ONE / A( K, J+1 )
! 332: CALL DCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
! 333: CALL DSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA )
! 334: ELSE
! 335: CALL DLASET( 'Full', 1, M-J-1, ZERO, ZERO,
! 336: $ A( K, J+2 ), LDA)
! 337: END IF
1.1 bertrand 338: END IF
339: END IF
340: J = J + 1
341: GO TO 10
342: 20 CONTINUE
343: *
344: ELSE
345: *
346: * .....................................................
347: * Factorize A as L*D*L**T using the lower triangle of A
348: * .....................................................
349: *
350: 30 CONTINUE
351: IF( J.GT.MIN( M, NB ) )
352: $ GO TO 40
353: *
354: * K is the column to be factorized
355: * when being called from DSYTRF_AA,
356: * > for the first block column, J1 is 1, hence J1+J-1 is J,
357: * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
358: *
359: K = J1+J-1
1.3 bertrand 360: IF( J.EQ.M ) THEN
361: *
362: * Only need to compute T(J, J)
363: *
364: MJ = 1
365: ELSE
366: MJ = M-J+1
367: END IF
1.1 bertrand 368: *
1.3 bertrand 369: * H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T,
370: * where H(J:M, J) has been initialized to be A(J:M, J)
1.1 bertrand 371: *
372: IF( K.GT.2 ) THEN
373: *
374: * K is the column to be factorized
375: * > for the first block column, K is J, skipping the first two
376: * columns
377: * > for the rest of the columns, K is J+1, skipping only the
378: * first column
379: *
1.3 bertrand 380: CALL DGEMV( 'No transpose', MJ, J-K1,
1.1 bertrand 381: $ -ONE, H( J, K1 ), LDH,
382: $ A( J, 1 ), LDA,
383: $ ONE, H( J, J ), 1 )
384: END IF
385: *
1.3 bertrand 386: * Copy H(J:M, J) into WORK
1.1 bertrand 387: *
1.3 bertrand 388: CALL DCOPY( MJ, H( J, J ), 1, WORK( 1 ), 1 )
1.1 bertrand 389: *
390: IF( J.GT.K1 ) THEN
391: *
1.3 bertrand 392: * Compute WORK := WORK - L(J:M, J-1) * T(J-1,J),
1.1 bertrand 393: * where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
394: *
395: ALPHA = -A( J, K-1 )
1.3 bertrand 396: CALL DAXPY( MJ, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 )
1.1 bertrand 397: END IF
398: *
399: * Set A(J, J) = T(J, J)
400: *
401: A( J, K ) = WORK( 1 )
402: *
403: IF( J.LT.M ) THEN
404: *
1.3 bertrand 405: * Compute WORK(2:M) = T(J, J) L((J+1):M, J)
406: * where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, J)
1.1 bertrand 407: *
408: IF( K.GT.1 ) THEN
409: ALPHA = -A( J, K )
410: CALL DAXPY( M-J, ALPHA, A( J+1, K-1 ), 1,
411: $ WORK( 2 ), 1 )
412: ENDIF
413: *
1.3 bertrand 414: * Find max(|WORK(2:M)|)
1.1 bertrand 415: *
416: I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1
417: PIV = WORK( I2 )
418: *
419: * Apply symmetric pivot
420: *
421: IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN
422: *
423: * Swap WORK(I1) and WORK(I2)
424: *
425: I1 = 2
426: WORK( I2 ) = WORK( I1 )
427: WORK( I1 ) = PIV
428: *
1.3 bertrand 429: * Swap A(I1+1:M, I1) with A(I2, I1+1:M)
1.1 bertrand 430: *
431: I1 = I1+J-1
432: I2 = I2+J-1
433: CALL DSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1,
434: $ A( I2, J1+I1 ), LDA )
435: *
1.3 bertrand 436: * Swap A(I2+1:M, I1) with A(I2+1:M, I2)
1.1 bertrand 437: *
1.5 ! bertrand 438: IF( I2.LT.M )
! 439: $ CALL DSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
! 440: $ A( I2+1, J1+I2-1 ), 1 )
1.1 bertrand 441: *
442: * Swap A(I1, I1) with A(I2, I2)
443: *
444: PIV = A( I1, J1+I1-1 )
445: A( I1, J1+I1-1 ) = A( I2, J1+I2-1 )
446: A( I2, J1+I2-1 ) = PIV
447: *
448: * Swap H(I1, I1:J1) with H(I2, I2:J1)
449: *
450: CALL DSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH )
451: IPIV( I1 ) = I2
452: *
453: IF( I1.GT.(K1-1) ) THEN
454: *
455: * Swap L(1:I1-1, I1) with L(1:I1-1, I2),
456: * skipping the first column
457: *
458: CALL DSWAP( I1-K1+1, A( I1, 1 ), LDA,
459: $ A( I2, 1 ), LDA )
460: END IF
461: ELSE
462: IPIV( J+1 ) = J+1
463: ENDIF
464: *
465: * Set A(J+1, J) = T(J+1, J)
466: *
467: A( J+1, K ) = WORK( 2 )
468: *
469: IF( J.LT.NB ) THEN
470: *
1.3 bertrand 471: * Copy A(J+1:M, J+1) into H(J+1:M, J),
1.1 bertrand 472: *
473: CALL DCOPY( M-J, A( J+1, K+1 ), 1,
474: $ H( J+1, J+1 ), 1 )
475: END IF
476: *
1.3 bertrand 477: * Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
478: * where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
1.1 bertrand 479: *
1.5 ! bertrand 480: IF( J.LT.(M-1) ) THEN
! 481: IF( A( J+1, K ).NE.ZERO ) THEN
! 482: ALPHA = ONE / A( J+1, K )
! 483: CALL DCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
! 484: CALL DSCAL( M-J-1, ALPHA, A( J+2, K ), 1 )
! 485: ELSE
! 486: CALL DLASET( 'Full', M-J-1, 1, ZERO, ZERO,
! 487: $ A( J+2, K ), LDA )
! 488: END IF
1.1 bertrand 489: END IF
490: END IF
491: J = J + 1
492: GO TO 30
493: 40 CONTINUE
494: END IF
495: RETURN
496: *
497: * End of DLASYF_AA
498: *
499: END
CVSweb interface <joel.bertrand@systella.fr>