| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=r64) | :: | ajac(10,10) | ||||
| integer | :: | n | ||||
| integer | :: | indx(10) | ||||
| real(kind=r64) | :: | d |
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
SUBROUTINE LUdecomposition(ajac,n,indx,d)
! SUBROUTINE INFORMATION:
! AUTHOR F. Winkelmann, adapted from Numerical Recipes
! DATE WRITTEN February 2000
! MODIFIED na
! RE-ENGINEERED na
! PURPOSE OF THIS SUBROUTINE:
! Performs LU decomposition of a matrix.
! METHODOLOGY EMPLOYED:
! na
! REFERENCES:
! na
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
! SUBROUTINE ARGUMENT DEFINITIONS:
REAL(r64) :: ajac(10,10) ! As input: matrix to be decomposed;
! as output: decomposed matrix
INTEGER :: n ! Dimension of matrix
INTEGER :: indx(10) ! Vector of row permutations
REAL(r64) :: d ! +1 if even number of row interchange is even, -1
! if odd
! SUBROUTINE PARAMETER DEFINITIONS:
! na
! INTERFACE BLOCK SPECIFICATIONS:
! na
! DERIVED TYPE DEFINITIONS:
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
INTEGER :: i,j,k ! Counters
INTEGER :: imax ! Temporary variable
! as output: decomposed matrix
REAL(r64) :: vv(500) ! Stores the implicit scaling of each row
REAL(r64) :: aamax ! Absolute value of largest element of matrix
REAL(r64) :: dum ! Temporary variable
REAL(r64) :: sum ! Sum of products of matrix elements
! FLOW
d=1.0d0
do i=1,n
aamax=0.0d0
do j=1,n
if (abs(ajac(i,j)).gt.aamax) aamax=abs(ajac(i,j))
end do
if (aamax.eq.0.0d0) CALL ShowFatalError('Singular matrix in LUdecomposition, window calculations')
vv(i)=1.0d0/aamax
end do
do j=1,n
do i=1,j-1
sum=ajac(i,j)
do k=1,i-1
sum=sum-ajac(i,k)*ajac(k,j)
end do
ajac(i,j)=sum
end do
aamax=0.0d0
do i=j,n
sum=ajac(i,j)
do k=1,j-1
sum=sum-ajac(i,k)*ajac(k,j)
end do
ajac(i,j)=sum
dum=vv(i)*abs(sum)
if (dum.ge.aamax) then
imax=i
aamax=dum
endif
end do
if (j.ne.imax)then
do k=1,n
dum=ajac(imax,k)
ajac(imax,k)=ajac(j,k)
ajac(j,k)=dum
end do
d=-d
vv(imax)=vv(j)
endif
indx(j)=imax
if(ajac(j,j).eq.0.0d0) ajac(j,j)=rTinyValue
if(j.ne.n)then
dum=1.0d0/ajac(j,j)
do i=j+1,n
ajac(i,j)=ajac(i,j)*dum
end do
endif
end do
return
END SUBROUTINE LUdecomposition