SUBROUTINE ludcmp(a,n,np,indx,d)
! Given a matrix a(1:n,l:n), with physical dimension np by np, this routine replaces it by
! the LU decomposition of a rowwise permutation of itself. a and n are input. a is output,
! arranged as in equation (2.3.14) above; indx(l:n) is an output vector that records the
! row permutation effected by the partial pivoting; d is output as +1 or -1 depending on whether
! the number of row interchanges was even or odd, respectively. This routine is used' in
! combination with lubksbto solve linear equations or invert a matrix.

INTEGER n,np,indx(n),NMAX
REAL d,a(np,np), TINY
PARAMETER(NMAX=500, TINY=1.0e-20)
INTEGER i,imax,j,k
REAL aamax, dum, sum, vv(NMAX)								! vv stores the implicit scaling of each row.

d=1.0d0														! No row interchanges yet. 

do i=1,n													! Loop over rows to get the implicit scaling information.
	aamax=0.0d0
	do j=1,n
		if (abs(a(i,j)).gt.aamax) aamax = abs(a(i,j))
	end do
	if (aamax.eq.0.0d0) pause' singular matrix in ludcmp'	! No nonzero largest element.
	vv(i)=1.0d0 /aamax										! Save the scaling.
end do

do j=1,n													! This is the loop over columns of Crout's method.
	do i=1,j-1												! This is equation (2.3.12) except for i = j.
		sum=a(i,j)
		do k=1,i-1
			sum = sum-a(i,k)*a(k,j)
		enddo
		a(i,j)=sum
	end do

	aamax=0.0d0												! Initialize for the search for largest pivot element.
	
	do i=j ,n												! This is i = j of equation (2.3.12) and i = j + 1. . .N
		sum=a(i ,j)
		do k=1,j-1
			sum=sum-a(i ,k)*a (k,j)
		end do

		a(i,j)=sum
		dum=vv(i)*abs(sum)									! Figure of merit for the pivot.
		
		if (dum.ge.aamax) then								! Is it better tha n the best so far?
			imax=i
			aamax=dum
		endif
	end do
	if (j.ne.imax)then										! Do we need to interchange rows?
		do k=1,n											! Yes, do so...
			dum=a(imax,k)
			a(imax,k)=a(j,k)
			a(j ,k)=dum
		end do
		d=-d												! ...and change the parity of d.
		vv(imax)=vv(j)										! Also interchange the scale factor.
	endif

	indx(j)=imax
	
	if(a(j,j).eq.0.0d0)a(j,j)=TINY							! if the pivot element is zero the matrix is singular 
															! (at least to the precision of the algorithm).
															! For some applications on singular matrices, 
															! it is desirable to substitute TINY !for zero.

	if(j.ne.n)then											! Now, finally, divide by the pivot element.
		dum=1./a(j,j)
		do i=j+1,n
			a(i,j)=a(i,j)*dum
		end do
	endif
end do														! Go back for the next column in the reduction.
return
END


SUBROUTINE lubksb(a,n,np,indx,b)
INTEGER n,np,indx(n)
REAL a(np,np),b(n)

!Solves the set of n linear equations A .X = B. Here a is input, not as the matrix A but
!rather as its LU decomposition, determined by the routine ludcmp. indx is input as the
!permutation vector returned by ludcmp. b(l: n) is input as the right-hand side vector B,
!and returns with the solution vector X. a, n, np, and indx are not modified by this routine
!and can be left in place for successive calls with different right-hand sides b. This routine
!takes into account the possibility that b will begin with many zero elements, so it is efficient
!for use in matrix inversion.

INTEGER i,ii,j,ll
REAL sum
ii=0									! When ii is set to a positive value, it will become the index
do i=1,n								! of the first nonvanishing element of b. We now do
	ll=indx(i)							! the forward substitution, equation (2.3.6). The only new
	sum=b(11)							! wrinkle is to unscramble the permutation as we go,
	b(ll)=b(i)
	if (ii.ne.0)then
		do j=ii,i-1
			sum=sum-a(i,j)*b(j)
		end do
		else if (sum.ne.0.0) then
			ii=i						! A nonzero element was encountered, so from now on we will
	endif								! have to do the sums in the loop above.
	b(i)=sum
end do

do i=n,1,-1								! Now we do the backsubstitution, equation (2.3.7)

	sum=b(i)

	do j=i+1,n
		sum=sum-a(i,j)*b(j)
	end do

	b(i)=sum/a(i,i)						! Store a component of the solution vector X.

end do
return
END