ludcmp Subroutine

public subroutine ludcmp(a, indx, d)

Calculates LU decomposition

Arguments

TypeIntentOptionalAttributesName
real(kind=double), intent(inout), dimension(:,:):: a
integer, intent(out), dimension(:):: indx
real(kind=double), intent(out) :: d

Calls

proc~~ludcmp~~CallsGraph proc~ludcmp ludcmp proc~swap swap proc~ludcmp->proc~swap proc~outerprod outerprod proc~ludcmp->proc~outerprod

Called by

proc~~ludcmp~~CalledByGraph proc~ludcmp ludcmp proc~matinv matinv proc~matinv->proc~ludcmp 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 ludcmp(a,indx,d)
    !+Calculates LU decomposition
    real(double), dimension(:,:),intent(INOUT):: a
    integer,dimension(:),  intent(OUT)        :: indx
    real(double),                intent(OUT)  :: d
    real(double), dimension(size(a,1))        :: vv
    integer,dimension(1)                      :: imaxloc
    integer :: j,n,imax
    n=size(indx)
    d=1.0
    vv=maxval(abs(a),dim=2)
    if(any(vv.eq.0.))stop 'singular matrix in ludcmp'
    vv=1.d0/vv
    do j=1,n
       imaxloc=maxloc(vv(j:n)*abs(a(j:n,j)))
       imax=(j-1)+imaxloc(1)
       if (j /= imax) then
          call swap(a(imax,:),a(j,:))
          d=-d
          vv(imax)=vv(j)
       endif
       indx(j)=imax
       if (a(j,j) == 0.0) a(j,j)=1.0d-20
       a(j+1:n,j)=a(j+1:n,j)/a(j,j)
       a(j+1:n,j+1:n)=a(j+1:n,j+1:n)-outerprod(a(j+1:n,j),a(j,j+1:n))
    enddo
  end subroutine ludcmp