LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ clahef_aa()

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 
)

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).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2017

Definition at line 146 of file clahef_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 A( LDA, * ), H( LDH, * ), WORK( * )
161 * ..
162 *
163 * =====================================================================
164 * .. Parameters ..
165  COMPLEX ZERO, ONE
166  parameter( zero = (0.0e+0, 0.0e+0), one = (1.0e+0, 0.0e+0) )
167 *
168 * .. Local Scalars ..
169  INTEGER J, K, K1, I1, I2, MJ
170  COMPLEX PIV, ALPHA
171 * ..
172 * .. External Functions ..
173  LOGICAL LSAME
174  INTEGER ICAMAX, ILAENV
175  EXTERNAL lsame, ilaenv, icamax
176 * ..
177 * .. External Subroutines ..
178  EXTERNAL clacgv, cgemv, cscal, caxpy, ccopy, cswap, claset,
179  $ xerbla
180 * ..
181 * .. Intrinsic Functions ..
182  INTRINSIC real, conjg, 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 CHETRF_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 clacgv( j-k1, a( 1, j ), 1 )
230  CALL cgemv( 'No transpose', mj, j-k1,
231  $ -one, h( j, k1 ), ldh,
232  $ a( 1, j ), 1,
233  $ one, h( j, j ), 1 )
234  CALL clacgv( j-k1, a( 1, j ), 1 )
235  END IF
236 *
237 * Copy H(i:n, i) into WORK
238 *
239  CALL ccopy( 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 = -conjg( a( k-1, j ) )
247  CALL caxpy( 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 ) = real( 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 caxpy( 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 = icamax( 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 cswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
285  $ a( j1+i1, i2 ), 1 )
286  CALL clacgv( i2-i1, a( j1+i1-1, i1+1 ), lda )
287  CALL clacgv( 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 cswap( 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 cswap( 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 cswap( 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 ccopy( 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 ccopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
337  CALL cscal( m-j-1, alpha, a( k, j+2 ), lda )
338  ELSE
339  CALL claset( '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 CHETRF_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 clacgv( j-k1, a( j, 1 ), lda )
385  CALL cgemv( 'No transpose', mj, j-k1,
386  $ -one, h( j, k1 ), ldh,
387  $ a( j, 1 ), lda,
388  $ one, h( j, j ), 1 )
389  CALL clacgv( j-k1, a( j, 1 ), lda )
390  END IF
391 *
392 * Copy H(J:N, J) into WORK
393 *
394  CALL ccopy( 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 = -conjg( a( j, k-1 ) )
402  CALL caxpy( 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 ) = real( 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 caxpy( 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 = icamax( 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 cswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
440  $ a( i2, j1+i1 ), lda )
441  CALL clacgv( i2-i1, a( i1+1, j1+i1-1 ), 1 )
442  CALL clacgv( 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 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 *
477  IF( j.LT.nb ) THEN
478 *
479 * Copy A(J+1:N, J+1) into H(J+1:N, J),
480 *
481  CALL ccopy( 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 ccopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
492  CALL cscal( m-j-1, alpha, a( j+2, k ), 1 )
493  ELSE
494  CALL claset( '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 CLAHEF_AA
506 *
Here is the call graph for this function:
Here is the caller graph for this function:
cgemv
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
clacgv
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:76
cscal
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:80
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:83
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:90
icamax
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:73