LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ zlahef_aa()

 subroutine zlahef_aa ( character UPLO, integer J1, integer M, integer NB, complex*16, dimension( lda, * ) A, integer LDA, integer, dimension( * ) IPIV, complex*16, dimension( ldh, * ) H, integer LDH, complex*16, dimension( * ) WORK )

ZLAHEF_AA

Purpose:
``` DLAHEF_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 ZHETRF_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*16 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*16 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*16 workspace, dimension (M).`

Definition at line 142 of file zlahef_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*16 A( LDA, * ), H( LDH, * ), WORK( * )
158* ..
159*
160* =====================================================================
161* .. Parameters ..
162 COMPLEX*16 ZERO, ONE
163 parameter( zero = (0.0d+0, 0.0d+0), one = (1.0d+0, 0.0d+0) )
164*
165* .. Local Scalars ..
166 INTEGER J, K, K1, I1, I2, MJ
167 COMPLEX*16 PIV, ALPHA
168* ..
169* .. External Functions ..
170 LOGICAL LSAME
171 INTEGER IZAMAX, ILAENV
172 EXTERNAL lsame, ilaenv, izamax
173* ..
174* .. External Subroutines ..
175 EXTERNAL zgemm, zgemv, zaxpy, zlacgv, zcopy, zscal, zswap,
176 \$ zlaset, xerbla
177* ..
178* .. Intrinsic Functions ..
179 INTRINSIC dble, dconjg, 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 ZHETRF_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 zlacgv( j-k1, a( 1, j ), 1 )
227 CALL zgemv( 'No transpose', mj, j-k1,
228 \$ -one, h( j, k1 ), ldh,
229 \$ a( 1, j ), 1,
230 \$ one, h( j, j ), 1 )
231 CALL zlacgv( j-k1, a( 1, j ), 1 )
232 END IF
233*
234* Copy H(i:n, i) into WORK
235*
236 CALL zcopy( 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 = -dconjg( a( k-1, j ) )
244 CALL zaxpy( 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 ) = dble( 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 zaxpy( 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 = izamax( 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 zswap( i2-i1-1, a( j1+i1-1, i1+1 ), lda,
282 \$ a( j1+i1, i2 ), 1 )
283 CALL zlacgv( i2-i1, a( j1+i1-1, i1+1 ), lda )
284 CALL zlacgv( 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 zswap( 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 zswap( 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 zswap( 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 zcopy( 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 zcopy( m-j-1, work( 3 ), 1, a( k, j+2 ), lda )
334 CALL zscal( m-j-1, alpha, a( k, j+2 ), lda )
335 ELSE
336 CALL zlaset( '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 ZHETRF_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 zlacgv( j-k1, a( j, 1 ), lda )
382 CALL zgemv( 'No transpose', mj, j-k1,
383 \$ -one, h( j, k1 ), ldh,
384 \$ a( j, 1 ), lda,
385 \$ one, h( j, j ), 1 )
386 CALL zlacgv( j-k1, a( j, 1 ), lda )
387 END IF
388*
389* Copy H(J:N, J) into WORK
390*
391 CALL zcopy( 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 = -dconjg( a( j, k-1 ) )
399 CALL zaxpy( 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 ) = dble( 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 zaxpy( 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 = izamax( 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 zswap( i2-i1-1, a( i1+1, j1+i1-1 ), 1,
437 \$ a( i2, j1+i1 ), lda )
438 CALL zlacgv( i2-i1, a( i1+1, j1+i1-1 ), 1 )
439 CALL zlacgv( 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 zswap( 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 zswap( 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 zswap( 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 zcopy( 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 zcopy( m-j-1, work( 3 ), 1, a( j+2, k ), 1 )
489 CALL zscal( m-j-1, alpha, a( j+2, k ), 1 )
490 ELSE
491 CALL zlaset( '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 ZLAHEF_AA
503*
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:74
Here is the call graph for this function:
Here is the caller graph for this function: