LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
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,
integer  INFO 
)

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).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization
                has been completed, but the block diagonal matrix D is
                exactly singular, and division by zero will occur if it
                is used to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 156 of file zlahef_aa.f.

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 *
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
zgemv
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
izamax
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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:53
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: