1: *> \brief \b ZLAHEF_AA
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZLAHEF_AA + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlahef_aa.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlahef_aa.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlahef_aa.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAHEF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
22: * H, LDH, WORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER UPLO
26: * INTEGER J1, M, NB, LDA, LDH, INFO
27: * ..
28: * .. Array Arguments ..
29: * INTEGER IPIV( * )
30: * COMPLEX*16 A( LDA, * ), H( LDH, * ), WORK( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DLAHEF_AA factorizes a panel of a complex hermitian matrix A using
40: *> the Aasen's algorithm. The panel consists of a set of NB rows of A
41: *> when UPLO is U, or a set of NB columns when UPLO is L.
42: *>
43: *> In order to factorize the panel, the Aasen's algorithm requires the
44: *> last row, or column, of the previous panel. The first row, or column,
45: *> of A is set to be the first row, or column, of an identity matrix,
46: *> which is used to factorize the first panel.
47: *>
48: *> The resulting J-th row of U, or J-th column of L, is stored in the
49: *> (J-1)-th row, or column, of A (without the unit diagonals), while
50: *> the diagonal and subdiagonal of A are overwritten by those of T.
51: *>
52: *> \endverbatim
53: *
54: * Arguments:
55: * ==========
56: *
57: *> \param[in] UPLO
58: *> \verbatim
59: *> UPLO is CHARACTER*1
60: *> = 'U': Upper triangle of A is stored;
61: *> = 'L': Lower triangle of A is stored.
62: *> \endverbatim
63: *>
64: *> \param[in] J1
65: *> \verbatim
66: *> J1 is INTEGER
67: *> The location of the first row, or column, of the panel
68: *> within the submatrix of A, passed to this routine, e.g.,
69: *> when called by ZHETRF_AA, for the first panel, J1 is 1,
70: *> while for the remaining panels, J1 is 2.
71: *> \endverbatim
72: *>
73: *> \param[in] M
74: *> \verbatim
75: *> M is INTEGER
76: *> The dimension of the submatrix. M >= 0.
77: *> \endverbatim
78: *>
79: *> \param[in] NB
80: *> \verbatim
81: *> NB is INTEGER
82: *> The dimension of the panel to be facotorized.
83: *> \endverbatim
84: *>
85: *> \param[in,out] A
86: *> \verbatim
87: *> A is COMPLEX*16 array, dimension (LDA,M) for
88: *> the first panel, while dimension (LDA,M+1) for the
89: *> remaining panels.
90: *>
91: *> On entry, A contains the last row, or column, of
92: *> the previous panel, and the trailing submatrix of A
93: *> to be factorized, except for the first panel, only
94: *> the panel is passed.
95: *>
96: *> On exit, the leading panel is factorized.
97: *> \endverbatim
98: *>
99: *> \param[in] LDA
100: *> \verbatim
101: *> LDA is INTEGER
102: *> The leading dimension of the array A. LDA >= max(1,N).
103: *> \endverbatim
104: *>
105: *> \param[out] IPIV
106: *> \verbatim
107: *> IPIV is INTEGER array, dimension (N)
108: *> Details of the row and column interchanges,
109: *> the row and column k were interchanged with the row and
110: *> column IPIV(k).
111: *> \endverbatim
112: *>
113: *> \param[in,out] H
114: *> \verbatim
115: *> H is COMPLEX*16 workspace, dimension (LDH,NB).
116: *>
117: *> \endverbatim
118: *>
119: *> \param[in] LDH
120: *> \verbatim
121: *> LDH is INTEGER
122: *> The leading dimension of the workspace H. LDH >= max(1,M).
123: *> \endverbatim
124: *>
125: *> \param[out] WORK
126: *> \verbatim
127: *> WORK is COMPLEX*16 workspace, dimension (M).
128: *> \endverbatim
129: *>
130: *> \param[out] INFO
131: *> \verbatim
132: *> INFO is INTEGER
133: *> = 0: successful exit
134: *> < 0: if INFO = -i, the i-th argument had an illegal value
135: *> > 0: if INFO = i, D(i,i) is exactly zero. The factorization
136: *> has been completed, but the block diagonal matrix D is
137: *> exactly singular, and division by zero will occur if it
138: *> is used to solve a system of equations.
139: *> \endverbatim
140: *
141: * Authors:
142: * ========
143: *
144: *> \author Univ. of Tennessee
145: *> \author Univ. of California Berkeley
146: *> \author Univ. of Colorado Denver
147: *> \author NAG Ltd.
148: *
149: *> \date December 2016
150: *
151: *> \ingroup complex16HEcomputational
152: *
153: * =====================================================================
154: SUBROUTINE ZLAHEF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
155: $ H, LDH, WORK, INFO )
156: *
157: * -- LAPACK computational routine (version 3.7.0) --
158: * -- LAPACK is a software package provided by Univ. of Tennessee, --
159: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160: * December 2016
161: *
162: IMPLICIT NONE
163: *
164: * .. Scalar Arguments ..
165: CHARACTER UPLO
166: INTEGER M, NB, J1, LDA, LDH, INFO
167: * ..
168: * .. Array Arguments ..
169: INTEGER IPIV( * )
170: COMPLEX*16 A( LDA, * ), H( LDH, * ), WORK( * )
171: * ..
172: *
173: * =====================================================================
174: * .. Parameters ..
175: COMPLEX*16 ZERO, ONE
176: PARAMETER ( ZERO = (0.0D+0, 0.0D+0), ONE = (1.0D+0, 0.0D+0) )
177: *
178: * .. Local Scalars ..
179: INTEGER J, K, K1, I1, I2
180: COMPLEX*16 PIV, ALPHA
181: * ..
182: * .. External Functions ..
183: LOGICAL LSAME
184: INTEGER IZAMAX, ILAENV
185: EXTERNAL LSAME, ILAENV, IZAMAX
186: * ..
187: * .. External Subroutines ..
188: EXTERNAL XERBLA
189: * ..
190: * .. Intrinsic Functions ..
191: INTRINSIC DBLE, DCONJG, MAX
192: * ..
193: * .. Executable Statements ..
194: *
195: INFO = 0
196: J = 1
197: *
198: * K1 is the first column of the panel to be factorized
199: * i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks
200: *
201: K1 = (2-J1)+1
202: *
203: IF( LSAME( UPLO, 'U' ) ) THEN
204: *
205: * .....................................................
206: * Factorize A as U**T*D*U using the upper triangle of A
207: * .....................................................
208: *
209: 10 CONTINUE
210: IF ( J.GT.MIN(M, NB) )
211: $ GO TO 20
212: *
213: * K is the column to be factorized
214: * when being called from ZHETRF_AA,
215: * > for the first block column, J1 is 1, hence J1+J-1 is J,
216: * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
217: *
218: K = J1+J-1
219: *
220: * H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
221: * where H(J:N, J) has been initialized to be A(J, J:N)
222: *
223: IF( K.GT.2 ) THEN
224: *
225: * K is the column to be factorized
226: * > for the first block column, K is J, skipping the first two
227: * columns
228: * > for the rest of the columns, K is J+1, skipping only the
229: * first column
230: *
231: CALL ZLACGV( J-K1, A( 1, J ), 1 )
232: CALL ZGEMV( 'No transpose', M-J+1, J-K1,
233: $ -ONE, H( J, K1 ), LDH,
234: $ A( 1, J ), 1,
235: $ ONE, H( J, J ), 1 )
236: CALL ZLACGV( J-K1, A( 1, J ), 1 )
237: END IF
238: *
239: * Copy H(i:n, i) into WORK
240: *
241: CALL ZCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
242: *
243: IF( J.GT.K1 ) THEN
244: *
245: * Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
246: * where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
247: *
248: ALPHA = -DCONJG( A( K-1, J ) )
249: CALL ZAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )
250: END IF
251: *
252: * Set A(J, J) = T(J, J)
253: *
254: A( K, J ) = DBLE( WORK( 1 ) )
255: *
256: IF( J.LT.M ) THEN
257: *
258: * Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
259: * where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
260: *
261: IF( K.GT.1 ) THEN
262: ALPHA = -A( K, J )
263: CALL ZAXPY( M-J, ALPHA, A( K-1, J+1 ), LDA,
264: $ WORK( 2 ), 1 )
265: ENDIF
266: *
267: * Find max(|WORK(2:n)|)
268: *
269: I2 = IZAMAX( M-J, WORK( 2 ), 1 ) + 1
270: PIV = WORK( I2 )
271: *
272: * Apply hermitian pivot
273: *
274: IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN
275: *
276: * Swap WORK(I1) and WORK(I2)
277: *
278: I1 = 2
279: WORK( I2 ) = WORK( I1 )
280: WORK( I1 ) = PIV
281: *
282: * Swap A(I1, I1+1:N) with A(I1+1:N, I2)
283: *
284: I1 = I1+J-1
285: I2 = I2+J-1
286: CALL ZSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA,
287: $ A( J1+I1, I2 ), 1 )
288: CALL ZLACGV( I2-I1, A( J1+I1-1, I1+1 ), LDA )
289: CALL ZLACGV( I2-I1-1, A( J1+I1, I2 ), 1 )
290: *
291: * Swap A(I1, I2+1:N) with A(I2, I2+1:N)
292: *
293: CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
294: $ A( J1+I2-1, I2+1 ), LDA )
295: *
296: * Swap A(I1, I1) with A(I2,I2)
297: *
298: PIV = A( I1+J1-1, I1 )
299: A( J1+I1-1, I1 ) = A( J1+I2-1, I2 )
300: A( J1+I2-1, I2 ) = PIV
301: *
302: * Swap H(I1, 1:J1) with H(I2, 1:J1)
303: *
304: CALL ZSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH )
305: IPIV( I1 ) = I2
306: *
307: IF( I1.GT.(K1-1) ) THEN
308: *
309: * Swap L(1:I1-1, I1) with L(1:I1-1, I2),
310: * skipping the first column
311: *
312: CALL ZSWAP( I1-K1+1, A( 1, I1 ), 1,
313: $ A( 1, I2 ), 1 )
314: END IF
315: ELSE
316: IPIV( J+1 ) = J+1
317: ENDIF
318: *
319: * Set A(J, J+1) = T(J, J+1)
320: *
321: A( K, J+1 ) = WORK( 2 )
322: IF( (A( K, J ).EQ.ZERO ) .AND.
323: $ ( (J.EQ.M) .OR. (A( K, J+1 ).EQ.ZERO))) THEN
324: IF(INFO .EQ. 0) THEN
325: INFO = J
326: END IF
327: END IF
328: *
329: IF( J.LT.NB ) THEN
330: *
331: * Copy A(J+1:N, J+1) into H(J:N, J),
332: *
333: CALL ZCOPY( M-J, A( K+1, J+1 ), LDA,
334: $ H( J+1, J+1 ), 1 )
335: END IF
336: *
337: * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
338: * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
339: *
340: IF( A( K, J+1 ).NE.ZERO ) THEN
341: ALPHA = ONE / A( K, J+1 )
342: CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
343: CALL ZSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA )
344: ELSE
345: CALL ZLASET( 'Full', 1, M-J-1, ZERO, ZERO,
346: $ A( K, J+2 ), LDA)
347: END IF
348: ELSE
349: IF( (A( K, J ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN
350: INFO = J
351: END IF
352: END IF
353: J = J + 1
354: GO TO 10
355: 20 CONTINUE
356: *
357: ELSE
358: *
359: * .....................................................
360: * Factorize A as L*D*L**T using the lower triangle of A
361: * .....................................................
362: *
363: 30 CONTINUE
364: IF( J.GT.MIN( M, NB ) )
365: $ GO TO 40
366: *
367: * K is the column to be factorized
368: * when being called from ZHETRF_AA,
369: * > for the first block column, J1 is 1, hence J1+J-1 is J,
370: * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
371: *
372: K = J1+J-1
373: *
374: * H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
375: * where H(J:N, J) has been initialized to be A(J:N, J)
376: *
377: IF( K.GT.2 ) THEN
378: *
379: * K is the column to be factorized
380: * > for the first block column, K is J, skipping the first two
381: * columns
382: * > for the rest of the columns, K is J+1, skipping only the
383: * first column
384: *
385: CALL ZLACGV( J-K1, A( J, 1 ), LDA )
386: CALL ZGEMV( 'No transpose', M-J+1, J-K1,
387: $ -ONE, H( J, K1 ), LDH,
388: $ A( J, 1 ), LDA,
389: $ ONE, H( J, J ), 1 )
390: CALL ZLACGV( J-K1, A( J, 1 ), LDA )
391: END IF
392: *
393: * Copy H(J:N, J) into WORK
394: *
395: CALL ZCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
396: *
397: IF( J.GT.K1 ) THEN
398: *
399: * Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
400: * where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
401: *
402: ALPHA = -DCONJG( A( J, K-1 ) )
403: CALL ZAXPY( M-J+1, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 )
404: END IF
405: *
406: * Set A(J, J) = T(J, J)
407: *
408: A( J, K ) = DBLE( WORK( 1 ) )
409: *
410: IF( J.LT.M ) THEN
411: *
412: * Compute WORK(2:N) = T(J, J) L((J+1):N, J)
413: * where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
414: *
415: IF( K.GT.1 ) THEN
416: ALPHA = -A( J, K )
417: CALL ZAXPY( M-J, ALPHA, A( J+1, K-1 ), 1,
418: $ WORK( 2 ), 1 )
419: ENDIF
420: *
421: * Find max(|WORK(2:n)|)
422: *
423: I2 = IZAMAX( M-J, WORK( 2 ), 1 ) + 1
424: PIV = WORK( I2 )
425: *
426: * Apply hermitian pivot
427: *
428: IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN
429: *
430: * Swap WORK(I1) and WORK(I2)
431: *
432: I1 = 2
433: WORK( I2 ) = WORK( I1 )
434: WORK( I1 ) = PIV
435: *
436: * Swap A(I1+1:N, I1) with A(I2, I1+1:N)
437: *
438: I1 = I1+J-1
439: I2 = I2+J-1
440: CALL ZSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1,
441: $ A( I2, J1+I1 ), LDA )
442: CALL ZLACGV( I2-I1, A( I1+1, J1+I1-1 ), 1 )
443: CALL ZLACGV( I2-I1-1, A( I2, J1+I1 ), LDA )
444: *
445: * Swap A(I2+1:N, I1) with A(I2+1:N, I2)
446: *
447: CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
448: $ A( I2+1, J1+I2-1 ), 1 )
449: *
450: * Swap A(I1, I1) with A(I2, I2)
451: *
452: PIV = A( I1, J1+I1-1 )
453: A( I1, J1+I1-1 ) = A( I2, J1+I2-1 )
454: A( I2, J1+I2-1 ) = PIV
455: *
456: * Swap H(I1, I1:J1) with H(I2, I2:J1)
457: *
458: CALL ZSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH )
459: IPIV( I1 ) = I2
460: *
461: IF( I1.GT.(K1-1) ) THEN
462: *
463: * Swap L(1:I1-1, I1) with L(1:I1-1, I2),
464: * skipping the first column
465: *
466: CALL ZSWAP( I1-K1+1, A( I1, 1 ), LDA,
467: $ A( I2, 1 ), LDA )
468: END IF
469: ELSE
470: IPIV( J+1 ) = J+1
471: ENDIF
472: *
473: * Set A(J+1, J) = T(J+1, J)
474: *
475: A( J+1, K ) = WORK( 2 )
476: IF( (A( J, K ).EQ.ZERO) .AND.
477: $ ( (J.EQ.M) .OR. (A( J+1, K ).EQ.ZERO)) ) THEN
478: IF (INFO .EQ. 0)
479: $ INFO = J
480: END IF
481: *
482: IF( J.LT.NB ) THEN
483: *
484: * Copy A(J+1:N, J+1) into H(J+1:N, J),
485: *
486: CALL ZCOPY( M-J, A( J+1, K+1 ), 1,
487: $ H( J+1, J+1 ), 1 )
488: END IF
489: *
490: * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
491: * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
492: *
493: IF( A( J+1, K ).NE.ZERO ) THEN
494: ALPHA = ONE / A( J+1, K )
495: CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
496: CALL ZSCAL( M-J-1, ALPHA, A( J+2, K ), 1 )
497: ELSE
498: CALL ZLASET( 'Full', M-J-1, 1, ZERO, ZERO,
499: $ A( J+2, K ), LDA )
500: END IF
501: ELSE
502: IF( (A( J, K ).EQ.ZERO) .AND. (J.EQ.M)
503: $ .AND. (INFO.EQ.0) ) INFO = J
504: END IF
505: J = J + 1
506: GO TO 30
507: 40 CONTINUE
508: END IF
509: RETURN
510: *
511: * End of ZLAHEF_AA
512: *
513: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlahef_aa.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack