LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
zlahef_aa.f
Go to the documentation of this file.
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 )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER UPLO
26 * INTEGER J1, M, NB, LDA, LDH
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 *
131 * Authors:
132 * ========
133 *
134 *> \author Univ. of Tennessee
135 *> \author Univ. of California Berkeley
136 *> \author Univ. of Colorado Denver
137 *> \author NAG Ltd.
138 *
139 *> \date November 2017
140 *
141 *> \ingroup complex16HEcomputational
142 *
143 * =====================================================================
144  SUBROUTINE zlahef_aa( UPLO, J1, M, NB, A, LDA, IPIV,
145  $ H, LDH, WORK )
146 *
147 * -- LAPACK computational routine (version 3.8.0) --
148 * -- LAPACK is a software package provided by Univ. of Tennessee, --
149 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150 * November 2017
151 *
152  IMPLICIT NONE
153 *
154 * .. Scalar Arguments ..
155  CHARACTER UPLO
156  INTEGER M, NB, J1, LDA, LDH
157 * ..
158 * .. Array Arguments ..
159  INTEGER IPIV( * )
160  COMPLEX*16 A( LDA, * ), H( LDH, * ), WORK( * )
161 * ..
162 *
163 * =====================================================================
164 * .. Parameters ..
165  COMPLEX*16 ZERO, ONE
166  parameter( zero = (0.0d+0, 0.0d+0), one = (1.0d+0, 0.0d+0) )
167 *
168 * .. Local Scalars ..
169  INTEGER J, K, K1, I1, I2, MJ
170  COMPLEX*16 PIV, ALPHA
171 * ..
172 * .. External Functions ..
173  LOGICAL LSAME
174  INTEGER IZAMAX, ILAENV
175  EXTERNAL lsame, ilaenv, izamax
176 * ..
177 * .. External Subroutines ..
178  EXTERNAL zgemm, zgemv, zaxpy, zlacgv, zcopy, zscal, zswap,
179  $ zlaset, xerbla
180 * ..
181 * .. Intrinsic Functions ..
182  INTRINSIC dble, dconjg, max
183 * ..
184 * .. Executable Statements ..
185 *
186  j = 1
187 *
188 * K1 is the first column of the panel to be factorized
189 * i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks
190 *
191  k1 = (2-j1)+1
192 *
193  IF( lsame( uplo, 'U' ) ) THEN
194 *
195 * .....................................................
196 * Factorize A as U**T*D*U using the upper triangle of A
197 * .....................................................
198 *
199  10 CONTINUE
200  IF ( j.GT.min(m, nb) )
201  $ GO TO 20
202 *
203 * K is the column to be factorized
204 * when being called from ZHETRF_AA,
205 * > for the first block column, J1 is 1, hence J1+J-1 is J,
206 * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
207 *
208  k = j1+j-1
209  IF( j.EQ.m ) THEN
210 *
211 * Only need to compute T(J, J)
212 *
213  mj = 1
214  ELSE
215  mj = m-j+1
216  END IF
217 *
218 * H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
219 * where H(J:N, J) has been initialized to be A(J, J:N)
220 *
221  IF( k.GT.2 ) THEN
222 *
223 * K is the column to be factorized
224 * > for the first block column, K is J, skipping the first two
225 * columns
226 * > for the rest of the columns, K is J+1, skipping only the
227 * first column
228 *
229  CALL zlacgv( j-k1, a( 1, j ), 1 )
230  CALL zgemv( 'No transpose', mj, j-k1,
231  $ -one, h( j, k1 ), ldh,
232  $ a( 1, j ), 1,
233  $ one, h( j, j ), 1 )
234  CALL zlacgv( j-k1, a( 1, j ), 1 )
235  END IF
236 *
237 * Copy H(i:n, i) into WORK
238 *
239  CALL zcopy( mj, h( j, j ), 1, work( 1 ), 1 )
240 *
241  IF( j.GT.k1 ) THEN
242 *
243 * Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
244 * where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
245 *
246  alpha = -dconjg( a( k-1, j ) )
247  CALL zaxpy( mj, alpha, a( k-2, j ), lda, work( 1 ), 1 )
248  END IF
249 *
250 * Set A(J, J) = T(J, J)
251 *
252  a( k, j ) = dble( work( 1 ) )
253 *
254  IF( j.LT.m ) THEN
255 *
256 * Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
257 * where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
258 *
259  IF( k.GT.1 ) THEN
260  alpha = -a( k, j )
261  CALL zaxpy( m-j, alpha, a( k-1, j+1 ), lda,
262  $ work( 2 ), 1 )
263  ENDIF
264 *
265 * Find max(|WORK(2:n)|)
266 *
267  i2 = izamax( m-j, work( 2 ), 1 ) + 1
268  piv = work( i2 )
269 *
270 * Apply hermitian pivot
271 *
272  IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
273 *
274 * Swap WORK(I1) and WORK(I2)
275 *
276  i1 = 2
277  work( i2 ) = work( i1 )
278  work( i1 ) = piv
279 *
280 * Swap A(I1, I1+1:N) with A(I1+1:N, I2)
281 *
282  i1 = i1+j-1
283  i2 = i2+j-1
284  CALL zswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
285  $ a( j1+i1, i2 ), 1 )
286  CALL zlacgv( i2-i1, a( j1+i1-1, i1+1 ), lda )
287  CALL zlacgv( i2-i1-1, a( j1+i1, i2 ), 1 )
288 *
289 * Swap A(I1, I2+1:N) with A(I2, I2+1:N)
290 *
291  IF( i2.LT.m )
292  $ CALL zswap( m-i2, a( j1+i1-1, i2+1 ), lda,
293  $ a( j1+i2-1, i2+1 ), lda )
294 *
295 * Swap A(I1, I1) with A(I2,I2)
296 *
297  piv = a( i1+j1-1, i1 )
298  a( j1+i1-1, i1 ) = a( j1+i2-1, i2 )
299  a( j1+i2-1, i2 ) = piv
300 *
301 * Swap H(I1, 1:J1) with H(I2, 1:J1)
302 *
303  CALL zswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
304  ipiv( i1 ) = i2
305 *
306  IF( i1.GT.(k1-1) ) THEN
307 *
308 * Swap L(1:I1-1, I1) with L(1:I1-1, I2),
309 * skipping the first column
310 *
311  CALL zswap( i1-k1+1, a( 1, i1 ), 1,
312  $ a( 1, i2 ), 1 )
313  END IF
314  ELSE
315  ipiv( j+1 ) = j+1
316  ENDIF
317 *
318 * Set A(J, J+1) = T(J, J+1)
319 *
320  a( k, j+1 ) = work( 2 )
321 *
322  IF( j.LT.nb ) THEN
323 *
324 * Copy A(J+1:N, J+1) into H(J:N, J),
325 *
326  CALL zcopy( m-j, a( k+1, j+1 ), lda,
327  $ h( j+1, j+1 ), 1 )
328  END IF
329 *
330 * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
331 * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
332 *
333  IF( j.LT.(m-1) ) THEN
334  IF( a( k, j+1 ).NE.zero ) THEN
335  alpha = one / a( k, j+1 )
336  CALL zcopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
337  CALL zscal( m-j-1, alpha, a( k, j+2 ), lda )
338  ELSE
339  CALL zlaset( 'Full', 1, m-j-1, zero, zero,
340  $ a( k, j+2 ), lda)
341  END IF
342  END IF
343  END IF
344  j = j + 1
345  GO TO 10
346  20 CONTINUE
347 *
348  ELSE
349 *
350 * .....................................................
351 * Factorize A as L*D*L**T using the lower triangle of A
352 * .....................................................
353 *
354  30 CONTINUE
355  IF( j.GT.min( m, nb ) )
356  $ GO TO 40
357 *
358 * K is the column to be factorized
359 * when being called from ZHETRF_AA,
360 * > for the first block column, J1 is 1, hence J1+J-1 is J,
361 * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
362 *
363  k = j1+j-1
364  IF( j.EQ.m ) THEN
365 *
366 * Only need to compute T(J, J)
367 *
368  mj = 1
369  ELSE
370  mj = m-j+1
371  END IF
372 *
373 * H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
374 * where H(J:N, J) has been initialized to be A(J:N, J)
375 *
376  IF( k.GT.2 ) THEN
377 *
378 * K is the column to be factorized
379 * > for the first block column, K is J, skipping the first two
380 * columns
381 * > for the rest of the columns, K is J+1, skipping only the
382 * first column
383 *
384  CALL zlacgv( j-k1, a( j, 1 ), lda )
385  CALL zgemv( 'No transpose', mj, j-k1,
386  $ -one, h( j, k1 ), ldh,
387  $ a( j, 1 ), lda,
388  $ one, h( j, j ), 1 )
389  CALL zlacgv( j-k1, a( j, 1 ), lda )
390  END IF
391 *
392 * Copy H(J:N, J) into WORK
393 *
394  CALL zcopy( mj, h( j, j ), 1, work( 1 ), 1 )
395 *
396  IF( j.GT.k1 ) THEN
397 *
398 * Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
399 * where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
400 *
401  alpha = -dconjg( a( j, k-1 ) )
402  CALL zaxpy( mj, alpha, a( j, k-2 ), 1, work( 1 ), 1 )
403  END IF
404 *
405 * Set A(J, J) = T(J, J)
406 *
407  a( j, k ) = dble( work( 1 ) )
408 *
409  IF( j.LT.m ) THEN
410 *
411 * Compute WORK(2:N) = T(J, J) L((J+1):N, J)
412 * where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
413 *
414  IF( k.GT.1 ) THEN
415  alpha = -a( j, k )
416  CALL zaxpy( m-j, alpha, a( j+1, k-1 ), 1,
417  $ work( 2 ), 1 )
418  ENDIF
419 *
420 * Find max(|WORK(2:n)|)
421 *
422  i2 = izamax( m-j, work( 2 ), 1 ) + 1
423  piv = work( i2 )
424 *
425 * Apply hermitian pivot
426 *
427  IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
428 *
429 * Swap WORK(I1) and WORK(I2)
430 *
431  i1 = 2
432  work( i2 ) = work( i1 )
433  work( i1 ) = piv
434 *
435 * Swap A(I1+1:N, I1) with A(I2, I1+1:N)
436 *
437  i1 = i1+j-1
438  i2 = i2+j-1
439  CALL zswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
440  $ a( i2, j1+i1 ), lda )
441  CALL zlacgv( i2-i1, a( i1+1, j1+i1-1 ), 1 )
442  CALL zlacgv( i2-i1-1, a( i2, j1+i1 ), lda )
443 *
444 * Swap A(I2+1:N, I1) with A(I2+1:N, I2)
445 *
446  IF( i2.LT.m )
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 *
477  IF( j.LT.nb ) THEN
478 *
479 * Copy A(J+1:N, J+1) into H(J+1:N, J),
480 *
481  CALL zcopy( m-j, a( j+1, k+1 ), 1,
482  $ h( j+1, j+1 ), 1 )
483  END IF
484 *
485 * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
486 * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
487 *
488  IF( j.LT.(m-1) ) THEN
489  IF( a( j+1, k ).NE.zero ) THEN
490  alpha = one / a( j+1, k )
491  CALL zcopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
492  CALL zscal( m-j-1, alpha, a( j+2, k ), 1 )
493  ELSE
494  CALL zlaset( 'Full', m-j-1, 1, zero, zero,
495  $ a( j+2, k ), lda )
496  END IF
497  END IF
498  END IF
499  j = j + 1
500  GO TO 30
501  40 CONTINUE
502  END IF
503  RETURN
504 *
505 * End of ZLAHEF_AA
506 *
507  END
zlahef_aa
subroutine zlahef_aa(UPLO, J1, M, NB, A, LDA, IPIV, H, LDH, WORK)
ZLAHEF_AA
Definition: zlahef_aa.f:146
zlacgv
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:76
zaxpy
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:90
zgemv
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:83
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:108
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:83
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80