1C Output from Public domain Ratfor, version 1.0
2      subroutine sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,flag)
3c
4C Purpose :  Computes Inner Products between columns of L^{-1}
5C            where L = abd is a Banded Matrix with 3 subdiagonals
6
7C In R, only called in one place, from sslrvg() [ ./sslvrg.f ] with flag = 0
8
9C The algorithm works in two passes:
10C
11C               Pass 1 computes (cj,ck) k=j,j-1,j-2,j-3 ;  j=nk, .. 1
12C               Pass 2 computes (cj,ck) k <= j-4  (If flag == 1 ).
13C
14C               A refinement of Elden's trick is used.
15c Args
16      integer ld4,nk,ldnk,flag
17      DOUBLE precision abd(ld4,nk), p1ip(ld4,nk), p2ip(ldnk,nk)
18c Locals
19      integer i,j,k
20      DOUBLE precision  wjm3(3),wjm2(2),wjm1(1),c0,c1,c2,c3
21c
22c     unnecessary initialization of c1 c2 c3 to keep g77 -Wall happy
23c
24      c1 = 0.0d0
25      c2 = 0.0d0
26      c3 = 0.0d0
27C
28C Pass 1
29      wjm3(1)=0d0
30      wjm3(2)=0d0
31      wjm3(3)=0d0
32      wjm2(1)=0d0
33      wjm2(2)=0d0
34      wjm1(1)=0d0
35      do 100 i=1,nk
36         j=nk-i+1
37         c0 = 1d0/abd(4,j)
38         if(j.le.nk-3)then
39            c1 = abd(1,j+3)*c0
40            c2 = abd(2,j+2)*c0
41            c3 = abd(3,j+1)*c0
42         else if(j.eq.nk-2)then
43            c1 = 0d0
44            c2 = abd(2,j+2)*c0
45            c3 = abd(3,j+1)*c0
46         else if(j.eq.nk-1)then
47            c1 = 0d0
48            c2 = 0d0
49            c3 = abd(3,j+1)*c0
50         else if(j.eq.nk)then
51            c1 = 0d0
52            c2 = 0d0
53            c3 = 0d0
54         endif
55         p1ip(1,j) = 0d0- (c1*wjm3(1)+c2*wjm3(2)+c3*wjm3(3))
56         p1ip(2,j) = 0d0- (c1*wjm3(2)+c2*wjm2(1)+c3*wjm2(2))
57         p1ip(3,j) = 0d0- (c1*wjm3(3)+c2*wjm2(2)+c3*wjm1(1))
58         p1ip(4,j) = c0**2 + c1**2*wjm3(1) + 2d0*c1*c2*wjm3(2)+
59     &        2d0*c1*c3*wjm3(3) + c2**2*wjm2(1) + 2d0*c2*c3*wjm2(2) +
60     &        c3**2*wjm1(1)
61         wjm3(1)=wjm2(1)
62         wjm3(2)=wjm2(2)
63         wjm3(3)=p1ip(2,j)
64         wjm2(1)=wjm1(1)
65         wjm2(2)=p1ip(3,j)
66         wjm1(1)=p1ip(4,j)
67 100  continue
68
69      if(flag.ne.0) then
70
71C     ____ Pass 2 _____  Compute p2ip  [never from R's code!]
72         do 120 i=1,nk
73            j=nk-i+1
74C           for(k=1;k<=4 & j+k-1<=nk;k=k+1) { p2ip(.) = .. }:
75            do k=1,4
76               if(j+k-1 .gt. nk) goto 120
77               p2ip(j,j+k-1) = p1ip(5-k,j)
78            end do
79 120     continue
80
81         do 170 i=1,nk
82            j=nk-i+1
83c           for(k=j-4;k>=1;k=k-1){
84            if(j-4 .ge. 1) then
85               do k= j-4,1, -1
86                  c0 = 1d0/abd(4,k)
87                  c1 = abd(1,k+3)*c0
88                  c2 = abd(2,k+2)*c0
89                  c3 = abd(3,k+1)*c0
90                  p2ip(k,j)= 0d0 - ( c1*p2ip(k+3,j) + c2*p2ip(k+2,j) +
91     &                 c3*p2ip(k+1,j) )
92               end do
93            endif
94 170     continue
95      endif
96      return
97      end
98
99