LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ zlahef_aa()

subroutine zlahef_aa ( character  UPLO,
integer  J1,
integer  M,
integer  NB,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex*16, dimension( ldh, * )  H,
integer  LDH,
complex*16, dimension( * )  WORK 
)

ZLAHEF_AA

Download ZLAHEF_AA + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DLAHEF_AA factorizes a panel of a complex hermitian matrix A using
 the Aasen's algorithm. The panel consists of a set of NB rows of A
 when UPLO is U, or a set of NB columns when UPLO is L.

 In order to factorize the panel, the Aasen's algorithm requires the
 last row, or column, of the previous panel. The first row, or column,
 of A is set to be the first row, or column, of an identity matrix,
 which is used to factorize the first panel.

 The resulting J-th row of U, or J-th column of L, is stored in the
 (J-1)-th row, or column, of A (without the unit diagonals), while
 the diagonal and subdiagonal of A are overwritten by those of T.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]J1
          J1 is INTEGER
          The location of the first row, or column, of the panel
          within the submatrix of A, passed to this routine, e.g.,
          when called by ZHETRF_AA, for the first panel, J1 is 1,
          while for the remaining panels, J1 is 2.
[in]M
          M is INTEGER
          The dimension of the submatrix. M >= 0.
[in]NB
          NB is INTEGER
          The dimension of the panel to be facotorized.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,M) for
          the first panel, while dimension (LDA,M+1) for the
          remaining panels.

          On entry, A contains the last row, or column, of
          the previous panel, and the trailing submatrix of A
          to be factorized, except for the first panel, only
          the panel is passed.

          On exit, the leading panel is factorized.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the row and column interchanges,
          the row and column k were interchanged with the row and
          column IPIV(k).
[in,out]H
          H is COMPLEX*16 workspace, dimension (LDH,NB).
[in]LDH
          LDH is INTEGER
          The leading dimension of the workspace H. LDH >= max(1,M).
[out]WORK
          WORK is COMPLEX*16 workspace, dimension (M).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2017

Definition at line 146 of file zlahef_aa.f.

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 *
Here is the call graph for this function:
Here is the caller graph for this function:
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
izamax
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:73
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
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:83
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80