LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ 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

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.```

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:100
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: