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