| 1 |
subroutine spoco(a,lda,n,rcond,z,info) |
| 2 |
integer lda,n,info |
| 3 |
real a(lda,1),z(1) |
| 4 |
real rcond |
| 5 |
c |
| 6 |
c spoco factors a real symmetric positive definite matrix |
| 7 |
c and estimates the condition of the matrix. |
| 8 |
c |
| 9 |
c if rcond is not needed, spofa is slightly faster. |
| 10 |
c to solve a*x = b , follow spoco by sposl. |
| 11 |
c to compute inverse(a)*c , follow spoco by sposl. |
| 12 |
c to compute determinant(a) , follow spoco by spodi. |
| 13 |
c to compute inverse(a) , follow spoco by spodi. |
| 14 |
c |
| 15 |
c on entry |
| 16 |
c |
| 17 |
c a real(lda, n) |
| 18 |
c the symmetric matrix to be factored. only the |
| 19 |
c diagonal and upper triangle are used. |
| 20 |
c |
| 21 |
c lda integer |
| 22 |
c the leading dimension of the array a . |
| 23 |
c |
| 24 |
c n integer |
| 25 |
c the order of the matrix a . |
| 26 |
c |
| 27 |
c on return |
| 28 |
c |
| 29 |
c a an upper triangular matrix r so that a = trans(r)*r |
| 30 |
c where trans(r) is the transpose. |
| 31 |
c the strict lower triangle is unaltered. |
| 32 |
c if info .ne. 0 , the factorization is not complete. |
| 33 |
c |
| 34 |
c rcond real |
| 35 |
c an estimate of the reciprocal condition of a . |
| 36 |
c for the system a*x = b , relative perturbations |
| 37 |
c in a and b of size epsilon may cause |
| 38 |
c relative perturbations in x of size epsilon/rcond . |
| 39 |
c if rcond is so small that the logical expression |
| 40 |
c 1.0 + rcond .eq. 1.0 |
| 41 |
c is true, then a may be singular to working |
| 42 |
c precision. in particular, rcond is zero if |
| 43 |
c exact singularity is detected or the estimate |
| 44 |
c underflows. if info .ne. 0 , rcond is unchanged. |
| 45 |
c |
| 46 |
c z real(n) |
| 47 |
c a work vector whose contents are usually unimportant. |
| 48 |
c if a is close to a singular matrix, then z is |
| 49 |
c an approximate null vector in the sense that |
| 50 |
c norm(a*z) = rcond*norm(a)*norm(z) . |
| 51 |
c if info .ne. 0 , z is unchanged. |
| 52 |
c |
| 53 |
c info integer |
| 54 |
c = 0 for normal return. |
| 55 |
c = k signals an error condition. the leading minor |
| 56 |
c of order k is not positive definite. |
| 57 |
c |
| 58 |
c linpack. this version dated 08/14/78 . |
| 59 |
c cleve moler, university of new mexico, argonne national lab. |
| 60 |
c |
| 61 |
c subroutines and functions |
| 62 |
c |
| 63 |
c linpack spofa |
| 64 |
c blas saxpy,sdot,sscal,sasum |
| 65 |
c fortran abs,amax1,real,sign |
| 66 |
c |
| 67 |
c internal variables |
| 68 |
c |
| 69 |
real sdot,ek,t,wk,wkm |
| 70 |
real anorm,s,sasum,sm,ynorm |
| 71 |
integer i,j,jm1,k,kb,kp1 |
| 72 |
c |
| 73 |
c |
| 74 |
c find norm of a using only upper half |
| 75 |
c |
| 76 |
do 30 j = 1, n |
| 77 |
z(j) = sasum(j,a(1,j),1) |
| 78 |
jm1 = j - 1 |
| 79 |
if (jm1 .lt. 1) go to 20 |
| 80 |
do 10 i = 1, jm1 |
| 81 |
z(i) = z(i) + abs(a(i,j)) |
| 82 |
10 continue |
| 83 |
20 continue |
| 84 |
30 continue |
| 85 |
anorm = 0.00 |
| 86 |
do 40 j = 1, n |
| 87 |
anorm = amax1(anorm,z(j)) |
| 88 |
40 continue |
| 89 |
c |
| 90 |
c factor |
| 91 |
c |
| 92 |
call spofa(a,lda,n,info) |
| 93 |
if (info .ne. 0) go to 180 |
| 94 |
c |
| 95 |
c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . |
| 96 |
c estimate = norm(z)/norm(y) where a*z = y and a*y = e . |
| 97 |
c the components of e are chosen to cause maximum local |
| 98 |
c growth in the elements of w where trans(r)*w = e . |
| 99 |
c the vectors are frequently rescaled to avoid overflow. |
| 100 |
c |
| 101 |
c solve trans(r)*w = e |
| 102 |
c |
| 103 |
ek = 1.00 |
| 104 |
do 50 j = 1, n |
| 105 |
z(j) = 0.00 |
| 106 |
50 continue |
| 107 |
do 110 k = 1, n |
| 108 |
if (z(k) .ne. 0.00) ek = sign(ek,-z(k)) |
| 109 |
if (abs(ek-z(k)) .le. a(k,k)) go to 60 |
| 110 |
s = a(k,k)/abs(ek-z(k)) |
| 111 |
call sscal(n,s,z,1) |
| 112 |
ek = s*ek |
| 113 |
60 continue |
| 114 |
wk = ek - z(k) |
| 115 |
wkm = -ek - z(k) |
| 116 |
s = abs(wk) |
| 117 |
sm = abs(wkm) |
| 118 |
wk = wk/a(k,k) |
| 119 |
wkm = wkm/a(k,k) |
| 120 |
kp1 = k + 1 |
| 121 |
if (kp1 .gt. n) go to 100 |
| 122 |
do 70 j = kp1, n |
| 123 |
sm = sm + abs(z(j)+wkm*a(k,j)) |
| 124 |
z(j) = z(j) + wk*a(k,j) |
| 125 |
s = s + abs(z(j)) |
| 126 |
70 continue |
| 127 |
if (s .ge. sm) go to 90 |
| 128 |
t = wkm - wk |
| 129 |
wk = wkm |
| 130 |
do 80 j = kp1, n |
| 131 |
z(j) = z(j) + t*a(k,j) |
| 132 |
80 continue |
| 133 |
90 continue |
| 134 |
100 continue |
| 135 |
z(k) = wk |
| 136 |
110 continue |
| 137 |
s = 1.00/sasum(n,z,1) |
| 138 |
call sscal(n,s,z,1) |
| 139 |
c |
| 140 |
c solve r*y = w |
| 141 |
c |
| 142 |
do 130 kb = 1, n |
| 143 |
k = n + 1 - kb |
| 144 |
if (abs(z(k)) .le. a(k,k)) go to 120 |
| 145 |
s = a(k,k)/abs(z(k)) |
| 146 |
call sscal(n,s,z,1) |
| 147 |
120 continue |
| 148 |
z(k) = z(k)/a(k,k) |
| 149 |
t = -z(k) |
| 150 |
call saxpy(k-1,t,a(1,k),1,z(1),1) |
| 151 |
130 continue |
| 152 |
s = 1.00/sasum(n,z,1) |
| 153 |
call sscal(n,s,z,1) |
| 154 |
c |
| 155 |
ynorm = 1.00 |
| 156 |
c |
| 157 |
c solve trans(r)*v = y |
| 158 |
c |
| 159 |
do 150 k = 1, n |
| 160 |
z(k) = z(k) - sdot(k-1,a(1,k),1,z(1),1) |
| 161 |
if (abs(z(k)) .le. a(k,k)) go to 140 |
| 162 |
s = a(k,k)/abs(z(k)) |
| 163 |
call sscal(n,s,z,1) |
| 164 |
ynorm = s*ynorm |
| 165 |
140 continue |
| 166 |
z(k) = z(k)/a(k,k) |
| 167 |
150 continue |
| 168 |
s = 1.00/sasum(n,z,1) |
| 169 |
call sscal(n,s,z,1) |
| 170 |
ynorm = s*ynorm |
| 171 |
c |
| 172 |
c solve r*z = v |
| 173 |
c |
| 174 |
do 170 kb = 1, n |
| 175 |
k = n + 1 - kb |
| 176 |
if (abs(z(k)) .le. a(k,k)) go to 160 |
| 177 |
s = a(k,k)/abs(z(k)) |
| 178 |
call sscal(n,s,z,1) |
| 179 |
ynorm = s*ynorm |
| 180 |
160 continue |
| 181 |
z(k) = z(k)/a(k,k) |
| 182 |
t = -z(k) |
| 183 |
call saxpy(k-1,t,a(1,k),1,z(1),1) |
| 184 |
170 continue |
| 185 |
c make znorm = 1.0 |
| 186 |
s = 1.00/sasum(n,z,1) |
| 187 |
call sscal(n,s,z,1) |
| 188 |
ynorm = s*ynorm |
| 189 |
c |
| 190 |
if (anorm .ne. 0.00) rcond = ynorm/anorm |
| 191 |
if (anorm .eq. 0.00) rcond = 0.00 |
| 192 |
180 continue |
| 193 |
return |
| 194 |
end |
Parent Directory
| 
Revision Log
| 
Revision Graph