LAPACK  3.8.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  CALL cswap( m-i2, a( j1+i1-1, i2+1 ), lda,
292  $ a( j1+i2-1, i2+1 ), lda )
293 *
294 * Swap A(I1, I1) with A(I2,I2)
295 *
296  piv = a( i1+j1-1, i1 )
297  a( j1+i1-1, i1 ) = a( j1+i2-1, i2 )
298  a( j1+i2-1, i2 ) = piv
299 *
300 * Swap H(I1, 1:J1) with H(I2, 1:J1)
301 *
302  CALL cswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
303  ipiv( i1 ) = i2
304 *
305  IF( i1.GT.(k1-1) ) THEN
306 *
307 * Swap L(1:I1-1, I1) with L(1:I1-1, I2),
308 * skipping the first column
309 *
310  CALL cswap( i1-k1+1, a( 1, i1 ), 1,
311  $ a( 1, i2 ), 1 )
312  END IF
313  ELSE
314  ipiv( j+1 ) = j+1
315  ENDIF
316 *
317 * Set A(J, J+1) = T(J, J+1)
318 *
319  a( k, j+1 ) = work( 2 )
320 *
321  IF( j.LT.nb ) THEN
322 *
323 * Copy A(J+1:N, J+1) into H(J:N, J),
324 *
325  CALL ccopy( m-j, a( k+1, j+1 ), lda,
326  $ h( j+1, j+1 ), 1 )
327  END IF
328 *
329 * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
330 * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
331 *
332  IF( a( k, j+1 ).NE.zero ) THEN
333  alpha = one / a( k, j+1 )
334  CALL ccopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
335  CALL cscal( m-j-1, alpha, a( k, j+2 ), lda )
336  ELSE
337  CALL claset( 'Full', 1, m-j-1, zero, zero,
338  $ a( k, j+2 ), lda)
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  CALL cswap( m-i2, a( i2+1, j1+i1-1 ), 1,
444  $ a( i2+1, j1+i2-1 ), 1 )
445 *
446 * Swap A(I1, I1) with A(I2, I2)
447 *
448  piv = a( i1, j1+i1-1 )
449  a( i1, j1+i1-1 ) = a( i2, j1+i2-1 )
450  a( i2, j1+i2-1 ) = piv
451 *
452 * Swap H(I1, I1:J1) with H(I2, I2:J1)
453 *
454  CALL cswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
455  ipiv( i1 ) = i2
456 *
457  IF( i1.GT.(k1-1) ) THEN
458 *
459 * Swap L(1:I1-1, I1) with L(1:I1-1, I2),
460 * skipping the first column
461 *
462  CALL cswap( i1-k1+1, a( i1, 1 ), lda,
463  $ a( i2, 1 ), lda )
464  END IF
465  ELSE
466  ipiv( j+1 ) = j+1
467  ENDIF
468 *
469 * Set A(J+1, J) = T(J+1, J)
470 *
471  a( j+1, k ) = work( 2 )
472 *
473  IF( j.LT.nb ) THEN
474 *
475 * Copy A(J+1:N, J+1) into H(J+1:N, J),
476 *
477  CALL ccopy( m-j, a( j+1, k+1 ), 1,
478  $ h( j+1, j+1 ), 1 )
479  END IF
480 *
481 * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
482 * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
483 *
484  IF( a( j+1, k ).NE.zero ) THEN
485  alpha = one / a( j+1, k )
486  CALL ccopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
487  CALL cscal( m-j-1, alpha, a( j+2, k ), 1 )
488  ELSE
489  CALL claset( 'Full', m-j-1, 1, zero, zero,
490  $ a( j+2, k ), lda )
491  END IF
492  END IF
493  j = j + 1
494  GO TO 30
495  40 CONTINUE
496  END IF
497  RETURN
498 *
499 * End of CLAHEF_AA
500 *
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:73
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:80
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:76
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:83
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:90
Here is the call graph for this function:
Here is the caller graph for this function: