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
79
80
81
82
83
84
85
86
87
|
subroutine bchfac ( w, nbands, nrow, diag )
c from * a practical guide to splines * by c. de boor
constructs cholesky factorization
c c = l * d * l-transpose
c with l unit lower triangular and d diagonal, for given matrix c of
c order n r o w , in case c is (symmetric) positive semidefinite
c and b a n d e d , having n b a n d s diagonals at and below the
c main diagonal.
c
c****** i n p u t ******
c nrow.....is the order of the matrix c .
c nbands.....indicates its bandwidth, i.e.,
c c(i,j) = 0 for abs(i-j) .gt. nbands .
c w.....workarray of size (nbands,nrow) containing the nbands diago-
c nals in its rows, with the main diagonal in row 1 . precisely,
c w(i,j) contains c(i+j-1,j), i=1,...,nbands, j=1,...,nrow.
c for example, the interesting entries of a seven diagonal sym-
c metric matrix c of order 9 would be stored in w as
c
c 11 22 33 44 55 66 77 88 99
c 21 32 43 54 65 76 87 98
c 31 42 53 64 75 86 97
c 41 52 63 74 85 96
c
c all other entries of w not identified in this way with an en-
c try of c are never referenced .
c diag.....is a work array of length nrow .
c
c****** o u t p u t ******
c w.....contains the cholesky factorization c = l*d*l-transp, with
c w(1,i) containing 1/d(i,i)
c and w(i,j) containing l(i-1+j,j), i=2,...,nbands.
c
c****** m e t h o d ******
c gauss elimination, adapted to the symmetry and bandedness of c , is
c used .
c near zero pivots are handled in a special way. the diagonal ele-
c ment c(n,n) = w(1,n) is saved initially in diag(n), all n. at the n-
c th elimination step, the current pivot element, viz. w(1,n), is com-
c pared with its original value, diag(n). if, as the result of prior
c elimination steps, this element has been reduced by about a word
c length, (i.e., if w(1,n)+diag(n) .le. diag(n)), then the pivot is de-
c clared to be zero, and the entire n-th row is declared to be linearly
c dependent on the preceding rows. this has the effect of producing
c x(n) = 0 when solving c*x = b for x, regardless of b. justific-
c ation for this is as follows. in contemplated applications of this
c program, the given equations are the normal equations for some least-
c squares approximation problem, diag(n) = c(n,n) gives the norm-square
c of the n-th basis function, and, at this point, w(1,n) contains the
c norm-square of the error in the least-squares approximation to the n-
c th basis function by linear combinations of the first n-1 . having
c w(1,n)+diag(n) .le. diag(n) signifies that the n-th function is lin-
c early dependent to machine accuracy on the first n-1 functions, there
c fore can safely be left out from the basis of approximating functions
c the solution of a linear system
c c*x = b
c is effected by the succession of the following t w o calls:
c call bchfac ( w, nbands, nrow, diag ) , to get factorization
c call bchslv ( w, nbands, nrow, b, x ) , to solve for x.
c
integer nbands,nrow, i,imax,j,jmax,n
real w(nbands,nrow),diag(nrow), ratio
if (nrow .gt. 1) go to 9
if (w(1,1) .gt. 0.) w(1,1) = 1./w(1,1)
return
c store diagonal of c in diag.
9 do 10 n=1,nrow
10 diag(n) = w(1,n)
c factorization .
do 20 n=1,nrow
if (w(1,n)+diag(n) .gt. diag(n)) go to 15
do 14 j=1,nbands
14 w(j,n) = 0.
go to 20
15 w(1,n) = 1./w(1,n)
imax = min0(nbands-1,nrow - n)
if (imax .lt. 1) go to 20
jmax = imax
do 18 i=1,imax
ratio = w(i+1,n)*w(1,n)
do 17 j=1,jmax
17 w(j,n+i) = w(j,n+i) - w(j+i,n)*ratio
jmax = jmax - 1
18 w(i+1,n) = ratio
20 continue
return
end
|
