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

◆ zhetrf_aa()

subroutine zhetrf_aa ( character uplo,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
complex*16, dimension( * ) work,
integer lwork,
integer info )

ZHETRF_AA

Download ZHETRF_AA + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZHETRF_AA computes the factorization of a complex hermitian matrix A
!> using the Aasen's algorithm.  The form of the factorization is
!>
!>    A = U**H*T*U  or  A = L*T*L**H
!>
!> where U (or L) is a product of permutation and unit upper (lower)
!> triangular matrices, and T is a hermitian tridiagonal matrix.
!>
!> This is the blocked version of the algorithm, calling Level 3 BLAS.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in,out]A
!>          A is COMPLEX*16 array, dimension (LDA,N)
!>          On entry, the hermitian matrix A.  If UPLO = 'U', the leading
!>          N-by-N upper triangular part of A contains the upper
!>          triangular part of the matrix A, and the strictly lower
!>          triangular part of A is not referenced.  If UPLO = 'L', the
!>          leading N-by-N lower triangular part of A contains the lower
!>          triangular part of the matrix A, and the strictly upper
!>          triangular part of A is not referenced.
!>
!>          On exit, the tridiagonal matrix is stored in the diagonals
!>          and the subdiagonals of A just below (or above) the diagonals,
!>          and L is stored below (or above) the subdiagonals, when UPLO
!>          is 'L' (or 'U').
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          On exit, it contains the details of the interchanges, i.e.,
!>          the row and column k of A were interchanged with the
!>          row and column IPIV(k).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The length of WORK.
!>          LWORK >= 1, if N >= 1, and LWORK >= 2*N, otherwise.
!>          For optimum performance LWORK >= N*(1+NB), where NB is
!>          the optimal blocksize, returned by ILAENV.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 131 of file zhetrf_aa.f.

