![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytrf_aa_2stage.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dsytrf_aa_2stage.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: *> \brief \b DSYTRF_AA_2STAGE
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DSYTRF_AA_2STAGE + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrf_aa_2stage.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrf_aa_2stage.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrf_aa_2stage.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV,
22: * IPIV2, WORK, LWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER UPLO
26: * INTEGER N, LDA, LTB, LWORK, INFO
27: * ..
28: * .. Array Arguments ..
29: * INTEGER IPIV( * ), IPIV2( * )
30: * DOUBLE PRECISION A( LDA, * ), TB( * ), WORK( * )
31: * ..
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DSYTRF_AA_2STAGE computes the factorization of a real symmetric matrix A
39: *> using the Aasen's algorithm. The form of the factorization is
40: *>
1.2 ! bertrand 41: *> A = U**T*T*U or A = L*T*L**T
1.1 bertrand 42: *>
43: *> where U (or L) is a product of permutation and unit upper (lower)
44: *> triangular matrices, and T is a symmetric band matrix with the
45: *> bandwidth of NB (NB is internally selected and stored in TB( 1 ), and T is
46: *> LU factorized with partial pivoting).
47: *>
48: *> This is the blocked version of the algorithm, calling Level 3 BLAS.
49: *> \endverbatim
50: *
51: * Arguments:
52: * ==========
53: *
54: *> \param[in] UPLO
55: *> \verbatim
56: *> UPLO is CHARACTER*1
57: *> = 'U': Upper triangle of A is stored;
58: *> = 'L': Lower triangle of A is stored.
59: *> \endverbatim
60: *>
61: *> \param[in] N
62: *> \verbatim
63: *> N is INTEGER
64: *> The order of the matrix A. N >= 0.
65: *> \endverbatim
66: *>
67: *> \param[in,out] A
68: *> \verbatim
69: *> A is DOUBLE PRECISION array, dimension (LDA,N)
70: *> On entry, the symmetric matrix A. If UPLO = 'U', the leading
71: *> N-by-N upper triangular part of A contains the upper
72: *> triangular part of the matrix A, and the strictly lower
73: *> triangular part of A is not referenced. If UPLO = 'L', the
74: *> leading N-by-N lower triangular part of A contains the lower
75: *> triangular part of the matrix A, and the strictly upper
76: *> triangular part of A is not referenced.
77: *>
78: *> On exit, L is stored below (or above) the subdiaonal blocks,
79: *> when UPLO is 'L' (or 'U').
80: *> \endverbatim
81: *>
82: *> \param[in] LDA
83: *> \verbatim
84: *> LDA is INTEGER
85: *> The leading dimension of the array A. LDA >= max(1,N).
86: *> \endverbatim
87: *>
88: *> \param[out] TB
89: *> \verbatim
90: *> TB is DOUBLE PRECISION array, dimension (LTB)
91: *> On exit, details of the LU factorization of the band matrix.
92: *> \endverbatim
93: *>
94: *> \param[in] LTB
95: *> \verbatim
1.2 ! bertrand 96: *> LTB is INTEGER
1.1 bertrand 97: *> The size of the array TB. LTB >= 4*N, internally
98: *> used to select NB such that LTB >= (3*NB+1)*N.
99: *>
100: *> If LTB = -1, then a workspace query is assumed; the
101: *> routine only calculates the optimal size of LTB,
102: *> returns this value as the first entry of TB, and
103: *> no error message related to LTB is issued by XERBLA.
104: *> \endverbatim
105: *>
1.2 ! bertrand 106: *> \param[out] IPIV
! 107: *> \verbatim
! 108: *> IPIV is INTEGER array, dimension (N)
! 109: *> On exit, it contains the details of the interchanges, i.e.,
! 110: *> the row and column k of A were interchanged with the
! 111: *> row and column IPIV(k).
! 112: *> \endverbatim
! 113: *>
! 114: *> \param[out] IPIV2
! 115: *> \verbatim
! 116: *> IPIV2 is INTEGER array, dimension (N)
! 117: *> On exit, it contains the details of the interchanges, i.e.,
! 118: *> the row and column k of T were interchanged with the
! 119: *> row and column IPIV2(k).
! 120: *> \endverbatim
! 121: *>
1.1 bertrand 122: *> \param[out] WORK
123: *> \verbatim
124: *> WORK is DOUBLE PRECISION workspace of size LWORK
125: *> \endverbatim
126: *>
127: *> \param[in] LWORK
128: *> \verbatim
1.2 ! bertrand 129: *> LWORK is INTEGER
1.1 bertrand 130: *> The size of WORK. LWORK >= N, internally used to select NB
131: *> such that LWORK >= N*NB.
132: *>
133: *> If LWORK = -1, then a workspace query is assumed; the
134: *> routine only calculates the optimal size of the WORK array,
135: *> returns this value as the first entry of the WORK array, and
136: *> no error message related to LWORK is issued by XERBLA.
137: *> \endverbatim
138: *>
139: *> \param[out] INFO
140: *> \verbatim
141: *> INFO is INTEGER
142: *> = 0: successful exit
143: *> < 0: if INFO = -i, the i-th argument had an illegal value.
144: *> > 0: if INFO = i, band LU factorization failed on i-th column
145: *> \endverbatim
146: *
147: * Authors:
148: * ========
149: *
150: *> \author Univ. of Tennessee
151: *> \author Univ. of California Berkeley
152: *> \author Univ. of Colorado Denver
153: *> \author NAG Ltd.
154: *
155: *> \date November 2017
156: *
157: *> \ingroup doubleSYcomputational
158: *
159: * =====================================================================
160: SUBROUTINE DSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV,
161: $ IPIV2, WORK, LWORK, INFO )
162: *
163: * -- LAPACK computational routine (version 3.8.0) --
164: * -- LAPACK is a software package provided by Univ. of Tennessee, --
165: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166: * November 2017
167: *
168: IMPLICIT NONE
169: *
170: * .. Scalar Arguments ..
171: CHARACTER UPLO
172: INTEGER N, LDA, LTB, LWORK, INFO
173: * ..
174: * .. Array Arguments ..
175: INTEGER IPIV( * ), IPIV2( * )
176: DOUBLE PRECISION A( LDA, * ), TB( * ), WORK( * )
177: * ..
178: *
179: * =====================================================================
180: * .. Parameters ..
181: DOUBLE PRECISION ZERO, ONE
182: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
183: *
184: * .. Local Scalars ..
185: LOGICAL UPPER, TQUERY, WQUERY
186: INTEGER I, J, K, I1, I2, TD
187: INTEGER LDTB, NB, KB, JB, NT, IINFO
188: DOUBLE PRECISION PIV
189: * ..
190: * .. External Functions ..
191: LOGICAL LSAME
192: INTEGER ILAENV
193: EXTERNAL LSAME, ILAENV
194: * ..
195: * .. External Subroutines ..
1.2 ! bertrand 196: EXTERNAL XERBLA, DCOPY, DLACPY,
1.1 bertrand 197: $ DLASET, DGBTRF, DGEMM, DGETRF,
198: $ DSYGST, DSWAP, DTRSM
199: * ..
200: * .. Intrinsic Functions ..
201: INTRINSIC MIN, MAX
202: * ..
203: * .. Executable Statements ..
204: *
205: * Test the input parameters.
206: *
207: INFO = 0
208: UPPER = LSAME( UPLO, 'U' )
209: WQUERY = ( LWORK.EQ.-1 )
210: TQUERY = ( LTB.EQ.-1 )
211: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
212: INFO = -1
213: ELSE IF( N.LT.0 ) THEN
214: INFO = -2
215: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
216: INFO = -4
217: ELSE IF ( LTB .LT. 4*N .AND. .NOT.TQUERY ) THEN
218: INFO = -6
219: ELSE IF ( LWORK .LT. N .AND. .NOT.WQUERY ) THEN
220: INFO = -10
221: END IF
222: *
223: IF( INFO.NE.0 ) THEN
224: CALL XERBLA( 'DSYTRF_AA_2STAGE', -INFO )
225: RETURN
226: END IF
227: *
228: * Answer the query
229: *
230: NB = ILAENV( 1, 'DSYTRF_AA_2STAGE', UPLO, N, -1, -1, -1 )
231: IF( INFO.EQ.0 ) THEN
232: IF( TQUERY ) THEN
233: TB( 1 ) = (3*NB+1)*N
234: END IF
235: IF( WQUERY ) THEN
236: WORK( 1 ) = N*NB
237: END IF
238: END IF
239: IF( TQUERY .OR. WQUERY ) THEN
240: RETURN
241: END IF
242: *
243: * Quick return
244: *
245: IF ( N.EQ.0 ) THEN
246: RETURN
247: ENDIF
248: *
249: * Determine the number of the block size
250: *
251: LDTB = LTB/N
252: IF( LDTB .LT. 3*NB+1 ) THEN
253: NB = (LDTB-1)/3
254: END IF
255: IF( LWORK .LT. NB*N ) THEN
256: NB = LWORK/N
257: END IF
258: *
259: * Determine the number of the block columns
260: *
261: NT = (N+NB-1)/NB
262: TD = 2*NB
263: KB = MIN(NB, N)
264: *
265: * Initialize vectors/matrices
266: *
267: DO J = 1, KB
268: IPIV( J ) = J
269: END DO
270: *
271: * Save NB
272: *
273: TB( 1 ) = NB
274: *
275: IF( UPPER ) THEN
276: *
277: * .....................................................
1.2 ! bertrand 278: * Factorize A as U**T*D*U using the upper triangle of A
1.1 bertrand 279: * .....................................................
280: *
281: DO J = 0, NT-1
282: *
283: * Generate Jth column of W and H
284: *
285: KB = MIN(NB, N-J*NB)
286: DO I = 1, J-1
287: IF( I .EQ. 1 ) THEN
288: * H(I,J) = T(I,I)*U(I,J) + T(I,I+1)*U(I+1,J)
289: IF( I .EQ. (J-1) ) THEN
290: JB = NB+KB
291: ELSE
292: JB = 2*NB
293: END IF
294: CALL DGEMM( 'NoTranspose', 'NoTranspose',
295: $ NB, KB, JB,
296: $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
297: $ A( (I-1)*NB+1, J*NB+1 ), LDA,
298: $ ZERO, WORK( I*NB+1 ), N )
299: ELSE
300: * H(I,J) = T(I,I-1)*U(I-1,J) + T(I,I)*U(I,J) + T(I,I+1)*U(I+1,J)
301: IF( I .EQ. J-1) THEN
302: JB = 2*NB+KB
303: ELSE
304: JB = 3*NB
305: END IF
306: CALL DGEMM( 'NoTranspose', 'NoTranspose',
307: $ NB, KB, JB,
308: $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
309: $ LDTB-1,
310: $ A( (I-2)*NB+1, J*NB+1 ), LDA,
311: $ ZERO, WORK( I*NB+1 ), N )
312: END IF
313: END DO
314: *
315: * Compute T(J,J)
316: *
317: CALL DLACPY( 'Upper', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
318: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
319: IF( J.GT.1 ) THEN
320: * T(J,J) = U(1:J,J)'*H(1:J)
321: CALL DGEMM( 'Transpose', 'NoTranspose',
322: $ KB, KB, (J-1)*NB,
323: $ -ONE, A( 1, J*NB+1 ), LDA,
324: $ WORK( NB+1 ), N,
325: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
326: * T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J)
327: CALL DGEMM( 'Transpose', 'NoTranspose',
328: $ KB, NB, KB,
329: $ ONE, A( (J-1)*NB+1, J*NB+1 ), LDA,
330: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
331: $ ZERO, WORK( 1 ), N )
332: CALL DGEMM( 'NoTranspose', 'NoTranspose',
333: $ KB, KB, NB,
334: $ -ONE, WORK( 1 ), N,
335: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
336: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
337: END IF
338: IF( J.GT.0 ) THEN
339: CALL DSYGST( 1, 'Upper', KB,
340: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
341: $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO )
342: END IF
343: *
344: * Expand T(J,J) into full format
345: *
346: DO I = 1, KB
347: DO K = I+1, KB
348: TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
349: $ = TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
350: END DO
351: END DO
352: *
353: IF( J.LT.NT-1 ) THEN
354: IF( J.GT.0 ) THEN
355: *
356: * Compute H(J,J)
357: *
358: IF( J.EQ.1 ) THEN
359: CALL DGEMM( 'NoTranspose', 'NoTranspose',
360: $ KB, KB, KB,
361: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
362: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
363: $ ZERO, WORK( J*NB+1 ), N )
364: ELSE
365: CALL DGEMM( 'NoTranspose', 'NoTranspose',
366: $ KB, KB, NB+KB,
367: $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
368: $ LDTB-1,
369: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
370: $ ZERO, WORK( J*NB+1 ), N )
371: END IF
372: *
373: * Update with the previous column
374: *
375: CALL DGEMM( 'Transpose', 'NoTranspose',
376: $ NB, N-(J+1)*NB, J*NB,
377: $ -ONE, WORK( NB+1 ), N,
378: $ A( 1, (J+1)*NB+1 ), LDA,
379: $ ONE, A( J*NB+1, (J+1)*NB+1 ), LDA )
380: END IF
381: *
382: * Copy panel to workspace to call DGETRF
383: *
384: DO K = 1, NB
385: CALL DCOPY( N-(J+1)*NB,
386: $ A( J*NB+K, (J+1)*NB+1 ), LDA,
387: $ WORK( 1+(K-1)*N ), 1 )
388: END DO
389: *
390: * Factorize panel
391: *
392: CALL DGETRF( N-(J+1)*NB, NB,
393: $ WORK, N,
394: $ IPIV( (J+1)*NB+1 ), IINFO )
395: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
396: c INFO = IINFO+(J+1)*NB
397: c END IF
398: *
399: * Copy panel back
400: *
401: DO K = 1, NB
402: CALL DCOPY( N-(J+1)*NB,
403: $ WORK( 1+(K-1)*N ), 1,
404: $ A( J*NB+K, (J+1)*NB+1 ), LDA )
405: END DO
406: *
407: * Compute T(J+1, J), zero out for GEMM update
408: *
409: KB = MIN(NB, N-(J+1)*NB)
410: CALL DLASET( 'Full', KB, NB, ZERO, ZERO,
411: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
412: CALL DLACPY( 'Upper', KB, NB,
413: $ WORK, N,
414: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
415: IF( J.GT.0 ) THEN
416: CALL DTRSM( 'R', 'U', 'N', 'U', KB, NB, ONE,
417: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
418: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
419: END IF
420: *
421: * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
422: * updates
423: *
424: DO K = 1, NB
425: DO I = 1, KB
426: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
427: $ = TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
428: END DO
429: END DO
430: CALL DLASET( 'Lower', KB, NB, ZERO, ONE,
431: $ A( J*NB+1, (J+1)*NB+1), LDA )
432: *
433: * Apply pivots to trailing submatrix of A
434: *
435: DO K = 1, KB
436: * > Adjust ipiv
437: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
438: *
439: I1 = (J+1)*NB+K
440: I2 = IPIV( (J+1)*NB+K )
441: IF( I1.NE.I2 ) THEN
442: * > Apply pivots to previous columns of L
443: CALL DSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
444: $ A( (J+1)*NB+1, I2 ), 1 )
1.2 ! bertrand 445: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
! 446: IF( I2.GT.(I1+1) )
! 447: $ CALL DSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
! 448: $ A( I1+1, I2 ), 1 )
1.1 bertrand 449: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
1.2 ! bertrand 450: IF( I2.LT.N )
! 451: $ CALL DSWAP( N-I2, A( I1, I2+1 ), LDA,
! 452: $ A( I2, I2+1 ), LDA )
1.1 bertrand 453: * > Swap A(I1, I1) with A(I2, I2)
454: PIV = A( I1, I1 )
455: A( I1, I1 ) = A( I2, I2 )
456: A( I2, I2 ) = PIV
457: * > Apply pivots to previous columns of L
458: IF( J.GT.0 ) THEN
459: CALL DSWAP( J*NB, A( 1, I1 ), 1,
460: $ A( 1, I2 ), 1 )
461: END IF
462: ENDIF
463: END DO
464: END IF
465: END DO
466: ELSE
467: *
468: * .....................................................
469: * Factorize A as L*D*L**T using the lower triangle of A
470: * .....................................................
471: *
472: DO J = 0, NT-1
473: *
474: * Generate Jth column of W and H
475: *
476: KB = MIN(NB, N-J*NB)
477: DO I = 1, J-1
478: IF( I.EQ.1 ) THEN
479: * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
480: IF( I .EQ. J-1) THEN
481: JB = NB+KB
482: ELSE
483: JB = 2*NB
484: END IF
485: CALL DGEMM( 'NoTranspose', 'Transpose',
486: $ NB, KB, JB,
487: $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
488: $ A( J*NB+1, (I-1)*NB+1 ), LDA,
489: $ ZERO, WORK( I*NB+1 ), N )
490: ELSE
491: * H(I,J) = T(I,I-1)*L(J,I-1)' + T(I,I)*L(J,I)' + T(I,I+1)*L(J,I+1)'
492: IF( I .EQ. J-1) THEN
493: JB = 2*NB+KB
494: ELSE
495: JB = 3*NB
496: END IF
497: CALL DGEMM( 'NoTranspose', 'Transpose',
498: $ NB, KB, JB,
499: $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
500: $ LDTB-1,
501: $ A( J*NB+1, (I-2)*NB+1 ), LDA,
502: $ ZERO, WORK( I*NB+1 ), N )
503: END IF
504: END DO
505: *
506: * Compute T(J,J)
507: *
508: CALL DLACPY( 'Lower', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
509: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
510: IF( J.GT.1 ) THEN
511: * T(J,J) = L(J,1:J)*H(1:J)
512: CALL DGEMM( 'NoTranspose', 'NoTranspose',
513: $ KB, KB, (J-1)*NB,
514: $ -ONE, A( J*NB+1, 1 ), LDA,
515: $ WORK( NB+1 ), N,
516: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
517: * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
518: CALL DGEMM( 'NoTranspose', 'NoTranspose',
519: $ KB, NB, KB,
520: $ ONE, A( J*NB+1, (J-1)*NB+1 ), LDA,
521: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
522: $ ZERO, WORK( 1 ), N )
523: CALL DGEMM( 'NoTranspose', 'Transpose',
524: $ KB, KB, NB,
525: $ -ONE, WORK( 1 ), N,
526: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
527: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
528: END IF
529: IF( J.GT.0 ) THEN
530: CALL DSYGST( 1, 'Lower', KB,
531: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
532: $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO )
533: END IF
534: *
535: * Expand T(J,J) into full format
536: *
537: DO I = 1, KB
538: DO K = I+1, KB
539: TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
540: $ = TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
541: END DO
542: END DO
543: *
544: IF( J.LT.NT-1 ) THEN
545: IF( J.GT.0 ) THEN
546: *
547: * Compute H(J,J)
548: *
549: IF( J.EQ.1 ) THEN
550: CALL DGEMM( 'NoTranspose', 'Transpose',
551: $ KB, KB, KB,
552: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
553: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
554: $ ZERO, WORK( J*NB+1 ), N )
555: ELSE
556: CALL DGEMM( 'NoTranspose', 'Transpose',
557: $ KB, KB, NB+KB,
558: $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
559: $ LDTB-1,
560: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
561: $ ZERO, WORK( J*NB+1 ), N )
562: END IF
563: *
564: * Update with the previous column
565: *
566: CALL DGEMM( 'NoTranspose', 'NoTranspose',
567: $ N-(J+1)*NB, NB, J*NB,
568: $ -ONE, A( (J+1)*NB+1, 1 ), LDA,
569: $ WORK( NB+1 ), N,
570: $ ONE, A( (J+1)*NB+1, J*NB+1 ), LDA )
571: END IF
572: *
573: * Factorize panel
574: *
575: CALL DGETRF( N-(J+1)*NB, NB,
576: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
577: $ IPIV( (J+1)*NB+1 ), IINFO )
578: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
579: c INFO = IINFO+(J+1)*NB
580: c END IF
581: *
582: * Compute T(J+1, J), zero out for GEMM update
583: *
584: KB = MIN(NB, N-(J+1)*NB)
585: CALL DLASET( 'Full', KB, NB, ZERO, ZERO,
586: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
587: CALL DLACPY( 'Upper', KB, NB,
588: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
589: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
590: IF( J.GT.0 ) THEN
591: CALL DTRSM( 'R', 'L', 'T', 'U', KB, NB, ONE,
592: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
593: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
594: END IF
595: *
596: * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
597: * updates
598: *
599: DO K = 1, NB
600: DO I = 1, KB
601: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
602: $ = TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
603: END DO
604: END DO
605: CALL DLASET( 'Upper', KB, NB, ZERO, ONE,
606: $ A( (J+1)*NB+1, J*NB+1), LDA )
607: *
608: * Apply pivots to trailing submatrix of A
609: *
610: DO K = 1, KB
611: * > Adjust ipiv
612: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
613: *
614: I1 = (J+1)*NB+K
615: I2 = IPIV( (J+1)*NB+K )
616: IF( I1.NE.I2 ) THEN
617: * > Apply pivots to previous columns of L
618: CALL DSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
619: $ A( I2, (J+1)*NB+1 ), LDA )
620: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
1.2 ! bertrand 621: IF( I2.GT.(I1+1) )
! 622: $ CALL DSWAP( I2-I1-1, A( I1+1, I1 ), 1,
! 623: $ A( I2, I1+1 ), LDA )
1.1 bertrand 624: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
1.2 ! bertrand 625: IF( I2.LT.N )
! 626: $ CALL DSWAP( N-I2, A( I2+1, I1 ), 1,
! 627: $ A( I2+1, I2 ), 1 )
1.1 bertrand 628: * > Swap A(I1, I1) with A(I2, I2)
629: PIV = A( I1, I1 )
630: A( I1, I1 ) = A( I2, I2 )
631: A( I2, I2 ) = PIV
632: * > Apply pivots to previous columns of L
633: IF( J.GT.0 ) THEN
634: CALL DSWAP( J*NB, A( I1, 1 ), LDA,
635: $ A( I2, 1 ), LDA )
636: END IF
637: ENDIF
638: END DO
639: *
640: * Apply pivots to previous columns of L
641: *
642: c CALL DLASWP( J*NB, A( 1, 1 ), LDA,
643: c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
644: END IF
645: END DO
646: END IF
647: *
648: * Factor the band matrix
649: CALL DGBTRF( N, N, NB, NB, TB, LDTB, IPIV2, INFO )
650: *
1.2 ! bertrand 651: RETURN
! 652: *
1.1 bertrand 653: * End of DSYTRF_AA_2STAGE
654: *
655: END