LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
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).`

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: