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

◆ dsytrf_aa()

subroutine dsytrf_aa ( character uplo,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
double precision, dimension( * ) work,
integer lwork,
integer info )

DSYTRF_AA

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

Purpose:
!>
!> DSYTRF_AA computes the factorization of a real symmetric matrix A
!> using the Aasen's algorithm.  The form of the factorization is
!>
!>    A = U**T*T*U  or  A = L*T*L**T
!>
!> where U (or L) is a product of permutation and unit upper (lower)
!> triangular matrices, and T is a symmetric 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 DOUBLE PRECISION array, dimension (LDA,N)
!>          On entry, the symmetric 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 DOUBLE PRECISION 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 dsytrf_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 DOUBLE PRECISION A( LDA, * ), WORK( * )
147* ..
148*
149* =====================================================================
150* .. Parameters ..
151 DOUBLE PRECISION ZERO, ONE
152 parameter( zero = 0.0d+0, one = 1.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 DOUBLE PRECISION ALPHA
159* ..
160* .. External Functions ..
161 LOGICAL LSAME
162 INTEGER ILAENV
163 EXTERNAL lsame, ilaenv
164* ..
165* .. External Subroutines ..
166 EXTERNAL dlasyf_aa, dgemm, dgemv, dscal, dcopy,
167 $ dswap,
168 $ xerbla
169* ..
170* .. Intrinsic Functions ..
171 INTRINSIC max
172* ..
173* .. Executable Statements ..
174*
175* Determine the block size
176*
177 nb = ilaenv( 1, 'DSYTRF_AA', uplo, n, -1, -1, -1 )
178*
179* Test the input parameters.
180*
181 info = 0
182 upper = lsame( uplo, 'U' )
183 lquery = ( lwork.EQ.-1 )
184 IF( n.LE.1 ) THEN
185 lwkmin = 1
186 lwkopt = 1
187 ELSE
188 lwkmin = 2*n
189 lwkopt = (nb+1)*n
190 END IF
191*
192 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
193 info = -1
194 ELSE IF( n.LT.0 ) THEN
195 info = -2
196 ELSE IF( lda.LT.max( 1, n ) ) THEN
197 info = -4
198 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
199 info = -7
200 END IF
201*
202 IF( info.EQ.0 ) THEN
203 work( 1 ) = lwkopt
204 END IF
205*
206 IF( info.NE.0 ) THEN
207 CALL xerbla( 'DSYTRF_AA', -info )
208 RETURN
209 ELSE IF( lquery ) THEN
210 RETURN
211 END IF
212*
213* Quick return
214*
215 IF( n.EQ.0 ) THEN
216 RETURN
217 ENDIF
218 ipiv( 1 ) = 1
219 IF( n.EQ.1 ) THEN
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**T*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 dcopy( 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 DLASYF;
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 dlasyf_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 dswap( 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 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 = a( j, j+1 )
289 a( j, j+1 ) = one
290 CALL dcopy( n-j, a( j-1, j+1 ), lda,
291 $ work( (j+1-j1+1)+jb*n ), 1 )
292 CALL dscal( 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=1 and K2= 0 for the first panel,
296* while K1=0 and K2=1 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 DGEMV
318*
319 j3 = j2
320 DO mj = nj-1, 1, -1
321 CALL dgemv( 'No transpose', mj, jb+1,
322 $ -one, work( j3-j1+1+k1*n ), n,
323 $ a( j1-k2, j3 ), 1,
324 $ one, a( j3, j3 ), lda )
325 j3 = j3 + 1
326 END DO
327*
328* Update off-diagonal block of J2-th block row with DGEMM
329*
330 CALL dgemm( 'Transpose', 'Transpose',
331 $ nj, n-j3+1, jb+1,
332 $ -one, a( j1-k2, j2 ), lda,
333 $ work( j3-j1+1+k1*n ), n,
334 $ one, a( j2, j3 ), lda )
335 END DO
336*
337* Recover T( J, J+1 )
338*
339 a( j, j+1 ) = alpha
340 END IF
341*
342* WORK(J+1, 1) stores H(J+1, 1)
343*
344 CALL dcopy( n-j, a( j+1, j+1 ), lda, work( 1 ), 1 )
345 END IF
346 GO TO 10
347 ELSE
348*
349* .....................................................
350* Factorize A as L*D*L**T using the lower triangle of A
351* .....................................................
352*
353* copy first column A(1:N, 1) into H(1:N, 1)
354* (stored in WORK(1:N))
355*
356 CALL dcopy( n, a( 1, 1 ), 1, work( 1 ), 1 )
357*
358* J is the main loop index, increasing from 1 to N in steps of
359* JB, where JB is the number of columns factorized by DLASYF;
360* JB is either NB, or N-J+1 for the last block
361*
362 j = 0
363 11 CONTINUE
364 IF( j.GE.n )
365 $ GO TO 20
366*
367* each step of the main loop
368* J is the last column of the previous panel
369* J1 is the first column of the current panel
370* K1 identifies if the previous column of the panel has been
371* explicitly stored, e.g., K1=1 for the first panel, and
372* K1=0 for the rest
373*
374 j1 = j+1
375 jb = min( n-j1+1, nb )
376 k1 = max(1, j)-j
377*
378* Panel factorization
379*
380 CALL dlasyf_aa( uplo, 2-k1, n-j, jb,
381 $ a( j+1, max(1, j) ), lda,
382 $ ipiv( j+1 ), work, n, work( n*nb+1 ) )
383*
384* Adjust IPIV and apply it back (J-th step picks (J+1)-th pivot)
385*
386 DO j2 = j+2, min(n, j+jb+1)
387 ipiv( j2 ) = ipiv( j2 ) + j
388 IF( (j2.NE.ipiv(j2)) .AND. ((j1-k1).GT.2) ) THEN
389 CALL dswap( j1-k1-2, a( j2, 1 ), lda,
390 $ a( ipiv(j2), 1 ), lda )
391 END IF
392 END DO
393 j = j + jb
394*
395* Trailing submatrix update, where
396* A(J2+1, J1-1) stores L(J2+1, J1) and
397* WORK(J2+1, 1) stores H(J2+1, 1)
398*
399 IF( j.LT.n ) THEN
400*
401* if first panel and JB=1 (NB=1), then nothing to do
402*
403 IF( j1.GT.1 .OR. jb.GT.1 ) THEN
404*
405* Merge rank-1 update with BLAS-3 update
406*
407 alpha = a( j+1, j )
408 a( j+1, j ) = one
409 CALL dcopy( n-j, a( j+1, j-1 ), 1,
410 $ work( (j+1-j1+1)+jb*n ), 1 )
411 CALL dscal( n-j, alpha, work( (j+1-j1+1)+jb*n ), 1 )
412*
413* K1 identifies if the previous column of the panel has been
414* explicitly stored, e.g., K1=1 and K2= 0 for the first panel,
415* while K1=0 and K2=1 for the rest
416*
417 IF( j1.GT.1 ) THEN
418*
419* Not first panel
420*
421 k2 = 1
422 ELSE
423*
424* First panel
425*
426 k2 = 0
427*
428* First update skips the first column
429*
430 jb = jb - 1
431 END IF
432*
433 DO j2 = j+1, n, nb
434 nj = min( nb, n-j2+1 )
435*
436* Update (J2, J2) diagonal block with DGEMV
437*
438 j3 = j2
439 DO mj = nj-1, 1, -1
440 CALL dgemv( 'No transpose', mj, jb+1,
441 $ -one, work( j3-j1+1+k1*n ), n,
442 $ a( j3, j1-k2 ), lda,
443 $ one, a( j3, j3 ), 1 )
444 j3 = j3 + 1
445 END DO
446*
447* Update off-diagonal block in J2-th block column with DGEMM
448*
449 CALL dgemm( 'No transpose', 'Transpose',
450 $ n-j3+1, nj, jb+1,
451 $ -one, work( j3-j1+1+k1*n ), n,
452 $ a( j2, j1-k2 ), lda,
453 $ one, a( j3, j2 ), lda )
454 END DO
455*
456* Recover T( J+1, J )
457*
458 a( j+1, j ) = alpha
459 END IF
460*
461* WORK(J+1, 1) stores H(J+1, 1)
462*
463 CALL dcopy( n-j, a( j+1, j+1 ), 1, work( 1 ), 1 )
464 END IF
465 GO TO 11
466 END IF
467*
468 20 CONTINUE
469 work( 1 ) = lwkopt
470 RETURN
471*
472* End of DSYTRF_AA
473*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
dlasyf_aa
subroutine dlasyf_aa(uplo, j1, m, nb, a, lda, ipiv, h, ldh, work)
DLASYF_AA
Definition dlasyf_aa.f:142
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: