![[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: *> \ingroup doubleSYcomputational
156: *
157: * =====================================================================
158: SUBROUTINE DSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV,
159: $ IPIV2, WORK, LWORK, INFO )
160: *
1.3 ! bertrand 161: * -- LAPACK computational routine --
1.1 bertrand 162: * -- LAPACK is a software package provided by Univ. of Tennessee, --
163: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
164: *
165: IMPLICIT NONE
166: *
167: * .. Scalar Arguments ..
168: CHARACTER UPLO
169: INTEGER N, LDA, LTB, LWORK, INFO
170: * ..
171: * .. Array Arguments ..
172: INTEGER IPIV( * ), IPIV2( * )
173: DOUBLE PRECISION A( LDA, * ), TB( * ), WORK( * )
174: * ..
175: *
176: * =====================================================================
177: * .. Parameters ..
178: DOUBLE PRECISION ZERO, ONE
179: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
180: *
181: * .. Local Scalars ..
182: LOGICAL UPPER, TQUERY, WQUERY
183: INTEGER I, J, K, I1, I2, TD
184: INTEGER LDTB, NB, KB, JB, NT, IINFO
185: DOUBLE PRECISION PIV
186: * ..
187: * .. External Functions ..
188: LOGICAL LSAME
189: INTEGER ILAENV
190: EXTERNAL LSAME, ILAENV
191: * ..
192: * .. External Subroutines ..
1.2 bertrand 193: EXTERNAL XERBLA, DCOPY, DLACPY,
1.1 bertrand 194: $ DLASET, DGBTRF, DGEMM, DGETRF,
195: $ DSYGST, DSWAP, DTRSM
196: * ..
197: * .. Intrinsic Functions ..
198: INTRINSIC MIN, MAX
199: * ..
200: * .. Executable Statements ..
201: *
202: * Test the input parameters.
203: *
204: INFO = 0
205: UPPER = LSAME( UPLO, 'U' )
206: WQUERY = ( LWORK.EQ.-1 )
207: TQUERY = ( LTB.EQ.-1 )
208: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
209: INFO = -1
210: ELSE IF( N.LT.0 ) THEN
211: INFO = -2
212: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
213: INFO = -4
214: ELSE IF ( LTB .LT. 4*N .AND. .NOT.TQUERY ) THEN
215: INFO = -6
216: ELSE IF ( LWORK .LT. N .AND. .NOT.WQUERY ) THEN
217: INFO = -10
218: END IF
219: *
220: IF( INFO.NE.0 ) THEN
221: CALL XERBLA( 'DSYTRF_AA_2STAGE', -INFO )
222: RETURN
223: END IF
224: *
225: * Answer the query
226: *
227: NB = ILAENV( 1, 'DSYTRF_AA_2STAGE', UPLO, N, -1, -1, -1 )
228: IF( INFO.EQ.0 ) THEN
229: IF( TQUERY ) THEN
230: TB( 1 ) = (3*NB+1)*N
231: END IF
232: IF( WQUERY ) THEN
233: WORK( 1 ) = N*NB
234: END IF
235: END IF
236: IF( TQUERY .OR. WQUERY ) THEN
237: RETURN
238: END IF
239: *
240: * Quick return
241: *
242: IF ( N.EQ.0 ) THEN
243: RETURN
244: ENDIF
245: *
246: * Determine the number of the block size
247: *
248: LDTB = LTB/N
249: IF( LDTB .LT. 3*NB+1 ) THEN
250: NB = (LDTB-1)/3
251: END IF
252: IF( LWORK .LT. NB*N ) THEN
253: NB = LWORK/N
254: END IF
255: *
256: * Determine the number of the block columns
257: *
258: NT = (N+NB-1)/NB
259: TD = 2*NB
260: KB = MIN(NB, N)
261: *
262: * Initialize vectors/matrices
263: *
264: DO J = 1, KB
265: IPIV( J ) = J
266: END DO
267: *
268: * Save NB
269: *
270: TB( 1 ) = NB
271: *
272: IF( UPPER ) THEN
273: *
274: * .....................................................
1.2 bertrand 275: * Factorize A as U**T*D*U using the upper triangle of A
1.1 bertrand 276: * .....................................................
277: *
278: DO J = 0, NT-1
279: *
280: * Generate Jth column of W and H
281: *
282: KB = MIN(NB, N-J*NB)
283: DO I = 1, J-1
284: IF( I .EQ. 1 ) THEN
285: * H(I,J) = T(I,I)*U(I,J) + T(I,I+1)*U(I+1,J)
286: IF( I .EQ. (J-1) ) THEN
287: JB = NB+KB
288: ELSE
289: JB = 2*NB
290: END IF
291: CALL DGEMM( 'NoTranspose', 'NoTranspose',
292: $ NB, KB, JB,
293: $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
294: $ A( (I-1)*NB+1, J*NB+1 ), LDA,
295: $ ZERO, WORK( I*NB+1 ), N )
296: ELSE
297: * 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)
298: IF( I .EQ. J-1) THEN
299: JB = 2*NB+KB
300: ELSE
301: JB = 3*NB
302: END IF
303: CALL DGEMM( 'NoTranspose', 'NoTranspose',
304: $ NB, KB, JB,
305: $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
306: $ LDTB-1,
307: $ A( (I-2)*NB+1, J*NB+1 ), LDA,
308: $ ZERO, WORK( I*NB+1 ), N )
309: END IF
310: END DO
311: *
312: * Compute T(J,J)
313: *
314: CALL DLACPY( 'Upper', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
315: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
316: IF( J.GT.1 ) THEN
317: * T(J,J) = U(1:J,J)'*H(1:J)
318: CALL DGEMM( 'Transpose', 'NoTranspose',
319: $ KB, KB, (J-1)*NB,
320: $ -ONE, A( 1, J*NB+1 ), LDA,
321: $ WORK( NB+1 ), N,
322: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
323: * T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J)
324: CALL DGEMM( 'Transpose', 'NoTranspose',
325: $ KB, NB, KB,
326: $ ONE, A( (J-1)*NB+1, J*NB+1 ), LDA,
327: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
328: $ ZERO, WORK( 1 ), N )
329: CALL DGEMM( 'NoTranspose', 'NoTranspose',
330: $ KB, KB, NB,
331: $ -ONE, WORK( 1 ), N,
332: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
333: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
334: END IF
335: IF( J.GT.0 ) THEN
336: CALL DSYGST( 1, 'Upper', KB,
337: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
338: $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO )
339: END IF
340: *
341: * Expand T(J,J) into full format
342: *
343: DO I = 1, KB
344: DO K = I+1, KB
345: TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
346: $ = TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
347: END DO
348: END DO
349: *
350: IF( J.LT.NT-1 ) THEN
351: IF( J.GT.0 ) THEN
352: *
353: * Compute H(J,J)
354: *
355: IF( J.EQ.1 ) THEN
356: CALL DGEMM( 'NoTranspose', 'NoTranspose',
357: $ KB, KB, KB,
358: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
359: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
360: $ ZERO, WORK( J*NB+1 ), N )
361: ELSE
362: CALL DGEMM( 'NoTranspose', 'NoTranspose',
363: $ KB, KB, NB+KB,
364: $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
365: $ LDTB-1,
366: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
367: $ ZERO, WORK( J*NB+1 ), N )
368: END IF
369: *
370: * Update with the previous column
371: *
372: CALL DGEMM( 'Transpose', 'NoTranspose',
373: $ NB, N-(J+1)*NB, J*NB,
374: $ -ONE, WORK( NB+1 ), N,
375: $ A( 1, (J+1)*NB+1 ), LDA,
376: $ ONE, A( J*NB+1, (J+1)*NB+1 ), LDA )
377: END IF
378: *
379: * Copy panel to workspace to call DGETRF
380: *
381: DO K = 1, NB
382: CALL DCOPY( N-(J+1)*NB,
383: $ A( J*NB+K, (J+1)*NB+1 ), LDA,
384: $ WORK( 1+(K-1)*N ), 1 )
385: END DO
386: *
387: * Factorize panel
388: *
389: CALL DGETRF( N-(J+1)*NB, NB,
390: $ WORK, N,
391: $ IPIV( (J+1)*NB+1 ), IINFO )
392: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
393: c INFO = IINFO+(J+1)*NB
394: c END IF
395: *
396: * Copy panel back
397: *
398: DO K = 1, NB
399: CALL DCOPY( N-(J+1)*NB,
400: $ WORK( 1+(K-1)*N ), 1,
401: $ A( J*NB+K, (J+1)*NB+1 ), LDA )
402: END DO
403: *
404: * Compute T(J+1, J), zero out for GEMM update
405: *
406: KB = MIN(NB, N-(J+1)*NB)
407: CALL DLASET( 'Full', KB, NB, ZERO, ZERO,
408: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
409: CALL DLACPY( 'Upper', KB, NB,
410: $ WORK, N,
411: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
412: IF( J.GT.0 ) THEN
413: CALL DTRSM( 'R', 'U', 'N', 'U', KB, NB, ONE,
414: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
415: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
416: END IF
417: *
418: * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
419: * updates
420: *
421: DO K = 1, NB
422: DO I = 1, KB
423: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
424: $ = TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
425: END DO
426: END DO
427: CALL DLASET( 'Lower', KB, NB, ZERO, ONE,
428: $ A( J*NB+1, (J+1)*NB+1), LDA )
429: *
430: * Apply pivots to trailing submatrix of A
431: *
432: DO K = 1, KB
433: * > Adjust ipiv
434: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
435: *
436: I1 = (J+1)*NB+K
437: I2 = IPIV( (J+1)*NB+K )
438: IF( I1.NE.I2 ) THEN
439: * > Apply pivots to previous columns of L
440: CALL DSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
441: $ A( (J+1)*NB+1, I2 ), 1 )
1.2 bertrand 442: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
443: IF( I2.GT.(I1+1) )
444: $ CALL DSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
445: $ A( I1+1, I2 ), 1 )
1.1 bertrand 446: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
1.2 bertrand 447: IF( I2.LT.N )
448: $ CALL DSWAP( N-I2, A( I1, I2+1 ), LDA,
449: $ A( I2, I2+1 ), LDA )
1.1 bertrand 450: * > Swap A(I1, I1) with A(I2, I2)
451: PIV = A( I1, I1 )
452: A( I1, I1 ) = A( I2, I2 )
453: A( I2, I2 ) = PIV
454: * > Apply pivots to previous columns of L
455: IF( J.GT.0 ) THEN
456: CALL DSWAP( J*NB, A( 1, I1 ), 1,
457: $ A( 1, I2 ), 1 )
458: END IF
459: ENDIF
460: END DO
461: END IF
462: END DO
463: ELSE
464: *
465: * .....................................................
466: * Factorize A as L*D*L**T using the lower triangle of A
467: * .....................................................
468: *
469: DO J = 0, NT-1
470: *
471: * Generate Jth column of W and H
472: *
473: KB = MIN(NB, N-J*NB)
474: DO I = 1, J-1
475: IF( I.EQ.1 ) THEN
476: * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
477: IF( I .EQ. J-1) THEN
478: JB = NB+KB
479: ELSE
480: JB = 2*NB
481: END IF
482: CALL DGEMM( 'NoTranspose', 'Transpose',
483: $ NB, KB, JB,
484: $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
485: $ A( J*NB+1, (I-1)*NB+1 ), LDA,
486: $ ZERO, WORK( I*NB+1 ), N )
487: ELSE
488: * 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)'
489: IF( I .EQ. J-1) THEN
490: JB = 2*NB+KB
491: ELSE
492: JB = 3*NB
493: END IF
494: CALL DGEMM( 'NoTranspose', 'Transpose',
495: $ NB, KB, JB,
496: $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
497: $ LDTB-1,
498: $ A( J*NB+1, (I-2)*NB+1 ), LDA,
499: $ ZERO, WORK( I*NB+1 ), N )
500: END IF
501: END DO
502: *
503: * Compute T(J,J)
504: *
505: CALL DLACPY( 'Lower', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
506: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
507: IF( J.GT.1 ) THEN
508: * T(J,J) = L(J,1:J)*H(1:J)
509: CALL DGEMM( 'NoTranspose', 'NoTranspose',
510: $ KB, KB, (J-1)*NB,
511: $ -ONE, A( J*NB+1, 1 ), LDA,
512: $ WORK( NB+1 ), N,
513: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
514: * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
515: CALL DGEMM( 'NoTranspose', 'NoTranspose',
516: $ KB, NB, KB,
517: $ ONE, A( J*NB+1, (J-1)*NB+1 ), LDA,
518: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
519: $ ZERO, WORK( 1 ), N )
520: CALL DGEMM( 'NoTranspose', 'Transpose',
521: $ KB, KB, NB,
522: $ -ONE, WORK( 1 ), N,
523: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
524: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
525: END IF
526: IF( J.GT.0 ) THEN
527: CALL DSYGST( 1, 'Lower', KB,
528: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
529: $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO )
530: END IF
531: *
532: * Expand T(J,J) into full format
533: *
534: DO I = 1, KB
535: DO K = I+1, KB
536: TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
537: $ = TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
538: END DO
539: END DO
540: *
541: IF( J.LT.NT-1 ) THEN
542: IF( J.GT.0 ) THEN
543: *
544: * Compute H(J,J)
545: *
546: IF( J.EQ.1 ) THEN
547: CALL DGEMM( 'NoTranspose', 'Transpose',
548: $ KB, KB, KB,
549: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
550: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
551: $ ZERO, WORK( J*NB+1 ), N )
552: ELSE
553: CALL DGEMM( 'NoTranspose', 'Transpose',
554: $ KB, KB, NB+KB,
555: $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
556: $ LDTB-1,
557: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
558: $ ZERO, WORK( J*NB+1 ), N )
559: END IF
560: *
561: * Update with the previous column
562: *
563: CALL DGEMM( 'NoTranspose', 'NoTranspose',
564: $ N-(J+1)*NB, NB, J*NB,
565: $ -ONE, A( (J+1)*NB+1, 1 ), LDA,
566: $ WORK( NB+1 ), N,
567: $ ONE, A( (J+1)*NB+1, J*NB+1 ), LDA )
568: END IF
569: *
570: * Factorize panel
571: *
572: CALL DGETRF( N-(J+1)*NB, NB,
573: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
574: $ IPIV( (J+1)*NB+1 ), IINFO )
575: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
576: c INFO = IINFO+(J+1)*NB
577: c END IF
578: *
579: * Compute T(J+1, J), zero out for GEMM update
580: *
581: KB = MIN(NB, N-(J+1)*NB)
582: CALL DLASET( 'Full', KB, NB, ZERO, ZERO,
583: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
584: CALL DLACPY( 'Upper', KB, NB,
585: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
586: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
587: IF( J.GT.0 ) THEN
588: CALL DTRSM( 'R', 'L', 'T', 'U', KB, NB, ONE,
589: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
590: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
591: END IF
592: *
593: * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
594: * updates
595: *
596: DO K = 1, NB
597: DO I = 1, KB
598: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
599: $ = TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
600: END DO
601: END DO
602: CALL DLASET( 'Upper', KB, NB, ZERO, ONE,
603: $ A( (J+1)*NB+1, J*NB+1), LDA )
604: *
605: * Apply pivots to trailing submatrix of A
606: *
607: DO K = 1, KB
608: * > Adjust ipiv
609: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
610: *
611: I1 = (J+1)*NB+K
612: I2 = IPIV( (J+1)*NB+K )
613: IF( I1.NE.I2 ) THEN
614: * > Apply pivots to previous columns of L
615: CALL DSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
616: $ A( I2, (J+1)*NB+1 ), LDA )
617: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
1.2 bertrand 618: IF( I2.GT.(I1+1) )
619: $ CALL DSWAP( I2-I1-1, A( I1+1, I1 ), 1,
620: $ A( I2, I1+1 ), LDA )
1.1 bertrand 621: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
1.2 bertrand 622: IF( I2.LT.N )
623: $ CALL DSWAP( N-I2, A( I2+1, I1 ), 1,
624: $ A( I2+1, I2 ), 1 )
1.1 bertrand 625: * > Swap A(I1, I1) with A(I2, I2)
626: PIV = A( I1, I1 )
627: A( I1, I1 ) = A( I2, I2 )
628: A( I2, I2 ) = PIV
629: * > Apply pivots to previous columns of L
630: IF( J.GT.0 ) THEN
631: CALL DSWAP( J*NB, A( I1, 1 ), LDA,
632: $ A( I2, 1 ), LDA )
633: END IF
634: ENDIF
635: END DO
636: *
637: * Apply pivots to previous columns of L
638: *
639: c CALL DLASWP( J*NB, A( 1, 1 ), LDA,
640: c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
641: END IF
642: END DO
643: END IF
644: *
645: * Factor the band matrix
646: CALL DGBTRF( N, N, NB, NB, TB, LDTB, IPIV2, INFO )
647: *
1.2 bertrand 648: RETURN
649: *
1.1 bertrand 650: * End of DSYTRF_AA_2STAGE
651: *
652: END