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