sinerp.f

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