sppdi.f
128 lines
| 3.3 KiB
| text/x-fortran
|
FortranFixedLexer
| r1601 | subroutine sppdi(ap,n,det,job) | |||
| integer n,job | ||||
| real ap(1) | ||||
| real det(2) | ||||
| c | ||||
| c sppdi computes the determinant and inverse | ||||
| c of a real symmetric positive definite matrix | ||||
| c using the factors computed by sppco or sppfa . | ||||
| c | ||||
| c on entry | ||||
| c | ||||
| c ap real (n*(n+1)/2) | ||||
| c the output from sppco or sppfa. | ||||
| c | ||||
| c n integer | ||||
| c the order of the matrix a . | ||||
| c | ||||
| c job integer | ||||
| c = 11 both determinant and inverse. | ||||
| c = 01 inverse only. | ||||
| c = 10 determinant only. | ||||
| c | ||||
| c on return | ||||
| c | ||||
| c ap the upper triangular half of the inverse . | ||||
| c the strict lower triangle is unaltered. | ||||
| c | ||||
| c det real(2) | ||||
| c determinant of original matrix if requested. | ||||
| c otherwise not referenced. | ||||
| c determinant = det(1) * 10.0**det(2) | ||||
| c with 1.0 .le. det(1) .lt. 10.0 | ||||
| c or det(1) .eq. 0.0 . | ||||
| c | ||||
| c error condition | ||||
| c | ||||
| c a division by zero will occur if the input factor contains | ||||
| c a zero on the diagonal and the inverse is requested. | ||||
| c it will not occur if the subroutines are called correctly | ||||
| c and if spoco or spofa has set info .eq. 0 . | ||||
| c | ||||
| c linpack. this version dated 08/14/78 . | ||||
| c cleve moler, university of new mexico, argonne national lab. | ||||
| c | ||||
| c subroutines and functions | ||||
| c | ||||
| c blas saxpy,sscal | ||||
| c fortran mod | ||||
| c | ||||
| c internal variables | ||||
| c | ||||
| real t | ||||
| real s | ||||
| integer i,ii,j,jj,jm1,j1,k,kj,kk,kp1,k1 | ||||
| c | ||||
| c compute determinant | ||||
| c | ||||
| if (job/10 .eq. 0) go to 70 | ||||
| det(1) = 1.0e0 | ||||
| det(2) = 0.0e0 | ||||
| s = 10.0e0 | ||||
| ii = 0 | ||||
| do 50 i = 1, n | ||||
| ii = ii + i | ||||
| det(1) = ap(ii)**2*det(1) | ||||
| c ...exit | ||||
| if (det(1) .eq. 0.0e0) go to 60 | ||||
| 10 if (det(1) .ge. 1.0e0) go to 20 | ||||
| det(1) = s*det(1) | ||||
| det(2) = det(2) - 1.0e0 | ||||
| go to 10 | ||||
| 20 continue | ||||
| 30 if (det(1) .lt. s) go to 40 | ||||
| det(1) = det(1)/s | ||||
| det(2) = det(2) + 1.0e0 | ||||
| go to 30 | ||||
| 40 continue | ||||
| 50 continue | ||||
| 60 continue | ||||
| 70 continue | ||||
| c | ||||
| c compute inverse(r) | ||||
| c | ||||
| if (mod(job,10) .eq. 0) go to 140 | ||||
| kk = 0 | ||||
| do 100 k = 1, n | ||||
| k1 = kk + 1 | ||||
| kk = kk + k | ||||
| ap(kk) = 1.0e0/ap(kk) | ||||
| t = -ap(kk) | ||||
| call sscal(k-1,t,ap(k1),1) | ||||
| kp1 = k + 1 | ||||
| j1 = kk + 1 | ||||
| kj = kk + k | ||||
| if (n .lt. kp1) go to 90 | ||||
| do 80 j = kp1, n | ||||
| t = ap(kj) | ||||
| ap(kj) = 0.0e0 | ||||
| call saxpy(k,t,ap(k1),1,ap(j1),1) | ||||
| j1 = j1 + j | ||||
| kj = kj + j | ||||
| 80 continue | ||||
| 90 continue | ||||
| 100 continue | ||||
| c | ||||
| c form inverse(r) * trans(inverse(r)) | ||||
| c | ||||
| jj = 0 | ||||
| do 130 j = 1, n | ||||
| j1 = jj + 1 | ||||
| jj = jj + j | ||||
| jm1 = j - 1 | ||||
| k1 = 1 | ||||
| kj = j1 | ||||
| if (jm1 .lt. 1) go to 120 | ||||
| do 110 k = 1, jm1 | ||||
| t = ap(kj) | ||||
| call saxpy(k,t,ap(j1),1,ap(k1),1) | ||||
| k1 = k1 + k | ||||
| kj = kj + 1 | ||||
| 110 continue | ||||
| 120 continue | ||||
| t = ap(jj) | ||||
| call sscal(j,t,ap(j1),1) | ||||
| 130 continue | ||||
| 140 continue | ||||
| return | ||||
| end | ||||
