1: *> \brief \b ZSYTRF_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 ZSYTRF_AA_2STAGE + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytrf_aa_2stage.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytrf_aa_2stage.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytrf_aa_2stage.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZSYTRF_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: * COMPLEX*16 A( LDA, * ), TB( * ), WORK( * )
31: * ..
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZSYTRF_AA_2STAGE computes the factorization of a complex symmetric matrix A
39: *> using the Aasen's algorithm. The form of the factorization is
40: *>
41: *> A = U**T*T*U or A = L*T*L**T
42: *>
43: *> where U (or L) is a product of permutation and unit upper (lower)
44: *> triangular matrices, and T is a complex 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 COMPLEX*16 array, dimension (LDA,N)
70: *> On entry, the hermitian 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 COMPLEX*16 array, dimension (LTB)
91: *> On exit, details of the LU factorization of the band matrix.
92: *> \endverbatim
93: *>
94: *> \param[in] LTB
95: *> \verbatim
96: *> LTB is INTEGER
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: *>
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 IPIV(k).
120: *> \endverbatim
121: *>
122: *> \param[out] WORK
123: *> \verbatim
124: *> WORK is COMPLEX*16 workspace of size LWORK
125: *> \endverbatim
126: *>
127: *> \param[in] LWORK
128: *> \verbatim
129: *> LWORK is INTEGER
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 complex16SYcomputational
156: *
157: * =====================================================================
158: SUBROUTINE ZSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV,
159: $ IPIV2, WORK, LWORK, INFO )
160: *
161: * -- LAPACK computational routine --
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: COMPLEX*16 A( LDA, * ), TB( * ), WORK( * )
174: * ..
175: *
176: * =====================================================================
177: * .. Parameters ..
178: COMPLEX*16 CZERO, CONE
179: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
180: $ CONE = ( 1.0D+0, 0.0D+0 ) )
181: *
182: * .. Local Scalars ..
183: LOGICAL UPPER, TQUERY, WQUERY
184: INTEGER I, J, K, I1, I2, TD
185: INTEGER LDTB, NB, KB, JB, NT, IINFO
186: COMPLEX*16 PIV
187: * ..
188: * .. External Functions ..
189: LOGICAL LSAME
190: INTEGER ILAENV
191: EXTERNAL LSAME, ILAENV
192: * ..
193: * .. External Subroutines ..
194: EXTERNAL XERBLA, ZCOPY, ZGBTRF, ZGEMM, ZGETRF,
195: $ ZLACPY, ZLASET, ZLASWP, ZTRSM, ZSWAP
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( 'ZSYTRF_AA_2STAGE', -INFO )
222: RETURN
223: END IF
224: *
225: * Answer the query
226: *
227: NB = ILAENV( 1, 'ZSYTRF_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: * .....................................................
275: * Factorize A as U**T*D*U using the upper triangle of A
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+1,I)*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 ZGEMM( 'NoTranspose', 'NoTranspose',
292: $ NB, KB, JB,
293: $ CONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
294: $ A( (I-1)*NB+1, J*NB+1 ), LDA,
295: $ CZERO, 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 ZGEMM( 'NoTranspose', 'NoTranspose',
304: $ NB, KB, JB,
305: $ CONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
306: $ LDTB-1,
307: $ A( (I-2)*NB+1, J*NB+1 ), LDA,
308: $ CZERO, WORK( I*NB+1 ), N )
309: END IF
310: END DO
311: *
312: * Compute T(J,J)
313: *
314: CALL ZLACPY( '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 ZGEMM( 'Transpose', 'NoTranspose',
319: $ KB, KB, (J-1)*NB,
320: $ -CONE, A( 1, J*NB+1 ), LDA,
321: $ WORK( NB+1 ), N,
322: $ CONE, 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 ZGEMM( 'Transpose', 'NoTranspose',
325: $ KB, NB, KB,
326: $ CONE, A( (J-1)*NB+1, J*NB+1 ), LDA,
327: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
328: $ CZERO, WORK( 1 ), N )
329: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
330: $ KB, KB, NB,
331: $ -CONE, WORK( 1 ), N,
332: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
333: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
334: END IF
335: *
336: * Expand T(J,J) into full format
337: *
338: DO I = 1, KB
339: DO K = I+1, KB
340: TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
341: $ = TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
342: END DO
343: END DO
344: IF( J.GT.0 ) THEN
345: c CALL CHEGST( 1, 'Upper', KB,
346: c $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
347: c $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO )
348: CALL ZTRSM( 'L', 'U', 'T', 'N', KB, KB, CONE,
349: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
350: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
351: CALL ZTRSM( 'R', 'U', 'N', 'N', KB, KB, CONE,
352: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
353: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
354: END IF
355: *
356: IF( J.LT.NT-1 ) THEN
357: IF( J.GT.0 ) THEN
358: *
359: * Compute H(J,J)
360: *
361: IF( J.EQ.1 ) THEN
362: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
363: $ KB, KB, KB,
364: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
365: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
366: $ CZERO, WORK( J*NB+1 ), N )
367: ELSE
368: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
369: $ KB, KB, NB+KB,
370: $ CONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
371: $ LDTB-1,
372: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
373: $ CZERO, WORK( J*NB+1 ), N )
374: END IF
375: *
376: * Update with the previous column
377: *
378: CALL ZGEMM( 'Transpose', 'NoTranspose',
379: $ NB, N-(J+1)*NB, J*NB,
380: $ -CONE, WORK( NB+1 ), N,
381: $ A( 1, (J+1)*NB+1 ), LDA,
382: $ CONE, A( J*NB+1, (J+1)*NB+1 ), LDA )
383: END IF
384: *
385: * Copy panel to workspace to call ZGETRF
386: *
387: DO K = 1, NB
388: CALL ZCOPY( N-(J+1)*NB,
389: $ A( J*NB+K, (J+1)*NB+1 ), LDA,
390: $ WORK( 1+(K-1)*N ), 1 )
391: END DO
392: *
393: * Factorize panel
394: *
395: CALL ZGETRF( N-(J+1)*NB, NB,
396: $ WORK, N,
397: $ IPIV( (J+1)*NB+1 ), IINFO )
398: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
399: c INFO = IINFO+(J+1)*NB
400: c END IF
401: *
402: * Copy panel back
403: *
404: DO K = 1, NB
405: CALL ZCOPY( N-(J+1)*NB,
406: $ WORK( 1+(K-1)*N ), 1,
407: $ A( J*NB+K, (J+1)*NB+1 ), LDA )
408: END DO
409: *
410: * Compute T(J+1, J), zero out for GEMM update
411: *
412: KB = MIN(NB, N-(J+1)*NB)
413: CALL ZLASET( 'Full', KB, NB, CZERO, CZERO,
414: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
415: CALL ZLACPY( 'Upper', KB, NB,
416: $ WORK, N,
417: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
418: IF( J.GT.0 ) THEN
419: CALL ZTRSM( 'R', 'U', 'N', 'U', KB, NB, CONE,
420: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
421: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
422: END IF
423: *
424: * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
425: * updates
426: *
427: DO K = 1, NB
428: DO I = 1, KB
429: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
430: $ = TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
431: END DO
432: END DO
433: CALL ZLASET( 'Lower', KB, NB, CZERO, CONE,
434: $ A( J*NB+1, (J+1)*NB+1), LDA )
435: *
436: * Apply pivots to trailing submatrix of A
437: *
438: DO K = 1, KB
439: * > Adjust ipiv
440: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
441: *
442: I1 = (J+1)*NB+K
443: I2 = IPIV( (J+1)*NB+K )
444: IF( I1.NE.I2 ) THEN
445: * > Apply pivots to previous columns of L
446: CALL ZSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
447: $ A( (J+1)*NB+1, I2 ), 1 )
448: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
449: IF( I2.GT.(I1+1) )
450: $ CALL ZSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
451: $ A( I1+1, I2 ), 1 )
452: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
453: IF( I2.LT.N )
454: $ CALL ZSWAP( N-I2, A( I1, I2+1 ), LDA,
455: $ A( I2, I2+1 ), LDA )
456: * > Swap A(I1, I1) with A(I2, I2)
457: PIV = A( I1, I1 )
458: A( I1, I1 ) = A( I2, I2 )
459: A( I2, I2 ) = PIV
460: * > Apply pivots to previous columns of L
461: IF( J.GT.0 ) THEN
462: CALL ZSWAP( J*NB, A( 1, I1 ), 1,
463: $ A( 1, I2 ), 1 )
464: END IF
465: ENDIF
466: END DO
467: END IF
468: END DO
469: ELSE
470: *
471: * .....................................................
472: * Factorize A as L*D*L**T using the lower triangle of A
473: * .....................................................
474: *
475: DO J = 0, NT-1
476: *
477: * Generate Jth column of W and H
478: *
479: KB = MIN(NB, N-J*NB)
480: DO I = 1, J-1
481: IF( I.EQ.1 ) THEN
482: * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
483: IF( I .EQ. (J-1) ) THEN
484: JB = NB+KB
485: ELSE
486: JB = 2*NB
487: END IF
488: CALL ZGEMM( 'NoTranspose', 'Transpose',
489: $ NB, KB, JB,
490: $ CONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
491: $ A( J*NB+1, (I-1)*NB+1 ), LDA,
492: $ CZERO, WORK( I*NB+1 ), N )
493: ELSE
494: * 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)'
495: IF( I .EQ. (J-1) ) THEN
496: JB = 2*NB+KB
497: ELSE
498: JB = 3*NB
499: END IF
500: CALL ZGEMM( 'NoTranspose', 'Transpose',
501: $ NB, KB, JB,
502: $ CONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
503: $ LDTB-1,
504: $ A( J*NB+1, (I-2)*NB+1 ), LDA,
505: $ CZERO, WORK( I*NB+1 ), N )
506: END IF
507: END DO
508: *
509: * Compute T(J,J)
510: *
511: CALL ZLACPY( 'Lower', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
512: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
513: IF( J.GT.1 ) THEN
514: * T(J,J) = L(J,1:J)*H(1:J)
515: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
516: $ KB, KB, (J-1)*NB,
517: $ -CONE, A( J*NB+1, 1 ), LDA,
518: $ WORK( NB+1 ), N,
519: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
520: * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
521: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
522: $ KB, NB, KB,
523: $ CONE, A( J*NB+1, (J-1)*NB+1 ), LDA,
524: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
525: $ CZERO, WORK( 1 ), N )
526: CALL ZGEMM( 'NoTranspose', 'Transpose',
527: $ KB, KB, NB,
528: $ -CONE, WORK( 1 ), N,
529: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
530: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
531: END IF
532: *
533: * Expand T(J,J) into full format
534: *
535: DO I = 1, KB
536: DO K = I+1, KB
537: TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
538: $ = TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
539: END DO
540: END DO
541: IF( J.GT.0 ) THEN
542: c CALL CHEGST( 1, 'Lower', KB,
543: c $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
544: c $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO )
545: CALL ZTRSM( 'L', 'L', 'N', 'N', KB, KB, CONE,
546: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
547: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
548: CALL ZTRSM( 'R', 'L', 'T', 'N', KB, KB, CONE,
549: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
550: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
551: END IF
552: *
553: * Symmetrize T(J,J)
554: *
555: DO I = 1, KB
556: DO K = I+1, KB
557: TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
558: $ = TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
559: END DO
560: END DO
561: *
562: IF( J.LT.NT-1 ) THEN
563: IF( J.GT.0 ) THEN
564: *
565: * Compute H(J,J)
566: *
567: IF( J.EQ.1 ) THEN
568: CALL ZGEMM( 'NoTranspose', 'Transpose',
569: $ KB, KB, KB,
570: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
571: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
572: $ CZERO, WORK( J*NB+1 ), N )
573: ELSE
574: CALL ZGEMM( 'NoTranspose', 'Transpose',
575: $ KB, KB, NB+KB,
576: $ CONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
577: $ LDTB-1,
578: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
579: $ CZERO, WORK( J*NB+1 ), N )
580: END IF
581: *
582: * Update with the previous column
583: *
584: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
585: $ N-(J+1)*NB, NB, J*NB,
586: $ -CONE, A( (J+1)*NB+1, 1 ), LDA,
587: $ WORK( NB+1 ), N,
588: $ CONE, A( (J+1)*NB+1, J*NB+1 ), LDA )
589: END IF
590: *
591: * Factorize panel
592: *
593: CALL ZGETRF( N-(J+1)*NB, NB,
594: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
595: $ IPIV( (J+1)*NB+1 ), IINFO )
596: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
597: c INFO = IINFO+(J+1)*NB
598: c END IF
599: *
600: * Compute T(J+1, J), zero out for GEMM update
601: *
602: KB = MIN(NB, N-(J+1)*NB)
603: CALL ZLASET( 'Full', KB, NB, CZERO, CZERO,
604: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
605: CALL ZLACPY( 'Upper', KB, NB,
606: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
607: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
608: IF( J.GT.0 ) THEN
609: CALL ZTRSM( 'R', 'L', 'T', 'U', KB, NB, CONE,
610: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
611: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
612: END IF
613: *
614: * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
615: * updates
616: *
617: DO K = 1, NB
618: DO I = 1, KB
619: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB ) =
620: $ TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
621: END DO
622: END DO
623: CALL ZLASET( 'Upper', KB, NB, CZERO, CONE,
624: $ A( (J+1)*NB+1, J*NB+1 ), LDA )
625: *
626: * Apply pivots to trailing submatrix of A
627: *
628: DO K = 1, KB
629: * > Adjust ipiv
630: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
631: *
632: I1 = (J+1)*NB+K
633: I2 = IPIV( (J+1)*NB+K )
634: IF( I1.NE.I2 ) THEN
635: * > Apply pivots to previous columns of L
636: CALL ZSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
637: $ A( I2, (J+1)*NB+1 ), LDA )
638: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
639: IF( I2.GT.(I1+1) )
640: $ CALL ZSWAP( I2-I1-1, A( I1+1, I1 ), 1,
641: $ A( I2, I1+1 ), LDA )
642: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
643: IF( I2.LT.N )
644: $ CALL ZSWAP( N-I2, A( I2+1, I1 ), 1,
645: $ A( I2+1, I2 ), 1 )
646: * > Swap A(I1, I1) with A(I2, I2)
647: PIV = A( I1, I1 )
648: A( I1, I1 ) = A( I2, I2 )
649: A( I2, I2 ) = PIV
650: * > Apply pivots to previous columns of L
651: IF( J.GT.0 ) THEN
652: CALL ZSWAP( J*NB, A( I1, 1 ), LDA,
653: $ A( I2, 1 ), LDA )
654: END IF
655: ENDIF
656: END DO
657: *
658: * Apply pivots to previous columns of L
659: *
660: c CALL ZLASWP( J*NB, A( 1, 1 ), LDA,
661: c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
662: END IF
663: END DO
664: END IF
665: *
666: * Factor the band matrix
667: CALL ZGBTRF( N, N, NB, NB, TB, LDTB, IPIV2, INFO )
668: *
669: RETURN
670: *
671: * End of ZSYTRF_AA_2STAGE
672: *
673: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsytrf_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