LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
dsytri2x.f
Go to the documentation of this file.
1*> \brief \b DSYTRI2X
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytri2x.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytri2x.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytri2x.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSYTRI2X( 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* DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DSYTRI2X 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*> DSYTRF.
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 DOUBLE PRECISION 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 DSYTRF.
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 DSYTRF.
86*> \endverbatim
87*>
88*> \param[out] WORK
89*> \verbatim
90*> WORK is DOUBLE PRECISION 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 hetri2x
117*
118* =====================================================================
119 SUBROUTINE dsytri2x( 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 DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* )
132* ..
133*
134* =====================================================================
135*
136* .. Parameters ..
137 DOUBLE PRECISION ONE, ZERO
138 parameter( one = 1.0d+0, zero = 0.0d+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 DOUBLE PRECISION AK, AKKP1, AKP1, D, T
147 DOUBLE PRECISION U01_I_J, U01_IP1_J
148 DOUBLE PRECISION U11_I_J, U11_IP1_J
149* ..
150* .. External Functions ..
151 LOGICAL LSAME
152 EXTERNAL lsame
153* ..
154* .. External Subroutines ..
155 EXTERNAL dsyconv, xerbla, dtrtri
156 EXTERNAL dgemm, dtrmm, dsyswapr
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( 'DSYTRI2X', -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 dsyconv( 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 dtrtri( 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 dtrmm('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 dgemm('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*
353* U11 = U11**T*invD1*U11 + U01**T*invD*U01
354*
355 DO i=1,nnb
356 DO j=i,nnb
357 a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
358 END DO
359 END DO
360*
361* U01 = U00**T*invD0*U01
362*
363 CALL dtrmm('L',uplo,'T','U',cut, nnb,
364 \$ one,a,lda,work,n+nb+1)
365
366*
367* Update U01
368*
369 DO i=1,cut
370 DO j=1,nnb
371 a(i,cut+j)=work(i,j)
372 END DO
373 END DO
374*
375* Next Block
376*
377 END DO
378*
379* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
380*
381 i=1
382 DO WHILE ( i .LE. n )
383 IF( ipiv(i) .GT. 0 ) THEN
384 ip=ipiv(i)
385 IF (i .LT. ip) CALL dsyswapr( uplo, n, a, lda, i ,ip )
386 IF (i .GT. ip) CALL dsyswapr( uplo, n, a, lda, ip ,i )
387 ELSE
388 ip=-ipiv(i)
389 i=i+1
390 IF ( (i-1) .LT. ip)
391 \$ CALL dsyswapr( uplo, n, a, lda, i-1 ,ip )
392 IF ( (i-1) .GT. ip)
393 \$ CALL dsyswapr( uplo, n, a, lda, ip ,i-1 )
394 ENDIF
395 i=i+1
396 END DO
397 ELSE
398*
399* LOWER...
400*
401* invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
402*
403 CALL dtrtri( uplo, 'U', n, a, lda, info )
404*
405* inv(D) and inv(D)*inv(U)
406*
407 k=n
408 DO WHILE ( k .GE. 1 )
409 IF( ipiv( k ).GT.0 ) THEN
410* 1 x 1 diagonal NNB
411 work(k,invd) = one / a( k, k )
412 work(k,invd+1) = 0
413 k=k-1
414 ELSE
415* 2 x 2 diagonal NNB
416 t = work(k-1,1)
417 ak = a( k-1, k-1 ) / t
418 akp1 = a( k, k ) / t
419 akkp1 = work(k-1,1) / t
420 d = t*( ak*akp1-one )
421 work(k-1,invd) = akp1 / d
422 work(k,invd) = ak / d
423 work(k,invd+1) = -akkp1 / d
424 work(k-1,invd+1) = -akkp1 / d
425 k=k-2
426 END IF
427 END DO
428*
429* inv(U**T) = (inv(U))**T
430*
431* inv(U**T)*inv(D)*inv(U)
432*
433 cut=0
434 DO WHILE (cut .LT. n)
435 nnb=nb
436 IF (cut + nnb .GT. n) THEN
437 nnb=n-cut
438 ELSE
439 count = 0
440* count negative elements,
441 DO i=cut+1,cut+nnb
442 IF (ipiv(i) .LT. 0) count=count+1
443 END DO
444* need a even number for a clear cut
445 IF (mod(count,2) .EQ. 1) nnb=nnb+1
446 END IF
447* L21 Block
448 DO i=1,n-cut-nnb
449 DO j=1,nnb
450 work(i,j)=a(cut+nnb+i,cut+j)
451 END DO
452 END DO
453* L11 Block
454 DO i=1,nnb
455 work(u11+i,i)=one
456 DO j=i+1,nnb
457 work(u11+i,j)=zero
458 END DO
459 DO j=1,i-1
460 work(u11+i,j)=a(cut+i,cut+j)
461 END DO
462 END DO
463*
464* invD*L21
465*
466 i=n-cut-nnb
467 DO WHILE (i .GE. 1)
468 IF (ipiv(cut+nnb+i) > 0) THEN
469 DO j=1,nnb
470 work(i,j)=work(cut+nnb+i,invd)*work(i,j)
471 END DO
472 i=i-1
473 ELSE
474 DO j=1,nnb
475 u01_i_j = work(i,j)
476 u01_ip1_j = work(i-1,j)
477 work(i,j)=work(cut+nnb+i,invd)*u01_i_j+
478 \$ work(cut+nnb+i,invd+1)*u01_ip1_j
479 work(i-1,j)=work(cut+nnb+i-1,invd+1)*u01_i_j+
480 \$ work(cut+nnb+i-1,invd)*u01_ip1_j
481 END DO
482 i=i-2
483 END IF
484 END DO
485*
486* invD1*L11
487*
488 i=nnb
489 DO WHILE (i .GE. 1)
490 IF (ipiv(cut+i) > 0) THEN
491 DO j=1,nnb
492 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j)
493 END DO
494 i=i-1
495 ELSE
496 DO j=1,nnb
497 u11_i_j = work(u11+i,j)
498 u11_ip1_j = work(u11+i-1,j)
499 work(u11+i,j)=work(cut+i,invd)*work(u11+i,j) +
500 \$ work(cut+i,invd+1)*u11_ip1_j
501 work(u11+i-1,j)=work(cut+i-1,invd+1)*u11_i_j+
502 \$ work(cut+i-1,invd)*u11_ip1_j
503 END DO
504 i=i-2
505 END IF
506 END DO
507*
508* L11**T*invD1*L11->L11
509*
510 CALL dtrmm('L',uplo,'T','U',nnb, nnb,
511 \$ one,a(cut+1,cut+1),lda,work(u11+1,1),n+nb+1)
512
513*
514 DO i=1,nnb
515 DO j=1,i
516 a(cut+i,cut+j)=work(u11+i,j)
517 END DO
518 END DO
519*
520 IF ( (cut+nnb) .LT. n ) THEN
521*
522* L21**T*invD2*L21->A(CUT+I,CUT+J)
523*
524 CALL dgemm('T','N',nnb,nnb,n-nnb-cut,one,a(cut+nnb+1,cut+1)
525 \$ ,lda,work,n+nb+1, zero, work(u11+1,1), n+nb+1)
526
527*
528* L11 = L11**T*invD1*L11 + U01**T*invD*U01
529*
530 DO i=1,nnb
531 DO j=1,i
532 a(cut+i,cut+j)=a(cut+i,cut+j)+work(u11+i,j)
533 END DO
534 END DO
535*
536* L01 = L22**T*invD2*L21
537*
538 CALL dtrmm('L',uplo,'T','U', n-nnb-cut, nnb,
539 \$ one,a(cut+nnb+1,cut+nnb+1),lda,work,n+nb+1)
540*
541* Update L21
542*
543 DO i=1,n-cut-nnb
544 DO j=1,nnb
545 a(cut+nnb+i,cut+j)=work(i,j)
546 END DO
547 END DO
548
549 ELSE
550*
551* L11 = L11**T*invD1*L11
552*
553 DO i=1,nnb
554 DO j=1,i
555 a(cut+i,cut+j)=work(u11+i,j)
556 END DO
557 END DO
558 END IF
559*
560* Next Block
561*
562 cut=cut+nnb
563 END DO
564*
565* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
566*
567 i=n
568 DO WHILE ( i .GE. 1 )
569 IF( ipiv(i) .GT. 0 ) THEN
570 ip=ipiv(i)
571 IF (i .LT. ip) CALL dsyswapr( uplo, n, a, lda, i ,ip )
572 IF (i .GT. ip) CALL dsyswapr( uplo, n, a, lda, ip ,i )
573 ELSE
574 ip=-ipiv(i)
575 IF ( i .LT. ip) CALL dsyswapr( uplo, n, a, lda, i ,ip )
576 IF ( i .GT. ip) CALL dsyswapr( uplo, n, a, lda, ip, i )
577 i=i-1
578 ENDIF
579 i=i-1
580 END DO
581 END IF
582*
583 RETURN
584*
585* End of DSYTRI2X
586*
587 END
588
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
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:100
subroutine dsytri2x(uplo, n, a, lda, ipiv, work, nb, info)
DSYTRI2X
Definition dsytri2x.f:120
subroutine dsyconv(uplo, way, n, a, lda, ipiv, e, info)
DSYCONV
Definition dsyconv.f:114
subroutine dtrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRMM
Definition dtrmm.f:177
subroutine dtrtri(uplo, diag, n, a, lda, info)
DTRTRI
Definition dtrtri.f:109