1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
|
subroutine dpofa(a,lda,n,info)
integer lda,n,info
double precision a(lda,*)
c!purpose
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!calling sequence
c
c subroutine dpofa(a,lda,n,info)
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!originator
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c!auxiliary routines
c
c blas ddot
c fortran sqrt
c
c!
c internal variables
c
double precision ddot,t
double precision s
integer j,jm1,k
c begin block with ...exits to 40
c
c
do 30 j = 1, n
info = j
s = 0.0d+0
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
if (s .le. 0.0d+0) go to 40
a(j,j) = sqrt(s)
30 continue
info = 0
40 continue
return
end
|
