forked from mothur/mothur
-
Notifications
You must be signed in to change notification settings - Fork 1
/
sgram.f
142 lines (129 loc) · 4.91 KB
/
sgram.f
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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
C Output from Public domain Ratfor, version 1.0
C PURPOSE
C Calculation of the cubic B-spline smoothness prior
C for "usual" interior knot setup.
C Uses BSPVD and INTRV in the CMLIB
C sgm[0-3](nb) Symmetric matrix
C whose (i,j)'th element contains the integral of
C B''(i,.) B''(j,.) , i=1,2 ... nb and j=i,...nb.
C Only the upper four diagonals are computed.
subroutine sgram(sg0,sg1,sg2,sg3,tb,nb)
c implicit none
C indices
integer nb
DOUBLE precision sg0(nb),sg1(nb),sg2(nb),sg3(nb), tb(nb+4)
c -------------
integer ileft,mflag, i,ii,jj, lentb
DOUBLE precision vnikx(4,3),work(16),yw1(4),yw2(4), wpt
c
c integer interv
c external interv
lentb=nb+4
C Initialise the sigma vectors
do 1 i=1,nb
sg0(i)=0.d0
sg1(i)=0.d0
sg2(i)=0.d0
sg3(i)=0.d0
1 continue
ileft = 1
do 2 i=1,nb
C Calculate a linear approximation to the
C second derivative of the non-zero B-splines
C over the interval [tb(i),tb(i+1)].
C call intrv(tb(1),(nb+1),tb(i),ilo,ileft,mflag)
call interv(tb(1), nb+1,tb(i), ileft, mflag)
C Left end second derivatives
C call bspvd (tb,4,3,tb(i),ileft,4,vnikx,work)
call bsplvd (tb,lentb,4,tb(i),ileft,work,vnikx,3)
C Put values into yw1
do 4 ii=1,4
yw1(ii) = vnikx(ii,3)
4 continue
C Right end second derivatives
C call bspvd (tb,4,3,tb(i+1),ileft,4,vnikx,work)
call bsplvd (tb,lentb,4,tb(i+1),ileft,work,vnikx,3)
C Slope*(length of interval) in Linear Approximation to B''
do 6 ii=1,4
yw2(ii) = vnikx(ii,3) - yw1(ii)
6 continue
wpt = tb(i+1) - tb(i)
C Calculate Contributions to the sigma vectors
if(ileft.ge.4) then
do 10 ii=1,4
jj=ii
sg0(ileft-4+ii) = sg0(ileft-4+ii) +
& wpt*(yw1(ii)*yw1(jj)+
& (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
& + yw2(ii)*yw2(jj)*0.3330d0)
jj=ii+1
if(jj.le.4)then
sg1(ileft+ii-4) = sg1(ileft+ii-4) +
& wpt* (yw1(ii)*yw1(jj) +
* (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
& +yw2(ii)*yw2(jj)*0.3330d0 )
endif
jj=ii+2
if(jj.le.4)then
sg2(ileft+ii-4) = sg2(ileft+ii-4) +
& wpt* (yw1(ii)*yw1(jj) +
* (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
& +yw2(ii)*yw2(jj)*0.3330d0 )
endif
jj=ii+3
if(jj.le.4)then
sg3(ileft+ii-4) = sg3(ileft+ii-4) +
& wpt* (yw1(ii)*yw1(jj) +
* (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
& +yw2(ii)*yw2(jj)*0.3330d0 )
endif
10 continue
else if(ileft.eq.3)then
do 20 ii=1,3
jj=ii
sg0(ileft-3+ii) = sg0(ileft-3+ii) +
& wpt* (yw1(ii)*yw1(jj) +
* (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
& +yw2(ii)*yw2(jj)*0.3330d0 )
jj=ii+1
if(jj.le.3)then
sg1(ileft+ii-3) = sg1(ileft+ii-3) +
& wpt* (yw1(ii)*yw1(jj) +
* (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
& +yw2(ii)*yw2(jj)*0.3330d0 )
endif
jj=ii+2
if(jj.le.3)then
sg2(ileft+ii-3) = sg2(ileft+ii-3) +
& wpt* (yw1(ii)*yw1(jj) +
* (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
& +yw2(ii)*yw2(jj)*0.3330d0 )
endif
20 continue
else if(ileft.eq.2)then
do 28 ii=1,2
jj=ii
sg0(ileft-2+ii) = sg0(ileft-2+ii) +
& wpt* (yw1(ii)*yw1(jj) +
* (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
& +yw2(ii)*yw2(jj)*0.3330d0 )
jj=ii+1
if(jj.le.2)then
sg1(ileft+ii-2) = sg1(ileft+ii-2) +
& wpt* (yw1(ii)*yw1(jj) +
* (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
& +yw2(ii)*yw2(jj)*0.3330d0 )
endif
28 continue
else if(ileft.eq.1)then
do 34 ii=1,1
jj=ii
sg0(ileft-1+ii) = sg0(ileft-1+ii) +
& wpt* (yw1(ii)*yw1(jj) +
* (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*0.5d0
& +yw2(ii)*yw2(jj)*0.3330d0 )
34 continue
endif
2 continue
return
end