| 1 |
afe |
1.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