1 C Output from Public domain Ratfor, version 1.0
2 subroutine sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,flag)
4 C Purpose : Computes Inner Products between columns of L^{-1}
5 C where L = abd is a Banded Matrix with 3 subdiagonals
7 C The algorithm works in two passes:
9 C Pass 1 computes (cj,ck) k=j,j-1,j-2,j-3 ; j=nk, .. 1
10 C Pass 2 computes (cj,ck) k <= j-4 (If flag == 1 ).
12 C A refinement of Elden's trick is used.
14 integer ld4,nk,ldnk,flag
15 DOUBLE precision abd(ld4,nk),p1ip(ld4,nk), p2ip(ldnk,nk)
18 DOUBLE precision wjm3(3),wjm2(2),wjm1(1),c0,c1,c2,c3
20 c unnecessary initialization of c1 c2 c3 to keep g77 -Wall happy
40 else if(j.eq.nk-2)then
44 else if(j.eq.nk-1)then
53 p1ip(1,j) = 0d0- (c1*wjm3(1)+c2*wjm3(2)+c3*wjm3(3))
54 p1ip(2,j) = 0d0- (c1*wjm3(2)+c2*wjm2(1)+c3*wjm2(2))
55 p1ip(3,j) = 0d0- (c1*wjm3(3)+c2*wjm2(2)+c3*wjm1(1))
56 p1ip(4,j) = c0**2 + c1**2*wjm3(1) + 2d0*c1*c2*wjm3(2)+
57 & 2d0*c1*c3*wjm3(3) + c2**2*wjm2(1) + 2d0*c2*c3*wjm2(2) +
74 C for(k=1;k<=4 & j+k-1<=nk;k=k+1) { p2ip(.) = .. }:
76 if(j+k-1 .gt. nk)goto 120
77 p2ip(j,j+k-1) = p1ip(5-k,j)
83 c for(k=j-4;k>=1;k=k-1){
90 p2ip(k,j)= 0d0 - ( c1*p2ip(k+3,j) + c2*p2ip(k+2,j) +