CalcMatrixInverse Subroutine

private subroutine CalcMatrixInverse(Matrix, InvMatrix)

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=r64),	intent(inout),		DIMENSION(:,:)	::	Matrix
real(kind=r64),	intent(inout),		DIMENSION(:,:)	::	InvMatrix
Called By

Help
Source Code

CalcMatrixInverse
Source Code

SUBROUTINE CalcMatrixInverse(Matrix,InvMatrix)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Jakob Asmundsson
          !       DATE WRITTEN   January 1999
          !       MODIFIED       September 2000 (RKS for EnergyPlus)
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! To find the inverse of Matrix, using partial pivoting.

          ! METHODOLOGY EMPLOYED:
          ! Inverse is found using partial pivoting and Gauss elimination

          ! REFERENCES:
          ! Any Linear Algebra book

          ! USE STATEMENTS:
          ! na

  IMPLICIT NONE    ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENTS:
  REAL(r64), DIMENSION(:,:), INTENT(INOUT) :: Matrix      ! Input Matrix
  REAL(r64), DIMENSION(:,:), INTENT(INOUT) :: InvMatrix   ! Inverse of Matrix

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  INTEGER      :: DimOfMatrix ! Matrix dimension
  REAL(r64), ALLOCATABLE, DIMENSION(:,:) :: Identity    ! Identity matrix
  INTEGER,   ALLOCATABLE, DIMENSION(:)   :: p           ! Vector containing the
                                                               ! pivot order
  INTEGER         :: temp        ! Temporary variable
  REAL(r64)  :: mm          ! Multiplier
  REAL(r64)  :: pivot       ! Pivot value
  INTEGER         :: piv         ! Pivot index
  INTEGER         :: i           ! Loop counter
  INTEGER         :: j           ! Loop counter
  INTEGER         :: k           ! Loop counter

          ! FLOW:
  DimOfMatrix = SIZE(Matrix, DIM=1)
  ALLOCATE(Identity(DimOfMatrix,DimOfMatrix))
  ALLOCATE(p(DimOfMatrix))

  Identity  = 0.0d0
  InvMatrix = 0.0d0

  DO i = 1, DimOfMatrix
    Identity(i,i) = 1.0d0
  END DO

  p=(/ (I, I=1,DimOfMatrix) /)

  DO j = 1, DimOfMatrix-1
    pivot = ABS(Matrix(p(j),j))
    piv   = j
    temp  = p(j)
    DO i = j+1, DimOfMatrix
      IF (ABS(Matrix(p(i),j))>pivot) THEN
        pivot = ABS(Matrix(p(i),j))
        piv   = i
      END IF
    END DO

    p(j)   = p(piv)
    p(piv) = temp
    DO i = j+1, DimOfMatrix
      mm               = Matrix(p(i),j)/Matrix(p(j),j)
      Matrix(p(i),j)   = 0.0d0
      Identity(p(i),:) = Identity(p(i),:) - mm*Identity(p(j),:)
      DO k = j+1, DimOfMatrix
        Matrix(p(i),k) = Matrix(p(i),k) - mm*Matrix(p(j),k)
      END DO
    END DO
  END DO

  DO i = DimOfMatrix, 1, -1
    InvMatrix(i,:) = Identity(p(i),:)
    DO j = i+1, DimOfMatrix
      InvMatrix(i,:) = InvMatrix(i,:) - Matrix(p(i),j)*InvMatrix(j,:)
    END DO
    InvMatrix(i,:) = InvMatrix(i,:)/Matrix(p(i),i)
  END DO

  DEALLOCATE(p)
  DEALLOCATE(Identity)

  RETURN

END SUBROUTINE CalcMatrixInverse
All Procedures
