LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dsytri2x()

subroutine dsytri2x ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
double precision, dimension( n+nb+1,* )  WORK,
integer  NB,
integer  INFO 
)

DSYTRI2X

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

Purpose:
 DSYTRI2X computes the inverse of a real symmetric indefinite matrix
 A using the factorization A = U*D*U**T or A = L*D*L**T computed by
 DSYTRF.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the details of the factorization are stored
          as an upper or lower triangular matrix.
          = 'U':  Upper triangular, form is A = U*D*U**T;
          = 'L':  Lower triangular, form is A = L*D*L**T.
[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 NNB diagonal matrix D and the multipliers
          used to obtain the factor U or L as computed by DSYTRF.

          On exit, if INFO = 0, the (symmetric) inverse of the original
          matrix.  If UPLO = 'U', the upper triangular part of the
          inverse is formed and the part of A below the diagonal is not
          referenced; if UPLO = 'L' the lower triangular part of the
          inverse is formed and the part of A above the diagonal is
          not referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]IPIV
          IPIV is INTEGER array, dimension (N)
          Details of the interchanges and the NNB structure of D
          as determined by DSYTRF.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (N+NB+1,NB+3)
[in]NB
          NB is INTEGER
          Block size
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
               inverse could not be computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2017

Definition at line 122 of file dsytri2x.f.

122 *
123 * -- LAPACK computational routine (version 3.7.1) --
124 * -- LAPACK is a software package provided by Univ. of Tennessee, --
125 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126 * June 2017
127 *
128 * .. Scalar Arguments ..
129  CHARACTER uplo
130  INTEGER info, lda, n, nb
131 * ..
132 * .. Array Arguments ..
133  INTEGER ipiv( * )
134  DOUBLE PRECISION a( lda, * ), work( n+nb+1,* )
135 * ..
136 *
137 * =====================================================================
138 *
139 * .. Parameters ..
140  DOUBLE PRECISION one, zero
141  parameter( one = 1.0d+0, zero = 0.0d+0 )
142 * ..
143 * .. Local Scalars ..
144  LOGICAL upper
145  INTEGER i, iinfo, ip, k, cut, nnb
146  INTEGER count
147  INTEGER j, u11, invd
148 
149  DOUBLE PRECISION ak, akkp1, akp1, d, t
150  DOUBLE PRECISION u01_i_j, u01_ip1_j
151  DOUBLE PRECISION u11_i_j, u11_ip1_j
152 * ..
153 * .. External Functions ..
154  LOGICAL lsame
155  EXTERNAL lsame
156 * ..
157 * .. External Subroutines ..
158  EXTERNAL dsyconv, xerbla, dtrtri
159  EXTERNAL dgemm, dtrmm, dsyswapr
160 * ..
161 * .. Intrinsic Functions ..
162  INTRINSIC max
163 * ..
164 * .. Executable Statements ..
165 *
166 * Test the input parameters.
167 *
168  info = 0
169  upper = lsame( uplo, 'U' )
170  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
171  info = -1
172  ELSE IF( n.LT.0 ) THEN
173  info = -2
174  ELSE IF( lda.LT.max( 1, n ) ) THEN
175  info = -4
176  END IF
177 *
178 * Quick return if possible
179 *
180 *
181  IF( info.NE.0 ) THEN
182  CALL xerbla( 'DSYTRI2X', -info )
183  RETURN
184  END IF
185  IF( n.EQ.0 )
186  $ RETURN
187 *
188 * Convert A
189 * Workspace got Non-diag elements of D
190 *
191  CALL dsyconv( uplo, 'C', n, a, lda, ipiv, work, iinfo )
192 *
193 * Check that the diagonal matrix D is nonsingular.
194 *
195  IF( upper ) THEN
196 *
197 * Upper triangular storage: examine D from bottom to top
198 *
199  DO info = n, 1, -1
200  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
201  $ RETURN
202  END DO
203  ELSE
204 *
205 * Lower triangular storage: examine D from top to bottom.
206 *
207  DO info = 1, n
208  IF( ipiv( info ).GT.0 .AND. a( info, info ).EQ.zero )
209  $ RETURN
210  END DO
211  END IF
212  info = 0
213 *
214 * Splitting Workspace
215 * U01 is a block (N,NB+1)
216 * The first element of U01 is in WORK(1,1)
217 * U11 is a block (NB+1,NB+1)
218 * The first element of U11 is in WORK(N+1,1)
219  u11 = n
220 * INVD is a block (N,2)
221 * The first element of INVD is in WORK(1,INVD)
222  invd = nb+2
223 
224  IF( upper ) THEN
225 *
226 * invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
227 *
228  CALL dtrtri( uplo, 'U', n, a, lda, info )
229 *
230 * inv(D) and inv(D)*inv(U)
231 *
232  k=1
233  DO WHILE ( k .LE. n )
234  IF( ipiv( k ).GT.0 ) THEN
235 * 1 x 1 diagonal NNB
236  work(k,invd) = one / a( k, k )
237  work(k,invd+1) = 0
238  k=k+1
239  ELSE
240 * 2 x 2 diagonal NNB
241  t = work(k+1,1)
242  ak = a( k, k ) / t
243  akp1 = a( k+1, k+1 ) / t
244  akkp1 = work(k+1,1) / t
245  d = t*( ak*akp1-one )
246  work(k,invd) = akp1 / d
247  work(k+1,invd+1) = ak / d
248  work(k,invd+1) = -akkp1 / d
249  work(k+1,invd) = -akkp1 / d
250  k=k+2
251  END IF
252  END DO
253 *
254 * inv(U**T) = (inv(U))**T
255 *
256 * inv(U**T)*inv(D)*inv(U)
257 *
258  cut=n
259  DO WHILE (cut .GT. 0)
260  nnb=nb
261  IF (cut .LE. nnb) THEN
262  nnb=cut
263  ELSE
264  count = 0
265 * count negative elements,
266  DO i=cut+1-nnb,cut
267  IF (ipiv(i) .LT. 0) count=count+1
268  END DO
269 * need a even number for a clear cut
270  IF (mod(count,2) .EQ. 1) nnb=nnb+1
271  END IF
272 
273  cut=cut-nnb
274 *
275 * U01 Block
276 *
277  DO i=1,cut
278  DO j=1,nnb
279  work(i,j)=a(i,cut+j)
280  END DO
281  END DO
282 *
283 * U11 Block
284 *
285  DO i=1,nnb
286  work(u11+i,i)=one
287  DO j=1,i-1
288  work(u11+i,j)=zero
289  END DO
290  DO j=i+1,nnb
291  work(u11+i,j)=a(cut+i,cut+j)
292  END DO
293  END DO
294 *
295 * invD*U01
296 *
297  i=1
298  DO WHILE (i .LE. cut)
299  IF (ipiv(i) > 0) THEN
300  DO j=1,nnb
301  work(i,j)=work(i,invd)*work(i,j)
302  END DO
303  i=i+1
304  ELSE
305  DO j=1,nnb
306  u01_i_j = work(i,j)
307  u01_ip1_j = work(i+1,j)
308  work(i,j)=work(i,invd)*u01_i_j+
309  $ work(i,invd+1)*u01_ip1_j
310  work(i+1,j)=work(i+1,invd)*u01_i_j+
311  $ work(i+1,invd+1)*u01_ip1_j
312  END DO
313  i=i+2
314  END IF
315  END DO
316 *
317 * invD1*U11
318 *
319  i=1
320  DO WHILE (i .LE. nnb)
321  IF (ipiv(cut+i) > 0) THEN
322  DO j=i,nnb
323  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
324  END DO
325  i=i+1
326  ELSE
327  DO j=i,nnb
328  u11_i_j = work(u11+i,j)
329  u11_ip1_j = work(u11+i+1,j)
330  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
331  $ work(cut+i,invd+1)*work(u11+i+1,j)
332  work(u11+i+1,j)=work(cut+i+1,invd)*u11_i_j+
333  $ work(cut+i+1,invd+1)*u11_ip1_j
334  END DO
335  i=i+2
336  END IF
337  END DO
338 *
339 * U11**T*invD1*U11->U11
340 *
341  CALL dtrmm('L','U','T','U',nnb, nnb,
342  $ one,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
343 *
344  DO i=1,nnb
345  DO j=i,nnb
346  a(cut+i,cut+j)=work(u11+i,j)
347  END DO
348  END DO
349 *
350 * U01**T*invD*U01->A(CUT+I,CUT+J)
351 *
352  CALL dgemm('T','N',nnb,nnb,cut,one,a(1,cut+1),lda,
353  $ work,n+nb+1, zero, work(u11+1,1), n+nb+1)
354 
355 *
356 * U11 = U11**T*invD1*U11 + U01**T*invD*U01
357 *
358  DO i=1,nnb
359  DO j=i,nnb
360  a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
361  END DO
362  END DO
363 *
364 * U01 = U00**T*invD0*U01
365 *
366  CALL dtrmm('L',uplo,'T','U',cut, nnb,
367  $ one,a,lda,work,n+nb+1)
368 
369 *
370 * Update U01
371 *
372  DO i=1,cut
373  DO j=1,nnb
374  a(i,cut+j)=work(i,j)
375  END DO
376  END DO
377 *
378 * Next Block
379 *
380  END DO
381 *
382 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
383 *
384  i=1
385  DO WHILE ( i .LE. n )
386  IF( ipiv(i) .GT. 0 ) THEN
387  ip=ipiv(i)
388  IF (i .LT. ip) CALL dsyswapr( uplo, n, a, lda, i ,ip )
389  IF (i .GT. ip) CALL dsyswapr( uplo, n, a, lda, ip ,i )
390  ELSE
391  ip=-ipiv(i)
392  i=i+1
393  IF ( (i-1) .LT. ip)
394  $ CALL dsyswapr( uplo, n, a, lda, i-1 ,ip )
395  IF ( (i-1) .GT. ip)
396  $ CALL dsyswapr( uplo, n, a, lda, ip ,i-1 )
397  ENDIF
398  i=i+1
399  END DO
400  ELSE
401 *
402 * LOWER...
403 *
404 * invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
405 *
406  CALL dtrtri( uplo, 'U', n, a, lda, info )
407 *
408 * inv(D) and inv(D)*inv(U)
409 *
410  k=n
411  DO WHILE ( k .GE. 1 )
412  IF( ipiv( k ).GT.0 ) THEN
413 * 1 x 1 diagonal NNB
414  work(k,invd) = one / a( k, k )
415  work(k,invd+1) = 0
416  k=k-1
417  ELSE
418 * 2 x 2 diagonal NNB
419  t = work(k-1,1)
420  ak = a( k-1, k-1 ) / t
421  akp1 = a( k, k ) / t
422  akkp1 = work(k-1,1) / t
423  d = t*( ak*akp1-one )
424  work(k-1,invd) = akp1 / d
425  work(k,invd) = ak / d
426  work(k,invd+1) = -akkp1 / d
427  work(k-1,invd+1) = -akkp1 / d
428  k=k-2
429  END IF
430  END DO
431 *
432 * inv(U**T) = (inv(U))**T
433 *
434 * inv(U**T)*inv(D)*inv(U)
435 *
436  cut=0
437  DO WHILE (cut .LT. n)
438  nnb=nb
439  IF (cut + nnb .GT. n) THEN
440  nnb=n-cut
441  ELSE
442  count = 0
443 * count negative elements,
444  DO i=cut+1,cut+nnb
445  IF (ipiv(i) .LT. 0) count=count+1
446  END DO
447 * need a even number for a clear cut
448  IF (mod(count,2) .EQ. 1) nnb=nnb+1
449  END IF
450 * L21 Block
451  DO i=1,n-cut-nnb
452  DO j=1,nnb
453  work(i,j)=a(cut+nnb+i,cut+j)
454  END DO
455  END DO
456 * L11 Block
457  DO i=1,nnb
458  work(u11+i,i)=one
459  DO j=i+1,nnb
460  work(u11+i,j)=zero
461  END DO
462  DO j=1,i-1
463  work(u11+i,j)=a(cut+i,cut+j)
464  END DO
465  END DO
466 *
467 * invD*L21
468 *
469  i=n-cut-nnb
470  DO WHILE (i .GE. 1)
471  IF (ipiv(cut+nnb+i) > 0) THEN
472  DO j=1,nnb
473  work(i,j)=work(cut+nnb+i,invd)*work(i,j)
474  END DO
475  i=i-1
476  ELSE
477  DO j=1,nnb
478  u01_i_j = work(i,j)
479  u01_ip1_j = work(i-1,j)
480  work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
481  $ work(cut+nnb+i,invd+1)*u01_ip1_j
482  work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
483  $ work(cut+nnb+i-1,invd)*u01_ip1_j
484  END DO
485  i=i-2
486  END IF
487  END DO
488 *
489 * invD1*L11
490 *
491  i=nnb
492  DO WHILE (i .GE. 1)
493  IF (ipiv(cut+i) > 0) THEN
494  DO j=1,nnb
495  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
496  END DO
497  i=i-1
498  ELSE
499  DO j=1,nnb
500  u11_i_j = work(u11+i,j)
501  u11_ip1_j = work(u11+i-1,j)
502  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
503  $ work(cut+i,invd+1)*u11_ip1_j
504  work(u11+i-1,j)=work(cut+i-1,invd+1)*u11_i_j+
505  $ work(cut+i-1,invd)*u11_ip1_j
506  END DO
507  i=i-2
508  END IF
509  END DO
510 *
511 * L11**T*invD1*L11->L11
512 *
513  CALL dtrmm('L',uplo,'T','U',nnb, nnb,
514  $ one,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
515 
516 *
517  DO i=1,nnb
518  DO j=1,i
519  a(cut+i,cut+j)=work(u11+i,j)
520  END DO
521  END DO
522 *
523  IF ( (cut+nnb) .LT. n ) THEN
524 *
525 * L21**T*invD2*L21->A(CUT+I,CUT+J)
526 *
527  CALL dgemm('T','N',nnb,nnb,n-nnb-cut,one,a(cut+nnb+1,cut+1)
528  $ ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
529 
530 *
531 * L11 = L11**T*invD1*L11 + U01**T*invD*U01
532 *
533  DO i=1,nnb
534  DO j=1,i
535  a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
536  END DO
537  END DO
538 *
539 * L01 = L22**T*invD2*L21
540 *
541  CALL dtrmm('L',uplo,'T','U', n-nnb-cut, nnb,
542  $ one,a(cut+nnb+1,cut+nnb+1),lda,work,n+nb+1)
543 *
544 * Update L21
545 *
546  DO i=1,n-cut-nnb
547  DO j=1,nnb
548  a(cut+nnb+i,cut+j)=work(i,j)
549  END DO
550  END DO
551 
552  ELSE
553 *
554 * L11 = L11**T*invD1*L11
555 *
556  DO i=1,nnb
557  DO j=1,i
558  a(cut+i,cut+j)=work(u11+i,j)
559  END DO
560  END DO
561  END IF
562 *
563 * Next Block
564 *
565  cut=cut+nnb
566  END DO
567 *
568 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
569 *
570  i=n
571  DO WHILE ( i .GE. 1 )
572  IF( ipiv(i) .GT. 0 ) THEN
573  ip=ipiv(i)
574  IF (i .LT. ip) CALL dsyswapr( uplo, n, a, lda, i ,ip )
575  IF (i .GT. ip) CALL dsyswapr( uplo, n, a, lda, ip ,i )
576  ELSE
577  ip=-ipiv(i)
578  IF ( i .LT. ip) CALL dsyswapr( uplo, n, a, lda, i ,ip )
579  IF ( i .GT. ip) CALL dsyswapr( uplo, n, a, lda, ip, i )
580  i=i-1
581  ENDIF
582  i=i-1
583  END DO
584  END IF
585 *
586  RETURN
587 *
588 * End of DSYTRI2X
589 *
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRMM
Definition: dtrmm.f:179
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine dsyswapr(UPLO, N, A, LDA, I1, I2)
DSYSWAPR applies an elementary permutation on the rows and columns of a symmetric matrix...
Definition: dsyswapr.f:104
subroutine dsyconv(UPLO, WAY, N, A, LDA, IPIV, E, INFO)
DSYCONV
Definition: dsyconv.f:116
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dtrtri(UPLO, DIAG, N, A, LDA, INFO)
DTRTRI
Definition: dtrtri.f:111
Here is the call graph for this function:
Here is the caller graph for this function: