1*59599516SKenneth E. Jansen subroutine lubksb(aa,n,np,indx,bb) 2*59599516SKenneth E. Jansenc--------------------------------------------------------------------- 3*59599516SKenneth E. Jansenc 4*59599516SKenneth E. Jansenc LU back substitution routine. 5*59599516SKenneth E. Jansenc 6*59599516SKenneth E. Jansenc--------------------------------------------------------------------- 7*59599516SKenneth E. Jansen include "common.h" 8*59599516SKenneth E. Jansen 9*59599516SKenneth E. Jansen dimension indx(n),aa(np,np),bb(n) 10*59599516SKenneth E. Jansen 11*59599516SKenneth E. Jansen ii=0 12*59599516SKenneth E. Jansen do i=1,n 13*59599516SKenneth E. Jansen ll=indx(i) 14*59599516SKenneth E. Jansen sum=bb(ll) 15*59599516SKenneth E. Jansen bb(ll)=bb(i) 16*59599516SKenneth E. Jansen if (ii.ne.0)then 17*59599516SKenneth E. Jansen do j=ii,i-1 18*59599516SKenneth E. Jansen sum=sum-aa(i,j)*bb(j) 19*59599516SKenneth E. Jansen enddo 20*59599516SKenneth E. Jansen else if (sum.ne.0.0) then 21*59599516SKenneth E. Jansen ii=i 22*59599516SKenneth E. Jansen endif 23*59599516SKenneth E. Jansen bb(i)=sum 24*59599516SKenneth E. Jansen enddo 25*59599516SKenneth E. Jansen do i=n,1,-1 26*59599516SKenneth E. Jansen sum=bb(i) 27*59599516SKenneth E. Jansen do j=i+1,n 28*59599516SKenneth E. Jansen sum=sum-aa(i,j)*bb(j) 29*59599516SKenneth E. Jansen enddo 30*59599516SKenneth E. Jansen bb(i)=sum/aa(i,i) 31*59599516SKenneth E. Jansen enddo 32*59599516SKenneth E. Jansen 33*59599516SKenneth E. Jansen return 34*59599516SKenneth E. Jansen end 35*59599516SKenneth E. Jansen 36*59599516SKenneth E. Jansen 37*59599516SKenneth E. Jansen 38