133*
134* -- LAPACK computational routine --
135* -- LAPACK is a software package provided by Univ. of Tennessee, --
136* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
137*
138 IMPLICIT NONE
139*
140* .. Scalar Arguments ..
141 CHARACTER UPLO
142 INTEGER N, LDA, LWORK, INFO
143* ..
144* .. Array Arguments ..
145 INTEGER IPIV( * )
146 COMPLEX*16 A( LDA, * ), WORK( * )
147* ..
148*
149* =====================================================================
150* .. Parameters ..
151 COMPLEX*16 ZERO, ONE
152 parameter( zero = (0.0d+0, 0.0d+0), one = (1.0d+0, 0.0d+0) )
153*
154* .. Local Scalars ..
155 LOGICAL LQUERY, UPPER
156 INTEGER J, LWKMIN, LWKOPT
157 INTEGER NB, MJ, NJ, K1, K2, J1, J2, J3, JB
158 COMPLEX*16 ALPHA
159* ..
160* .. External Functions ..
161 LOGICAL LSAME
162 INTEGER ILAENV
163 EXTERNAL lsame, ilaenv
164* ..
165* .. External Subroutines ..
166 EXTERNAL zlahef_aa, zgemm, zgemv, zcopy, zscal, zswap,
167 $ xerbla
168* ..
169* .. Intrinsic Functions ..
170 INTRINSIC dble, dconjg, max
171* ..
172* .. Executable Statements ..
173*
174* Determine the block size
175*
176 nb = ilaenv( 1, 'ZHETRF_AA', uplo, n, -1, -1, -1 )
177*
178* Test the input parameters.
179*
180 info = 0
181 upper = lsame( uplo, 'U' )
182 lquery = ( lwork.EQ.-1 )
183 IF( n.LE.1 ) THEN
184 lwkmin = 1
185 lwkopt = 1
186 ELSE
187 lwkmin = 2*n
188 lwkopt = (nb+1)*n
189 END IF
190*
191 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
192 info = -1
193 ELSE IF( n.LT.0 ) THEN
194 info = -2
195 ELSE IF( lda.LT.max( 1, n ) ) THEN
196 info = -4
197 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
198 info = -7
199 END IF
200*
201 IF( info.EQ.0 ) THEN
202 work( 1 ) = lwkopt
203 END IF
204*
205 IF( info.NE.0 ) THEN
206 CALL xerbla( 'ZHETRF_AA', -info )
207 RETURN
208 ELSE IF( lquery ) THEN
209 RETURN
210 END IF
211*
212* Quick return
213*
214 IF( n.EQ.0 ) THEN
215 RETURN
216 ENDIF
217 ipiv( 1 ) = 1
218 IF( n.EQ.1 ) THEN
219 a( 1, 1 ) = dble( a( 1, 1 ) )
220 RETURN
221 END IF
222*
223* Adjust block size based on the workspace size
224*
225 IF( lwork.LT.((1+nb)*n) ) THEN
226 nb = ( lwork-n ) / n
227 END IF
228*
229 IF( upper ) THEN
230*
231* .....................................................
232* Factorize A as U**H*D*U using the upper triangle of A
233* .....................................................
234*
235* copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N))
236*
237 CALL zcopy( n, a( 1, 1 ), lda, work( 1 ), 1 )
238*
239* J is the main loop index, increasing from 1 to N in steps of
240* JB, where JB is the number of columns factorized by ZLAHEF;
241* JB is either NB, or N-J+1 for the last block
242*
243 j = 0
244 10 CONTINUE
245 IF( j.GE.n )
246 $ GO TO 20
247*
248* each step of the main loop
249* J is the last column of the previous panel
250* J1 is the first column of the current panel
251* K1 identifies if the previous column of the panel has been
252* explicitly stored, e.g., K1=1 for the first panel, and
253* K1=0 for the rest
254*
255 j1 = j + 1
256 jb = min( n-j1+1, nb )
257 k1 = max(1, j)-j
258*
259* Panel factorization
260*
261 CALL zlahef_aa( uplo, 2-k1, n-j, jb,
262 $ a( max(1, j), j+1 ), lda,
263 $ ipiv( j+1 ), work, n, work( n*nb+1 ) )
264*
265* Adjust IPIV and apply it back (J-th step picks (J+1)-th pivot)
266*
267 DO j2 = j+2, min(n, j+jb+1)
268 ipiv( j2 ) = ipiv( j2 ) + j
269 IF( (j2.NE.ipiv(j2)) .AND. ((j1-k1).GT.2) ) THEN
270 CALL zswap( j1-k1-2, a( 1, j2 ), 1,
271 $ a( 1, ipiv(j2) ), 1 )
272 END IF
273 END DO
274 j = j + jb
275*
276* Trailing submatrix update, where
277* the row A(J1-1, J2-1:N) stores U(J1, J2+1:N) and
278* WORK stores the current block of the auxiriarly matrix H
279*
280 IF( j.LT.n ) THEN
281*
282* if the first panel and JB=1 (NB=1), then nothing to do
283*
284 IF( j1.GT.1 .OR. jb.GT.1 ) THEN
285*
286* Merge rank-1 update with BLAS-3 update
287*
288 alpha = dconjg( a( j, j+1 ) )
289 a( j, j+1 ) = one
290 CALL zcopy( n-j, a( j-1, j+1 ), lda,
291 $ work( (j+1-j1+1)+jb*n ), 1 )
292 CALL zscal( n-j, alpha, work( (j+1-j1+1)+jb*n ), 1 )
293*
294* K1 identifies if the previous column of the panel has been
295* explicitly stored, e.g., K1=0 and K2=1 for the first panel,
296* and K1=1 and K2=0 for the rest
297*
298 IF( j1.GT.1 ) THEN
299*
300* Not first panel
301*
302 k2 = 1
303 ELSE
304*
305* First panel
306*
307 k2 = 0
308*
309* First update skips the first column
310*
311 jb = jb - 1
312 END IF
313*
314 DO j2 = j+1, n, nb
315 nj = min( nb, n-j2+1 )
316*
317* Update (J2, J2) diagonal block with ZGEMV
318*
319 j3 = j2
320 DO mj = nj-1, 1, -1
321 CALL zgemm( 'Conjugate transpose', 'Transpose',
322 $ 1, mj, jb+1,
323 $ -one, a( j1-k2, j3 ), lda,
324 $ work( (j3-j1+1)+k1*n ), n,
325 $ one, a( j3, j3 ), lda )
326 j3 = j3 + 1
327 END DO
328*
329* Update off-diagonal block of J2-th block row with ZGEMM
330*
331 CALL zgemm( 'Conjugate transpose', 'Transpose',
332 $ nj, n-j3+1, jb+1,
333 $ -one, a( j1-k2, j2 ), lda,
334 $ work( (j3-j1+1)+k1*n ), n,
335 $ one, a( j2, j3 ), lda )
336 END DO
337*
338* Recover T( J, J+1 )
339*
340 a( j, j+1 ) = dconjg( alpha )
341 END IF
342*
343* WORK(J+1, 1) stores H(J+1, 1)
344*
345 CALL zcopy( n-j, a( j+1, j+1 ), lda, work( 1 ), 1 )
346 END IF
347 GO TO 10
348 ELSE
349*
350* .....................................................
351* Factorize A as L*D*L**H using the lower triangle of A
352* .....................................................
353*
354* copy first column A(1:N, 1) into H(1:N, 1)
355* (stored in WORK(1:N))
356*
357 CALL zcopy( n, a( 1, 1 ), 1, work( 1 ), 1 )
358*
359* J is the main loop index, increasing from 1 to N in steps of
360* JB, where JB is the number of columns factorized by ZLAHEF;
361* JB is either NB, or N-J+1 for the last block
362*
363 j = 0
364 11 CONTINUE
365 IF( j.GE.n )
366 $ GO TO 20
367*
368* each step of the main loop
369* J is the last column of the previous panel
370* J1 is the first column of the current panel
371* K1 identifies if the previous column of the panel has been
372* explicitly stored, e.g., K1=1 for the first panel, and
373* K1=0 for the rest
374*
375 j1 = j+1
376 jb = min( n-j1+1, nb )
377 k1 = max(1, j)-j
378*
379* Panel factorization
380*
381 CALL zlahef_aa( uplo, 2-k1, n-j, jb,
382 $ a( j+1, max(1, j) ), lda,
383 $ ipiv( j+1 ), work, n, work( n*nb+1 ) )
384*
385* Adjust IPIV and apply it back (J-th step picks (J+1)-th pivot)
386*
387 DO j2 = j+2, min(n, j+jb+1)
388 ipiv( j2 ) = ipiv( j2 ) + j
389 IF( (j2.NE.ipiv(j2)) .AND. ((j1-k1).GT.2) ) THEN
390 CALL zswap( j1-k1-2, a( j2, 1 ), lda,
391 $ a( ipiv(j2), 1 ), lda )
392 END IF
393 END DO
394 j = j + jb
395*
396* Trailing submatrix update, where
397* A(J2+1, J1-1) stores L(J2+1, J1) and
398* WORK(J2+1, 1) stores H(J2+1, 1)
399*
400 IF( j.LT.n ) THEN
401*
402* if the first panel and JB=1 (NB=1), then nothing to do
403*
404 IF( j1.GT.1 .OR. jb.GT.1 ) THEN
405*
406* Merge rank-1 update with BLAS-3 update
407*
408 alpha = dconjg( a( j+1, j ) )
409 a( j+1, j ) = one
410 CALL zcopy( n-j, a( j+1, j-1 ), 1,
411 $ work( (j+1-j1+1)+jb*n ), 1 )
412 CALL zscal( n-j, alpha, work( (j+1-j1+1)+jb*n ), 1 )
413*
414* K1 identifies if the previous column of the panel has been
415* explicitly stored, e.g., K1=0 and K2=1 for the first panel,
416* and K1=1 and K2=0 for the rest
417*
418 IF( j1.GT.1 ) THEN
419*
420* Not first panel
421*
422 k2 = 1
423 ELSE
424*
425* First panel
426*
427 k2 = 0
428*
429* First update skips the first column
430*
431 jb = jb - 1
432 END IF
433*
434 DO j2 = j+1, n, nb
435 nj = min( nb, n-j2+1 )
436*
437* Update (J2, J2) diagonal block with ZGEMV
438*
439 j3 = j2
440 DO mj = nj-1, 1, -1
441 CALL zgemm( 'No transpose',
442 $ 'Conjugate transpose',
443 $ mj, 1, jb+1,
444 $ -one, work( (j3-j1+1)+k1*n ), n,
445 $ a( j3, j1-k2 ), lda,
446 $ one, a( j3, j3 ), lda )
447 j3 = j3 + 1
448 END DO
449*
450* Update off-diagonal block of J2-th block column with ZGEMM
451*
452 CALL zgemm( 'No transpose', 'Conjugate transpose',
453 $ n-j3+1, nj, jb+1,
454 $ -one, work( (j3-j1+1)+k1*n ), n,
455 $ a( j2, j1-k2 ), lda,
456 $ one, a( j3, j2 ), lda )
457 END DO
458*
459* Recover T( J+1, J )
460*
461 a( j+1, j ) = dconjg( alpha )
462 END IF
463*
464* WORK(J+1, 1) stores H(J+1, 1)
465*
466 CALL zcopy( n-j, a( j+1, j+1 ), 1, work( 1 ), 1 )
467 END IF
468 GO TO 11
469 END IF
470*
471 20 CONTINUE
472 work( 1 ) = lwkopt
473 RETURN
474*
475* End of ZHETRF_AA
476*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
zlahef_aa
subroutine zlahef_aa(uplo, j1, m, nb, a, lda, ipiv, h, ldh, work)
ZLAHEF_AA
Definition zlahef_aa.f:142
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
zscal
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: