Source Code for Project Mixer

Functions in Class Base.Container.Reactor     [Go to Top]

! Base_Container_Reactor.f
! - DLL routines for class <Component>Base.Container.Reactor
! DATE: Monday, September 10, 2001  TIME: 03:57:13 PM
! The skeleton of this file is generated by SansGUI(tm)

! Attribute indices in class version [1.0.alpha.7]
! 1: fConcentration - Concentration
! 2: rReactor - Reactor Table
! 3: rConstant - Constant Matrix
! 4: iPartIndex - Part Index (1-Based)

! ======================================================================
! SG_xInitSize - Resize for Init
! ----------------------------------------------------------------------
      integer function SG_xInitSize_Base_Container_Reactor(self,        &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xInitSize_Base_Container_Reactor
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"

      integer, dimension(*) :: iNumReact
      POINTER(PTR_iNumReact, iNumReact)

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xInitSize_Base_Container_Reactor = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! register the reactor part for solver to resize matrices and tables
      PTR_zValues = simCtrl%pzValues
      PTR_iNumReact = zValues(SG_NDX_CTRL_INUMREACT)%vData
      if (iNumReact(1) .lt. 0) then
          iNumReact(1) = 1
      else
          iNumReact(1) = iNumReact(1) + 1;
      end if

      SG_xInitSize_Base_Container_Reactor = SG_R_OK
      return
      end

! ======================================================================
! SG_xInit - Initialization
! ----------------------------------------------------------------------
      integer function SG_xInit_Base_Container_Reactor(self,            &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xInit_Base_Container_Reactor
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"
      integer :: iPartNdx
      integer :: iPartSN
      integer, dimension(*) :: iPartIndex
      integer, dimension(*) :: iNumReact
      POINTER(PTR_iPartIndex, iPartIndex)
      POINTER(PTR_iNumReact, iNumReact)
      integer, parameter :: SG_NDX_OBJ_REACTORTABLE = 1
      integer, parameter :: SG_NDX_IPARTINDEX = 4

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xInit_Base_Container_Reactor = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! set user type for future identification - to prevent class name
      ! comparison during component class run-time type checking
      call MIX_SET_USER_TYPE(self, MIX_USER_TYPE_REACTOR)

      ! check to see if there are two reference objects
      if (iRefObjs .ne. 2) then
          cMessage = 'Need a Reactor Table and a Constant Matrix.'C
          !important - reset the number of reactors in the simControl object
          PTR_zValues = simCtrl%pzValues
          PTR_iNumReact = zValues(SG_NDX_CTRL_INUMREACT)%vData
          iNumReact(1) = 0
          SG_xInit_Base_Container_Reactor = SG_R_STOP
          return
      end if

      ! register the reactor in the reactor table and find the 1-based index
      ! registerReactor is implemented in Table_Reactor.f and is called via
      ! the function declared in Mixer_1_0F.h
      iPartSN = self%nCmpnNo
      PTR_refObject = pRefObjs(SG_NDX_OBJ_REACTORTABLE)
      iPartNdx = registerReactor(refObject, iPartSN)
      if (iPartNdx .lt. 1) then
          cMessage = 'Cannot register this part in the reactor table.'C
          PTR_zValues = simCtrl%pzValues
          PTR_iNumReact = zValues(SG_NDX_CTRL_INUMREACT)%vData
          iNumReact(1) = 0
          SG_xInit_Base_Container_Reactor = SG_R_STOP
          return
      else
          ! record the 1-based index for bi-directional reference
          PTR_zValues = self%pzValues
          PTR_iPartIndex = zValues(SG_NDX_IPARTINDEX)%vData
          iPartIndex(1) = iPartNdx
      end if

      SG_xInit_Base_Container_Reactor = SG_R_OK
      return
      end

! ======================================================================
! SG_xPreEval - Pre-Evaluation
! ----------------------------------------------------------------------
      integer function SG_xPreEval_Base_Container_Reactor(self,         &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xPreEval_Base_Container_Reactor
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"
      integer :: i
      integer, dimension(*) :: iLinkInfo
      real*4, dimension(*) :: fConcentration
      real*4, dimension(*) :: fLinkConc
      type (SG_VALU), dimension(*) :: lnkValues
      POINTER(PTR_fConcentration, fConcentration)
      POINTER(PTR_lnkValues, lnkValues)
      POINTER(PTR_fLinkConc, fLinkConc)
      POINTER(PTR_iLinkInfo, iLinkInfo)
      integer, parameter :: SG_NDX_FCONCENTRATION = 1

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xPreEval_Base_Container_Reactor = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! deposit Reactor Concentration to all the output links
      PTR_zValues = self%pzValues
      PTR_fConcentration = zValues(SG_NDX_FCONCENTRATION)%vData
      do i = 1, iLnkObjs
          PTR_lnkObject = pLnkObjs(i)
          PTR_lnkValues = lnkObject%pzValues
          PTR_iLinkInfo = lnkValues(SG_NDX_LINK_ILINKINFO)%vData
          if (iLinkInfo(1) .eq. SG_LINK_OUT) then
              PTR_fLinkConc = lnkValues(SG_NDX_LINK_FCONCENT)%vData
              fLinkConc(1) = fConcentration(1)
          end if
      end do

      SG_xPreEval_Base_Container_Reactor = SG_R_OK
      return
      end

! ======================================================================
! SG_xEval - Evaluation
! ----------------------------------------------------------------------
      integer function SG_xEval_Base_Container_Reactor(self,            &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xEval_Base_Container_Reactor
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xEval_Base_Container_Reactor = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! start loading the constant matrix and RHS vector,
      ! using the law of conservation
      if (iRefObjs .lt. 2) then
          cMessage =
     &    'Constant Matrix and Reactor Table objects are required.'C
          ! important - reset the number of reactors in the simControl object
          PTR_zValues = simCtrl%pzValues
          PTR_iNumReact = zValues(SG_NDX_CTRL_INUMREACT)%vData
          iNumReact(1) = 0
          SG_xEval_Base_Container_Reactor = SG_R_STOP
          return
      end if

      PTR_reactorTable = pRefObjs(SG_NDX_OBJ_REACTORTABLE)
      PTR_constMatrix = pRefObjs(SG_NDX_OBJ_CONSTMATRIX)

      do i = 1, iLnkObjs
          ! going through all the links and adjacent objects
          PTR_lnkObject = pLnkObjs(i)
          PTR_lnkValues = lnkObject%pzValues
          PTR_iLinkInfo = lnkValues(SG_NDX_LINK_ILINKINFO)%vData
          bInput = (iLinkInfo(1) .eq. SG_LINK_IN)
          if (bInput) then
              PTR_adjObject = pAdjObjs(i)
              iType = MIX_GET_USER_TYPE(adjObject)

              select case (iType)
              case (MIX_USER_TYPE_REACTOR)
                PTR_zValues = self%pzValues
                PTR_iPartIndex = zValues(SG_NDX_IPARTINDEX)%vData
                PTR_adjObject = pAdjObjs(i)
                PTR_adjValues = adjObject%pzValues
                PTR_iAdjPartNdx = adjValues(SG_NDX_IPARTINDEX)%vData
                PTR_fFlowRate = lnkValues(SG_NDX_LINK_FFLOWRATE)%vData
                call loadMatrixConstant(constMatrix, iPartIndex(1),
     &               iAdjPartNdx(1), fFlowRate(1), bInput )

              case (MIX_USER_TYPE_SOURCE)
                PTR_fFlowRate = lnkValues(SG_NDX_LINK_FFLOWRATE)%vData
                PTR_fConcent = lnkValues(SG_NDX_LINK_FCONCENT)%vData
                if (fFlowRate(1).lt. 0. .or. fConcent(1).lt. 0.) then
                    cMessage = 'Source to the reactor has not been initialized.'C
                    ! important - reset the number of reactors in SimControl
                    PTR_zValues = simCtrl%pzValues
                    PTR_iNumReact = zValues(SG_NDX_CTRL_INUMREACT)%vData
                    iNumReact(1) = 0
                    SG_xEval_Base_Container_Reactor = SG_R_STOP
                    return
                endif
                PTR_zValues = self%pzValues
                PTR_iPartIndex = zValues(SG_NDX_IPARTINDEX)%vData
                call loadTableConstant(reactorTable, iPartIndex(1),
     &               fConcent(1) * fFlowRate(1), bInput )

              case (MIX_USER_TYPE_SINK)
                ! do nothing

              case DEFAULT
                ! do nothing
              end select
          else    ! output
              PTR_adjObject = pAdjObjs(i)
              iType = MIX_GET_USER_TYPE(adjObject)
              select case (iType)
              case (MIX_USER_TYPE_REACTOR)
                ! load constant to the main diagnal cell
                ! same behavior for both reactor and sink in the output
                PTR_zValues = self%pzValues
                PTR_iPartIndex = zValues(SG_NDX_IPARTINDEX)%vData
                PTR_fFlowRate = lnkValues(SG_NDX_LINK_FFLOWRATE)%vData
                call loadMatrixConstant(constMatrix, iPartIndex(1),
     &               iPartIndex(1), fFlowRate(1), bInput )

              case (MIX_USER_TYPE_SINK)
                ! load constant to the main diagnal cell
                ! same behavior for both reactor and sink in the output
                PTR_zValues = self%pzValues
                PTR_iPartIndex = zValues(SG_NDX_IPARTINDEX)%vData
                PTR_fFlowRate = lnkValues(SG_NDX_LINK_FFLOWRATE)%vData
                call loadMatrixConstant(constMatrix, iPartIndex(1),
     &               iPartIndex(1), fFlowRate(1), bInput )

              case (MIX_USER_TYPE_SOURCE)
                ! do nothing

              case DEFAULT
                ! do nothing
              end select
          end if
      end do

      SG_xEval_Base_Container_Reactor = SG_R_OK
      return
      end

! ======================================================================
! SG_xPostEval - Post-Evaluation
! ----------------------------------------------------------------------
      integer function SG_xPostEval_Base_Container_Reactor(self,        &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xPostEval_Base_Container_Reactor
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"
      integer :: i
      integer, dimension(*) :: iNumReact
      integer, dimension(*) :: iPartIndex
      real*8, dimension(*) :: dSolution
      real*4, dimension(*) :: fConcentration
      real*4, dimension(*) :: fLinkConc
      integer, dimension(*) :: iLinkInfo
      type (SG_OBJ) :: reactorTable
      type (SG_VALU), dimension(*) :: tblValues
      type (SG_VALU), dimension(*) :: lnkValues
      POINTER(PTR_reactorTable, reactorTable)
      POINTER(PTR_tblValues, tblValues)
      POINTER(PTR_lnkValues, lnkValues)
      POINTER(PTR_dSolution, dSolution)
      POINTER(PTR_fConcentration, fConcentration)
      POINTER(PTR_fLinkConc, fLinkConc)
      POINTER(PTR_iLinkInfo, iLinkInfo)
      POINTER(PTR_iNumReact, iNumReact)
      POINTER(PTR_iPartIndex, iPartIndex)
      integer, parameter :: SG_NDX_IPARTINDEX = 4
      integer, parameter :: SG_NDX_FCONCENTRATION = 1
      integer, parameter :: SG_NDX_OBJ_REACTORTABLE = 1

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xPostEval_Base_Container_Reactor = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! check the existence of the reactor table
      if (iRefObjs .lt. 2) then
          cMessage =
     &    'Constant Matrix and Reactor Table objects are required.'C
          ! important - reset the number of reactors in the simControl object
          PTR_zValues = simCtrl%pzValues
          PTR_iNumReact = zValues(SG_NDX_CTRL_INUMREACT)%vData
          iNumReact(1) = 0
          SG_xPostEval_Base_Container_Reactor = SG_R_STOP
          return
      end if

      ! fetch the solution from the reactor table
      ! notice that iPartIndex is 1-based
      PTR_reactorTable = pRefObjs(SG_NDX_OBJ_REACTORTABLE)
      PTR_tblValues = reactorTable%pzValues
      PTR_dSolution = tblValues(SG_NDX_TBL_DSOLUTION)%vData
      PTR_zValues = self%pzValues
      PTR_fConcentration = zValues(SG_NDX_FCONCENTRATION)%vData
      PTR_iPartIndex = zValues(SG_NDX_IPARTINDEX)%vData
      fConcentration(1) = SNGL(dSolution(iPartIndex(1)))

      ! deposit Reactor Concentration to all the output links
      do i = 1, iLnkObjs
          PTR_lnkObject = pLnkObjs(i)
          PTR_lnkValues = lnkObject%pzValues
          PTR_iLinkInfo = lnkValues(SG_NDX_LINK_ILINKINFO)%vData
          if (iLinkInfo(1) .eq. SG_LINK_OUT) then
              PTR_fLinkConc = lnkValues(SG_NDX_LINK_FCONCENT)%vData
              fLinkConc(1) = fConcentration(1)
          end if
      end do

! Base_Container_Sink.f
! - DLL routines for class <Component>Base.Container.Sink
! DATE: Monday, September 10, 2001  TIME: 03:57:13 PM
! The skeleton of this file is generated by SansGUI(tm)

! Attribute indices in class version [1.0.alpha.7]
! 1: fConcentration - Concentration

! ======================================================================
! SG_xInit - Initialization
! ----------------------------------------------------------------------
      integer function SG_xInit_Base_Container_Sink(self,               &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xInit_Base_Container_Sink
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xInit_Base_Container_Sink = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! set user type for future identification - to prevent class name
      ! comparison during component class run-time type checking
      call MIX_SET_USER_TYPE(self, MIX_USER_TYPE_SINK)

      SG_xInit_Base_Container_Sink = SG_R_OK
      return
      end

! ======================================================================
! SG_xPostEval - Post-Evaluation
! ----------------------------------------------------------------------
      integer function SG_xPostEval_Base_Container_Sink(self,           &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xPostEval_Base_Container_Sink
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"

      real*4, dimension(*) :: fConcentration
      real*4, dimension(*) :: fLinkConc
      type (SG_VALU), dimension(*) :: lnkValues
      POINTER(PTR_fConcentration, fConcentration)
      POINTER(PTR_lnkValues, lnkValues)
      POINTER(PTR_fLinkConc, fLinkConc)
      integer, parameter :: SG_NDX_FCONCENTRATION = 1

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xPostEval_Base_Container_Sink = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! fetch concentration from input link - one big requirement here is
      ! that this routine in all sinks shall be executed AFTER all the
      ! PostEval routines in the reactors have been executed.  Use the
      ! name order or other execution sequence control in the simControl
      ! object.
      if (iLnkObjs .gt. 0) then
          PTR_zValues = self%pzValues
          PTR_fConcentration = zValues(SG_NDX_FCONCENTRATION)%vData
          PTR_lnkObject = pLnkObjs(1)
          PTR_lnkValues = lnkObject%pzValues
          PTR_fLinkConc = lnkValues(SG_NDX_LINK_FCONCENT)%vData
          fConcentration(1) = fLinkConc(1)
      end if

! Base_Source.f
! - DLL routines for class <Component>Base.Source
! DATE: Monday, September 10, 2001  TIME: 03:57:13 PM
! The skeleton of this file is generated by SansGUI(tm)

! Attribute indices in class version [1.0.alpha.3]
! 1: fConcentration - Initial Concentration

! ======================================================================
! SG_xInit - Initialization
! ----------------------------------------------------------------------
      integer function SG_xInit_Base_Source(self,                       &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xInit_Base_Source
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xInit_Base_Source = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! set user type for future identification - to prevent class name
      ! comparison during component class run-time type checking
      call MIX_SET_USER_TYPE(self, MIX_USER_TYPE_SOURCE)

      SG_xInit_Base_Source = SG_R_OK
      return
      end

! ======================================================================
! SG_xPreEval - Pre-Evaluation
! ----------------------------------------------------------------------
      integer function SG_xPreEval_Base_Source(self,                    &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xPreEval_Base_Source
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"
      real*4, dimension(*) :: fConcentration
      real*4, dimension(*) :: fLinkConc
      type (SG_VALU), dimension(*) :: lnkValues
      POINTER(PTR_fConcentration, fConcentration)
      POINTER(PTR_lnkValues, lnkValues)
      POINTER(PTR_fLinkConc, fLinkConc)
      integer, parameter :: SG_NDX_FCONCENTRATION = 1

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xPreEval_Base_Source = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! deposit Source Concentration to the output link - should be only one
      if (iLnkObjs .gt. 0) then
          PTR_zValues = self%pzValues
          PTR_fConcentration = zValues(SG_NDX_FCONCENTRATION)%vData
          PTR_lnkObject = pLnkObjs(1)
          PTR_lnkValues = lnkObject%pzValues
          PTR_fLinkConc = lnkValues(SG_NDX_LINK_FCONCENT)%vData
          fLinkConc(1) = fConcentration(1)
      end if

! Base_Source_Variable.f
! - DLL routines for class <Component>Base.Source.Variable
! DATE: Monday, September 10, 2001  TIME: 03:57:13 PM
! The skeleton of this file is generated by SansGUI(tm)

! Attribute indices in class version [1.0.alpha.4]
! 1: fConcentration - Initial Concentration
! 2: fSteady - Steady State Concentration
! 3: fCurrent - Current Concentration

! ======================================================================
! SG_xInit - Initialization
! ----------------------------------------------------------------------
      integer function SG_xInit_Base_Source_Variable(self,              &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xInit_Base_Source_Variable
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      integer :: SG_xInit_Base_Source            ! base class function

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xInit_Base_Source_Variable = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      return
      end

! ======================================================================
! SG_xPreEval - Pre-Evaluation
! ----------------------------------------------------------------------
      integer function SG_xPreEval_Base_Source_Variable(self,           &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xPreEval_Base_Source_Variable
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"
      integer :: SG_xPreEval_Base_Source         ! base class function
      real*4 :: EXP
      real*4, dimension(*) :: fCurTime
      real*4, dimension(*) :: fInit
      real*4, dimension(*) :: fSteady
      real*4, dimension(*) :: fConcentration
      real*4, dimension(*) :: fLinkConc
      type (SG_VALU), dimension(*) :: lnkValues
      POINTER(PTR_fCurTime, fCurTime)
      POINTER(PTR_fInit, fInit)
      POINTER(PTR_fSteady, fSteady)
      POINTER(PTR_fConcentration, fConcentration)
      POINTER(PTR_fLinkConc, fLinkConc)
      POINTER(PTR_lnkValues, lnkValues)
      integer, parameter :: SG_NDX_FCONCENTRATION = 1
      integer, parameter :: SG_NDX_FSTEADY = 2
      integer, parameter :: SG_NDX_FCURRENT = 3

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xPreEval_Base_Source_Variable = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! use simple (exponential) lag for the first-order step response
      ! the time constant is hard coded with 1 minute (not shown in terms)
      PTR_zValues = simCtrl%pzValues
      PTR_fCurTime = zValues(SG_NDX_CTRL_FCURTIME)%vData
      PTR_zValues = self%pzValues
      PTR_fInit = zValues(SG_NDX_FCONCENTRATION)%vData
      PTR_fSteady = zValues(SG_NDX_FSTEADY)%vData
      PTR_fConcentration = zValues(SG_NDX_FCURRENT)%vData
      fConcentration(1) = (fSteady(1) - fInit(1)) *
     &                    (1.0 - EXP(-fCurTime(1))) + fInit(1)

      ! deposit Source Concentration to the output link - should be only one
      if (iLnkObjs .gt. 0) then
          PTR_lnkObject = pLnkObjs(1)
          PTR_lnkValues = lnkObject%pzValues
          PTR_fLinkConc = lnkValues(SG_NDX_LINK_FCONCENT)%vData
          fLinkConc(1) = fConcentration(1)
      end if

! Collection_Solver.f
! - DLL routines for class <Reference>Collection.Solver
! DATE: Monday, September 10, 2001  TIME: 03:57:13 PM
! The skeleton of this file is generated by SansGUI(tm)

! Attribute indices in class version [1.0.alpha.7]
! 1: rReactor - Reactor Table
! 2: rConstant - Constant Matrix
! 3: rInverse - Inverse Matrix
! 4: rScratch - Scratch Matrix for Temporaries

! Define MIX_WITH_MATLAB in the compilation option to activate MATLAB solution.
! If MIX_WITH_MATLAB is defined, it requires the MATLAB Engine from Mathworks.
! The MATLAB Engine include and library files will be needed to build the DLL.

!DEC$ IF DEFINED (MIX_WITH_MATLAB)
      module Mixer_ML
          integer :: pMatlabEng      ! MATLAB Engine
          integer :: pMatlabC        ! constant matrix
          integer :: pMatlabI        ! inverse matrix
          integer :: pMatlabR        ! RHS vector
          integer :: pMatlabS        ! solution vector
      end module Mixer_ML
!DEC$ END IF

      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"
      integer, dimension(*) :: iNumVars
      integer, dimension(*) :: iSolver
      POINTER(PTR_iNumVars, iNumVars)
      POINTER(PTR_iSolver, iSolver)
!DEC$ IF DEFINED (MIX_WITH_MATLAB)
      integer :: ENGOPEN
      integer :: MXCREATEDOUBLEMATRIX
!DEC$ END IF

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xBgnRun_Collection_Solver = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      PTR_zValues = simCtrl%pzValues
      PTR_iNumVars = zValues(SG_NDX_CTRL_INUMREACT)%vData
      PTR_iSolver = zValues(SG_NDX_CTRL_ISOLVER)%vData

      if (iSolver(1) .eq. MIX_SOLVE_MATLAB) then
!DEC$ IF DEFINED (MIX_WITH_MATLAB)
          ! the following statement may not be thread safe
          ! open connection to the local MATLAB Engine
          pMatlabEng = ENGOPEN('matlab ')
          if (pMatlabEng .eq. 0) then
              cMessage = 'Cannot open MATLAB Engine.'C
              SG_xBgnRun_Collection_Solver = SG_R_STOP
              return
          end if
          ! create all matrices for MATLAB Engine
          pMatlabC = MXCREATEDOUBLEMATRIX(iNumVars(1), iNumVars(1), 0)
          pMatlabR = MXCREATEDOUBLEMATRIX(iNumVars(1), 1, 0)
          pMatlabS = MXCREATEDOUBLEMATRIX(iNumVars(1), 1, 0)
          pMatlabI = MXCREATEDOUBLEMATRIX(iNumVars(1), iNumVars(1), 0)
          if (pMatlabC .eq. 0 .or. pMatlabR .eq. 0 .or.
     &        pMatlabS .eq. 0 .or. pMatlabI .eq. 0 ) then
              call ENGCLOSE(pMatlabEng)
              cMessage = 'Cannot create MATLAB matrices.'C
              SG_xBgnRun_Collection_Solver = SG_R_STOP
              return
          end if
!DEC$ ELSE
          ! call MATLAB option is set, but no MATLAB access code
          cMessage = 'This version does not have MATLAB support.'C
          SG_xBgnRun_Collection_Solver = SG_R_STOP
          return
!DEC$ END IF
      end if

      SG_xBgnRun_Collection_Solver = SG_R_OK
      return
      end

! ======================================================================
! SG_xEval - Evaluation
! ----------------------------------------------------------------------
      integer function SG_xEval_Collection_Solver(self,                 &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xEval_Collection_Solver
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"
      real*8, dimension(*) :: dConstantMatrix
      real*8, dimension(*) :: dInverseMatrix
      real*8, dimension(*) :: dScratchMatrix
      real*8, dimension(*) :: dConstantRHS
      real*8, dimension(*) :: dSolutionVector
      real*8, dimension(*) :: dScratchVector1
      real*8, dimension(*) :: dScratchVector2
      integer, dimension(*) :: iNumReact
      integer, dimension(*) :: iSolver
      logical :: solveLinearEquationsGauss
      logical :: solveLinearEquationsLUD
      logical :: solveLinearEquationsMatlab
      logical :: bResult
      POINTER(PTR_dConstantMatrix, dConstantMatrix)
      POINTER(PTR_dInverseMatrix, dInverseMatrix)
      POINTER(PTR_dScratchMatrix, dScratchMatrix)
      POINTER(PTR_dConstantRHS, dConstantRHS)
      POINTER(PTR_dSolutionVector, dSolutionVector)
      POINTER(PTR_dScratchVector1, dScratchVector1)
      POINTER(PTR_dScratchVector2, dScratchVector2)
      POINTER(PTR_iNumReact, iNumReact)
      POINTER(PTR_iSolver, iSolver)
      integer, parameter :: SG_NDX_RREACTOR = 1
      integer, parameter :: SG_NDX_RCONSTANT = 2
      integer, parameter :: SG_NDX_RINVERSE = 3
      integer, parameter :: SG_NDX_RSCRATCH = 4

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xEval_Collection_Solver = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! Here we start to solve the simultaneous equations because all the
      ! parts' evaluation routines have been called to load the constant
      ! matrix and reactor table (RHS constants)
      PTR_zValues = simCtrl%pzValues
      PTR_iNumReact = zValues(SG_NDX_CTRL_INUMREACT)%vData
      PTR_iSolver = zValues(SG_NDX_CTRL_ISOLVER)%vData

      if (iRefObjs .lt. 4) then
          cMessage = 'Either Constant, Inverse, Scratch Matrix or Reacto
     &r Table is missing.'C
          !important - reset the number of reactors in the simControl object
          iNumReact(1) = 0
          SG_xEval_Collection_Solver = SG_R_STOP
          return
      end if

      ! The matrices are in column major order
      PTR_refObject = pRefObjs(SG_NDX_RCONSTANT)
      PTR_zValues = refObject%pzValues
      PTR_dConstantMatrix = zValues(SG_NDX_MTX_DELEMENT)%vData
      PTR_refObject = pRefObjs(SG_NDX_RINVERSE)
      PTR_zValues = refObject%pzValues
      PTR_dInverseMatrix = zValues(SG_NDX_MTX_DELEMENT)%vData

      PTR_refObject = pRefObjs(SG_NDX_RSCRATCH)
      PTR_zValues = refObject%pzValues
      PTR_dScratchMatrix = zValues(SG_NDX_MTX_DELEMENT)%vData

      PTR_refObject = pRefObjs(SG_NDX_RREACTOR)
      PTR_zValues = refObject%pzValues
      PTR_dConstantRHS = zValues(SG_NDX_TBL_DCONSTANT)%vData
      PTR_dSolutionVector = zValues(SG_NDX_TBL_DSOLUTION)%vData
      PTR_dScratchVector1 = zValues(SG_NDX_TBL_DSCRATCH1)%vData
      PTR_dScratchVector2 = zValues(SG_NDX_TBL_DSCRATCH2)%vData

      ! call the solver routine
      select case (iSolver(1))
      case (MIX_SOLVE_GAUSS)
        bResult = solveLinearEquationsGauss(dConstantMatrix,
     &            dConstantRHS, dSolutionVector, dScratchMatrix,
     &            dScratchVector1, iNumReact )
        if (.not. bResult) then
          cMessage='No solution can be found using Gauss Elimination.'C
          ! important - reset the number of reactors in simControl object
          iNumReact(1) = 0
          SG_xEval_Collection_Solver = SG_R_STOP
          return
        end if
      case (MIX_SOLVE_LUDECOMP)
        bResult = solveLinearEquationsLUD(dConstantMatrix,
     &            dConstantRHS, dSolutionVector, dInverseMatrix,
     &            dScratchMatrix, dScratchVector1, dScratchVector2,
     &            iNumReact )
        if (.not. bResult) then
          cMessage='No solution can be found using LU Decomposition.'C
          ! important - reset the number of reactors in simControl object
          iNumReact(1) = 0
          SG_xEval_Collection_Solver = SG_R_STOP
          return
        end if
      case (MIX_SOLVE_MATLAB)
        bResult = solveLinearEquationsMatlab(dConstantMatrix,
     &            dConstantRHS, dSolutionVector, dInverseMatrix,
     &            iNumReact )
        if (.not. bResult) then
          cMessage='Cannot locate/use the MATLAB Engine to solve it.'C
          ! important - reset the number of reactors in simControl object
          iNumReact(1) = 0
          SG_xEval_Collection_Solver = SG_R_STOP
          return
        end if
      case DEFAULT
        cMessage = 'Unknown solver type.  Check the SimControl object.'C
        ! important - reset the number of reactors in simControl object
        iNumReact(1) = 0
        SG_xEval_Collection_Solver = SG_R_STOP
        return
      end select

      SG_xEval_Collection_Solver = SG_R_OK
      return
      end

! ======================================================================
! SG_xPostEval - Post-Evaluation
! ----------------------------------------------------------------------
      integer function SG_xPostEval_Collection_Solver(self,             &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xPostEval_Collection_Solver
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"
      real*4, dimension(*) :: fCurTime
      real*4, dimension(*) :: fTimeInc
      POINTER(PTR_fCurTime, fCurTime)
      POINTER(PTR_fTimeInc, fTimeInc)

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xPostEval_Collection_Solver = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! Advance system clock - access with care in other reference objects
      PTR_zValues = simCtrl%pzValues
      PTR_fCurTime = zValues(SG_NDX_CTRL_FCURTIME)%vData
      PTR_fTimeInc = zValues(SG_NDX_CTRL_FTIMEINC)%vData
      fCurTime(1) = fCurTime(1) + fTimeInc(1)

      SG_xPostEval_Collection_Solver = SG_R_OK
      return
      end

! ======================================================================
! SG_xEndRun - End Run
! ----------------------------------------------------------------------
      integer function SG_xEndRun_Collection_Solver(self,               &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xEndRun_Collection_Solver
!DEC$ END IF

      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"
      real*4, dimension(*) :: fCurTime
      integer, dimension(*) :: iNumReact
!DEC$ IF DEFINED (MIX_WITH_MATLAB)
      integer, dimension(*) :: iSolver
      POINTER(PTR_iSolver, iSolver)
!DEC$ END IF
      POINTER(PTR_fCurTime, fCurTime)
      POINTER(PTR_iNumReact, iNumReact)

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xEndRun_Collection_Solver = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! Solver's duty to reset the reactor count for the next run
      PTR_zValues = simCtrl%pzValues
      PTR_iNumReact = zValues(SG_NDX_CTRL_INUMREACT)%vData
      iNumReact(1) = 0

      ! Reset the current time to 0 if it is preferred
      ! PTR_fCurTime = zValues(SG_NDX_CTRL_FCURTIME)%vData
      ! fCurTime(1) = 0.

!DEC$ IF DEFINED (MIX_WITH_MATLAB)
      PTR_iSolver = zValues(SG_NDX_CTRL_ISOLVER)%vData
      if (iSolver(1) .eq. MIX_SOLVE_MATLAB) then
          if (pMatlabEng .ne. 0) call matlabCleanUp()
      end if
!DEC$ END IF

      SG_xEndRun_Collection_Solver = SG_R_OK
      return
      end

! ======================================================================
! solveLinearEquationsGauss - solve the linear simultaneous equations
!     using naive Gauss Elimination
! ----------------------------------------------------------------------
! There is no inverse matrix calculated when using this method.
! ----------------------------------------------------------------------
      logical function solveLinearEquationsGauss(dConstantMatrix,
     &    dConstantRHS, dSolutionVector, dScratchMatrix, dScratchVector,
     &    iNumVars )
      implicit none
      integer :: iNumVars
      intent (in) iNumVars
      real*8, dimension(iNumVars, iNumVars) :: dConstantMatrix
      real*8, dimension(iNumVars) :: dConstantRHS
      real*8, dimension(iNumVars) :: dSolutionVector
      real*8, dimension(iNumVars, iNumVars) :: dScratchMatrix
      real*8, dimension(iNumVars) :: dScratchVector
      integer :: i, j, k
      real*8 :: dFactor, dDivisor, dSum

      ! remember the constant matrix and RHS constant vector
      dScratchMatrix = dConstantMatrix
      dScratchVector = dConstantRHS

      do k = 1, iNumVars - 1
          do i = k + 1, iNumVars
              dDivisor = dScratchMatrix(k, k)
              if (dDivisor .eq. 0.D0) then
                  solveLinearEquationsGauss = .false.
                  return
              end if
              dFactor = dScratchMatrix(i, k) / dDivisor

              do j = k + 1, iNumVars
                  dScratchMatrix(i, j) = dScratchMatrix(i, j) -
     &                                   dFactor * dScratchMatrix(k, j)
              end do
              dScratchVector(i) = dScratchVector(i) -
     &                            dFactor * dScratchVector(k)
          end do
      end do

      dDivisor = dScratchMatrix(iNumVars, iNumVars)
      if (dDivisor .eq. 0.D0) then
          solveLinearEquationsGauss = .false.
          return
      end if
      dSolutionVector(iNumVars) = dScratchVector(iNumVars) / dDivisor

      do i = iNumVars - 1, 1, -1
          dSum = 0.D0
          do j = i + 1, iNumVars
              dSum = dSum + dScratchMatrix(i, j) * dSolutionVector(j)
          end do
          dDivisor = dScratchMatrix(i, i)
          if (dDivisor .eq. 0.D0) then
              solveLinearEquationsGauss = .false.
              return
          end if
          dSolutionVector(i) = (dScratchVector(i) - dSum) / dDivisor
      end do

      solveLinearEquationsGauss = .true.
      return
      end

! ======================================================================
! solveLinearEquationsLUD - solve the linear simultaneous equations
!     using LU Decomposition
! ----------------------------------------------------------------------
      logical function solveLinearEquationsLUD(dConstantMatrix,
     &    dConstantRHS, dSolutionVector, dInverseMatrix, dScratchMatrix,
     &    dScratchVector1, dScratchVector2, iNumVars )
      implicit none
      integer :: iNumVars
      intent (in) iNumVars
      real*8, dimension(iNumVars, iNumVars) :: dConstantMatrix
      real*8, dimension(iNumVars) :: dConstantRHS
      real*8, dimension(iNumVars) :: dSolutionVector
      real*8, dimension(iNumVars, iNumVars) :: dInverseMatrix
      real*8, dimension(iNumVars, iNumVars) :: dScratchMatrix
      real*8, dimension(iNumVars) :: dScratchVector1
      real*8, dimension(iNumVars) :: dScratchVector2
      logical :: bResult
      logical :: decompose, substitute, inverse

      ! remember the constant matrix
      dScratchMatrix = dConstantMatrix
      bResult = decompose(dScratchMatrix, iNumVars)
      if (.not. bResult) then
          solveLinearEquationsLUD = .false.
          return
      end if

      ! remember the LU matrix in dScratchMatrix for inverse matrix calculation
      dInverseMatrix = dScratchMatrix

      ! copy the RHS vector for inverse matrix calculation
      dScratchVector1 = dConstantRHS
      bResult = substitute(dScratchMatrix, dScratchVector1,
     &                     dSolutionVector, iNumVars )
      if (.not. bResult) then
          solveLinearEquationsLUD = .false.
          return
      end if

      ! inverse matrix should contain the LU decomposed result
      dScratchMatrix = dInverseMatrix

      solveLinearEquationsLUD = inverse(dScratchMatrix, dInverseMatrix,
     &    dConstantRHS, dScratchVector1, dScratchVector2, iNumVars )
      return
      end

! ----------------------------------------------------------------------
      logical function decompose(dConstantMatrix, iNumVars)
      implicit none
      integer :: iNumVars
      intent (in) iNumVars
      real*8, dimension(iNumVars, iNumVars) :: dConstantMatrix
      integer :: i, j, k
      real*8 :: dFactor, dDivisor

      do k = 1, iNumVars - 1
          do i = k + 1, iNumVars
              dDivisor = dConstantMatrix(k, k)
              if (dDivisor .eq. 0.D0) then
                  decompose = .false.
                  return
              end if
              dFactor = dConstantMatrix(i, k) / dDivisor
              dConstantMatrix(i, k) = dFactor
              do j = k + 1, iNumVars
                  dConstantMatrix(i, j) = dConstantMatrix(i, j) -
     &                     dFactor * dConstantMatrix(k, j)
              end do
          end do
      end do
      decompose = .true.
      return
      end

! ----------------------------------------------------------------------
      logical function substitute(dMatrixLU, dVectorRHS,
     &    dVectorSolution, iNumVars )
      implicit none
      integer :: iNumVars
      intent (in) iNumVars
      real*8, dimension(iNumVars, iNumVars) :: dMatrixLU
      real*8, dimension(iNumVars) :: dVectorRHS
      real*8, dimension(iNumVars) :: dVectorSolution
      integer :: i, j
      real*8 :: dSum, dDivisor

      ! forward substitution
      do i = 2, iNumVars
          dSum = dVectorRHS(i)
          do j = 1, i - 1
              dSum = dSum - dMatrixLU(i, j) * dVectorRHS(j)
          end do
          dVectorRHS(i) = dSum
      end do

      ! backward substitution
      dDivisor = dMatrixLU(iNumVars, iNumVars)
      if (dDivisor .eq. 0.D0) then
          substitute = .false.
          return
      end if
      dVectorSolution(iNumVars) = dVectorRHS(iNumVars) / dDivisor
      do i = iNumVars - 1, 1, -1
          dSum = 0.D0
          do j = i + 1, iNumVars
              dSum = dSum + dMatrixLU(i, j) * dVectorSolution(j)
          end do
          dDivisor = dMatrixLU(i, i)
          if (dDivisor .eq. 0.D0) then
              substitute = .false.
              return
          end if
          dVectorSolution(i) = (dVectorRHS(i) - dSum) / dDivisor
      end do
      substitute = .true.
      return
      end

! ----------------------------------------------------------------------
      logical function inverse(dMatrixLU, dMatrixInverse, dVectorRHS,
     &    dVectorScratch1, dVectorScratch2, iNumVars )
      implicit none
      integer :: iNumVars
      intent (in) iNumVars
      real*8, dimension(iNumVars, iNumVars) :: dMatrixLU
      real*8, dimension(iNumVars, iNumVars) :: dMatrixInverse
      real*8, dimension(iNumVars) :: dVectorRHS
      real*8, dimension(iNumVars) :: dVectorScratch1
      real*8, dimension(iNumVars) :: dVectorScratch2
      integer :: i, j
      logical :: bResult
      logical :: substitute

      ! calling decompose() is not necessary because the input is LU result
      dVectorScratch1 = dVectorRHS

      do i = 1, iNumVars
          do j = 1, iNumVars
              if (i .eq. j) then
                  dVectorScratch1(j) = 1.D0
              else
                  dVectorScratch1(j) = 0.D0
              end if
          end do
          ! use dVectorScratch2 to receive solution vector X
          bResult = substitute(dMatrixLU, dVectorScratch1,
     &                         dVectorScratch2, iNumVars )
          if (.not. bResult) then
              inverse = .false.
              return
          end if
          do j = 1, iNumVars
              dMatrixInverse(j, i) = dVectorScratch2(j)
          end do
      end do
      inverse = .true.
      return
      end

! ======================================================================
! solveLinearEquationsMatlab - solve the linear simultaneous equations
!     using the MATLAB Engine
! ----------------------------------------------------------------------
      logical function solveLinearEquationsMatlab(dConstantMatrix,
     &    dConstantRHS, dSolutionVector, dInverseMatrix, iNumVars )
!DEC$ IF DEFINED (MIX_WITH_MATLAB)
      use Mixer_ML
!DEC$ END IF
      implicit none
      integer :: iNumVars
      intent (in) iNumVars
      integer :: iSizeM
      real*8, dimension(iNumVars, iNumVars) :: dConstantMatrix
      real*8, dimension(iNumVars) :: dConstantRHS
      real*8, dimension(iNumVars) :: dSolutionVector
      real*8, dimension(iNumVars, iNumVars) :: dInverseMatrix

      integer :: MXGETPR
      integer :: ENGGETVARIABLE

!DEC$ IF DEFINED (MIX_WITH_MATLAB)
      iSizeM = iNumVars * iNumVars
      dSolutionVector = 0.D0
      dInverseMatrix = 0.D0
      ! prepare input matrices
      call MXCOPYREAL8TOPTR(dConstantMatrix, MXGETPR(pMatlabC), iSizeM)
      call MXCOPYREAL8TOPTR(dConstantRHS, MXGETPR(pMatlabR), iNumVars)

      call ENGPUTVARIABLE(pMatlabEng, 'C', pMatlabC)
      call ENGPUTVARIABLE(pMatlabEng, 'R', pMatlabR)
      call ENGPUTVARIALBE(pMatlabEng, 'S', pMatlabS)
      call ENGPUTVARIABLE(pMatlabEng, 'I', pMatlabI)

      ! now execute MATLAB commands - both must be executed
      call ENGEVALSTRING(pMatlabEng, 'S = C \ R;')
      call ENGEVALSTRING(pMatlabEng, 'I = inv(C);')

      ! now fetch the solutions
      pMatlabS = ENGGETVARIABLE(pMatlabEng, 'S')
      if (pMatlabS .eq. 0) then
          call matlabCleanUp()
          solveLinearEquationsMatlab = .false.
          return
      end if
      call MXCOPYPTRTOREAL8(MXGETPR(pMatlabS),dSolutionVector,iNumVars)

      pMatlabI = ENGGETVARIABLE(pMatlabEng, 'I')
      if (pMatlabI .eq. 0) then
          call matlabCleanUp()
          solveLinearEquationsMatlab = .false.
          return
      end if
      call MXCOPYPTRTOREAL8(MXGETPR(pMatlabI),dInverseMatrix,iSizeM)

      solveLinearEquationsMatlab = .true.
!DEC$ ELSE
      ! call MATLAB option is set, but no MATLAB access code
      solveLinearEquationsMatlab = .false.
!DEC$ END IF
      return
      end

! ======================================================================
! matlabCleanUp - subroutine to clean up MATLAB access
! ----------------------------------------------------------------------
!DEC$ IF DEFINED (MIX_WITH_MATLAB)
      subroutine matlabCleanUp()
      use Mixer_ML
      implicit none
      call MXDESTROYARRAY(pMatlabC)
      call MXDESTROYARRAY(pMatlabR)
      call MXDESTROYARRAY(pMatlabS)
      call MXDESTROYARRAY(pMatlabI)
      call ENGCLOSE(pMatlabEng)
      return
      end
!DEC$ END IF

! Matrix_Calculation.f
! - DLL routines for class <Reference>Matrix.Calculation
! DATE: Monday, September 10, 2001  TIME: 03:57:13 PM
! The skeleton of this file is generated by SansGUI(tm)

! Attribute indices in class version [1.0.alpha.5]
! 1: iCols - Number of Columns
! 2: iRows - Number of Rows
! 3: iSheets - Number of Sheets
! 4: dElement - Elements in Matrix

! ======================================================================
! SG_xInitSize - Resize for Init
! ----------------------------------------------------------------------
      integer function SG_xInitSize_Matrix_Calculation(self,            &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xInitSize_Matrix_Calculation
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"
      integer, dimension(*) :: iNumReact
      integer, dimension(*) :: iCols
      integer, dimension(*) :: iRows
      integer, dimension(*) :: iSheets
      POINTER(PTR_iNumReact, iNumReact)
      POINTER(PTR_iCols, iCols)
      POINTER(PTR_iRows, iRows)
      POINTER(PTR_iSheets, iSheets)
      integer, parameter :: SG_NDX_ICOLS = 1
      integer, parameter :: SG_NDX_IROWS = 2
      integer, parameter :: SG_NDX_ISHEETS = 3

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xInitSize_Matrix_Calculation = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! resize constant matrix according to the number of reactors registered
      PTR_zValues = simCtrl%pzValues
      PTR_iNumReact = zValues(SG_NDX_CTRL_INUMREACT)%vData
      if (iNumReact(1) .lt. 1) then
          cMessage = 'At least one reactor part is required.'C
          iNumReact(1) = 0
          SG_xInitSize_Matrix_Calculation = SG_R_STOP
          return
      end if

      PTR_zValues = self%pzValues
      PTR_iCols = zValues(SG_NDX_ICOLS)%vData
      PTR_iRows = zValues(SG_NDX_IROWS)%vData
      PTR_iSheets = zValues(SG_NDX_ISHEETS)%vData

      iCols(1) = iNumReact(1)
      iRows(1) = iNumReact(1)
      iSheets(1) = 1          ! always

      SG_xInitSize_Matrix_Calculation = SG_R_OK
      return
      end

! ======================================================================
! SG_xPreEval - Pre-Evaluation
! ----------------------------------------------------------------------
      integer function SG_xPreEval_Matrix_Calculation(self,             &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xPreEval_Matrix_Calculation
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      integer :: iSize
      real*8, dimension(*) :: dElement
      POINTER(PTR_dElement, dElement)
      integer, parameter :: SG_NDX_DELEMENT = 4

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xPreEval_Matrix_Calculation = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! initialize the matrix to 0 to prepare for loading constants
      PTR_zValues = self%pzValues
      PTR_dElement = zValues(SG_NDX_DELEMENT)%vData
      iSize = zValues(SG_NDX_DELEMENT)%iSize
      call MIX_ZERO_OUT_DARRAY(dElement, iSize)

      SG_xPreEval_Matrix_Calculation = SG_R_OK
      return
      end

! ======================================================================
! loadMatrixConstant - load a simultaneous equation constant into the
!                      constant matrix
! ----------------------------------------------------------------------
! ARGUMENT:
!     self       - the Matrix.Constant object to load constants to
!     ciPartNdx  - the registered 1-based index of the reactor (row)
!     ciLoadNdx  - the 1-based index of the reactor associated (column)
!     cfValue    - the constant value to be loaded
!     cbInput    - TRUE, if input to the part; FALSE, if output
!
! RETURN:
!     1-based index of the reactor (row); 0 when failed.
!
! NOTES:
!     The convention for loading equation constants (LHS) is that
!     the input values are negative and output values are positive.
!     The matrix is stored in column major order.
! ----------------------------------------------------------------------
      integer function loadMatrixConstant(self, ciPartNdx, ciLoadNdx,
     &                                    cfValue, cbInput )
      include "../Mixer_1_0F/Mixer_T.h"
      type (SG_OBJ) :: self
      integer :: ciPartNdx, ciLoadNdx
      real*4 :: cfValue
      logical :: cbInput
      intent (in) ciPartNdx, ciLoadNdx, cfValue, cbInput

      integer :: iNdx                ! to store linearized index -- column major
      integer, dimension(*) :: iSize ! size of rows or columns (square matrix)
      real*8, dimension(*) :: dElement
      type (SG_VALU), dimension(*) :: zValues
      POINTER(PTR_zValues, zValues)
      POINTER(PTR_iSize, iSize)
      POINTER(PTR_dElement, dElement)
      integer, parameter :: SG_NDX_IROWS = 2
      integer, parameter :: SG_NDX_DELEMENT = 4

      ! should be a square matrix
      PTR_zValues = self%pzValues
      PTR_iSize = zValues(SG_NDX_IROWS)%vData

      if (ciPartNdx .lt. 1 .or. ciLoadNdx .lt. 1 .or.
     &    ciPartNdx .gt. iSize(1) .or. ciLoadNdx .gt. iSize(1) ) then
          ! out of bound error
          loadMatrixConstant = 0
          return
      end if

      ! calculate linear index
      iNdx = iSize(1) * (ciLoadNdx - 1) + ciPartNdx
      PTR_dElement = zValues(SG_NDX_DELEMENT)%vData
      if (cbInput) then
          dElement(iNdx) = dElement(iNdx) - DBLE(cfValue)
      else
          dElement(iNdx) = dElement(iNdx) + DBLE(cfValue)
      end if
      

      loadMatrixConstant = ciPartNdx
      return
      end

! Table_Reactor.f
! - DLL routines for class <Reference>Table.Reactor
! DATE: Monday, September 10, 2001  TIME: 03:57:13 PM
! The skeleton of this file is generated by SansGUI(tm)

! Attribute indices in class version [1.0.alpha.9]
! 1: iSize - Table Size
! 2: iSheets - Number of Sheets
! 3: iPartSN - Part Internal Serial Number
! 4: dConstant - Constants in Mass Balance Eqns
! 5: dSolution - Solutions in Mass Balance Eqns
! 6: dScratch1 - Vector Buffer 1 for Temporaries
! 7: dScratch2 - Vector Buffer 2 for Temporaries

! ======================================================================
! SG_xInitSize - Resize for Init
! ----------------------------------------------------------------------
      integer function SG_xInitSize_Table_Reactor(self,                 &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xInitSize_Table_Reactor
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      include "../Mixer_1_0F/Mixer_1_0F.h"
      integer, dimension(*) :: iNumReact
      integer, dimension(*) :: iSize
      integer, dimension(*) :: iSheets
      POINTER(PTR_iNumReact, iNumReact)
      POINTER(PTR_iSize, iSize)
      POINTER(PTR_iSheets, iSheets)
      integer, parameter :: SG_NDX_ISIZE = 1
      integer, parameter :: SG_NDX_ISHEETS = 2

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xInitSize_Table_Reactor = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! resize reactor table according to the number of reactors
      PTR_zValues = simCtrl%pzValues
      PTR_iNumReact = zValues(SG_NDX_CTRL_INUMREACT)%vData
      if (iNumReact(1) .lt. 1) then
          cMessage = 'At least one reactor part is required.'C
          iNumReact(1) = 0
          SG_xInitSize_Table_Reactor = SG_R_STOP
          return
      end if

      PTR_zValues = self%pzValues
      PTR_iSize = zValues(SG_NDX_ISIZE)%vData
      PTR_iSheets = zValues(SG_NDX_ISHEETS)%vData
      iSize(1) = iNumReact(1)
      iSheets(1) = 1             ! always

      SG_xInitSize_Table_Reactor = SG_R_OK
      return
      end

! ======================================================================
! SG_xInit - Initialization
! ----------------------------------------------------------------------
      integer function SG_xInit_Table_Reactor(self,                     &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xInit_Table_Reactor
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      integer :: iSize
      integer, dimension(*) :: iPartSN
      POINTER(PTR_iPartSN, iPartSN)
      integer, parameter :: SG_NDX_IPARTSN = 3

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xInit_Table_Reactor = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      ! zero out the reactor parts registration
      ! all the other vectors need to be initialized in each Pre-Eval cycle
      PTR_zValues = self%pzValues
      PTR_iPartSN = zValues(SG_NDX_IPARTSN)%vData
      iSize = zValues(SG_NDX_IPARTSN)%iSize
      call MIX_ZERO_OUT_IARRAY(iPartSN, iSize)

      SG_xInit_Table_Reactor = SG_R_OK
      return
      end

! ======================================================================
! SG_xPreEval - Pre-Evaluation
! ----------------------------------------------------------------------
      integer function SG_xPreEval_Table_Reactor(self,                  &
     &                        simCtrl, chgChild,                        &
     &                        pRefObjs, iRefObjs,                       &
     &                        pAdjObjs, iAdjObjs,                       &
     &                        pLnkObjs, iLnkObjs,                       &
     &                        cMessage, cCommand, pOutFile )
!DEC$ IF DEFINED (_DLL)
!DEC$ ATTRIBUTES DLLEXPORT :: SG_xPreEval_Table_Reactor
!DEC$ END IF
      include "SGdllf.h"

      ! TODO: declare your local variables here

      integer :: iSize
      real*8, dimension(*) :: dConstant
      real*8, dimension(*) :: dSolution
      real*8, dimension(*) :: dScratch1
      real*8, dimension(*) :: dScratch2
      POINTER(PTR_dConstant, dConstant)
      POINTER(PTR_dSolution, dSolution)
      POINTER(PTR_dScratch1, dScratch1)
      POINTER(PTR_dScratch2, dScratch2)
      integer, parameter :: SG_NDX_DCONSTANT = 4
      integer, parameter :: SG_NDX_DSOLUTION = 5
      integer, parameter :: SG_NDX_DSCRATCH1 = 6
      integer, parameter :: SG_NDX_DSCRATCH2 = 7

      if (self%nSGobjSchema .ne. SG_OBJ_SCHEMA) then
          SG_xPreEval_Table_Reactor = SG_R_SCHM
          return
      end if

      ! TODO: put your simulator code here

      PTR_zValues = self%pzValues
      iSize = zValues(SG_NDX_DCONSTANT)%iSize
      PTR_dConstant = zValues(SG_NDX_DCONSTANT)%vData
      PTR_dSolution = zValues(SG_NDX_DSOLUTION)%vData
      PTR_dScratch1 = zValues(SG_NDX_DSCRATCH1)%vData
      PTR_dScratch2 = zValues(SG_NDX_DSCRATCH2)%vData

      call MIX_ZERO_OUT_DARRAY(dConstant, iSize)
      call MIX_ZERO_OUT_DARRAY(dSolution, iSize)
      call MIX_ZERO_OUT_DARRAY(dScratch1, iSize)
      call MIX_ZERO_OUT_DARRAY(dScratch2, iSize)

      SG_xPreEval_Table_Reactor = SG_R_OK
      return
      end

! ======================================================================
! registerReactor - register reactor to the table
! ----------------------------------------------------------------------
! ARGUMENTS:
!     self     - the Table.Reactor object to provide registration
!     ciPartSN - the unique serial number of the reactor
!
! RETURN VALUE:
!     1-based index of the reactor; 0 when failed
! ----------------------------------------------------------------------
      integer function registerReactor(self, ciPartSN)
      include "../Mixer_1_0F/Mixer_T.h"
      type (SG_OBJ) :: self
      integer :: ciPartSN
      intent (in) ciPartSN

      integer :: i
      integer :: iSize
      integer, dimension(*) :: iPartSN
      type (SG_VALU), dimension(*) :: zValues
      POINTER(PTR_iPartSN, iPartSN)
      POINTER(PTR_zValues, zValues)
      integer, parameter :: SG_NDX_IPARTSN = 3

      PTR_zValues = self%pzValues
      PTR_iPartSN = zValues(SG_NDX_IPARTSN)%vData
      iSize = zValues(SG_NDX_IPARTSN)%iSize

      do i = 1, iSize
          if (iPartSN(i) .lt. 1) then
              !it is not occupied, register it in the table
              iPartSN(i) = ciPartSN
              registerReactor = i
              return
          else if (iPartSN(i) .eq. ciPartSN) then
              ! it has been registered before
              registerReactor = i
              return
          end if
      end do

      registerReactor = 0   ! failed
      return
      end

! ======================================================================
! loadTableConstant - load a constant into the RHS vector
! ----------------------------------------------------------------------
! ARGUMENTS:
!     self      - the Table.Reactor object to load RHS constants to
!     ciPartNdx - the registered 1-based index of the reactor
!     cfValue   - the constant value to be loaded
!     cbInput   - .TRUE., if input to the part; .FALSE., if output
!
! RETURN VALUE:
!     1-based index of the reactor; 0 when failed
!
! NOTES:
!     The convention for loading RHS constants is that the input values
!     are positive and output values are negative.
! ----------------------------------------------------------------------
      integer function loadTableConstant(self, ciPartNdx,
     &                                   cfValue, cbInput )
      include "../Mixer_1_0F/Mixer_T.h"
      type (SG_OBJ) :: self
            integer :: ciPartNdx
      real*4 :: cfValue
      logical :: cbInput
      intent (in) ciPartNdx, cfValue, cbInput

      integer :: iSize
      type (SG_VALU), dimension(*) :: zValues
      real*8, dimension(*) :: dConstant
      POINTER(PTR_zValues, zValues)
      POINTER(PTR_dConstant, dConstant)
      integer, parameter :: SG_NDX_DCONSTANT = 4

      PTR_zValues = self%pzValues
      PTR_dConstant = zValues(SG_NDX_DCONSTANT)%vData
      iSize = zValues(SG_NDX_DCONSTANT)%iSize

      if (ciPartNdx .lt. 1 .or. ciPartNdx .gt. iSize) then
          loadTableConstant = 0    ! out of bound error
          return
      end if

      ! See NOTES above for the sign of the value
      if (cbInput) then
          dConstant(ciPartNdx) = dConstant(ciPartNdx) + DBLE(cfValue)
      else
          dConstant(ciPartNdx) = dConstant(ciPartNdx) - DBLE(cfValue)
      end if

      loadTableConstant = ciPartNdx
      return
      end

! This header file contains just the SG_OBJ and SG_VALU data structures
! to be included by the routines that need to access the structure members.
! ======================================================================
      implicit none

! data structure declaration
      type SG_VALU
          sequence
          integer :: nType
          integer :: iSize
          integer :: iCols
          integer :: iRows
          integer(4) :: vData
      end type SG_VALU

      type SG_OBJ
          sequence
          integer :: nSGobjSchema
          integer :: iStatus
          integer :: iUserData
          integer :: iNumVars
          integer(4) :: pzValues
          integer :: iVerMajor
          integer :: iVerMinor
          integer :: iVerPatch
          integer :: iVerBuild
          integer :: nCmpnNo
          integer(4) :: pcObjName
          integer(4) :: pcCmpnName
          integer(4) :: pcClassPath
          integer(4) :: pcCmpnPath
          integer(4) :: psVarNames
      end type SG_OBJ

! This implementation file contains common routines for component classes.
! ======================================================================
! MIX_SET_USER_TYPE - set user type for run-time class identification
! ----------------------------------------------------------------------
      integer function MIX_SET_USER_TYPE(self, iType)
      include "Mixer_T.h"
      type (SG_OBJ) :: self
      integer :: iType

      self%iUserData = iType
      MIX_SET_USER_TYPE = iType
      return
      end

! ======================================================================
! MIX_GET_USER_TYPE - get user type for run-time class identification
! ----------------------------------------------------------------------
      integer function MIX_GET_USER_TYPE(self)
      include "Mixer_T.h"
      type (SG_OBJ) :: self

      MIX_GET_USER_TYPE = self%iUserData
      return
      end

! ======================================================================
! MIX_ZERO_OUT_IARRAY - zero out an integer array
! ----------------------------------------------------------------------
      subroutine MIX_ZERO_OUT_IARRAY(iArray, iSize)
      implicit none
      integer :: iSize
      intent (in) iSize
      integer, dimension(iSize) :: iArray

      iArray = 0
      return
      end

! ======================================================================
! MIX_ZERO_OUT_DARRAY - zero out a double precision array
! ----------------------------------------------------------------------
      subroutine MIX_ZERO_OUT_DARRAY(dArray, iSize)
      implicit none
      integer :: iSize
      intent (in) iSize
      real*8, dimension(iSize) :: dArray

      dArray = 0.D0
      return
      end

! Mixer_1_0F.h - manually created header file for macros and
!                function prototype, called from various classes

      ! from class Table.Reactor
      integer, parameter :: SG_NDX_TBL_IPARTSN = 3
      integer, parameter :: SG_NDX_TBL_DCONSTANT = 4
      integer, parameter :: SG_NDX_TBL_DSOLUTION = 5
      integer, parameter :: SG_NDX_TBL_DSCRATCH1 = 6
      integer, parameter :: SG_NDX_TBL_DSCRATCH2 = 7

      ! from class Matrix.Calculation
      integer, parameter :: SG_NDX_MTX_DELEMENT = 4

      ! from class simControl.Mixer
      integer, parameter :: SG_NDX_CTRL_ISOLVER = 10
      integer, parameter :: SG_NDX_CTRL_FTIMEINC = 11
      integer, parameter :: SG_NDX_CTRL_FCURTIME = 12
      integer, parameter :: SG_NDX_CTRL_INUMREACT = 13

      ! from class Link.Pipe
      integer, parameter :: SG_NDX_LINK_ILINKINFO = 1
      integer, parameter :: SG_NDX_LINK_FFLOWRATE = 2
      integer, parameter :: SG_NDX_LINK_FCONCENT = 3

      ! for class type info in iUserData
      integer, parameter :: MIX_USER_TYPE_UNKNOWN = 0
      integer, parameter :: MIX_USER_TYPE_REACTOR = 1
      integer, parameter :: MIX_USER_TYPE_SOURCE = 2
      integer, parameter :: MIX_USER_TYPE_SINK = 3

      ! solution methods
      integer, parameter :: MIX_SOLVE_GAUSS = 0
      integer, parameter :: MIX_SOLVE_LUDECOMP = 1
      integer, parameter :: MIX_SOLVE_MATLAB = 2

      ! external functions
      integer MIX_SET_USER_TYPE
      integer MIX_GET_USER_TYPE
      integer registerReactor
      integer loadTableConstant
      integer loadMatrixConstant