LAPACK  3.10.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.

Definition at line 142 of file zlahef_aa.f.

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