Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
tree: 2911cd2a6e
Fetching contributors…

Cannot retrieve contributors at this time

79 lines (77 sloc) 2.18 kB
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,*)
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
Jump to Line
Something went wrong with that request. Please try again.