LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
ssytri2x.f
Go to the documentation of this file.
1 *> \brief \b SSYTRI2X
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SSYTRI2X + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytri2x.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytri2x.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytri2x.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, LDA, N, NB
26 * ..
27 * .. Array Arguments ..
28 * INTEGER IPIV( * )
29 * REAL A( LDA, * ), WORK( N+NB+1,* )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> SSYTRI2X computes the inverse of a real symmetric indefinite matrix
39 *> A using the factorization A = U*D*U**T or A = L*D*L**T computed by
40 *> SSYTRF.
41 *> \endverbatim
42 *
43 * Arguments:
44 * ==========
45 *
46 *> \param[in] UPLO
47 *> \verbatim
48 *> UPLO is CHARACTER*1
49 *> Specifies whether the details of the factorization are stored
50 *> as an upper or lower triangular matrix.
51 *> = 'U': Upper triangular, form is A = U*D*U**T;
52 *> = 'L': Lower triangular, form is A = L*D*L**T.
53 *> \endverbatim
54 *>
55 *> \param[in] N
56 *> \verbatim
57 *> N is INTEGER
58 *> The order of the matrix A. N >= 0.
59 *> \endverbatim
60 *>
61 *> \param[in,out] A
62 *> \verbatim
63 *> A is REAL array, dimension (LDA,N)
64 *> On entry, the NNB diagonal matrix D and the multipliers
65 *> used to obtain the factor U or L as computed by SSYTRF.
66 *>
67 *> On exit, if INFO = 0, the (symmetric) inverse of the original
68 *> matrix. If UPLO = 'U', the upper triangular part of the
69 *> inverse is formed and the part of A below the diagonal is not
70 *> referenced; if UPLO = 'L' the lower triangular part of the
71 *> inverse is formed and the part of A above the diagonal is
72 *> not referenced.
73 *> \endverbatim
74 *>
75 *> \param[in] LDA
76 *> \verbatim
77 *> LDA is INTEGER
78 *> The leading dimension of the array A. LDA >= max(1,N).
79 *> \endverbatim
80 *>
81 *> \param[in] IPIV
82 *> \verbatim
83 *> IPIV is INTEGER array, dimension (N)
84 *> Details of the interchanges and the NNB structure of D
85 *> as determined by SSYTRF.
86 *> \endverbatim
87 *>
88 *> \param[out] WORK
89 *> \verbatim
90 *> WORK is REAL array, dimension (N+NB+1,NB+3)
91 *> \endverbatim
92 *>
93 *> \param[in] NB
94 *> \verbatim
95 *> NB is INTEGER
96 *> Block size
97 *> \endverbatim
98 *>
99 *> \param[out] INFO
100 *> \verbatim
101 *> INFO is INTEGER
102 *> = 0: successful exit
103 *> < 0: if INFO = -i, the i-th argument had an illegal value
104 *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
105 *> inverse could not be computed.
106 *> \endverbatim
107 *
108 * Authors:
109 * ========
110 *
111 *> \author Univ. of Tennessee
112 *> \author Univ. of California Berkeley
113 *> \author Univ. of Colorado Denver
114 *> \author NAG Ltd.
115 *
116 *> \ingroup realSYcomputational
117 *
118 * =====================================================================
119  SUBROUTINE ssytri2x( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
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 *
586  END
587 
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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 ssytri2x(UPLO, N, A, LDA, IPIV, WORK, NB, INFO)
SSYTRI2X
Definition: ssytri2x.f:120
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