| C Modified 2002-05-20 for R to add a tolerance of positive definiteness. |
| C |
| c |
| c dpofa factors a double precision symmetric positive definite |
| c matrix. |
| c |
| c dpofa is usually called by dpoco, but it can be called |
| c directly with a saving in time if rcond is not needed. |
| c (time for dpoco) = (1 + 18/n)*(time for dpofa) . |
| c |
| c on entry |
| c |
| c a double precision(lda, n) |
| c the symmetric matrix to be factored. only the |
| c diagonal and upper triangle are used. |
| c |
| c lda integer |
| c the leading dimension of the array a . |
| c |
| c n integer |
| c the order of the matrix a . |
| c |
| c on return |
| c |
| c a an upper triangular matrix r so that a = trans(r)*r |
| c where trans(r) is the transpose. |
| c the strict lower triangle is unaltered. |
| c if info .ne. 0 , the factorization is not complete. |
| c |
| c info integer |
| c = 0 for normal return. |
| c = k signals an error condition. the leading minor |
| c of order k is not positive definite. |
| 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 ddot |
| c fortran dsqrt |
| c |
| subroutine dpofa(a,lda,n,info) |
| integer lda,n,info |
| double precision a(lda,n) |
| c |
| c internal variables |
| c |
| double precision ddot,t,eps |
| double precision s |
| integer j,jm1,k |
| data eps/1.d-14/ |
| |
| c begin block with ...exits to 40 |
| c |
| c |
| do 30 j = 1, n |
| info = j |
| s = 0.0d0 |
| jm1 = j - 1 |
| if (jm1 .lt. 1) go to 20 |
| do 10 k = 1, jm1 |
| t = a(k,j) - ddot(k-1,a(1,k),1,a(1,j),1) |
| t = t/a(k,k) |
| a(k,j) = t |
| s = s + t*t |
| 10 continue |
| 20 continue |
| s = a(j,j) - s |
| c ......exit |
| c if (s .le. 0.0d0) go to 40 |
| if (s .le. eps * abs(a(j,j))) go to 40 |
| a(j,j) = dsqrt(s) |
| 30 continue |
| info = 0 |
| 40 continue |
| return |
| end |
