lubksb Subroutine

public subroutine lubksb(a, indx, b)

Does LU back substitution

Arguments

TypeIntentOptionalAttributesName
real(kind=double), intent(in), dimension(:,:):: a
integer, intent(in), dimension(:):: indx
real(kind=double), intent(inout), dimension(:):: b

Called by

proc~~lubksb~~CalledByGraph proc~lubksb lubksb proc~matinv matinv proc~matinv->proc~lubksb proc~linsolve linsolve proc~linsolve->proc~matinv proc~colrad colrad proc~colrad->proc~linsolve proc~ndmc ndmc proc~ndmc->proc~colrad proc~fida_weights_los fida_weights_los proc~fida_weights_los->proc~colrad proc~dcx dcx proc~dcx->proc~colrad proc~fida_f fida_f proc~fida_f->proc~colrad proc~pfida_f pfida_f proc~pfida_f->proc~colrad proc~read_equilibrium read_equilibrium proc~read_equilibrium->proc~colrad proc~attenuate attenuate proc~attenuate->proc~colrad proc~halo halo proc~halo->proc~colrad proc~pfida_mc pfida_mc proc~pfida_mc->proc~colrad proc~fida_weights_mc fida_weights_mc proc~fida_weights_mc->proc~colrad proc~fida_mc fida_mc proc~fida_mc->proc~colrad program~fidasim fidasim program~fidasim->proc~ndmc program~fidasim->proc~fida_weights_los program~fidasim->proc~dcx program~fidasim->proc~fida_f program~fidasim->proc~pfida_f program~fidasim->proc~read_equilibrium program~fidasim->proc~halo program~fidasim->proc~pfida_mc program~fidasim->proc~fida_weights_mc program~fidasim->proc~fida_mc proc~pnpa_mc pnpa_mc program~fidasim->proc~pnpa_mc proc~pnpa_f pnpa_f program~fidasim->proc~pnpa_f proc~npa_mc npa_mc program~fidasim->proc~npa_mc proc~npa_f npa_f program~fidasim->proc~npa_f proc~npa_weights npa_weights program~fidasim->proc~npa_weights proc~pnpa_mc->proc~attenuate proc~pnpa_f->proc~attenuate proc~npa_mc->proc~attenuate proc~npa_f->proc~attenuate proc~npa_weights->proc~attenuate

Contents

Source Code


Source Code

  subroutine lubksb(a,indx,b)
    !+ Does LU back substitution
    real(double), dimension(:,:),intent(IN)   :: a
    integer,dimension(:),  intent(IN)         :: indx
    real(double), dimension(:),  intent(INOUT):: b
    integer       :: i,n,ii,ll
    real(double)  :: summ
    n=size(indx)
    ii=0
    do i=1,n
       ll=indx(i)
       summ=b(ll)
       b(ll)=b(i)
       if (ii /= 0) then
          summ=summ-dot_product(a(i,ii:i-1),b(ii:i-1))
       else if (summ /= 0.0) then
          ii=i
       endif
       b(i)=summ
    enddo
    do i=n,1,-1
       b(i) = (b(i)-dot_product(a(i,i+1:n),b(i+1:n)))/a(i,i)
    enddo
  end subroutine lubksb