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

Definition at line 142 of file clahef_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 A( LDA, * ), H( LDH, * ), WORK( * )
158 * ..
159 *
160 * =====================================================================
161 * .. Parameters ..
162  COMPLEX ZERO, ONE
163  parameter( zero = (0.0e+0, 0.0e+0), one = (1.0e+0, 0.0e+0) )
164 *
165 * .. Local Scalars ..
166  INTEGER J, K, K1, I1, I2, MJ
167  COMPLEX PIV, ALPHA
168 * ..
169 * .. External Functions ..
170  LOGICAL LSAME
171  INTEGER ICAMAX, ILAENV
172  EXTERNAL lsame, ilaenv, icamax
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL clacgv, cgemv, cscal, caxpy, ccopy, cswap, claset,
176  $ xerbla
177 * ..
178 * .. Intrinsic Functions ..
179  INTRINSIC real, conjg, 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 CHETRF_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 clacgv( j-k1, a( 1, j ), 1 )
227  CALL cgemv( 'No transpose', mj, j-k1,
228  $ -one, h( j, k1 ), ldh,
229  $ a( 1, j ), 1,
230  $ one, h( j, j ), 1 )
231  CALL clacgv( j-k1, a( 1, j ), 1 )
232  END IF
233 *
234 * Copy H(i:n, i) into WORK
235 *
236  CALL ccopy( 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 = -conjg( a( k-1, j ) )
244  CALL caxpy( 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 ) = real( 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 caxpy( 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 = icamax( 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 cswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
282  $ a( j1+i1, i2 ), 1 )
283  CALL clacgv( i2-i1, a( j1+i1-1, i1+1 ), lda )
284  CALL clacgv( 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 cswap( 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 cswap( 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 cswap( 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 ccopy( 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 ccopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
334  CALL cscal( m-j-1, alpha, a( k, j+2 ), lda )
335  ELSE
336  CALL claset( '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 CHETRF_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 clacgv( j-k1, a( j, 1 ), lda )
382  CALL cgemv( 'No transpose', mj, j-k1,
383  $ -one, h( j, k1 ), ldh,
384  $ a( j, 1 ), lda,
385  $ one, h( j, j ), 1 )
386  CALL clacgv( j-k1, a( j, 1 ), lda )
387  END IF
388 *
389 * Copy H(J:N, J) into WORK
390 *
391  CALL ccopy( 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 = -conjg( a( j, k-1 ) )
399  CALL caxpy( 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 ) = real( 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 caxpy( 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 = icamax( 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 cswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
437  $ a( i2, j1+i1 ), lda )
438  CALL clacgv( i2-i1, a( i1+1, j1+i1-1 ), 1 )
439  CALL clacgv( 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 cswap( 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 cswap( 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 cswap( 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 ccopy( 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 ccopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
489  CALL cscal( m-j-1, alpha, a( j+2, k ), 1 )
490  ELSE
491  CALL claset( '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 CLAHEF_AA
503 *
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:71
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:74
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:106
Here is the call graph for this function:
Here is the caller graph for this function: