LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ ssytri2x()

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

SSYTRI2X

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

Purpose:
 SSYTRI2X 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
 SSYTRF.
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 REAL 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 SSYTRF.

          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 SSYTRF.
[out]WORK
          WORK is REAL 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 ssytri2x.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  REAL a( lda, * ), work( n+nb+1,* )
135 * ..
136 *
137 * =====================================================================
138 *
139 * .. Parameters ..
140  REAL one, zero
141  parameter( one = 1.0e+0, zero = 0.0e+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  REAL ak, akkp1, akp1, d, t
150  REAL u01_i_j, u01_ip1_j
151  REAL u11_i_j, u11_ip1_j
152 * ..
153 * .. External Functions ..
154  LOGICAL lsame
155  EXTERNAL lsame
156 * ..
157 * .. External Subroutines ..
158  EXTERNAL ssyconv, xerbla, strtri
159  EXTERNAL sgemm, strmm, ssyswapr
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( 'SSYTRI2X', -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 ssyconv( 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 strtri( 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 strmm('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 sgemm('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 * U11 = U11**T*invD1*U11 + U01**T*invD*U01
356 *
357  DO i=1,nnb
358  DO j=i,nnb
359  a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
360  END DO
361  END DO
362 *
363 * U01 = U00**T*invD0*U01
364 *
365  CALL strmm('L',uplo,'T','U',cut, nnb,
366  $ one,a,lda,work,n+nb+1)
367 
368 *
369 * Update U01
370 *
371  DO i=1,cut
372  DO j=1,nnb
373  a(i,cut+j)=work(i,j)
374  END DO
375  END DO
376 *
377 * Next Block
378 *
379  END DO
380 *
381 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
382 *
383  i=1
384  DO WHILE ( i .LE. n )
385  IF( ipiv(i) .GT. 0 ) THEN
386  ip=ipiv(i)
387  IF (i .LT. ip) CALL ssyswapr( uplo, n, a, lda, i ,ip )
388  IF (i .GT. ip) CALL ssyswapr( uplo, n, a, lda, ip ,i )
389  ELSE
390  ip=-ipiv(i)
391  i=i+1
392  IF ( (i-1) .LT. ip)
393  $ CALL ssyswapr( uplo, n, a, lda, i-1 ,ip )
394  IF ( (i-1) .GT. ip)
395  $ CALL ssyswapr( uplo, n, a, lda, ip ,i-1 )
396  ENDIF
397  i=i+1
398  END DO
399  ELSE
400 *
401 * LOWER...
402 *
403 * invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
404 *
405  CALL strtri( uplo, 'U', n, a, lda, info )
406 *
407 * inv(D) and inv(D)*inv(U)
408 *
409  k=n
410  DO WHILE ( k .GE. 1 )
411  IF( ipiv( k ).GT.0 ) THEN
412 * 1 x 1 diagonal NNB
413  work(k,invd) = one / a( k, k )
414  work(k,invd+1) = 0
415  k=k-1
416  ELSE
417 * 2 x 2 diagonal NNB
418  t = work(k-1,1)
419  ak = a( k-1, k-1 ) / t
420  akp1 = a( k, k ) / t
421  akkp1 = work(k-1,1) / t
422  d = t*( ak*akp1-one )
423  work(k-1,invd) = akp1 / d
424  work(k,invd) = ak / d
425  work(k,invd+1) = -akkp1 / d
426  work(k-1,invd+1) = -akkp1 / d
427  k=k-2
428  END IF
429  END DO
430 *
431 * inv(U**T) = (inv(U))**T
432 *
433 * inv(U**T)*inv(D)*inv(U)
434 *
435  cut=0
436  DO WHILE (cut .LT. n)
437  nnb=nb
438  IF (cut + nnb .GT. n) THEN
439  nnb=n-cut
440  ELSE
441  count = 0
442 * count negative elements,
443  DO i=cut+1,cut+nnb
444  IF (ipiv(i) .LT. 0) count=count+1
445  END DO
446 * need a even number for a clear cut
447  IF (mod(count,2) .EQ. 1) nnb=nnb+1
448  END IF
449 * L21 Block
450  DO i=1,n-cut-nnb
451  DO j=1,nnb
452  work(i,j)=a(cut+nnb+i,cut+j)
453  END DO
454  END DO
455 * L11 Block
456  DO i=1,nnb
457  work(u11+i,i)=one
458  DO j=i+1,nnb
459  work(u11+i,j)=zero
460  END DO
461  DO j=1,i-1
462  work(u11+i,j)=a(cut+i,cut+j)
463  END DO
464  END DO
465 *
466 * invD*L21
467 *
468  i=n-cut-nnb
469  DO WHILE (i .GE. 1)
470  IF (ipiv(cut+nnb+i) > 0) THEN
471  DO j=1,nnb
472  work(i,j)=work(cut+nnb+i,invd)*work(i,j)
473  END DO
474  i=i-1
475  ELSE
476  DO j=1,nnb
477  u01_i_j = work(i,j)
478  u01_ip1_j = work(i-1,j)
479  work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
480  $ work(cut+nnb+i,invd+1)*u01_ip1_j
481  work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
482  $ work(cut+nnb+i-1,invd)*u01_ip1_j
483  END DO
484  i=i-2
485  END IF
486  END DO
487 *
488 * invD1*L11
489 *
490  i=nnb
491  DO WHILE (i .GE. 1)
492  IF (ipiv(cut+i) > 0) THEN
493  DO j=1,nnb
494  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
495  END DO
496  i=i-1
497  ELSE
498  DO j=1,nnb
499  u11_i_j = work(u11+i,j)
500  u11_ip1_j = work(u11+i-1,j)
501  work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
502  $ work(cut+i,invd+1)*u11_ip1_j
503  work(u11+i-1,j)=work(cut+i-1,invd+1)*u11_i_j+
504  $ work(cut+i-1,invd)*u11_ip1_j
505  END DO
506  i=i-2
507  END IF
508  END DO
509 *
510 * L11**T*invD1*L11->L11
511 *
512  CALL strmm('L',uplo,'T','U',nnb, nnb,
513  $ one,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
514 
515 *
516  DO i=1,nnb
517  DO j=1,i
518  a(cut+i,cut+j)=work(u11+i,j)
519  END DO
520  END DO
521 *
522  IF ( (cut+nnb) .LT. n ) THEN
523 *
524 * L21**T*invD2*L21->A(CUT+I,CUT+J)
525 *
526  CALL sgemm('T','N',nnb,nnb,n-nnb-cut,one,a(cut+nnb+1,cut+1)
527  $ ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
528 
529 *
530 * L11 = L11**T*invD1*L11 + U01**T*invD*U01
531 *
532  DO i=1,nnb
533  DO j=1,i
534  a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
535  END DO
536  END DO
537 *
538 * L01 = L22**T*invD2*L21
539 *
540  CALL strmm('L',uplo,'T','U', n-nnb-cut, nnb,
541  $ one,a(cut+nnb+1,cut+nnb+1),lda,work,n+nb+1)
542 *
543 * Update L21
544 *
545  DO i=1,n-cut-nnb
546  DO j=1,nnb
547  a(cut+nnb+i,cut+j)=work(i,j)
548  END DO
549  END DO
550 
551  ELSE
552 *
553 * L11 = L11**T*invD1*L11
554 *
555  DO i=1,nnb
556  DO j=1,i
557  a(cut+i,cut+j)=work(u11+i,j)
558  END DO
559  END DO
560  END IF
561 *
562 * Next Block
563 *
564  cut=cut+nnb
565  END DO
566 *
567 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
568 *
569  i=n
570  DO WHILE ( i .GE. 1 )
571  IF( ipiv(i) .GT. 0 ) THEN
572  ip=ipiv(i)
573  IF (i .LT. ip) CALL ssyswapr( uplo, n, a, lda, i ,ip )
574  IF (i .GT. ip) CALL ssyswapr( uplo, n, a, lda, ip ,i )
575  ELSE
576  ip=-ipiv(i)
577  IF ( i .LT. ip) CALL ssyswapr( uplo, n, a, lda, i ,ip )
578  IF ( i .GT. ip) CALL ssyswapr( uplo, n, a, lda, ip ,i )
579  i=i-1
580  ENDIF
581  i=i-1
582  END DO
583  END IF
584 *
585  RETURN
586 *
587 * End of SSYTRI2X
588 *
subroutine ssyconv(UPLO, WAY, N, A, LDA, IPIV, E, INFO)
SSYCONV
Definition: ssyconv.f:116
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
Definition: strmm.f:179
subroutine strtri(UPLO, DIAG, N, A, LDA, INFO)
STRTRI
Definition: strtri.f:111
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine ssyswapr(UPLO, N, A, LDA, I1, I2)
SSYSWAPR applies an elementary permutation on the rows and columns of a symmetric matrix...
Definition: ssyswapr.f:104
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
Here is the call graph for this function:
Here is the caller graph for this function: