LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 140 of file clahef_aa.f.

142*
143* -- LAPACK computational routine --
144* -- LAPACK is a software package provided by Univ. of Tennessee, --
145* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146*
147 IMPLICIT NONE
148*
149* .. Scalar Arguments ..
150 CHARACTER UPLO
151 INTEGER M, NB, J1, LDA, LDH
152* ..
153* .. Array Arguments ..
154 INTEGER IPIV( * )
155 COMPLEX A( LDA, * ), H( LDH, * ), WORK( * )
156* ..
157*
158* =====================================================================
159* .. Parameters ..
160 COMPLEX ZERO, ONE
161 parameter( zero = (0.0e+0, 0.0e+0), one = (1.0e+0, 0.0e+0) )
162*
163* .. Local Scalars ..
164 INTEGER J, K, K1, I1, I2, MJ
165 COMPLEX PIV, ALPHA
166* ..
167* .. External Functions ..
168 LOGICAL LSAME
169 INTEGER ICAMAX, ILAENV
170 EXTERNAL lsame, ilaenv, icamax
171* ..
172* .. External Subroutines ..
173 EXTERNAL clacgv, cgemv, cscal, caxpy, ccopy, cswap,
174 $ claset, xerbla
175* ..
176* .. Intrinsic Functions ..
177 INTRINSIC real, conjg, max
178* ..
179* .. Executable Statements ..
180*
181 j = 1
182*
183* K1 is the first column of the panel to be factorized
184* i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks
185*
186 k1 = (2-j1)+1
187*
188 IF( lsame( uplo, 'U' ) ) THEN
189*
190* .....................................................
191* Factorize A as U**T*D*U using the upper triangle of A
192* .....................................................
193*
194 10 CONTINUE
195 IF ( j.GT.min(m, nb) )
196 $ GO TO 20
197*
198* K is the column to be factorized
199* when being called from CHETRF_AA,
200* > for the first block column, J1 is 1, hence J1+J-1 is J,
201* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
202*
203 k = j1+j-1
204 IF( j.EQ.m ) THEN
205*
206* Only need to compute T(J, J)
207*
208 mj = 1
209 ELSE
210 mj = m-j+1
211 END IF
212*
213* H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
214* where H(J:N, J) has been initialized to be A(J, J:N)
215*
216 IF( k.GT.2 ) THEN
217*
218* K is the column to be factorized
219* > for the first block column, K is J, skipping the first two
220* columns
221* > for the rest of the columns, K is J+1, skipping only the
222* first column
223*
224 CALL clacgv( j-k1, a( 1, j ), 1 )
225 CALL cgemv( 'No transpose', mj, j-k1,
226 $ -one, h( j, k1 ), ldh,
227 $ a( 1, j ), 1,
228 $ one, h( j, j ), 1 )
229 CALL clacgv( j-k1, a( 1, j ), 1 )
230 END IF
231*
232* Copy H(i:n, i) into WORK
233*
234 CALL ccopy( mj, h( j, j ), 1, work( 1 ), 1 )
235*
236 IF( j.GT.k1 ) THEN
237*
238* Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
239* where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
240*
241 alpha = -conjg( a( k-1, j ) )
242 CALL caxpy( mj, alpha, a( k-2, j ), lda, work( 1 ), 1 )
243 END IF
244*
245* Set A(J, J) = T(J, J)
246*
247 a( k, j ) = real( work( 1 ) )
248*
249 IF( j.LT.m ) THEN
250*
251* Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
252* where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
253*
254 IF( k.GT.1 ) THEN
255 alpha = -a( k, j )
256 CALL caxpy( m-j, alpha, a( k-1, j+1 ), lda,
257 $ work( 2 ), 1 )
258 ENDIF
259*
260* Find max(|WORK(2:n)|)
261*
262 i2 = icamax( m-j, work( 2 ), 1 ) + 1
263 piv = work( i2 )
264*
265* Apply hermitian pivot
266*
267 IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
268*
269* Swap WORK(I1) and WORK(I2)
270*
271 i1 = 2
272 work( i2 ) = work( i1 )
273 work( i1 ) = piv
274*
275* Swap A(I1, I1+1:N) with A(I1+1:N, I2)
276*
277 i1 = i1+j-1
278 i2 = i2+j-1
279 CALL cswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
280 $ a( j1+i1, i2 ), 1 )
281 CALL clacgv( i2-i1, a( j1+i1-1, i1+1 ), lda )
282 CALL clacgv( i2-i1-1, a( j1+i1, i2 ), 1 )
283*
284* Swap A(I1, I2+1:N) with A(I2, I2+1:N)
285*
286 IF( i2.LT.m )
287 $ CALL cswap( m-i2, a( j1+i1-1, i2+1 ), lda,
288 $ a( j1+i2-1, i2+1 ), lda )
289*
290* Swap A(I1, I1) with A(I2,I2)
291*
292 piv = a( i1+j1-1, i1 )
293 a( j1+i1-1, i1 ) = a( j1+i2-1, i2 )
294 a( j1+i2-1, i2 ) = piv
295*
296* Swap H(I1, 1:J1) with H(I2, 1:J1)
297*
298 CALL cswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
299 ipiv( i1 ) = i2
300*
301 IF( i1.GT.(k1-1) ) THEN
302*
303* Swap L(1:I1-1, I1) with L(1:I1-1, I2),
304* skipping the first column
305*
306 CALL cswap( i1-k1+1, a( 1, i1 ), 1,
307 $ a( 1, i2 ), 1 )
308 END IF
309 ELSE
310 ipiv( j+1 ) = j+1
311 ENDIF
312*
313* Set A(J, J+1) = T(J, J+1)
314*
315 a( k, j+1 ) = work( 2 )
316*
317 IF( j.LT.nb ) THEN
318*
319* Copy A(J+1:N, J+1) into H(J:N, J),
320*
321 CALL ccopy( m-j, a( k+1, j+1 ), lda,
322 $ h( j+1, j+1 ), 1 )
323 END IF
324*
325* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
326* where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
327*
328 IF( j.LT.(m-1) ) THEN
329 IF( a( k, j+1 ).NE.zero ) THEN
330 alpha = one / a( k, j+1 )
331 CALL ccopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
332 CALL cscal( m-j-1, alpha, a( k, j+2 ), lda )
333 ELSE
334 CALL claset( 'Full', 1, m-j-1, zero, zero,
335 $ a( k, j+2 ), lda)
336 END IF
337 END IF
338 END IF
339 j = j + 1
340 GO TO 10
341 20 CONTINUE
342*
343 ELSE
344*
345* .....................................................
346* Factorize A as L*D*L**T using the lower triangle of A
347* .....................................................
348*
349 30 CONTINUE
350 IF( j.GT.min( m, nb ) )
351 $ GO TO 40
352*
353* K is the column to be factorized
354* when being called from CHETRF_AA,
355* > for the first block column, J1 is 1, hence J1+J-1 is J,
356* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
357*
358 k = j1+j-1
359 IF( j.EQ.m ) THEN
360*
361* Only need to compute T(J, J)
362*
363 mj = 1
364 ELSE
365 mj = m-j+1
366 END IF
367*
368* H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
369* where H(J:N, J) has been initialized to be A(J:N, J)
370*
371 IF( k.GT.2 ) THEN
372*
373* K is the column to be factorized
374* > for the first block column, K is J, skipping the first two
375* columns
376* > for the rest of the columns, K is J+1, skipping only the
377* first column
378*
379 CALL clacgv( j-k1, a( j, 1 ), lda )
380 CALL cgemv( 'No transpose', mj, j-k1,
381 $ -one, h( j, k1 ), ldh,
382 $ a( j, 1 ), lda,
383 $ one, h( j, j ), 1 )
384 CALL clacgv( j-k1, a( j, 1 ), lda )
385 END IF
386*
387* Copy H(J:N, J) into WORK
388*
389 CALL ccopy( mj, h( j, j ), 1, work( 1 ), 1 )
390*
391 IF( j.GT.k1 ) THEN
392*
393* Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
394* where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
395*
396 alpha = -conjg( a( j, k-1 ) )
397 CALL caxpy( mj, alpha, a( j, k-2 ), 1, work( 1 ), 1 )
398 END IF
399*
400* Set A(J, J) = T(J, J)
401*
402 a( j, k ) = real( work( 1 ) )
403*
404 IF( j.LT.m ) THEN
405*
406* Compute WORK(2:N) = T(J, J) L((J+1):N, J)
407* where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
408*
409 IF( k.GT.1 ) THEN
410 alpha = -a( j, k )
411 CALL caxpy( m-j, alpha, a( j+1, k-1 ), 1,
412 $ work( 2 ), 1 )
413 ENDIF
414*
415* Find max(|WORK(2:n)|)
416*
417 i2 = icamax( m-j, work( 2 ), 1 ) + 1
418 piv = work( i2 )
419*
420* Apply hermitian pivot
421*
422 IF( (i2.NE.2) .AND. (piv.NE.0) ) THEN
423*
424* Swap WORK(I1) and WORK(I2)
425*
426 i1 = 2
427 work( i2 ) = work( i1 )
428 work( i1 ) = piv
429*
430* Swap A(I1+1:N, I1) with A(I2, I1+1:N)
431*
432 i1 = i1+j-1
433 i2 = i2+j-1
434 CALL cswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
435 $ a( i2, j1+i1 ), lda )
436 CALL clacgv( i2-i1, a( i1+1, j1+i1-1 ), 1 )
437 CALL clacgv( i2-i1-1, a( i2, j1+i1 ), lda )
438*
439* Swap A(I2+1:N, I1) with A(I2+1:N, I2)
440*
441 IF( i2.LT.m )
442 $ CALL cswap( m-i2, a( i2+1, j1+i1-1 ), 1,
443 $ a( i2+1, j1+i2-1 ), 1 )
444*
445* Swap A(I1, I1) with A(I2, I2)
446*
447 piv = a( i1, j1+i1-1 )
448 a( i1, j1+i1-1 ) = a( i2, j1+i2-1 )
449 a( i2, j1+i2-1 ) = piv
450*
451* Swap H(I1, I1:J1) with H(I2, I2:J1)
452*
453 CALL cswap( i1-1, h( i1, 1 ), ldh, h( i2, 1 ), ldh )
454 ipiv( i1 ) = i2
455*
456 IF( i1.GT.(k1-1) ) THEN
457*
458* Swap L(1:I1-1, I1) with L(1:I1-1, I2),
459* skipping the first column
460*
461 CALL cswap( i1-k1+1, a( i1, 1 ), lda,
462 $ a( i2, 1 ), lda )
463 END IF
464 ELSE
465 ipiv( j+1 ) = j+1
466 ENDIF
467*
468* Set A(J+1, J) = T(J+1, J)
469*
470 a( j+1, k ) = work( 2 )
471*
472 IF( j.LT.nb ) THEN
473*
474* Copy A(J+1:N, J+1) into H(J+1:N, J),
475*
476 CALL ccopy( m-j, a( j+1, k+1 ), 1,
477 $ h( j+1, j+1 ), 1 )
478 END IF
479*
480* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
481* where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
482*
483 IF( j.LT.(m-1) ) THEN
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 END IF
494 j = j + 1
495 GO TO 30
496 40 CONTINUE
497 END IF
498 RETURN
499*
500* End of CLAHEF_AA
501*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
icamax
integer function icamax(n, cx, incx)
ICAMAX
Definition icamax.f:71
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
clacgv
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:72
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:104
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: