LAPACK  3.10.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.

Definition at line 119 of file ssytri2x.f.

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