CalculateInverseMatrix Subroutine

private subroutine CalculateInverseMatrix()

Arguments

None
Called By

Source Code

SUBROUTINE CalculateInverseMatrix

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Rick Strand
          !       DATE WRITTEN   Dec 1995
          !       MODIFIED       June 2000 RKS (made routine generic to allow for 2-D solutions)
          !       RE-ENGINEERED  June 1996, February 1997 RKS

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine computes the inverse of AMat for use
          ! in the calculation of the CTFs.

          ! METHODOLOGY EMPLOYED:
          ! Uses row elimination to zero the off-diagonal terms of
          ! AMat while performing the same operations on another
          ! matrix which starts as the identity matrix.  Once AMat
          ! has been converted to an identity matrix(I), the other
          ! matrix which started as the I will then be the inverse
          ! of A.  This algorithm has been customized for a
          ! tri-diagonal matrix.

          ! REFERENCES:
          ! Any linear algebra test (this is a generic routine).

          ! USE STATEMENTS:
          ! none

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

          ! SUBROUTINE ARGUMENT DEFINITIONS:
          ! na

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  REAL(r64), ALLOCATABLE, DIMENSION(:,:) :: AMat1  ! Intermediate calculation matrix equivalent at first to AMat
  INTEGER                                       :: ic     ! Loop counter
  INTEGER                                       :: ir     ! Loop counter
  INTEGER                                       :: irr    ! Loop counter

          ! FLOW:

          ! Subroutine initializations ...
  ALLOCATE(AMat1(rcmax,rcmax))

  AMat1 = AMat          ! Set AMat1 = AMat to avoid AMat changes
  AInv  = IdenMatrix    ! Set AInv to Identity Matrix

          ! Use Gaussian elimination to zero all of the elements of AMat left
          ! of the diagonal.
          ! This DO loop will cycle through each of the rows in AMat except the
          ! last row which is handled later because it does not have to be used
          ! to eliminate any other rows.  The index ir is the current row
          ! number and also the column of the current diagonal element.

  DO ir = 1, rcmax-1    ! Begin forward elimination loop ...

          ! Factor all of the elements of the row being used to zero the next
          ! row in both AMat and AInv by the diagonal element of this row.
          ! We should only need to factor the elements to the right of the
          ! diagonal since those to the right of it should be zero.
    DO ic = ir+1, rcmax
      AMat1(ir,ic) = AMat1(ir,ic)/AMat1(ir,ir)
    END DO

          ! In the forward elimination process, all the elements in AInv to the
          ! right of the diagonal are zero so they do not need to be factored.
    DO ic = 1, ir
      AInv(ir,ic) = AInv(ir,ic)/AMat1(ir,ir)
    END DO

    AMat1(ir,ir) = 1.0d0      ! By definition, the diagonal of AMat is now 1.

          ! Use this factored row to eliminate the off-diagonal element of the
          ! rows below the current one (ir)...

    DO irr = ir+1, rcmax    ! Start of row reduction loop...

      DO ic = ir+1, rcmax
        AMat1(irr,ic) = AMat1(irr,ic) - AMat1(irr,ir)*AMat1(ir,ic)
      END DO

          ! Now, determine the effect on the next row of AInv.  Again, all of
          ! the elements in AInv to the right of the diagonal are zero, so they
          ! can be ignored.

      DO ic = 1,ir
        AInv(irr,ic) = AInv(irr,ic) - AMat1(irr,ir)*AInv(ir,ic)
      END DO

      AMat1(irr,ir) = 0.0d0 ! By definition, the element to the left of the
                            ! diagonal in the next row of AMat is now zero.

    END DO                  ! ...end of row reduction loop

  END DO            ! ... end of the forward elimination loop.

          ! Factor the last row of AInv by the current value of the last
          ! diagonal element of AMat. After this is done, all of the diagonal
          ! elements of AMat are unity and all of the elements in AMat left of
          ! the diagonal are zero.

  DO ic = 1,rcmax
    AInv(rcmax,ic) = AInv(rcmax,ic)/AMat1(rcmax,rcmax)
  END DO
  AMat1(rcmax,rcmax) = 1.0D0

          ! Now, use back substitution to eliminate the elements to the right
          ! of the diagonal in AMat.  The procedure is similar to the forward
          ! elimination process except that we only have to operate on AInv,
          ! though now all of the columns of AInv may be non-zero.

          ! This DO loop will cycle through the remaining rows which are not
          ! yet diagonalized in reverse order.  Note that the only effect on
          ! AMat is that the off-diagonal element is zeroed.  The diagonal
          ! (which has already been set to unity) is not effected by this row
          ! elimination process.
          ! In the following code ir is the column being zeroed and irr is the
          ! row being worked on

  DO ir = rcmax, 2, -1  ! Begin reverse elimination loop ...
    DO irr = 1, ir-1
      DO ic = 1,rcmax
        AInv(irr,ic) = AInv(irr,ic)-AMat1(irr,ir)*AInv(ir,ic)
      END DO
      AMat1(irr,ir) = 0.0d0
    END DO
  END DO                ! ... end of reverse elimination loop.

          ! At this point, AMat1 is equal to the identity matrix (I)
          ! and AInv is equal to the inverse of AMat.

  DEALLOCATE(AMat1)

  RETURN

END SUBROUTINE CalculateInverseMatrix
All Procedures
