1: *> \brief \b ZHETRF_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 ZHETRF_AA_2STAGE + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHETRF_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: *> ZHETRF_AA_2STAGE computes the factorization of a double hermitian matrix A
39: *> using the Aasen's algorithm. The form of the factorization is
40: *>
41: *> A = U**H*T*U or A = L*T*L**H
42: *>
43: *> where U (or L) is a product of permutation and unit upper (lower)
44: *> triangular matrices, and T is a hermitian 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 ZHETRF_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 ZERO, ONE
179: PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ),
180: $ ONE = ( 1.0E+0, 0.0E+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, ZLACGV, ZLACPY,
195: $ ZLASET, ZGBTRF, ZGEMM, ZGETRF,
196: $ ZHEGST, ZSWAP, ZTRSM
197: * ..
198: * .. Intrinsic Functions ..
199: INTRINSIC DCONJG, 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( 'ZHETRF_AA_2STAGE', -INFO )
223: RETURN
224: END IF
225: *
226: * Answer the query
227: *
228: NB = ILAENV( 1, 'ZHETRF_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 U**H*D*U 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: $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
295: $ A( (I-1)*NB+1, J*NB+1 ), LDA,
296: $ ZERO, 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: $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
307: $ LDTB-1,
308: $ A( (I-2)*NB+1, J*NB+1 ), LDA,
309: $ ZERO, 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( 'Conjugate transpose', 'NoTranspose',
320: $ KB, KB, (J-1)*NB,
321: $ -ONE, A( 1, J*NB+1 ), LDA,
322: $ WORK( NB+1 ), N,
323: $ ONE, 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( 'Conjugate transpose', 'NoTranspose',
326: $ KB, NB, KB,
327: $ ONE, A( (J-1)*NB+1, J*NB+1 ), LDA,
328: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
329: $ ZERO, WORK( 1 ), N )
330: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
331: $ KB, KB, NB,
332: $ -ONE, WORK( 1 ), N,
333: $ A( (J-2)*NB+1, J*NB+1 ), LDA,
334: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
335: END IF
336: IF( J.GT.0 ) THEN
337: CALL ZHEGST( 1, 'Upper', KB,
338: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
339: $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO )
340: END IF
341: *
342: * Expand T(J,J) into full format
343: *
344: DO I = 1, KB
345: TB( TD+1 + (J*NB+I-1)*LDTB )
346: $ = REAL( TB( TD+1 + (J*NB+I-1)*LDTB ) )
347: DO K = I+1, KB
348: TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB )
349: $ = DCONJG( 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 ZGEMM( '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 ZGEMM( '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 ZGEMM( 'Conjugate 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 ZGETRF
383: *
384: DO K = 1, NB
385: CALL ZCOPY( 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 ZGETRF( 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: *
403: * Copy only L-factor
404: *
405: CALL ZCOPY( N-K-(J+1)*NB,
406: $ WORK( K+1+(K-1)*N ), 1,
407: $ A( J*NB+K, (J+1)*NB+K+1 ), LDA )
408: *
409: * Transpose U-factor to be copied back into T(J+1, J)
410: *
411: CALL ZLACGV( K, WORK( 1+(K-1)*N ), 1 )
412: END DO
413: *
414: * Compute T(J+1, J), zero out for GEMM update
415: *
416: KB = MIN(NB, N-(J+1)*NB)
417: CALL ZLASET( 'Full', KB, NB, ZERO, ZERO,
418: $ TB( TD+NB+1 + (J*NB)*LDTB) , LDTB-1 )
419: CALL ZLACPY( 'Upper', KB, NB,
420: $ WORK, N,
421: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
422: IF( J.GT.0 ) THEN
423: CALL ZTRSM( 'R', 'U', 'N', 'U', KB, NB, ONE,
424: $ A( (J-1)*NB+1, J*NB+1 ), LDA,
425: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
426: END IF
427: *
428: * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
429: * updates
430: *
431: DO K = 1, NB
432: DO I = 1, KB
433: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
434: $ = DCONJG( TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB ) )
435: END DO
436: END DO
437: CALL ZLASET( 'Lower', KB, NB, ZERO, ONE,
438: $ A( J*NB+1, (J+1)*NB+1), LDA )
439: *
440: * Apply pivots to trailing submatrix of A
441: *
442: DO K = 1, KB
443: * > Adjust ipiv
444: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
445: *
446: I1 = (J+1)*NB+K
447: I2 = IPIV( (J+1)*NB+K )
448: IF( I1.NE.I2 ) THEN
449: * > Apply pivots to previous columns of L
450: CALL ZSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
451: $ A( (J+1)*NB+1, I2 ), 1 )
452: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
453: IF( I2.GT.(I1+1) ) THEN
454: CALL ZSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
455: $ A( I1+1, I2 ), 1 )
456: CALL ZLACGV( I2-I1-1, A( I1+1, I2 ), 1 )
457: END IF
458: CALL ZLACGV( I2-I1, A( I1, I1+1 ), LDA )
459: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
460: IF( I2.LT.N )
461: $ CALL ZSWAP( N-I2, A( I1, I2+1 ), LDA,
462: $ A( I2, I2+1 ), LDA )
463: * > Swap A(I1, I1) with A(I2, I2)
464: PIV = A( I1, I1 )
465: A( I1, I1 ) = A( I2, I2 )
466: A( I2, I2 ) = PIV
467: * > Apply pivots to previous columns of L
468: IF( J.GT.0 ) THEN
469: CALL ZSWAP( J*NB, A( 1, I1 ), 1,
470: $ A( 1, I2 ), 1 )
471: END IF
472: ENDIF
473: END DO
474: END IF
475: END DO
476: ELSE
477: *
478: * .....................................................
479: * Factorize A as L*D*L**H using the lower triangle of A
480: * .....................................................
481: *
482: DO J = 0, NT-1
483: *
484: * Generate Jth column of W and H
485: *
486: KB = MIN(NB, N-J*NB)
487: DO I = 1, J-1
488: IF( I.EQ.1 ) THEN
489: * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
490: IF( I .EQ. (J-1) ) THEN
491: JB = NB+KB
492: ELSE
493: JB = 2*NB
494: END IF
495: CALL ZGEMM( 'NoTranspose', 'Conjugate transpose',
496: $ NB, KB, JB,
497: $ ONE, TB( TD+1 + (I*NB)*LDTB ), LDTB-1,
498: $ A( J*NB+1, (I-1)*NB+1 ), LDA,
499: $ ZERO, WORK( I*NB+1 ), N )
500: ELSE
501: * 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)'
502: IF( I .EQ. (J-1) ) THEN
503: JB = 2*NB+KB
504: ELSE
505: JB = 3*NB
506: END IF
507: CALL ZGEMM( 'NoTranspose', 'Conjugate transpose',
508: $ NB, KB, JB,
509: $ ONE, TB( TD+NB+1 + ((I-1)*NB)*LDTB ),
510: $ LDTB-1,
511: $ A( J*NB+1, (I-2)*NB+1 ), LDA,
512: $ ZERO, WORK( I*NB+1 ), N )
513: END IF
514: END DO
515: *
516: * Compute T(J,J)
517: *
518: CALL ZLACPY( 'Lower', KB, KB, A( J*NB+1, J*NB+1 ), LDA,
519: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
520: IF( J.GT.1 ) THEN
521: * T(J,J) = L(J,1:J)*H(1:J)
522: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
523: $ KB, KB, (J-1)*NB,
524: $ -ONE, A( J*NB+1, 1 ), LDA,
525: $ WORK( NB+1 ), N,
526: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
527: * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
528: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
529: $ KB, NB, KB,
530: $ ONE, A( J*NB+1, (J-1)*NB+1 ), LDA,
531: $ TB( TD+NB+1 + ((J-1)*NB)*LDTB ), LDTB-1,
532: $ ZERO, WORK( 1 ), N )
533: CALL ZGEMM( 'NoTranspose', 'Conjugate transpose',
534: $ KB, KB, NB,
535: $ -ONE, WORK( 1 ), N,
536: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
537: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1 )
538: END IF
539: IF( J.GT.0 ) THEN
540: CALL ZHEGST( 1, 'Lower', KB,
541: $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
542: $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO )
543: END IF
544: *
545: * Expand T(J,J) into full format
546: *
547: DO I = 1, KB
548: TB( TD+1 + (J*NB+I-1)*LDTB )
549: $ = REAL( TB( TD+1 + (J*NB+I-1)*LDTB ) )
550: DO K = I+1, KB
551: TB( TD-(K-(I+1)) + (J*NB+K-1)*LDTB )
552: $ = DCONJG( TB( TD+(K-I)+1 + (J*NB+I-1)*LDTB ) )
553: END DO
554: END DO
555: *
556: IF( J.LT.NT-1 ) THEN
557: IF( J.GT.0 ) THEN
558: *
559: * Compute H(J,J)
560: *
561: IF( J.EQ.1 ) THEN
562: CALL ZGEMM( 'NoTranspose', 'Conjugate transpose',
563: $ KB, KB, KB,
564: $ ONE, TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
565: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
566: $ ZERO, WORK( J*NB+1 ), N )
567: ELSE
568: CALL ZGEMM( 'NoTranspose', 'Conjugate transpose',
569: $ KB, KB, NB+KB,
570: $ ONE, TB( TD+NB+1 + ((J-1)*NB)*LDTB ),
571: $ LDTB-1,
572: $ A( J*NB+1, (J-2)*NB+1 ), LDA,
573: $ ZERO, WORK( J*NB+1 ), N )
574: END IF
575: *
576: * Update with the previous column
577: *
578: CALL ZGEMM( 'NoTranspose', 'NoTranspose',
579: $ N-(J+1)*NB, NB, J*NB,
580: $ -ONE, A( (J+1)*NB+1, 1 ), LDA,
581: $ WORK( NB+1 ), N,
582: $ ONE, A( (J+1)*NB+1, J*NB+1 ), LDA )
583: END IF
584: *
585: * Factorize panel
586: *
587: CALL ZGETRF( N-(J+1)*NB, NB,
588: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
589: $ IPIV( (J+1)*NB+1 ), IINFO )
590: c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
591: c INFO = IINFO+(J+1)*NB
592: c END IF
593: *
594: * Compute T(J+1, J), zero out for GEMM update
595: *
596: KB = MIN(NB, N-(J+1)*NB)
597: CALL ZLASET( 'Full', KB, NB, ZERO, ZERO,
598: $ TB( TD+NB+1 + (J*NB)*LDTB) , LDTB-1 )
599: CALL ZLACPY( 'Upper', KB, NB,
600: $ A( (J+1)*NB+1, J*NB+1 ), LDA,
601: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
602: IF( J.GT.0 ) THEN
603: CALL ZTRSM( 'R', 'L', 'C', 'U', KB, NB, ONE,
604: $ A( J*NB+1, (J-1)*NB+1 ), LDA,
605: $ TB( TD+NB+1 + (J*NB)*LDTB ), LDTB-1 )
606: END IF
607: *
608: * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
609: * updates
610: *
611: DO K = 1, NB
612: DO I = 1, KB
613: TB( TD-NB+K-I+1 + (J*NB+NB+I-1)*LDTB )
614: $ = DCONJG( TB( TD+NB+I-K+1 + (J*NB+K-1)*LDTB ) )
615: END DO
616: END DO
617: CALL ZLASET( 'Upper', KB, NB, ZERO, ONE,
618: $ A( (J+1)*NB+1, J*NB+1), LDA )
619: *
620: * Apply pivots to trailing submatrix of A
621: *
622: DO K = 1, KB
623: * > Adjust ipiv
624: IPIV( (J+1)*NB+K ) = IPIV( (J+1)*NB+K ) + (J+1)*NB
625: *
626: I1 = (J+1)*NB+K
627: I2 = IPIV( (J+1)*NB+K )
628: IF( I1.NE.I2 ) THEN
629: * > Apply pivots to previous columns of L
630: CALL ZSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
631: $ A( I2, (J+1)*NB+1 ), LDA )
632: * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
633: IF( I2.GT.(I1+1) ) THEN
634: CALL ZSWAP( I2-I1-1, A( I1+1, I1 ), 1,
635: $ A( I2, I1+1 ), LDA )
636: CALL ZLACGV( I2-I1-1, A( I2, I1+1 ), LDA )
637: END IF
638: CALL ZLACGV( I2-I1, A( I1+1, I1 ), 1 )
639: * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
640: IF( I2.LT.N )
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: RETURN
667: *
668: * End of ZHETRF_AA_2STAGE
669: *
670: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetrf_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