LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
subroutine clahef_aa ( character  UPLO,
integer  J1,
integer  M,
integer  NB,
complex, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
complex, dimension( ldh, * )  H,
integer  LDH,
complex, dimension( * )  WORK,
integer  INFO 
)

CLAHEF_AA

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

Purpose: 
 CLAHEF_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 CHETRF_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 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 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 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 clahef_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 a( lda, * ), h( ldh, * ), work( * )
171 * ..
172 *
173 * =====================================================================
174 * .. Parameters ..
175  COMPLEX zero, one
176  parameter ( zero = (0.0e+0, 0.0e+0), one = (1.0e+0, 0.0e+0) )
177 *
178 * .. Local Scalars ..
179  INTEGER j, k, k1, i1, i2
180  COMPLEX piv, alpha
181 * ..
182 * .. External Functions ..
183  LOGICAL lsame
184  INTEGER icamax, ilaenv
185  EXTERNAL lsame, ilaenv, icamax
186 * ..
187 * .. External Subroutines ..
188  EXTERNAL xerbla
189 * ..
190 * .. Intrinsic Functions ..
191  INTRINSIC REAL, conjg, 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 CHETRF_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 clacgv( j-k1, a( 1, j ), 1 )
232  CALL cgemv( '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 clacgv( j-k1, a( 1, j ), 1 )
237  END IF
238 *
239 * Copy H(i:n, i) into WORK
240 *
241  CALL ccopy( 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 = -conjg( a( k-1, j ) )
249  CALL caxpy( 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 ) = REAL( 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 caxpy( 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 = icamax( 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 cswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
287  $ a( j1+i1, i2 ), 1 )
288  CALL clacgv( i2-i1, a( j1+i1-1, i1+1 ), lda )
289  CALL clacgv( 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 cswap( 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 cswap( 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 cswap( 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 ccopy( 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 ccopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
343  CALL cscal( m-j-1, alpha, a( k, j+2 ), lda )
344  ELSE
345  CALL claset( '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 CHETRF_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 clacgv( j-k1, a( j, 1 ), lda )
386  CALL cgemv( '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 clacgv( j-k1, a( j, 1 ), lda )
391  END IF
392 *
393 * Copy H(J:N, J) into WORK
394 *
395  CALL ccopy( 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 = -conjg( a( j, k-1 ) )
403  CALL caxpy( 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 ) = REAL( 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 caxpy( 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 = icamax( 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 cswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
441  $ a( i2, j1+i1 ), lda )
442  CALL clacgv( i2-i1, a( i1+1, j1+i1-1 ), 1 )
443  CALL clacgv( 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 cswap( 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 cswap( 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 cswap( 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 ccopy( 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 ccopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
496  CALL cscal( m-j-1, alpha, a( j+2, k ), 1 )
497  ELSE
498  CALL claset( '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 CLAHEF_AA
512 *
icamax
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:53
claset
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
cscal
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
clacgv
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:76
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: