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 = ENGG