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