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/csytrf_aa_2stage.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrf_aa_2stage.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrf_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*U**T 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: *> The size of the array TB. LTB >= 4*N, internally
97: *> used to select NB such that LTB >= (3*NB+1)*N.
98: *>
99: *> If LTB = -1, then a workspace query is assumed; the
100: *> routine only calculates the optimal size of LTB,
101: *> returns this value as the first entry of TB, and
102: *> no error message related to LTB is issued by XERBLA.
103: *> \endverbatim
104: *>
105: *> \param[out] IPIV
106: *> \verbatim
107: *> IPIV is INTEGER array, dimension (N)
108: *> On exit, it contains the details of the interchanges, i.e.,
109: *> the row and column k of A were interchanged with the
110: *> row and column IPIV(k).
111: *> \endverbatim
112: *>
113: *> \param[out] IPIV2
114: *> \verbatim
115: *> IPIV is INTEGER array, dimension (N)
116: *> On exit, it contains the details of the interchanges, i.e.,
117: *> the row and column k of T were interchanged with the
118: *> row and column IPIV(k).
119: *> \endverbatim
120: *>
121: *> \param[out] WORK
122: *> \verbatim
123: *> WORK is COMPLEX*16 workspace of size LWORK
124: *> \endverbatim
125: *>
126: *> \param[in] LWORK
127: *> \verbatim
128: *> The size of WORK. LWORK >= N, internally used to select NB
129: *> such that LWORK >= N*NB.
130: *>
131: *> If LWORK = -1, then a workspace query is assumed; the
132: *> routine only calculates the optimal size of the WORK array,
133: *> returns this value as the first entry of the WORK array, and
134: *> no error message related to LWORK is issued by XERBLA.
135: *> \endverbatim
136: *>
137: *> \param[out] INFO
138: *> \verbatim
139: *> INFO is INTEGER
140: *> = 0: successful exit
141: *> < 0: if INFO = -i, the i-th argument had an illegal value.
142: *> > 0: if INFO = i, band LU factorization failed on i-th column
143: *> \endverbatim
144: *
145: * Authors:
146: * ========
147: *
148: *> \author Univ. of Tennessee
149: *> \author Univ. of California Berkeley
150: *> \author Univ. of Colorado Denver
151: *> \author NAG Ltd.
152: *
153: *> \date November 2017
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 (version 3.8.0) --
162: * -- LAPACK is a software package provided by Univ. of Tennessee, --
163: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
164: * November 2017
165: *
166: IMPLICIT NONE
167: *
168: * .. Scalar Arguments ..
169: CHARACTER UPLO
170: INTEGER N, LDA, LTB, LWORK, INFO
171: * ..
172: * .. Array Arguments ..
173: INTEGER IPIV( * ), IPIV2( * )
174: COMPLEX*16 A( LDA, * ), TB( * ), WORK( * )
175: * ..
176: *
177: * =====================================================================
178: * .. Parameters ..
179: COMPLEX*16 CZERO, CONE
180: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
181: $ CONE = ( 1.0D+0, 0.0D+0 ) )
182: *
183: * .. Local Scalars ..
184: LOGICAL UPPER, TQUERY, WQUERY
185: INTEGER I, J, K, I1, I2, TD
186: INTEGER LDTB, NB, KB, JB, NT, IINFO
187: COMPLEX*16 PIV
188: * ..
189: * .. External Functions ..
190: LOGICAL LSAME
191: INTEGER ILAENV
192: EXTERNAL LSAME, ILAENV
193: * ..
194: * .. External Subroutines ..
195: EXTERNAL XERBLA, ZCOPY, ZGBTRF, ZGEMM, ZGETRF,
196: $ ZLACPY, ZLASET, ZLASWP, ZTRSM, ZSWAP
197: * ..
198: * .. Intrinsic Functions ..
199: INTRINSIC MIN, MAX
200: * ..
201: * .. Executable Statements ..
202: *
203: * Test the input parameters.
204: *
205: INFO = 0
206: UPPER = LSAME( UPLO, 'U' )
207: WQUERY = ( LWORK.EQ.-1 )
208: TQUERY = ( LTB.EQ.-1 )
209: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
210: INFO = -1
211: ELSE IF( N.LT.0 ) THEN
212: INFO = -2
213: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
214: INFO = -4
215: ELSE IF ( LTB .LT. 4*N .AND. .NOT.TQUERY ) THEN
216: INFO = -6
217: ELSE IF ( LWORK .LT. N .AND. .NOT.WQUERY ) THEN
218: INFO = -10
219: END IF
220: *
221: IF( INFO.NE.0 ) THEN
222: CALL XERBLA( 'ZSYTRF_AA_2STAGE', -INFO )
223: RETURN
224: END IF
225: *
226: * Answer the query
227: *
228: NB = ILAENV( 1, 'ZSYTRF_AA_2STAGE', UPLO, N, -1, -1, -1 )
229: IF( INFO.EQ.0 ) THEN
230: IF( TQUERY ) THEN
231: TB( 1 ) = (3*NB+1)*N
232: END IF
233: IF( WQUERY ) THEN
234: WORK( 1 ) = N*NB
235: END IF
236: END IF
237: IF( TQUERY .OR. WQUERY ) THEN
238: RETURN
239: END IF
240: *
241: * Quick return
242: *
243: IF ( N.EQ.0 ) THEN
244: RETURN
245: ENDIF
246: *
247: * Determine the number of the block size
248: *
249: LDTB = LTB/N
250: IF( LDTB .LT. 3*NB+1 ) THEN
251: NB = (LDTB-1)/3
252: END IF
253: IF( LWORK .LT. NB*N ) THEN
254: NB = LWORK/N
255: END IF
256: *
257: * Determine the number of the block columns
258: *
259: NT = (N+NB-1)/NB
260: TD = 2*NB
261: KB = MIN(NB, N)
262: *
263: * Initialize vectors/matrices
264: *
265: DO J = 1, KB
266: IPIV( J ) = J
267: END DO
268: *
269: * Save NB
270: *
271: TB( 1 ) = NB
272: *
273: IF( UPPER ) THEN
274: *
275: * .....................................................
276: * Factorize A as L*D*L**T using the upper triangle of A
277: * .....................................................
278: *
279: DO J = 0, NT-1
280: *
281: * Generate Jth column of W and H
282: *
283: KB = MIN(NB, N-J*NB)
284: DO I = 1, J-1
285: IF( I.EQ.1 ) THEN
286: * H(I,J) = T(I,I)*U(I,J) + T(I+1,I)*U(I+1,J)
287: IF( I .EQ. (J-1) ) THEN
288: JB = NB+KB
289: ELSE
290: JB = 2*NB
291: END IF
292: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
293: $ NB, KB, JB,
294: $ CONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
295: $ A( (I-1)*NB+1, J*NB+1 ), LDA,
296: $ CZERO, WORK( I*NB+1 ), N )
297: ELSE
298: * 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)
299: IF( I .EQ. (J-1) ) THEN
300: JB = 2*NB+KB
301: ELSE
302: JB = 3*NB
303: END IF
304: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
305: $ NB, KB, JB,
306: $ CONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
307: $ LDTB-1,
308: $ A( (I-2)*NB+1, J*NB+1 ), LDA,
309: $ CZERO, WORK( I*NB+1 ), N )
310: END IF
311: END DO
312: *
313: * Compute T(J,J)
314: *
315: CALL ZLACPY( 'Upper', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
316: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
317: IF( J.GT.1 ) THEN
318: * T(J,J) = U(1:J,J)'*H(1:J)
319: CALL ZGEMM( 'Transpose', 'NoTranspose',
320: $ KB, KB, (J-1)*NB,
321: $ -CONE, A( 1, J*NB+1 ), LDA,
322: $ WORK( NB+1 ), N,
323: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
324: * T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J)
325: CALL ZGEMM( 'Transpose', 'NoTranspose',
326: $ KB, NB, KB,
327: $ CONE, A( (J-1)*NB+1, J*NB+1 ), LDA,
328: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
329: $ CZERO, WORK( 1 ), N )
330: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
331: $ KB, KB, NB,
332: $ -CONE, WORK( 1 ), N,
333: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
334: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
335: END IF
336: *
337: * Expand T(J,J) into full format
338: *
339: DO I = 1, KB
340: DO K = I+1, KB
341: TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
342: $ = TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
343: END DO
344: END DO
345: IF( J.GT.0 ) THEN
346: c CALL CHEGST( 1, 'Upper', KB,
347: c $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
348: c $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO )
349: CALL ZTRSM( 'L', 'U', 'T', 'N', KB, KB, CONE,
350: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
351: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
352: CALL ZTRSM( 'R', 'U', 'N', 'N', KB, KB, CONE,
353: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
354: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
355: END IF
356: *
357: IF( J.LT.NT-1 ) THEN
358: IF( J.GT.0 ) THEN
359: *
360: * Compute H(J,J)
361: *
362: IF( J.EQ.1 ) THEN
363: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
364: $ KB, KB, KB,
365: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
366: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
367: $ CZERO, WORK( J*NB+1 ), N )
368: ELSE
369: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
370: $ KB, KB, NB+KB,
371: $ CONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
372: $ LDTB-1,
373: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
374: $ CZERO, WORK( J*NB+1 ), N )
375: END IF
376: *
377: * Update with the previous column
378: *
379: CALL ZGEMM( 'Transpose', 'NoTranspose',
380: $ NB, N-(J+1)*NB, J*NB,
381: $ -CONE, WORK( NB+1 ), N,
382: $ A( 1, (J+1)*NB+1 ), LDA,
383: $ CONE, A( J*NB+1, (J+1)*NB+1 ), LDA )
384: END IF
385: *
386: * Copy panel to workspace to call ZGETRF
387: *
388: DO K = 1, NB
389: CALL ZCOPY( N-(J+1)*NB,
390: $ A( J*NB+K, (J+1)*NB+1 ), LDA,
391: $ WORK( 1+(K-1)*N ), 1 )
392: END DO
393: *
394: * Factorize panel
395: *
396: CALL ZGETRF( N-(J+1)*NB, NB,
397: $ WORK, N,
398: $ IPIV( (J+1)*NB+1 ), IINFO )
399: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
400: c INFO = IINFO+(J+1)*NB
401: c END IF
402: *
403: * Copy panel back
404: *
405: DO K = 1, NB
406: CALL ZCOPY( N-(J+1)*NB,
407: $ WORK( 1+(K-1)*N ), 1,
408: $ A( J*NB+K, (J+1)*NB+1 ), LDA )
409: END DO
410: *
411: * Compute T(J+1, J), zero out for GEMM update
412: *
413: KB = MIN(NB, N-(J+1)*NB)
414: CALL ZLASET( 'Full', KB, NB, CZERO, CZERO,
415: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
416: CALL ZLACPY( 'Upper', KB, NB,
417: $ WORK, N,
418: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
419: IF( J.GT.0 ) THEN
420: CALL ZTRSM( 'R', 'U', 'N', 'U', KB, NB, CONE,
421: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
422: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
423: END IF
424: *
425: * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
426: * updates
427: *
428: DO K = 1, NB
429: DO I = 1, KB
430: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
431: $ = TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
432: END DO
433: END DO
434: CALL ZLASET( 'Lower', KB, NB, CZERO, CONE,
435: $ A( J*NB+1, (J+1)*NB+1), LDA )
436: *
437: * Apply pivots to trailing submatrix of A
438: *
439: DO K = 1, KB
440: * > Adjust ipiv
441: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
442: *
443: I1 = (J+1)*NB+K
444: I2 = IPIV( (J+1)*NB+K )
445: IF( I1.NE.I2 ) THEN
446: * > Apply pivots to previous columns of L
447: CALL ZSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
448: $ A( (J+1)*NB+1, I2 ), 1 )
449: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
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: CALL ZSWAP( N-I2, A( I1, I2+1 ), LDA,
454: $ A( I2, I2+1 ), LDA )
455: * > Swap A(I1, I1) with A(I2, I2)
456: PIV = A( I1, I1 )
457: A( I1, I1 ) = A( I2, I2 )
458: A( I2, I2 ) = PIV
459: * > Apply pivots to previous columns of L
460: IF( J.GT.0 ) THEN
461: CALL ZSWAP( J*NB, A( 1, I1 ), 1,
462: $ A( 1, I2 ), 1 )
463: END IF
464: ENDIF
465: END DO
466: END IF
467: END DO
468: ELSE
469: *
470: * .....................................................
471: * Factorize A as L*D*L**T using the lower triangle of A
472: * .....................................................
473: *
474: DO J = 0, NT-1
475: *
476: * Generate Jth column of W and H
477: *
478: KB = MIN(NB, N-J*NB)
479: DO I = 1, J-1
480: IF( I.EQ.1 ) THEN
481: * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
482: IF( I .EQ. (J-1) ) THEN
483: JB = NB+KB
484: ELSE
485: JB = 2*NB
486: END IF
487: CALL ZGEMM( 'NoTranspose', 'Transpose',
488: $ NB, KB, JB,
489: $ CONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
490: $ A( J*NB+1, (I-1)*NB+1 ), LDA,
491: $ CZERO, WORK( I*NB+1 ), N )
492: ELSE
493: * 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)'
494: IF( I .EQ. (J-1) ) THEN
495: JB = 2*NB+KB
496: ELSE
497: JB = 3*NB
498: END IF
499: CALL ZGEMM( 'NoTranspose', 'Transpose',
500: $ NB, KB, JB,
501: $ CONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
502: $ LDTB-1,
503: $ A( J*NB+1, (I-2)*NB+1 ), LDA,
504: $ CZERO, WORK( I*NB+1 ), N )
505: END IF
506: END DO
507: *
508: * Compute T(J,J)
509: *
510: CALL ZLACPY( 'Lower', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
511: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
512: IF( J.GT.1 ) THEN
513: * T(J,J) = L(J,1:J)*H(1:J)
514: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
515: $ KB, KB, (J-1)*NB,
516: $ -CONE, A( J*NB+1, 1 ), LDA,
517: $ WORK( NB+1 ), N,
518: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
519: * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
520: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
521: $ KB, NB, KB,
522: $ CONE, A( J*NB+1, (J-1)*NB+1 ), LDA,
523: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
524: $ CZERO, WORK( 1 ), N )
525: CALL ZGEMM( 'NoTranspose', 'Transpose',
526: $ KB, KB, NB,
527: $ -CONE, WORK( 1 ), N,
528: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
529: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
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: IF( J.GT.0 ) THEN
541: c CALL CHEGST( 1, 'Lower', KB,
542: c $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
543: c $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO )
544: CALL ZTRSM( 'L', 'L', 'N', 'N', KB, KB, CONE,
545: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
546: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
547: CALL ZTRSM( 'R', 'L', 'T', 'N', KB, KB, CONE,
548: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
549: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
550: END IF
551: *
552: * Symmetrize T(J,J)
553: *
554: DO I = 1, KB
555: DO K = I+1, KB
556: TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
557: $ = TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
558: END DO
559: END DO
560: *
561: IF( J.LT.NT-1 ) THEN
562: IF( J.GT.0 ) THEN
563: *
564: * Compute H(J,J)
565: *
566: IF( J.EQ.1 ) THEN
567: CALL ZGEMM( 'NoTranspose', 'Transpose',
568: $ KB, KB, KB,
569: $ CONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
570: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
571: $ CZERO, WORK( J*NB+1 ), N )
572: ELSE
573: CALL ZGEMM( 'NoTranspose', 'Transpose',
574: $ KB, KB, NB+KB,
575: $ CONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
576: $ LDTB-1,
577: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
578: $ CZERO, WORK( J*NB+1 ), N )
579: END IF
580: *
581: * Update with the previous column
582: *
583: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
584: $ N-(J+1)*NB, NB, J*NB,
585: $ -CONE, A( (J+1)*NB+1, 1 ), LDA,
586: $ WORK( NB+1 ), N,
587: $ CONE, A( (J+1)*NB+1, J*NB+1 ), LDA )
588: END IF
589: *
590: * Factorize panel
591: *
592: CALL ZGETRF( N-(J+1)*NB, NB,
593: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
594: $ IPIV( (J+1)*NB+1 ), IINFO )
595: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
596: c INFO = IINFO+(J+1)*NB
597: c END IF
598: *
599: * Compute T(J+1, J), zero out for GEMM update
600: *
601: KB = MIN(NB, N-(J+1)*NB)
602: CALL ZLASET( 'Full', KB, NB, CZERO, CZERO,
603: $ TB( TD+NB+1 + (J*NB)*LDTB), LDTB-1 )
604: CALL ZLACPY( 'Upper', KB, NB,
605: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
606: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
607: IF( J.GT.0 ) THEN
608: CALL ZTRSM( 'R', 'L', 'T', 'U', KB, NB, CONE,
609: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
610: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
611: END IF
612: *
613: * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
614: * updates
615: *
616: DO K = 1, NB
617: DO I = 1, KB
618: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB ) =
619: $ TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB )
620: END DO
621: END DO
622: CALL ZLASET( 'Upper', KB, NB, CZERO, CONE,
623: $ A( (J+1)*NB+1, J*NB+1 ), LDA )
624: *
625: * Apply pivots to trailing submatrix of A
626: *
627: DO K = 1, KB
628: * > Adjust ipiv
629: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
630: *
631: I1 = (J+1)*NB+K
632: I2 = IPIV( (J+1)*NB+K )
633: IF( I1.NE.I2 ) THEN
634: * > Apply pivots to previous columns of L
635: CALL ZSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
636: $ A( I2, (J+1)*NB+1 ), LDA )
637: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
638: CALL ZSWAP( I2-I1-1, A( I1+1, I1 ), 1,
639: $ A( I2, I1+1 ), LDA )
640: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
641: CALL ZSWAP( N-I2, A( I2+1, I1 ), 1,
642: $ A( I2+1, I2 ), 1 )
643: * > Swap A(I1, I1) with A(I2, I2)
644: PIV = A( I1, I1 )
645: A( I1, I1 ) = A( I2, I2 )
646: A( I2, I2 ) = PIV
647: * > Apply pivots to previous columns of L
648: IF( J.GT.0 ) THEN
649: CALL ZSWAP( J*NB, A( I1, 1 ), LDA,
650: $ A( I2, 1 ), LDA )
651: END IF
652: ENDIF
653: END DO
654: *
655: * Apply pivots to previous columns of L
656: *
657: c CALL ZLASWP( J*NB, A( 1, 1 ), LDA,
658: c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
659: END IF
660: END DO
661: END IF
662: *
663: * Factor the band matrix
664: CALL ZGBTRF( N, N, NB, NB, TB, LDTB, IPIV2, INFO )
665: *
666: * End of ZSYTRF_AA_2STAGE
667: *
668: 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