CalculateExponentialMatrix Subroutine

private subroutine CalculateExponentialMatrix(delt)

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=r64)				::	delt
Called By

Help
Source Code

CalculateExponentialMatrix
Source Code

SUBROUTINE CalculateExponentialMatrix(delt)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         Russ Taylor
          !       DATE WRITTEN   June 1990
          !       MODIFIED       Dec 1995, Apr 1996, RKS; June 2000 RKS
          !       RE-ENGINEERED  June 1996, RKS; Nov 1999, LKL;

          ! PURPOSE OF THIS SUBROUTINE:
          ! This subroutine computes the exponential matrix exp(AMat*delt) for
          ! use in the state space method for the calculation of CTFs.

          ! METHODOLOGY EMPLOYED:
          ! Uses the method of Taylor expansion combined with scaling and
          ! squaring to most efficiently compute the exponential matrix.  The
          ! steps in the procedure are outlined in Seem's dissertation in
          ! Appendix A, page 128.  Exponential matrix multiplication modified
          ! to take advantage of the characteristic form of AMat.  AMat starts
          ! out as a tri-diagonal matrix.  Each time AMat is raised to a higher
          ! power two extra non-zero diagonals are added.  ExponMatrix now
          ! recognizes this.  This should speed up the calcs somewhat.  Also, a
          ! new cut-off criteria based on the significant figures of double-
          ! precision variables has been added.  The main loop for higher powers
          ! of AMat is now stopped whenever these powers of AMat will no longer
          ! add to the summation (AExp) instead ofstopping potentially at the
          ! artifical limit of AMat**100.

          ! REFERENCES:
          ! Seem, J.E.  "Modeling of Heat Transfer in Buildings",
          !  Department of Mechanical Engineering, University of
          !  Wisconsin-Madison, 1987.
          ! Strand, R.K. "Testing Design Description for the CTF
          !  Calculation Code in BEST", BSO internal document,
          !  May/June 1996.

          ! USE STATEMENTS:
          ! none

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


          ! SUBROUTINE ARGUMENT DEFINITIONS:
  REAL(r64) :: delt  ! Time step of the resulting CTFs

          ! SUBROUTINE PARAMETER DEFINITIONS:

  REAL(r64), PARAMETER :: DPLimit = 1.0d-20
          ! This argument is nice, but not sure it's accurate -- LKL Nov 1999.
          ! Parameter set to the significant figures limit of double
          ! precision variables plus a safety factor.- The argument for setting this parameter to 1E-20 involves the
          ! number of significant figures for REAL(r64) variables which is 16 and the largest power to which
          ! AMat will be raised which is 100.  This would be a factor of 1E-18.  A factor of "safety" of another 100
          ! arrives at the value chosen.  It is argued that if one number is 1E-16 larger than a second number, then
          ! adding the second to the first will not effect the first.  However, on the conservative side, there could
          ! be up to 100 numbers which might, added together, still could effect the original number.  Each
          ! successive power of AMat will have terms smaller than the previous power.  Thus, when the ratio between
          ! the terms of the latest power of AMat and the total (AExp) is less than DPLim, all further powers of
          ! AMat will have absolutely no effect on the REAL(r64) value of AExp.  Thus, there is no need to
          ! continue the calculation.  In effect, AExp has "converged".  In REAL(r64)ity, 1E-16 would probably guarantee
          ! convergence since AMat terms drop off quickly, but the extra powers allows for differences between
          ! computer platforms.

          ! INTERFACE BLOCK SPECIFICATIONS
          ! na

          ! DERIVED TYPE DEFINITIONS
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  REAL(r64)                         :: AMatRowNorm      ! Row norm for AMat
  REAL(r64)                         :: AMatRowNormMax   ! Largest row norm for AMat
  REAL(r64), ALLOCATABLE, DIMENSION(:,:) :: AMat1    ! AMat factored by (delt/2^k)
  REAL(r64), ALLOCATABLE, DIMENSION(:,:) :: AMato    ! AMat raised to the previous power (power of AMat1-1)
  REAL(r64), ALLOCATABLE, DIMENSION(:,:) :: AMatN    ! Current value of AMat raised to power n (n = 1,2...)
  LOGICAL                                       :: Backup   ! Used when numerics get to small in Exponentiation
  REAL(r64)                                     :: CheckVal ! Used to avoid possible overflow from Double->REAL(r64)->Integer
  REAL(r64)                         :: fact     ! Intermediate calculation variable (delt/2^k)
  INTEGER                                       :: i        ! Loop counter
  INTEGER                                       :: ic       ! Loop counter
  INTEGER                                       :: ict      ! Loop counter
  INTEGER                                       :: idm      ! Loop counter
  INTEGER                                       :: ir       ! Loop counter
  INTEGER                                       :: isq      ! Loop counter
  INTEGER                                       :: j        ! Loop counter
  INTEGER                                       :: k        ! Power of 2 which is used to factor AMat
  INTEGER                                       :: l        ! Theoretical power to which the A matrix must be
                                                            ! raised to accurately calculate the exponential matrix
  LOGICAL                                       :: SigFigLimit      ! Significant figure limit logical, true
                                                                    ! when exponential calculation loop can be exited (i.e.
                                                                    ! the significant figure limit for REAL(r64)
                                                                    ! variables reached)

          ! FLOW:
  ALLOCATE(AMat1(rcmax,rcmax))
  ALLOCATE(AMato(rcmax,rcmax))
  ALLOCATE(AMatN(rcmax,rcmax))

  ! Subroutine initializations.  AMat is assigned to local variable AMat1 to
  ! avoid the corruption of the original AMat 2-d array.
  AMat1=AMat

  !  Other arrays are initialized to zero.
  AExp=0.0D0
  AMato=0.0D0
  AMatN=0.0D0


          ! Step 1, page 128 (Seem's thesis):  Compute the matrix row norm.
          ! See equation (A.3) which states that the matrix row norm is the
          ! maximum summation of the elements in a row of AMat multiplied by
          ! the time step.

  AMatRowNormMax = 0.0D0            ! Start of Step 1 ...

  DO i = 1,rcmax

    AMatRowNorm = 0.0D0
    DO j = 1,rcmax
      AMatRowNorm = AMatRowNorm+ABS(AMat1(i,j))
    END DO

    AMatRowNorm = AMatRowNorm*delt

    AMatRowNormMax=MAX(AMatRowNormMax,AMatRowNorm)

  END DO                            ! ... end of Step 1.

          ! Step 2, page 128:  Find smallest integer k such that
          ! AMatRowNormMax< = 2^k

  k = INT(LOG(AMatRowNormMax)/LOG(2.0D0))+1 !Objexx:Num Handle AMatRowNormMax=0

          ! Step 3, page 128:  Divide (AMat*delt) by 2^k.  This section of code
          ! takes advantage of the fact that AMat is tridiagonal.  Thus, it
          ! only factors the elements of the AMat that are known to be non-zero.

  fact  = delt/(2.0d0**k)   ! Start of Step 3 ...
  AMat1 = AMat1*fact        ! ... end of Step 3.

          ! Step 4, page 128:  Calculate l, the highest power to which AMat
          ! must be taken theoretically to accurately calculate its exponential.
          ! This is based on a paper by Cadzow and Martens ("Discrete-Time and
          ! Computer Control Systems",Prentice-Hall, pp. 389-390, 1970).  This
          ! number is now used as the maximum power to which AMat must be
          ! raised in order to calculate the exponential matrix.  A new cut-off
          ! criteria based on the number of significant figures in a double-
          ! precision variable is used as a more practical limit on the
          ! exponentiation algorithm.

  CheckVal = MIN(3.0d0*AMatRowNormMax+6.d0,100.0d0)
  l        = INT(CheckVal)

          ! Step 5, page 128:  Calculate the exponential.  First, add the
          ! linear term to the identity matrix.
   AExp = AMat1 + IdenMatrix        ! Start of Step 5 ...

          ! Now, add successive terms to the expansion as per the standard
          ! exponential formula.  AMato contains the last "power" of AMat
          ! which saves the program from having to remultiply the entire power
          ! of AMat each time.  Since this is still the linear power of AMat,
          ! AMat1 is still tridiagonal in nature.
  AMato = AMat1

  i = 1   ! Initialize the counter for the following DO loop

          ! The following DO WHILE loop continues to raise AMat to successive
          ! powers and add it to the exponential matrix (AExp).
  DO WHILE (i < l) ! Begin power raising loop ...

    i = i+1               ! Increment the loop counter
    SigFigLimit = .true.  ! Set the significant factor limit flag

    DO ir = 1, rcmax    ! Begin matrix multiplication loop ...
          ! The following matrix multiplication could be "optimized" since
          ! for one-dimensional heat transfer AMat is 3-diagonal, AMat squared
          ! is 5-diagonal, etc.  However, the code can be much simpler if we
          ! ignore this fact and just do a generic matrix multiplication.
          ! For 2-D heat transfer, the number of off-diagonal non-zero terms
          ! is slightly more complicated as well.
      DO ic = 1, rcmax
        AMatN(ir,ic) = 0.0d0
        DO ict = 1, rcmax
          ! Make sure the next term won't cause an underflow.  If it will end up being
          ! so small as to go below TinyLimit, then ignore it since it won't add anything
          ! to AMatN anyway.
          IF (ABS(AMat1(ict,ic)) > TinyLimit) THEN
            IF (ABS(AMato(ir,ict)) > ABS(REAL(i,r64)*TinyLimit/AMat1(ict,ic))) &
              AMatN(ir,ic) = AMatN(ir,ic)+AMato(ir,ict)*AMat1(ict,ic)/REAL(i,r64)
          END IF
        END DO
      END DO
    END DO          ! ... end of matrix multiplication loop.

          ! Update AMato and AExp matrices
    AMato = AMatN
    AExp  = AExp + AMato

          ! The next DO loop tests the significant figures limit criteria to
          ! see if any values in AExp are still changing appreciably.
    DO ir = 1, rcmax
      DO ic = 1, rcmax
          ! Test of limit criteria:
        IF (ABS(AExp(ir,ic)) > TinyLimit) THEN  ! Next line divides by AExp entry so it
                                                ! must be checked to avoid dividing by zero.
          ! If the ratio between any current element in the power
          ! of AMat and its corresponding element in AExp is
          ! greater than the number which might effect the overall
          ! exponential matrix based on stability criteria, then
          ! continue raising AMat to another power (SigFigLimit = false).

          IF (ABS(AMato(ir,ic)/AExp(ir,ic)) > DPLimit) THEN
            SigFigLimit = .FALSE.
            EXIT ! DO loop (anytime SigFigLimit is false, AMat must continue to be raised another power)
          END IF

        ELSE        ! There are still elements of AExp which are zero, so
                    ! the raising of AMat to higher powers should continue.

          SigFigLimit = .FALSE.
          EXIT ! DO loop (anytime SigFigLimit is false, AMat must continue to be raised another power)

        END IF
      END DO
      IF (.NOT.SigFigLimit) EXIT ! DO loop (anytime SigFigLimit is false, AMat must continue to be raised another power)
    END DO

          ! Compute next term, only if necessary.  If SigFigLimit is still true,
          ! then all of the new terms being added to AExp are too small to
          ! affect it.  Thus, there is no need to continue this do loop further.

    IF (SigFigLimit) i = 100  ! SigFigLimit is still true, set i to maximum possible
                              ! value of l (100).

  END DO            ! ... end of power raising loop and Step 5.

          ! Step 6, page 128:
          ! Square AExp "k times" to obtain the actual exponential matrix
          ! (remember that AExp was scaled earlier in this routine).

  DO isq = 1,k      ! Begin squaring DO loop and Step 6 ...

          ! Use AMato to store the old values of AExp
    AMato  = AExp
    Backup = .true.
    AExp   = 0.0D0

          ! Multiply the old value of AExp (AMato) by itself and store in AExp.
    DO ir = 1,rcmax
      DO ic = 1,rcmax
        DO idm = 1,rcmax
          IF (ABS(AMato(ir,idm)*AMato(idm,ic)) > TinyLimit) THEN
            AExp(ir,ic) = AExp(ir,ic)+AMato(ir,idm)*AMato(idm,ic)
            Backup=.false.
          END IF
        END DO
      END DO
    END DO
    ! Backup is true when every item of AExp didnt pass the TinyLimit test
    IF (Backup) THEN
      AExp=AMato
      EXIT
    END IF

  END DO            ! ... end of squaring loop and Step 6.

  DEALLOCATE(AMat1)
  DEALLOCATE(AMato)
  DEALLOCATE(AMatN)

  RETURN            ! AExp is now the exponential of AMat.

END SUBROUTINE CalculateExponentialMatrix
All Procedures
