pywafo/wafo/source/rind2007/rindmod.f

! Programs available in module RINDMOD :
!
!   1) setConstants
!   2) RINDD
!
! SETCONSTANTS set member variables controlling the performance of RINDD
!
! CALL setConstants(method,xcscale,abseps,releps,coveps,maxpts,minpts,nit,xcutoff,Nc1c2)
!
! METHOD  = INTEGER defining the SCIS integration method
!             1 Integrate by SADAPT for Ndim<9 and by KRBVRC otherwise
!             2 Integrate by SADAPT for Ndim<20 and by KRBVRC otherwise
!             3 Integrate by KRBVRC by Genz (1993) (Fast Ndim<101) (default)
!             4 Integrate by KROBOV by Genz (1992) (Fast Ndim<101)
!             5 Integrate by RCRUDE by Genz (1992) (Slow Ndim<1001)
!             6 Integrate by SOBNIED               (Fast Ndim<1041)
!             7 Integrate by DKBVRC by Genz (2003) (Fast Ndim<1001)
!
!   XCSCALE = REAL to scale the conditinal probability density, i.e.,
!             f_{Xc} = exp(-0.5*Xc*inv(Sxc)*Xc + XcScale) (default XcScale =0)
!   ABSEPS  = REAL absolute error tolerance.       (default 0)
!   RELEPS  = REAL relative error tolerance.       (default 1e-3)
!   COVEPS  = REAL error tolerance in Cholesky factorization (default 1e-13)
!   MAXPTS  = INTEGER, maximum number of function values allowed. This
!             parameter can be used to limit the time. A sensible
!             strategy is to start with MAXPTS = 1000*N, and then
!             increase MAXPTS if ERROR is too large.
!             (Only for METHOD~=0) (default 40000)
!   MINPTS  = INTEGER, minimum number of function values allowed.
!             (Only for METHOD~=0) (default 0)
!   NIT     = INTEGER, maximum number of Xt variables to integrate
!             This parameter can be used to limit the time.
!             If NIT is less than the rank of the covariance matrix,
!             the returned result is a upper bound for the true value
!             of the integral.  (default 1000)
!   XCUTOFF = REAL cut off value where the marginal normal
!             distribution is truncated. (Depends on requested
!             accuracy. A value between 4 and 5 is reasonable.)
!  NC1C2    = number of times to use the regression equation to restrict
!            integration area. Nc1c2 = 1,2 is recommended. (default 2)
!
!
!RIND computes  E[Jacobian*Indicator|Condition]*f_{Xc}(xc(:,ix))
!
! where
!     "Indicator" = I{ H_lo(i) < X(i) < H_up(i), i=1:Nt+Nd }
!     "Jacobian"  = J(X(Nt+1),...,X(Nt+Nd+Nc)), special case is
!     "Jacobian"  = |X(Nt+1)*...*X(Nt+Nd)|=|Xd(1)*Xd(2)..Xd(Nd)|
!     "condition" = Xc=xc(:,ix),  ix=1,...,Nx.
!     X = [Xt; Xd ;Xc], a stochastic vector of Multivariate Gaussian
!         variables where Xt,Xd and Xc have the length Nt, Nd and Nc,
!         respectively.
!         (Recommended limitations Nx, Nt<101, Nd<7 and NIT,Nc<11)
! (RIND = Random Integration N Dimensions)
!
!CALL RINDD(E,err,terr,S,m,xc,Nt,indI,Blo,Bup,INFIN);
!
!        E = expectation/density as explained above size 1 x Nx (out)
!      ERR = estimated sampling error  size 1 x Nx (out)
!     TERR = estimated truncation error size 1 x Nx (out)
!        S = Covariance matrix of X=[Xt;Xd;Xc] size N x N (N=Nt+Nd+Nc) (in)
!        m = the expectation of X=[Xt;Xd;Xc]   size N x 1              (in)
!       xc = values to condition on            size Nc x Nx            (in)
!     indI = vector of indices to the different barriers in the        (in)
!            indicator function,  length NI, where   NI = Nb+1
!            (NB! restriction  indI(1)=0, indI(NI)=Nt+Nd )
!  Blo,Bup = Lower and upper barrier coefficients used to compute the   (in)
!            integration limits A and B, respectively.
!            size  Mb x Nb. If  Mb<Nc+1 then
!            Blo(Mb+1:Nc+1,:) is assumed to be zero.
!    INFIN = INTEGER, array of integration limits flags:  size 1 x Nb   (in)
!            if INFIN(I) < 0, Ith limits are (-infinity, infinity);
!            if INFIN(I) = 0, Ith limits are (-infinity, Hup(I)];
!            if INFIN(I) = 1, Ith limits are [Hlo(I), infinity);
!            if INFIN(I) = 2, Ith limits are [Hlo(I), Hup(I)].
!
! The relation to the integration limits Hlo and Hup are as follows
!    IF INFIN(j)>=0,
!      IF INFIN(j)~=0,  A(i)=Blo(1,j)+Blo(2:Mb,j).'*xc(1:Mb-1,ix),
!      IF INFIN(j)~=1,  B(i)=Bup(1,j)+Bup(2:Mb,j).'*xc(1:Mb-1,ix),
!
!            where i=indI(j-1)+1:indI(j), j=1:NI-1, ix=1:Nx
!            Thus the integration limits may change with the conditional
!            variables.
!Example:
! The indices, indI=[0 3 5 6], and coefficients Blo=[0 0 -1],
! Bup=[0 0 5], INFIN=[0 1 2]
! means that   A = [-inf -inf -inf 0 0 -1]  B = [0 0 0 inf inf 5]
!
!
! (Recommended limitations Nx,Nt<101, Nd<7 and Nc<11)
! Also note that the size information have to be transferred to RINDD
! through the input arguments E,S,m,Nt,IndI,Blo,Bup and INFIN
!
! For further description see the modules
!
! References
! Podgorski et al. (2000)
! "Exact distributions for apparent waves in irregular seas"
! Ocean Engineering,  Vol 27, no 1, pp979-1016.                        (RINDD)
!
! R. Ambartzumian, A. Der Kiureghian, V. Ohanian and H.
! Sukiasian (1998)
! "Multinormal probabilities by sequential conditioned
!  importance sampling: theory and application"                        (MVNFUN)
! Probabilistic Engineering Mechanics, Vol. 13, No 4. pp 299-308
!
! Alan Genz (1992)
! 'Numerical Computation of Multivariate Normal Probabilites'          (MVNFUN)
! J. computational Graphical Statistics, Vol.1, pp 141--149
!
! Alan Genz and Koon-Shing Kwong (2000?)
! 'Numerical Evaluation of Singular Multivariate Normal Distributions' (MVNFUN,COVSRT)
! Computational Statistics and Data analysis
!
!
! P. A. Brodtkorb (2004),                                 (RINDD, MVNFUN, COVSRT)
! Numerical evaluation of multinormal expectations
! In Lund university report series
! and in the Dr.Ing thesis:
! The probability of Occurrence of dangerous Wave Situations at Sea.
! Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU,
! Trondheim, Norway.

! Tested on:  DIGITAL UNIX Fortran90 compiler
!             PC pentium II with Lahey Fortran90 compiler
!             Solaris with SunSoft F90 compiler Version 1.0.1.0  (21229283)
! History:
! Revised pab aug. 2009
! -renamed from rind2007 to rindmod
! Revised pab July 2007
! - separated the absolute error into ERR and TERR.
! - renamed from alanpab24 -> rind2007
! revised pab 23may2004
! RIND module totally rewritten according to the last reference.


      MODULE GLOBALCONST        ! global constants
      IMPLICIT NONE
      DOUBLE PRECISION, PARAMETER :: gSQTWPI1= 0.39894228040143D0  !=1/sqrt(2*pi)
      DOUBLE PRECISION, PARAMETER :: gSQPI1  = 0.56418958354776D0  !=1/sqrt(pi)
      DOUBLE PRECISION, PARAMETER :: gSQPI   = 1.77245385090552D0  !=sqrt(pi)
      DOUBLE PRECISION, PARAMETER :: gSQTW   = 1.41421356237310D0  !=sqrt(2)
      DOUBLE PRECISION, PARAMETER :: gSQTW1  = 0.70710678118655D0  !=1/sqrt(2)
      DOUBLE PRECISION, PARAMETER :: gPI1    = 0.31830988618379D0  !=1/pi
      DOUBLE PRECISION, PARAMETER :: gPI     = 3.14159265358979D0  !=pi
      DOUBLE PRECISION, PARAMETER :: gTWPI   = 6.28318530717958D0  !=2*pi
      DOUBLE PRECISION, PARAMETER :: gSQTWPI = 2.50662827463100D0  !=sqrt(2*pi)
      DOUBLE PRECISION, PARAMETER :: gONE    = 1.D0
      DOUBLE PRECISION, PARAMETER :: gTWO    = 2.D0
      DOUBLE PRECISION, PARAMETER :: gHALF   = 0.5D0
      DOUBLE PRECISION, PARAMETER :: gZERO   = 0.D0
      DOUBLE PRECISION, PARAMETER :: gINFINITY = 37.D0 ! SQRT(-gTWO*LOG(1.D+12*TINY(gONE)))
!     Set gINFINITY (infinity).
!     Such that EXP(-2.x^2) > 10^(12) times TINY
!     SAVE gINFINITY
      END MODULE GLOBALCONST

      MODULE RINDMOD
      USE GLOBALCONST
!      USE PRINTMOD  ! used for debugging only
      IMPLICIT NONE
      PRIVATE
      PUBLIC :: RINDD, SetConstants
	PUBLIC :: mCovEps, mAbsEps,mRelEps, mXcutOff, mXcScale
      PUBLIC :: mNc1c2, mNIT, mMaxPts,mMinPts, mMethod, mSmall
      private :: preInit
      private :: initIntegrand
      private :: initfun,mvnfun,cvsrtxc,covsrt1,covsrt,rcscale,rcswap
      private :: cleanUp

	INTERFACE RINDD
      MODULE PROCEDURE RINDD
      END INTERFACE

	INTERFACE SetConstants
      MODULE PROCEDURE SetConstants
      END INTERFACE

! mInfinity = what is considered as infinite value in FI
! mFxcEpss  = if fxc is less, do not compute E(...|Xc)
! mXcEps2   =  if any Var(Xc(j)|Xc(1),...,Xc(j-1)) <= XCEPS2 then return NAN
      double precision, parameter :: mInfinity = 8.25d0 ! 37.0d0
      double precision, parameter :: mFxcEpss = 1.0D-20
      double precision, save      :: mXcEps2   = 2.3d-16
!     Constants defining accuracy of integration:
!     mCovEps = termination criteria for Cholesky decomposition
!     mAbsEps = requested absolute tolerance
!     mRelEps = requested relative tolerance
!     mXcutOff = truncation value to c1c2
!     mXcScale = scale factor in the exponential (in order to avoid overflow)
!     mNc1c2  = number of times to use function c1c2, i.e.,regression
!               equation to restrict integration area.
!     mNIT    = maximum number of Xt variables to integrate
!     mMethod = integration method:
!            1 Integrate all by SADAPT if NDIM<9 otherwise by KRBVRC (default)
!            2 Integrate all by SADAPT if NDIM<19 otherwise by KRBVRC
!            3 Integrate all by KRBVRC by Genz (1998) (Fast and reliable)
!            4 Integrate all by KROBOV by Genz (1992) (Fast and reliable)
!            5 Integrate all by RCRUDE by Genz (1992) (Reliable)
!            6 Integrate all by SOBNIED by Hong and Hickernell
!            7 Integrate all by DKBVRC by Genz (2003) (Fast Ndim<1001)
      double precision, save :: mCovEps  = 1.0d-10
      double precision, save :: mAbsEps  = 0.01d0
      double precision, save :: mRelEps  = 0.01d0
      double precision, save :: mXcutOff = 5.d0
      double precision, save :: mXcScale = 0.0d0
      integer, save :: mNc1c2  = 2
      integer, save :: mNIT    = 1000
      integer, save :: mMaxPts = 40000
      integer, save :: mMinPts = 0
      integer, save :: mMethod = 3


!     Integrand variables:
!     mBIG    = Cholesky Factor/Covariance matrix:
!               Upper triangular part is the cholesky factor
!               Lower triangular part contains the conditional
!               standarddeviations
!               (mBIG2 is only used if mNx>1)
!     mCDI    = Cholesky DIagonal elements
!     mA,mB   = Integration limits
!     mINFI   = integrationi limit flags
!     mCm     = conditional mean
!     mINFIXt,
!     mINFIXd = # redundant variables of Xt and Xd,
!              respectively
!     mIndex1,
!     mIndex2 = indices to the variables original place.   Size  Ntdc
!     xedni   = indices to the variables new place.        Size  Ntdc
!     mNt     = # Xt variables
!     mNd     = # Xd variables
!     mNc     = # Xc variables
!     mNtd    = mNt + mNd
!     mNtdc   = mNt + mNd + mNc
!     mNx     = # different integration limits

      double precision,allocatable, dimension(:,:) :: mBIG,mBIG2
      double precision,allocatable, dimension(:)   :: mA,mB,mCDI,mCm
      INTEGER, DIMENSION(:),ALLOCATABLE :: mInfi,mIndex1,mIndex2,mXedni
      INTEGER,SAVE :: mNt,mNd,mNc,mNtdc, mNtd, mNx ! Size information
      INTEGER,SAVE :: mInfiXt,mInfiXd
      logical,save :: mInitIntegrandCalled = .FALSE.

      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: mCDIXd, mCmXd
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: mXd, mXc, mY
      double precision, save :: mSmall = 2.3d-16

!     variables set in initfun and used in mvnfun:
      INTEGER, PRIVATE :: mI0,mNdleftN0
      DOUBLE PRECISION, PRIVATE :: mE1,mD1, mVAL0

      contains
      subroutine setConstants(method,xcscale,abseps,releps,coveps,
     &     maxpts,minpts,nit,xcutoff,Nc1c2)
      double precision, optional, intent(in) :: xcscale,abseps,releps
     $     ,coveps, xcutoff
      integer, optional,intent(in) :: method,nit,maxpts,minpts,Nc1c2
      double precision, parameter :: one = 1.0d0
      mSmall = spacing(one)
      if (present(method))  mMethod = method
      if (present(xcscale)) mXcScale = xcscale
      if (present(abseps)) mAbsEps = max(abseps,mSmall)
      if (present(releps)) mRelEps = max(releps,0.0d0)
      if (present(coveps)) mCovEps = max(coveps,1d-12)
      if (present(maxpts)) mMaxPts = maxpts
      if (present(minpts)) mMinPts = minpts
      if (present(nit))    mNit    = nit
      if (present(xcutOff)) mXcutOff = xCutOff
      if (present(Nc1c2))  mNc1c2   = max(Nc1c2,1)
!      print *, 'method=', mMethod
!      print *, 'xcscale=', mXcScale
!      print *, 'abseps=', mAbsEps
!      print *, 'releps=', mRelEps
!      print *, 'coveps=', mCovEps
!      print *, 'maxpts=', mMaxPts
!      print *, 'minpts=', mMinPts
!      print *, 'nit=',    mNit
!      print *, 'xcutOff=', mXcutOff
!      print *, 'Nc1c2=',  mNc1c2
      end subroutine setConstants

      subroutine preInit(BIG,Xc,Nt,inform)
      double precision,dimension(:,:), intent(in) :: BIG
      double precision,dimension(:,:), intent(in) :: Xc
      integer, intent(in)  :: Nt
      integer, intent(out) :: inform
!     Local variables
      integer :: I,J
      inform = 0
      mInitIntegrandCalled = .FALSE.
!     Find the size information
!~~~~~~~~~~~~~~~~~~~~~~~~~~
      mNt   = Nt
      mNc   = SIZE( Xc, dim = 1 )
      mNx   = MAX( SIZE( Xc, dim = 2), 1 )
      mNtdc = SIZE( BIG, dim = 1 )
      ! make sure it does not exceed Ntdc-Nc
      IF (mNt+mNc.GT.mNtdc) mNt = mNtdc - mNc
      mNd  = mNtdc-mNt-mNc
      mNtd = mNt+mNd
      IF (mNd < 0) THEN
!         PRINT *,'RIND Nt,Nd,Nc,Ntdc=',Nt,Nd,Nc,Ntdc
         ! Size information inconsistent
         inform = 3
         return
      ENDIF

 !     PRINT *,'Nt Nd Nc Ntd Ntdc,',Nt, Nd, Nc, Ntd, Ntdc

! ALLOCATION
!~~~~~~~~~~~~
      IF (mNd>0) THEN
         ALLOCATE(mXd(mNd),mCmXd(mNd),mCDIXd(mNd))
         mCmXd(:)  = gZERO
         mCDIXd(:) = gZERO
         mxd(:)    = gZERO
      END IF
      ALLOCATE(mBIG(mNtdc,mNtdc),mCm(mNtdc),mY(mNtd))
      ALLOCATE(mIndex1(mNtdc),mA(mNtd),mB(mNtd),mINFI(mNtd),mXc(mNc))
      ALLOCATE(mCDI(mNtd),mXedni(mNtdc),mIndex2(mNtdc))

!     Initialization
!~~~~~~~~~~~~~~~~~~~~~
!     Copy upper triangular of input matrix, only.
      do i = 1,mNtdc
         mBIG(1:i,i)    = BIG(1:i,i)
      end do

      mIndex2 = (/(J,J=1,mNtdc)/)

!      CALL mexprintf('BIG Before CovsrtXc'//CHAR(10))
!      CALL ECHO(BIG)
!     sort BIG by decreasing cond. variance for Xc
      CALL CVSRTXC(mNt,mNd,mBIG,mIndex2,INFORM)
!      CALL mexprintf('BIG after CovsrtXc'//CHAR(10))
!      CALL ECHO(BIG)

      IF (INFORM.GT.0) return ! degenerate case exit VALS=0 for all
                                ! (should perhaps return NaN instead??)


      DO I=mNtdc,1,-1
         J = mIndex2(I)            ! covariance matrix according to index2
         mXedni(J) = I
      END DO

      IF (mNx>1) THEN
         ALLOCATE(mBIG2(mNtdc,mNtdc))
         do i = 1,mNtdc
            mBIG2(1:i,i) = mBIG(1:i,i) !Copy input matrix
         end do
      ENDIF
      return
      end subroutine preInit
      subroutine initIntegrand(ix,Xc,Ex,indI,Blo,Bup,INFIN,
     &     fxc,value,abserr,NDIM,inform)
      integer, intent(in) :: ix ! integrand number
      double precision, dimension(:),intent(in) :: Ex
      double precision, dimension(:,:), intent(in) :: Xc,Blo,Bup
      integer, dimension(:), intent(in) :: indI,INFIN
      double precision, intent(out) :: fxc,value,abserr
      integer, intent(out) :: NDIM, inform
!     Locals
      DOUBLE PRECISION   :: SQ0,xx,quant
      integer :: I,J
      inform = 0
      NDIM   = 0
      VALUE   = gZERO
      fxc     = gONE
      abserr  = mSmall

      IF (mInitIntegrandCalled)  then
         do i = 1,mNtdc
            mBIG(1:i,i)    = mBIG2(1:i,i) !Copy input matrix
         end do
      else
         mInitIntegrandCalled = .TRUE.
      endif

                                ! Set the original means of the variables
      mCm(:)  = Ex(mIndex2(1:mNtdc)) !   Cm(1:Ntdc)  =Ex (index1(1:Ntdc))
      IF (mNc>0) THEN
         mXc(:) = Xc(:,ix)
                                !mXc(1:Nc)    = Xc(1:Nc,ix)
         QUANT = DBLE(mNc)*LOG(gSQTWPI1)
         I = mNtdc
         DO J = 1, mNc
!     Iterative conditioning on the last Nc variables
            SQ0 = mBIG(I,I)     ! SQRT(Var(X(i)|X(i+1),X(i+2),...,X(Ntdc)))
            xx = (mXc(mIndex2(I) - mNtd) - mCm(I))/SQ0
                                !Trick to calculate
                                !fxc = fxc*SQTWPI1*EXP(-0.5*(XX**2))/SQ0
            QUANT = QUANT - gHALF*xx*xx  - LOG(SQ0)
                                ! conditional mean (expectation)
                                ! E(X(1:i-1)|X(i),X(i+1),...,X(Ntdc))
            mCm(1:I-1) = mCm(1:I-1) + xx*mBIG(1:I-1,I)
            I = I-1
         ENDDO
                                ! Calculating the
                                ! fxc probability density for i=Ntdc-J+1,
                                ! fXc=f(X(i)|X(i+1),X(i+2)...X(Ntdc))*
                                !     f(X(i+1)|X(i+2)...X(Ntdc))*..*f(X(Ntdc))
         fxc = EXP(QUANT+mXcScale)

                                ! if fxc small:  don't bother
                                ! calculating it, goto end
         IF (fxc  < mFxcEpss) then
            abserr = gONE
            inform = 1
            return
         endif
      END IF
!     Set integration limits mA,mB and mINFI
!     NOTE: mA and mB are integration limits with mCm subtracted
      CALL setIntLimits(mXc,indI,Blo,Bup,INFIN,inform)
      if (inform>0) return
      mIndex1(:) = mIndex2(:)
      CALL COVSRT(.FALSE., mNt,mNd,mBIG,mCm,mA,mB,mINFI,
     &        mINDEX1,mINFIXt,mINFIXd,NDIM,mY,mCDI)

      CALL INITFUN(VALUE,abserr,INFORM)
!     IF INFORM>0 : degenerate case:
!     Integral can be calculated excactly, ie.
!     mean of deterministic variables outside the barriers,
!     or NDIM = 1
      return
      end subroutine initIntegrand
      subroutine cleanUp
!     Deallocate all work arrays and vectors
      IF (ALLOCATED(mXc))     DEALLOCATE(mXc)
      IF (ALLOCATED(mXd))     DEALLOCATE(mXd)
      IF (ALLOCATED(mCm))     DEALLOCATE(mCm)
      IF (ALLOCATED(mBIG2))   DEALLOCATE(mBIG2)
      IF (ALLOCATED(mBIG))    DEALLOCATE(mBIG)
      IF (ALLOCATED(mIndex2)) DEALLOCATE(mIndex2)
      IF (ALLOCATED(mIndex1)) DEALLOCATE(mIndex1)
      IF (ALLOCATED(mXedni))  DEALLOCATE(mXedni)
      IF (ALLOCATED(mA))      DEALLOCATE(mA)
      IF (ALLOCATED(mB))      DEALLOCATE(mB)
      IF (ALLOCATED(mY))      DEALLOCATE(mY)
      IF (ALLOCATED(mCDI))    DEALLOCATE(mCDI)
      IF (ALLOCATED(mCDIXd))  DEALLOCATE(mCDIXd)
      IF (ALLOCATED(mCmXd))   DEALLOCATE(mCmXd)
      IF (ALLOCATED(mINFI))   DEALLOCATE(mINFI)
      end subroutine cleanUp
      function integrandBound(I0,N,Y,FINY) result (bound1)
      use FIMOD
      integer, intent(in) :: I0,N,FINY
      double precision, intent(in) :: Y
      double precision :: bound1
! locals
      integer :: I,IK,FINA, FINB
      double precision :: AI,BI,D1,E1
      double precision :: TMP
!     Computes the upper bound for the intgrand
      bound1 = gzero
      if (FINY<1) return
      FINA = 0
      FINB = 0
      IK = 2
      DO I = I0, N
             ! E(Y(I) | Y(1))/STD(Y(IK)|Y(1))
         TMP = mBIG(IK-1,I)*Y
         IF (mINFI(I) > -1) then
!     May have infinite int. Limits if Nd>0
            IF ( mINFI(I) .NE. 0 ) THEN
               IF ( FINA .EQ. 1 ) THEN
                  AI = MAX( AI, mA(I) - tmp )
               ELSE
                  AI   = mA(I) - tmp
                  FINA = 1
               END IF
            END IF
            IF ( mINFI(I) .NE. 1 ) THEN
               IF ( FINB .EQ. 1 ) THEN
                  BI = MIN( BI, mB(I) - tmp)
               ELSE
                  BI   = mB(I) - tmp
                  FINB = 1
               END IF
            END IF
         endif

         IF (I.EQ.N.OR.mBIG(IK+1,I+1)>gZERO) THEN
            CALL MVNLMS( AI, BI,2*FINA+FINB-1, D1, E1 )
            IF (D1<E1) bound1 = E1-D1
            return
         ENDIF
      ENDDO
      RETURN
      end function integrandBound
      SUBROUTINE INITFUN(VALUE,abserr,INFORM)
      USE JACOBMOD
      use SWAPMOD
      USE FIMOD
!      USE GLOBALDATA, ONLY: NIT,EPS2,EPS,xCutOff,NC1C2,ABSEPS
      IMPLICIT NONE
      DOUBLE PRECISION, INTENT(OUT) :: VALUE,abserr
      INTEGER, INTENT(out) :: INFORM
! local variables:
      INTEGER ::  N,NdleftO
      INTEGER :: I, J,  FINA, FINB
      DOUBLE PRECISION :: AI, BI, D0,E0,R12, R13, R23
      DOUBLE PRECISION :: xCut , upError,loError,maxTruncError
      LOGICAL :: useC1C2,isXd

!
!     Integrand subroutine
!
!INITFUN initialize the Multivariate Normal integrand function
! COF  - conditional sorted ChOlesky Factor of the covariance matrix (IN)
! CDI  - Cholesky DIagonal elements used to calculate the mean
! Cm   - conditional mean of Xd and Xt given Xc, E(Xd,Xt|Xc)
! xd   - variables to the jacobian variable, need no initialization size Nd
! xc   - conditional variables (IN)
! INDEX1 - if INDEX1(I)>Nt then variable no. I is one of the Xd
!          variables otherwise it is one of Xt.

      !PRINT *,'Mvnfun,ndim',Ndim
      INFORM = 0
      VALUE  = gZERO
      abserr = max(mCovEps , 6.0d0*mSmall)
      mVAL0   = gONE


      mNdleftN0 = mNd               ! Counter for number of Xd variables left

      mI0  = 0
      FINA = 0
      FINB = 0
      N = mNt + mNd - mINFIXt - mINFIXd-1
      IF (mINFIXt+mINFIXd > 0) THEN
!     CHCKLIM Check if the conditional mean Cm = E(Xt,Xd|Xc) for the
!     deterministic variables are between the barriers, i.e.,
!     A=Hlo-Cm< 0 <B=Hup-Cm
!     INFIN  INTEGER, array of integration limits flags:
!            if INFIN(I) < 0, Ith limits are (-infinity, infinity);
!            if INFIN(I) = 0, Ith limits are (-infinity, B(I)];
!            if INFIN(I) = 1, Ith limits are [A(I), infinity);
!            if INFIN(I) = 2, Ith limits are [A(I), B(I)].
         I = N+1
         DO J=1, mINFIXt + mINFIXd
            I = I + 1
            IF (mINFI(I)>-1) THEN
               IF ((mINFI(I).NE.0).AND.(mAbsEps  < mA(I))) GOTO 200
               IF ((mINFI(I).NE.1).AND.(mB(I) < -mAbsEps )) GOTO 200
            ENDIF
         ENDDO

         IF (mINFIXd>0) THEN
            ! Redundant variables of Xd: replace Xd with the mean
            I = mNt + mNd !-INFIS
            J = mNdleftN0-mINFIXd

            DO WHILE (mNdleftN0>J)
               isXd = (mNt < mIndex1(I))
               IF (isXd) THEN
                  mXd (mNdleftN0) =  mCm (I)
                  mNdleftN0 = mNdleftN0-1
               END IF
               I = I-1
            ENDDO
         ENDIF
         IF (N+1 < 1) THEN
!     Degenerate case, No relevant variables left to integrate
!     Print *,'rind ndim1',Ndim1
            IF (mNd>0) THEN
               VALUE = jacob (mXd,mXc) ! jacobian of xd,xc
            ELSE
               VALUE = gONE
            END IF
            GOTO 200
         ENDIF
      ENDIF
      IF (mNIT<=100) THEN
         xCut = mXcutOff

         J = 1
         DO I = 2, N+1
            IF (mBIG(J+1,I)>gZERO) THEN
               J = J + 1
            ELSE
               ! Add xCut std to deterministic variables to get an upper
               ! bound for integral
              mA(I) =  mA(I) - xCut * mBIG(I,J)
              mB(I) =  mB(I) + xCut * mBIG(I,J)
            ENDIF
         END DO
      ELSE
         xCut = gZERO
      ENDIF

      NdleftO = mNdleftN0
      useC1C2 = (1<=mNc1c2)
      DO I = 1, N+1
         IF (mINFI(I) > -1) then
!     May have infinite int. Limits if Nd>0
            IF ( mINFI(I) .NE. 0 ) THEN
               IF ( FINA .EQ. 1 ) THEN
                  AI = MAX( AI, mA(I) )
               ELSE
                  AI   = mA(I)
                  FINA = 1
               END IF
            END IF
            IF ( mINFI(I) .NE. 1 ) THEN
               IF ( FINB .EQ. 1 ) THEN
                  BI = MIN( BI, mB(I) )
               ELSE
                  BI   = mB(I)
                  FINB = 1
               END IF
            END IF
         endif
         isXd = (mINDEX1(I)>mNt)
         IF (isXd) THEN         ! Save the mean for Xd
            mCmXd(mNdleftN0)  = mCm(I)
            mCDIXd(mNdleftN0) = mCDI(I)
            mNdleftN0        = mNdleftN0-1
         END IF

         IF (I.EQ.N+1.OR.mBIG(2,I+1)>gZERO) THEN
            IF (useC1C2.AND.I<N) THEN
               mY(:) = gZERO


               CALL MVNLMS( AI, BI,2*FINA+FINB-1, D0, E0 )
               IF (D0>=E0) GOTO 200

               CALL C1C2(I+1,N+1,1,mA,mB,mINFI,mY,mBIG,AI,BI,FINA,FINB)
               CALL MVNLMS( AI, BI,2*FINA+FINB-1, mD1, mE1 )
               IF (mD1>=mE1) GOTO 200
               maxTruncError = FI(-ABS(mXcutOff))*dble(mNc1c2)
               upError = abs(E0-mE1)
               loError = abs(D0-mD1)
               if (upError>mSmall) then
                  upError = upError*integrandBound(I+1,N+1,BI,FINB)
               endif
               if (loError>mSmall) then
                  loError = loError*integrandBound(I+1,N+1,AI,FINA)
               endif
               abserr  = abserr + min(upError + loError,maxTruncError)
               !CALL printvar(log10(loError+upError+msmall),'lo+up-err')
            ELSE
               CALL MVNLMS( AI, BI,2*FINA+FINB-1, mD1, mE1 )
               IF (mD1>=mE1) GOTO 200
            ENDIF
            !CALL MVNLMS( AI, BI,2*FINA+FINB-1, mD1, mE1 )
            !IF (mD1>=mE1) GOTO 200
            IF ( NdleftO<=0) THEN
               IF (mNd>0) mVAL0 = JACOB(mXd,mXc)
               SELECT CASE (I-N)
               CASE (1)   !IF (I.EQ.N+1) THEN
                  VALUE  = (mE1-mD1)*mVAL0
                  abserr = abserr*mVAL0
                  GO TO 200
               CASE (0)     !ELSEIF (I.EQ.N) THEN
                                !D1=1/sqrt(1-rho^2)=1/STD(X(I+1)|X(1))
                  mD1 = SQRT( gONE + mBIG(1,I+1)*mBIG(1,I+1) )
                  mINFI(2) = mINFI(I+1)
                  mA(1) = AI
                  mB(1) = BI
                  mINFI(1) = 2*FINA+FINB-1
                  IF ( mINFI(2) .NE. 0 ) mA(2) = mA(I+1)/mD1
                  IF ( mINFI(2) .NE. 1 ) mB(2) = mB(I+1)/mD1
                  VALUE = BVNMVN( mA, mB,mINFI,mBIG(1,I+1)/mD1 )*mVAL0
                  abserr = (abserr+1.0d-14)*mVAL0
                  GO TO 200
               CASE ( -1 )  !ELSEIF (I.EQ.N-1) THEN
                  IF (.FALSE.) THEN
! TODO :this needs further checking! (it should work though)
                  !1/D1= sqrt(1-r12^2) = STD(X(I+1)|X(1))
                  !1/E1=  STD(X(I+2)|X(1)X(I+1))
                  !D1  = BIG(I+1,1)
                  !E1  = BIG(I+2,2)

                  mD1 = gONE/SQRT( gONE + mBIG(1,I+1)*mBIG(1,I+1) )
                  R12 = mBIG( 1, I+1 ) * mD1
                  if (mBIG(3,I+2)>gZERO) then
                     mE1 = gONE/SQRT( gONE + mBIG(1,I+2)*mBIG(1,I+2) +
     &                    mBIG(2,I+2)*mBIG(2,I+2) )
                     R13 = mBIG( 1, I+2 ) * mE1
                     R23 = mBIG( 2, I+2 ) * (mE1 * mD1) + R12 * R13
                  else
                     mE1  = mCDI(I+2)
                     R13 = mBIG( 1, I+2 ) * mE1
                     R23 = mE1*mD1 + R12 * R13
                     IF ((mE1  <  gZERO).AND. mINFI(I+2)>-1) THEN
                        CALL SWAP(mA(I+2),mB(I+2))
                        IF (mINFI(I+2).NE. 2) mINFI(I+2) = 1-mINFI(I+2)
                     END IF
                     !R23 = BIG( 2, I+2 ) * (E1 * D1) + R12 * R13
                  endif
                  mINFI(2) = mINFI(I+1)
                  mINFI(3) = mINFI(I+2)
                  mA(1) = AI
                  mB(1) = BI
                  mINFI(1) = 2*FINA+FINB-1
                  IF ( mINFI(2) .NE. 0 ) mA(2) = mA(I+1) * mD1
                  IF ( mINFI(2) .NE. 1 ) mB(2) = mB(I+1) * mD1
                  IF ( mINFI(3) .NE. 0 ) mA(3) = mA(I+2) * mE1
                  IF ( mINFI(3) .NE. 1 ) mB(3) = mB(I+2) * mE1
                  if(.false.) then
                     CALL PRINTVECD((/R12, R13, R23 /),'R12 = ')
                     CALL PRINTVECD((/mD1, mE1 /),'D1 = ')
                     CALL PRINTVECD(mBIG(1,1:3),'BIG(1,1:3) = ')
                     CALL PRINTVECD(mBIG(2,2:3),'BIG(2,2:3) = ')
                     CALL PRINTVECD(mBIG(1:3,1),'BIG(1:3,1) = ')
                     CALL PRINTVECD(mBIG(2:3,2),'BIG(2:3,2) = ')
                     CALL PRINTVECD(mA(1:I+2),'A = ')
                     CALL PRINTVECD(mB(1:I+2),'B = ')
                     CALL PRINTVECI(mINFI(1:I+2),'INFI = ')
                     CALL PRINTVECI(mINDEX1(1:I+2),'index1 = ')
                  endif
                  VALUE = TVNMVN( mA, mB,mINFI,
     &                 (/R12, R13, R23 /),1.0d-13) * mVAL0
                  ABSERR = (ABSERR + 1.0d-13)*mVAL0
                  GOTO 200
                  ENDIF
               END SELECT !ENDIF
            ENDIF
            ABSERR = mVAL0*ABSERR
            mVAL0 = mVAL0 * (mE1-mD1)
            mI0   = I
            RETURN
         ENDIF
      ENDDO
      RETURN
 200  INFORM = 1
      RETURN
      END SUBROUTINE INITFUN
!
!     Integrand subroutine
!
      FUNCTION MVNFUN( Ndim, W ) RESULT (VAL)
      USE JACOBMOD
      USE FIMOD
      IMPLICIT NONE
      INTEGER, INTENT (IN) :: Ndim
      DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: W
      DOUBLE PRECISION :: VAL
! local variables:
      INTEGER ::  N,I, J, FINA, FINB
      INTEGER ::  NdleftN, NdleftO ,IK
      DOUBLE PRECISION :: TMP, AI, BI, DI, EI
      LOGICAL :: useC1C2, isXd
!MVNFUN Multivariate Normal integrand function
! where the integrand is transformed from an integral
! having integration limits A and B to an
! integral having  constant integration limits i.e.
!   B                                  1
! int jacob(xd,xc)*f(xd,xt)dxt dxd = int F2(W) dW
!  A                                  0
!
! W    - new transformed integration variables, valid range 0..1
!        The vector must have the length Ndim returned from Covsrt
! mBIG - conditional sorted ChOlesky Factor of the covariance matrix (IN)
! mCDI - Cholesky DIagonal elements used to calculate the mean
! mCm  - conditional mean of Xd and Xt given Xc, E(Xd,Xt|Xc)
! mXd  - variables to the jacobian variable, need no initialization size Nd
! mXc  - conditional variables (IN)
! mINDEX1 - if mINDEX1(I)>Nt then variable No. I is one of the Xd
!           variables otherwise it is one of Xt

      !PRINT *,'Mvnfun,ndim',Ndim

!      xCut = gZERO ! xCutOff

      N    = mNt+mNd-mINFIXt-mINFIXd-1
      IK   = 1                    ! Counter for Ndim
      FINA = 0
      FINB = 0

      NdleftN = mNdleftN0          ! Counter for number of Xd variables left
      VAL     = mVAL0
      NdleftO = mNd - mINFIXd
      mY(IK) = FIINV( mD1 + W(IK)*( mE1 - mD1 ) )
      useC1C2 = (IK+1.LE.mNc1c2)
      IF (useC1C2) THEN
         ! Calculate the conditional mean
         ! E(Y(I) | Y(1),...Y(I0))/STD(Y(I)|Y(1),,,,Y(I0))
         mY(mI0+1:N+1) = mBIG(IK, mI0+1:N+1)*mY(IK)
      ENDIF
      IF (NdleftO.GT.NdleftN ) THEN
         mXd(NdleftN+1:NdleftO) = mCmXd(NdleftN+1:NdleftO)+
     &        mY(IK) * mCDIXd(NdleftN+1:NdleftO)
      ENDIF
      NdleftO = NdleftN
      IK = 2                    !=IK+1


      DO I = mI0+1, N+1
         IF (useC1C2) THEN
             TMP = mY(I)
          ELSE
            TMP = 0.d0
            DO J = 1, IK-1
               ! E(Y(I) | Y(1),...Y(IK-1))/STD(Y(IK)|Y(1),,,,Y(IK-1))
               TMP = TMP + mBIG(J,I)*mY(J)
            END DO
         ENDIF
         IF (mINFI(I) < 0) GO TO 100
            ! May have infinite int. Limits if Nd>0
         IF ( mINFI(I) .NE. 0 ) THEN
            IF ( FINA .EQ. 1 ) THEN
               AI = MAX( AI, mA(I) - TMP)
            ELSE
               AI = mA(I) - TMP
               FINA = 1
            END IF
            IF (FINB.EQ.1.AND.BI<=AI) GOTO 200
         END IF
         IF ( mINFI(I) .NE. 1 ) THEN
            IF ( FINB .EQ. 1 ) THEN
               BI = MIN( BI, mB(I) - TMP)
            ELSE
               BI = mB(I) - TMP
               FINB = 1
            END IF
            IF (FINA.EQ.1.AND.BI<=AI) GOTO 200
         END IF
 100     isXd = (mNt<mINDEX1(I))
         IF (isXd) THEN
!     Save the mean of xd and Covariance diagonal element
            ! Conditional mean E(Xi|X1,..X)
            mCmXd(NdleftN)  = mCm(I) + TMP * mCDI(I)
             ! Covariance diagonal
            mCDIXd(NdleftN) = mCDI(I)
            NdleftN        = NdleftN - 1
         END IF
         IF (I == N+1 .OR. mBIG(IK+1,I+1) > gZERO ) THEN
            IF (useC1C2) THEN
!     Note: for J =I+1:N+1:  Y(J) = conditional expectation, E(Yj|Y1,...Yk)
               CALL C1C2(I+1,N+1,IK,mA,mB,mINFI,mY,mBIG,AI,BI,FINA,FINB)
            ENDIF
            CALL MVNLMS( AI, BI, 2*FINA+FINB-1, DI, EI )
            IF ( DI >= EI ) GO TO 200
            VAL = VAL * ( EI - DI )

            IF ( I <= N .OR. (NdleftN < NdleftO)) THEN
               mY(IK) = FIINV( DI + W(IK)*( EI - DI ) )
               IF (NdleftN < NdleftO ) THEN
                  mXd(NdleftN+1:NdleftO) = mCmXd(NdleftN+1:NdleftO)+
     &                 mY(IK) * mCDIXd(NdleftN+1:NdleftO)
                  NdleftO = NdleftN
               ENDIF
               useC1C2 = (IK+1<=mNc1c2)
               IF (useC1C2) THEN

                  ! E(Y(J) | Y(1),...Y(I))/STD(Y(J)|Y(1),,,,Y(I))
                  mY(I+1:N+1) = mY(I+1:N+1) + mBIG(IK, I+1:N+1)*mY(IK)
               ENDIF
            ENDIF
            IK   = IK + 1
            FINA = 0
            FINB = 0
         END IF
      END DO
      IF (mNd>0) VAL = VAL * jacob(mXd,mXc)
      RETURN
 200  VAL = gZERO
      RETURN
      END FUNCTION MVNFUN


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!******************* RINDD - the main program *********************!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE RINDD(VALS,ERR,TERR,Big,Ex,Xc,Nt,
     &               indI,Blo,Bup,INFIN)
      USE RCRUDEMOD
      USE KRBVRCMOD
      USE ADAPTMOD
      USE KROBOVMOD
      USE DKBVRCMOD
      USE SSOBOLMOD
      IMPLICIT NONE
      DOUBLE PRECISION, DIMENSION(:  ), INTENT(out):: VALS, ERR ,TERR
      DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: BIG
      DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: Xc
      DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: Ex
      DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: Blo, Bup
      INTEGER,          DIMENSION(:), INTENT(in) :: indI,INFIN
      INTEGER,                          INTENT(in) :: Nt
!      DOUBLE PRECISION,                 INTENT(in) :: XcScale
! local variables
      INTEGER :: ix, INFORM, NDIM, MAXPTS, MINPTS
      DOUBLE PRECISION :: VALUE,fxc,absERR,absERR2
      double precision :: LABSEPS,LRELEPS


      VALS(:) = gZERO
      ERR(:)  = gONE
      TERR(:) = gONE

      call preInit(BIG,Xc,Nt,inform)
      IF (INFORM.GT.0) GOTO 110 ! degenerate case exit VALS=0 for all
                                ! (should perhaps return NaN instead??)

!     Now the loop over all different values of
!     variables Xc (the one one is conditioning on)
!     is started. The density f_{Xc}(xc(:,ix))
!     will be computed and denoted by  fxc.
      DO  ix = 1, mNx
         call initIntegrand(ix,Xc,Ex,indI,Blo,Bup,infin,
     &        fxc,value,abserr,NDIM,inform)


         IF (INFORM.GT.0) GO TO 100

         MAXPTS  = mMAXPTS
         MINPTS  = mMINPTS
         LABSEPS = max(mABSEPS-abserr,0.2D0*mABSEPS)       !*fxc
         LRELEPS = mRELEPS
         ABSERR2 = mSmall

         SELECT CASE (mMethod)
         CASE (:1)
            IF (NDIM < 9) THEN
               CALL SADAPT(NDIM,MAXPTS,MVNFUN,LABSEPS,
     &              LRELEPS,ABSERR2,VALUE,INFORM)
               VALUE = MAX(VALUE,gZERO)
            ELSE
               CALL KRBVRC(NDIM, MINPTS, MAXPTS, MVNFUN,LABSEPS,LRELEPS,
     &              ABSERR2, VALUE, INFORM )
            ENDIF
         CASE (2)
!        Call the subregion adaptive integration subroutine
            IF ( NDIM .GT. 19.) THEN
!     print *, 'Ndim too large for SADMVN => Calling KRBVRC'
               CALL KRBVRC( NDIM, MINPTS, MAXPTS, MVNFUN, LABSEPS,
     &              LRELEPS, ABSERR2, VALUE, INFORM )
            ELSE
               CALL SADAPT(NDIM,MAXPTS,MVNFUN,LABSEPS,
     &              LRELEPS,ABSERR2,VALUE,INFORM)
               VALUE = MAX(VALUE,gZERO)
            ENDIF
         CASE (3)               !       Call the Lattice rule integration procedure
            CALL KRBVRC( NDIM, MINPTS, MAXPTS, MVNFUN, LABSEPS,
     &           LRELEPS, ABSERR2, VALUE, INFORM )
         CASE (4)               !       Call the Lattice rule
                                !       integration procedure
            CALL KROBOV( NDIM, MINPTS, MAXPTS, MVNFUN, LABSEPS,
     &           LRELEPS,ABSERR2, VALUE, INFORM )
         CASE (5)    ! Call Crude Monte Carlo integration procedure
            CALL RANMC( NDIM, MAXPTS, MVNFUN, LABSEPS,
     &           LRELEPS, ABSERR2, VALUE, INFORM )
         CASE (6)       !       Call the scrambled Sobol sequence rule integration procedure
           CALL SOBNIED( NDIM, MINPTS, MAXPTS, MVNFUN, LABSEPS, LRELEPS,
     &           ABSERR2, VALUE, INFORM )
         CASE (7:)
           CALL DKBVRC( NDIM, MINPTS, MAXPTS, MVNFUN, LABSEPS, LRELEPS,
     &           ABSERR2, VALUE, INFORM )
         END SELECT

!     IF (INFORM.gt.0) print *,'RIND, INFORM,error =',inform,error
 100     VALS(ix) = VALUE*fxc
         IF (SIZE(ERR, DIM = 1).EQ.mNx) ERR(ix)   = abserr2*fxc
         IF (SIZE(TERR, DIM = 1).EQ.mNx) TERR(ix) = abserr*fxc
      ENDDO                     !ix

 110  CONTINUE
      call cleanUp
      RETURN
      END SUBROUTINE RINDD

      SUBROUTINE setIntLimits(xc,indI,Blo,Bup,INFIN,inform)
      IMPLICIT NONE
      DOUBLE PRECISION, DIMENSION(:  ), INTENT(in)  :: xc
      INTEGER,          DIMENSION(:  ), INTENT(in)  :: indI,INFIN
      DOUBLE PRECISION, DIMENSION(:,:), INTENT(in)  :: Blo,Bup
      integer, intent(out) :: inform
!Local variables
      INTEGER   :: I, J, K, L,Mb1,Nb,NI,Nc
      DOUBLE PRECISION :: xCut, SQ0
!this procedure set mA,mB and mInfi according to Blo/Bup and INFIN
!
!     INFIN  INTEGER, array of integration limits flags:
!            if INFIN(I) < 0, Ith limits are (-infinity, infinity);
!            if INFIN(I) = 0, Ith limits are (-infinity, mB(I)];
!            if INFIN(I) = 1, Ith limits are [mA(I), infinity);
!            if INFIN(I) = 2, Ith limits are [mA(I), mB(I)].
! Note on member variables:
! mXedni     = indices to the variables new place after cvsrtXc.  Size Ntdc
! mCm        = E(Xt,Xd|Xc), i.e., conditional mean given Xc
! mBIG(:,1:Ntd) = Cov(Xt,Xd|Xc)

      xCut = ABS(mInfinity)
      Mb1 = size(Blo,DIM=1)-1
      Nb = size(Blo,DIM=2)
      NI = size(indI,DIM=1)
      Nc = size(xc,DIM=1)
      if (Mb1>Nc .or. Nb.NE.NI-1) then
!     size of variables inconsistent
         inform = 4
         return
      endif

!      IF (Mb.GT.Nc+1) print *,'barrier: Mb,Nc =',Mb,Nc
!      IF (Nb.NE.NI-1) print *,'barrier: Nb,NI =',Nb,NI
      DO J = 2, NI
         DO I = indI (J - 1) + 1 , indI (J)
            L        = mXedni(I)
            mINFI(L) = INFIN(J-1)
            SQ0      = SQRT(mBIG(L,L))
            mA(L)    = -xCut*SQ0
            mB(L)    =  xCut*SQ0
            IF (mINFI(L).GE.0) THEN
               IF  (mINFI(L).NE.0) THEN
                  mA(L) = Blo (1, J - 1)-mCm(L)
                  DO K = 1, Mb1
                     mA(L) = mA(L)+Blo(K+1,J-1)*xc(K)
                  ENDDO         ! K
                  ! This can only be done if
                  if (mA(L)< -xCut*SQ0) mINFI(L) = mINFI(L)-2
               ENDIF
               IF  (mINFI(L).NE.1) THEN
                  mB(L) = Bup (1, J - 1)-mCm(L)
                  DO K = 1, Mb1
                     mB(L) = mB(L)+Bup(K+1,J-1)*xc(K)
                  ENDDO
                  if (xCut*SQ0<mB(L)) mINFI(L) = mINFI(L)-1
               ENDIF            !
            ENDIF
         ENDDO                  ! I
      ENDDO                     ! J
!     print * ,'barrier hup:',size(Hup),Hup(xedni(1:indI(NI)))
!     print * ,'barrier hlo:',size(Hlo),Hlo(xedni(1:indI(NI)))
      RETURN
      END SUBROUTINE setIntLimits


      FUNCTION GETTMEAN(A,B,INFJ,PRB) RESULT (MEAN1)
      USE GLOBALCONST
      IMPLICIT NONE
!     GETTMEAN Returns the expected mean, E(I(a<x<b)*X)
      DOUBLE PRECISION, INTENT(IN) :: A,B,PRB
      INTEGER , INTENT(IN) :: INFJ
      DOUBLE PRECISION :: MEAN1
      DOUBLE PRECISION :: YL,YU
!     DOUBLE PRECISION, PARAMETER:: ZERO = 0.0D0, HALF = 0.5D0

      IF ( PRB .GT. gZERO) THEN
         YL = gZERO
         YU = gZERO
         IF (INFJ.GE.0) THEN
            IF (INFJ .NE. 0) YL =-EXP(-gHALF*(A*A))*gSQTWPI1
            IF (INFJ .NE. 1) YU =-EXP(-gHALF*(B*B))*gSQTWPI1
         ENDIF
         MEAN1 = ( YU - YL )/PRB
      ELSE
         SELECT CASE (INFJ)
         CASE (:-1)
            MEAN1 = gZERO
         CASE (0)
            MEAN1 = B
         CASE (1)
            MEAN1 = A
         CASE (2:)
            MEAN1 = ( A + B ) * gHALF
         END SELECT
      END IF
      RETURN
      END FUNCTION
      SUBROUTINE ADJLIMITS(A,B, infi)
!      USE GLOBALDATA, ONLY : xCutOff
      IMPLICIT NONE
!     Adjust INFI when integration limits A and/or B is too far out in the tail
      DOUBLE PRECISION, INTENT(IN)     :: A,B
      INTEGER,          INTENT(IN OUT) :: infi
!      DOUBLE PRECISION, PARAMETER :: xCutOff = 8.D0
      IF (infi>-1) THEN
         IF (infi.NE.0)THEN
            IF (A  <  -mXcutOff) THEN
               infi = infi-2
!               CALL mexprintf('ADJ A')
            ENDIF
         ENDIF
         IF (infi.NE.1) THEN
            IF (mXCutOff  <  B) THEN
               infi = infi-1
!               CALL mexprintf('ADJ B')
            ENDIF
         END IF
      END IF
      RETURN
      END SUBROUTINE ADJLIMITS
      SUBROUTINE C1C2(I0,I1,IK,A,B,INFIN, Cm, BIG, AJ, BJ, FINA,FINB)
! The regression equation for the conditional distr. of Y given X=x
! is equal  to the conditional expectation of Y given X=x, i.e.,
!
!       E(Y|X=x) = E(Y) + Cov(Y,X)/Var(X)[x-E(X)]
!
!  Let
!     x1 = (x-E(X))/SQRT(Var(X)) be zero mean,
!     C1< x1 <C2,
!     B1(I) = COV(Y(I),X)/SQRT(Var(X)).
!  Then the process  Y(I) with mean Cm(I) can be written as
!
!       y(I) = Cm(I) + x1*B1(I) + Delta(I) for  I=I0,...,I1.
!
!  where SQ(I) = sqrt(Var(Y|X)) is the standard deviation of Delta(I).
!
!  Since we are truncating all Gaussian  variables to
!  the interval [-C,C], then if for any I
!
!  a) Cm(I)+x1*B1(I)-C*SQ(I)>B(I)  or
!
!  b) Cm(I)+x1*B1(I)+C*SQ(I)<A(I)  then
!
!  the (XIND|Xn=xn) = 0 !!!!!!!!!
!
!  Consequently, for increasing the accuracy (by excluding possible
!  discontinuouities) we shall exclude such values for which (XIND|X1=x1) = 0.
!  Hence we assume that if Aj<x<Bj any of the previous conditions are
!  satisfied
!
!  OBSERVE!!, Aj, Bj has to be set to (the normalized) lower and upper bounds
!  of possible values for x1,respectively, i.e.,
!           Aj=max((A-E(X))/SQRT(Var(X)),-C), Bj=min((B-E(X))/SQRT(Var(X)),C)
!  before calling C1C2 subroutine.
!
      USE GLOBALCONST
!      USE GLOBALDATA, ONLY : EPS2,EPS ,xCutOff
      IMPLICIT NONE
      DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: Cm, A,B !, B1, SQ
      DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: BIG
      INTEGER,          DIMENSION(:), INTENT(in) :: INFIN
      DOUBLE PRECISION,            INTENT(inout) :: AJ,BJ
      INTEGER,                     INTENT(inout) :: FINA, FINB
      INTEGER, INTENT(IN) :: I0, I1, IK

!     Local variables
      DOUBLE PRECISION   :: xCut
!      DOUBLE PRECISION, PARAMETER :: TOL = 1.0D-16
      DOUBLE PRECISION :: CSQ, BdSQ0, LTOL
      INTEGER :: INFI, I

      xCut = MIN(ABS(mXcutOff),mInfinity)
      LTOL = mSmall ! EPS2
!      AJ = MAX(AJ,-xCut)
!      BJ = MIN(BJ,xCut)
!      IF (AJ.GE.BJ) GO TO 112
!      CALL PRINTVAR(AJ,TXT='BC1C2: AJ')
!      CALL PRINTVAR(BJ,TXT='BC1C2: BJ')

      IF (I1 < I0)  RETURN       !Not able to change integration limits
      DO I = I0,I1
!     C = xCutOff
         INFI = INFIN(I)
         IF (INFI>-1) THEN
            !BdSQ0 = B1(I)
            !CSQ   = xCut * SQ(I)
            BdSQ0  = BIG(IK,I)
            CSQ    = xCut * BIG(I,IK)
            IF (BdSQ0 > LTOL) THEN
               IF ( INFI .NE. 0 ) THEN
                  IF (FINA.EQ.1) THEN
                     AJ = MAX(AJ,(A(I) - Cm(I) - CSQ)/BdSQ0)
                  ELSE
                     AJ = (A(I) - Cm(I) - CSQ)/BdSQ0
                     FINA = 1
                  ENDIF
                  IF (FINB.GT.0) AJ = MIN(AJ,BJ)
               END IF
               IF ( INFI .NE. 1 ) THEN
                  IF (FINB.EQ.1) THEN
                     BJ = MIN(BJ,(B(I) - Cm(I) + CSQ)/BdSQ0)
                  ELSE
                     BJ = (B(I) - Cm(I) + CSQ)/BdSQ0
                     FINB = 1
                  ENDIF
                  IF (FINA.GT.0) BJ = MAX(AJ,BJ)
               END IF
             ELSEIF (BdSQ0  <  -LTOL) THEN
                IF ( INFI .NE. 0 ) THEN
                  IF (FINB.EQ.1) THEN
                     BJ = MIN(BJ,(A(I) - Cm(I) - CSQ)/BdSQ0)
                  ELSE
                     BJ = (A(I) - Cm(I) - CSQ)/BdSQ0
                     FINB = 1
                  ENDIF
                  IF (FINA.GT.0) BJ = MAX(AJ,BJ)
               END IF
               IF ( INFI .NE. 1 ) THEN
                  IF (FINA.EQ.1) THEN
                     AJ = MAX(AJ,(B(I) - Cm(I) + CSQ)/BdSQ0)
                  ELSE
                     AJ = (B(I) - Cm(I) + CSQ)/BdSQ0
                     FINA = 1
                  ENDIF
                  IF (FINB.GT.0) AJ = MIN(AJ,BJ)
               END IF
            END IF
         ENDIF
      END DO
!      IF (FINA>0 .AND. FINB>0) THEN
!         IF (AJ<BJ) THEN
!            IF (AJ   <= -xCut) FINA = 0
!            IF (xCut <= BJ   ) FINB = 0
!         ENDIF
!      ENDIF
!      CALL PRINTVAR(AJ,TXT='AC1C2: AJ')
!      CALL PRINTVAR(BJ,TXT='AC1C2: BJ')
      RETURN
      END SUBROUTINE C1C2
      SUBROUTINE CVSRTXC (Nt,Nd,R,index1,INFORM)
!      USE GLOBALDATA, ONLY :  XCEPS2
!      USE GLOBALCONST
      IMPLICIT NONE
      INTEGER, INTENT(in) :: Nt,Nd
      DOUBLE PRECISION, DIMENSION(:,:), INTENT(inout) :: R
      INTEGER,          DIMENSION(:  ), INTENT(inout) :: index1
      INTEGER, INTENT(out) :: INFORM
! local variables
      DOUBLE PRECISION, DIMENSION(:), allocatable   :: SQ
      INTEGER,          DIMENSION(1)                :: m
      INTEGER :: M1,K,I,J,Ntdc,Ntd,Nc, LO
      DOUBLE PRECISION :: LTOL, maxSQ
!     if any Var(Xc(j)|Xc(1),...,Xc(j-1)) <= XCEPS2 then return NAN
      double precision :: XCEPS2
!CVSRTXC calculate the conditional covariance matrix of Xt and Xd given Xc
! as well as the cholesky factorization for the Xc variable(s)
! The Xc variables are sorted by the largest conditional covariance
!
! R         = In : Cov(X) where X=[Xt Xd Xc] is stochastic vector
!             Out: sorted Conditional Covar. matrix, i.e.,
!                    [ Cov([Xt,Xd] | Xc)    Shape N X N  (N=Ntdc=Nt+Nd+Nc)
! index1    = In/Out : permutation vector giving the indices to the variables
!             original place.   Size  Ntdc
! INFORM    = Out, Returns
!             0 If Normal termination.
!             1 If R is degenerate, i.e., Cov(Xc) is singular.
!
! R=Cov([Xt,Xd,Xc]) is a covariance matrix of the stochastic
! vector X=[Xt Xd Xc] where the variables Xt, Xd and Xc have the size
! Nt, Nd and Nc, respectively.
! Xc are the conditional variables.
! Xd and Xt are the variables to integrate.
!(Xd,Xt = variables in the jacobian and indicator respectively)
!
! Note: CVSRTXC only works on the upper triangular part of R

      INFORM = 0
      Ntdc = size(R,DIM=1)
      Ntd  = Nt   + Nd
      Nc   = Ntdc - Ntd

      IF (Nc < 1) RETURN



      ALLOCATE(SQ(1:Ntdc))
      maxSQ = gZERO
      DO I = 1, Ntdc
         SQ(I)  = R(I,I)
         if (SQ(I)>maxSQ) maxSQ = SQ(I)
      ENDDO

      XCEPS2 = Ntdc*mSmall*maxSQ
      mXcEps2 = XCEPS2
      LTOL   = mSmall
      LO = 1
      K = Ntdc
      DO I = 1, Nc             ! Condsort Xc
         m  = K+1-MAXLOC(SQ(K:Ntd+1:-1))
         M1 = m(1)
         IF (SQ(m1)<=XCEPS2) THEN
!     PRINT *,'CVSRTXC: Degenerate case of Xc(Nc-J+1) for J=',ix
            !CALL mexprintf('CVSRTXC: Degenerate case of Xc(Nc-J+1)')
            INFORM = 1
            GOTO 200   ! RETURN    !degenerate case
         ENDIF
         IF (M1.NE.K) THEN
            ! Symmetric row column permuations
            ! Swap row and columns, but only upper triangular part
            CALL RCSWAP( M1, K, Ntdc,Ntd, R,INDEX1,SQ)
         END IF
         R(K,K) = SQRT(SQ(K))
         IF (K .EQ. LO) GOTO 200
         R(LO:K-1,K) = R(LO:K-1,K)/R(K,K)
! Cov(Xi,Xj|Xk,Xk+1,..,Xn)  = ....
!         Cov(Xi,Xj|Xk+1,..,Xn) - Cov(Xi,Xk|Xk+1,..Xn)*Cov(Xj,Xk|Xk+1,..Xn)
         DO J = LO,K-1
                                ! Var(Xj | Xk,Xk+1,...,Xn)
            SQ(J)  =  R(J,J) - R(J,K)*R(J,K)
            IF (SQ(J)<=LTOL.AND.J<=Ntd) THEN
               IF (LO < J) THEN
                  CALL RCSWAP(LO, J, Ntdc,Ntd, R,INDEX1,SQ)
               ENDIF
               R(LO,LO:K-1) = gZERO
               IF (SQ(LO) < -10.0D0*SQRT(LTOL)) THEN
                  ! inform = 2
                  !R(LO,K) = gZERO
                 ! CALL mexprintf('Negative definit BIG!'//CHAR(10))
               ENDIF
               SQ(LO) = gZERO
               LO = LO + 1
            ELSE
               R(J,J) = SQ(J)
               R(LO:J-1,J) = R(LO:J-1,J) - R(LO:J-1,K)*R(J,K)
            ENDIF
         END DO
         K = K - 1
      ENDDO
 200  DEALLOCATE(SQ)
      RETURN
      END SUBROUTINE CVSRTXC

      SUBROUTINE  RCSCALE(chkLim,K,K0,N1,N,K1,INFIS,CDI,Cm,
     &     R,A,B,INFI,INDEX1,Y)
      USE GLOBALCONST
      USE SWAPMOD
      IMPLICIT NONE
!RCSCALE:  Scale  covariance matrix and limits
!
!   CALL  RCSCALE( k, k0, N1, N,K1, CDI,Cm,R,A, B, INFIN,index1,Y)
!
!  chkLim  = TRUE  if check if variable K is redundant
!            FALSE
!    K     = index to variable which is deterministic,i.e.,
!            STD(Xk|X1,...Xr) = 0
!    N1    = Number of significant variables of [Xt,Xd]
!    N     = length(Xt)+length(Xd)
!    K1    = index to current variable we are conditioning on.
!   CDI    = Cholesky diagonal elements which contains either
!               CDI(J) = STD(Xj | X1,...,Xj-1,Xc) if Xj is stochastic given
!                        X1,...Xj, Xc
!              or
!               CDI(J) = COV(Xj,Xk | X1,..,Xk-1,Xc  )/STD(Xk | X1,..,Xk-1,Xc)
!               if Xj is determinstically determined given X1,..,Xk,Xc
!               for some k<j.
!   Cm     = conditional mean
!    R     = Matrix containing cholesky factor for
!            X = [Xt,Xd,Xc] in upper triangular part. Lower triangular
!            part contains conditional stdevs. size Ntdc X Ntdc
! INDEX1   = permutation index vector. (index to variables original place).
!    A,B   = lower and upper integration limit, respectively.
!    INFIN = if INFIN(I) < 0, Ith limits are (-infinity, infinity);
!            if INFIN(I) = 0, Ith limits are (-infinity, B(I)];
!            if INFIN(I) = 1, Ith limits are [A(I), infinity);
!            if INFIN(I) = 2, Ith limits are [A(I), B(I)].
!    Y     = work vector
!
! NOTE: RCSWAP works only on the upper triangular part of C
! + check if variable k is redundant
!  If the conditional covariance matrix diagonal entry is zero,
!  permute limits and/or rows, if necessary.
      LOGICAL, INTENT(IN) :: chkLim
      INTEGER, INTENT(IN) :: K, K0,N
      INTEGER, INTENT(INOUT) :: N1,K1,INFIS
      DOUBLE PRECISION, DIMENSION(:),  INTENT(INOUT) :: CDI,A,B,Cm
      DOUBLE PRECISION, DIMENSION(:,:),INTENT(INOUT) :: R
      INTEGER,          DIMENSION(:),  INTENT(INOUT) :: INFI,INDEX1
      DOUBLE PRECISION, DIMENSION(:),OPTIONAL,INTENT(INOUT) :: Y
!Local variables
      DOUBLE PRECISION, PARAMETER :: LTOL = 1.0D-16
      double precision :: xCut
      DOUBLE PRECISION :: D,AK,BK,CVDIAG
      INTEGER :: KK,K00, KKold, I, J, Ntdc, INFK
      K00 = K0
      DO WHILE( (0 < K00).AND. (ABS(R(K00,K)).LE.LTOL) )
         R(K00,K) = gZERO
         K00      = K00 - 1
      ENDDO
      IF (K00.GT.0) THEN
      !  CDI(K) = COV(Xk Xj| X1,..,Xj-1,Xc  )/STD(Xj | X1,..,Xj-1,Xc)
         CDI(K) = R(K00,K)
         A(K)   = A(K)/CDI(K)
         B(K)   = B(K)/CDI(K)

         IF ((CDI(K)  <  gZERO).AND. INFI(K).GE. 0) THEN
            CALL SWAP(A(K),B(K))
            IF (INFI(K).NE. 2) INFI(K) = 1-INFI(K)
         END IF


                                !Scale conditional covariances
         R(1:K00,K) = R(1:K00,K)/CDI(K)
                                !Scale conditional standard dev.s used in regression eq.
         R(K,1:K00) = R(K,1:K00)/ABS(CDI(K))


         R(K00+1:K,K)   = gZERO
         !R(K,K00+1:K-1) = gZERO ! original
         R(K,K00:K-1) = gZERO    ! check this

         !
         if (chkLim.AND.K00>1) then
            ! Check if variable is redundant
            ! TODO:  this chklim-block does not work correctly yet
            xCut = mInfinity
            I = 1
            Ak = R(I,K)*xCut
            Bk = - (R(I,K))*xCut
            if (INFI(I)>=0) then
               if (INFI(I).ne.0) then
                  Ak = -(R(I,K))*MAX(A(I),-xCut)
               endif
               if (INFI(I).ne.1) then
                  Bk = - (R(I,K))*MIN(B(I),xCut)
               endif
            endif

            if (R(I,K)<gZERO) THEN
               CALL SWAP(Ak,Bk)
            endif
            !call printvar(infi(k),'infi(k)')
            !call printvar(A(k),'AK')
            !call printvar(B(k),'BK')
            !call printvar(Ak,'AK')
            !call printvar(Bk,'BK')
            INFK = INFI(K)
            Ak   = A(K)+Ak
            Bk   = B(K)+Bk
            D = gZERO
            DO I = 2, K00-1
               D = D + ABS(R(I,K))
            END DO
            CVDIAG = abs(R(k,k00))
            !call printvar(cvdiag,'cvdiag')
            Ak = (Ak + (D+cvdiag)*xCut)
            Bk = (Bk - (D+cvdiag)*xCut)
            !call printvar(Ak,'AK')
            !call printvar(Bk,'BK')
            ! If Ak<-xCut and xCut<Bk then variable Xk is redundant
            CALL ADJLIMITS(Ak,Bk,INFK)

! Should change this to check against A(k00) and B(k00)
            IF (INFK < 0) THEN
               !variable is redundnant
               !                     CALL mexPrintf('AdjLim'//CHAR(10))
               IF ( K < N1 ) THEN
                  CALL RCSWAP( K, N1, N1,N, R,INDEX1,Cm, A, B, INFI)

                  ! move conditional standarddeviations
                  R(K,1:K0) = R(N1,1:K0)
                  CDI(K)    = CDI(N1)

                  IF (PRESENT(Y)) THEN
                     Y(K) = Y(N1)
                  ENDIF
               ENDIF
               CDI(N1)    = gZERO
               R(1:N1,N1) = gZERO
               R(N1,1:N1) = gZERO

               INFIS = INFIS+1
               N1    = N1-1
              ! CALL printvar(index1(N1),'index1(n1)')
              ! CALL mexPrintf('RCSCALE: updated N1')
              ! CALL printvar(INFIS,'INFIS ')
               return
            END IF
         endif
         KKold = K
         KK = K-1
         DO I = K0, K00+1, -1
            DO WHILE ((I.LE.KK) .AND. ABS(R(I,KK)).GT.LTOL)
               DO J = 1,I       !K0
                  ! SWAP Covariance matrix
                  CALL SWAP(R(J,KK),R(J,KKold))
                  ! SWAP conditional standarddeviations
                  CALL SWAP(R(KK,J),R(KKold,J))
               END DO
               CALL SWAP(CDI(KK),CDI(KKold))
               CALL SWAP(Cm(KK),Cm(KKold))
               CALL SWAP(INDEX1(KK),INDEX1(KKold))
               CALL SWAP(A(KK),A(KKold))
               CALL SWAP(B(KK),B(KKold))
               CALL SWAP(INFI(KK),INFI(KKold))
               IF (PRESENT(Y)) THEN
                   CALL SWAP(Y(KK),Y(KKold))
               ENDIF
               Ntdc = SIZE(R,DIM=1)
               IF (N < Ntdc) THEN
                  ! SWAP Xc entries, i.e, Cov(Xt,Xc) and Cov(Xd,Xc)
                  DO J = N+1, Ntdc
                     CALL SWAP( R(KK,J), R(KKold,J) )
                  END DO
               ENDIF
               KKold = KK
               KK    = KK - 1
            ENDDO
         END DO
         IF (KK < K1) THEN
            K1 = K1 + 1
!            CALL mexPrintf('RCSCALE: updated K1'//CHAR(10))
         END IF
!         CALL PRINTVAR(K,TXT='K')
!         CALL PRINTVAR(KK,TXT='KK')
!         CALL PRINTVAR(K1,TXT='K1')
!         CALL PRINTVAR(K00,TXT='K00')
!         CALL PRINTVAR(K0,TXT='K0')
!         CALL PRINTCOF(N,A,B,INFI,R,INDEX1)
      ELSE
!     Remove variable if it is conditional independent of all other variables
!         CALL mexPrintf('RCSCALE ERROR*********************'//char(10))
!         call PRINTCOF(N,A,B,INFI,R,INDEX1)
!         CALL mexPrintf('RCSCALE ERROR*********************'//char(10))
      ENDIF
!      if (chkLim) then
!         call PRINTCOF(N,A,B,INFI,R,INDEX1)
!      endif
      END SUBROUTINE RCSCALE

      SUBROUTINE COVSRT(BCVSRT, Nt,Nd,R,Cm,A,B,INFI,INDEX1,
     &     INFIS,INFISD, NDIM, Y, CDI )
      USE FIMOD
      USE SWAPMOD
      USE GLOBALCONST
!      USE GLOBALDATA, ONLY : EPS2,NIT,xCutOff
      IMPLICIT NONE
!COVSRT  sort integration limits and determine Cholesky factor.
!
!     Nt, Nd = size info about Xt and Xd variables.
!     R      = Covariance/Cholesky factored matrix for [Xt,Xd,Xc] (in)
!              On input:
!              Note: Only upper triangular part is needed/used.
!               1a) the first upper triangular the Nt + Nd times Nt + Nd
!                   block contains COV([Xt,Xd]|Xc)
!                 (conditional covariance matrix for Xt and Xd given Xc)
!               2a) The upper triangular part of the Nt+Nd+Nc times Nc
!                   last block contains the cholesky matrix for Xc,
!                   i.e.,
!
!              On output:
!               1b) part 2a) mentioned above is unchanged, only necessary
!                 permutations according to INDEX1 is done.
!               2b) part 1a) mentioned above is changed to a special
!                  form of cholesky matrix: (N = Nt+Nd-INFIS-INFISD)
!   C = COVARIANCE
!   R(1,1) = 1
!   R(1,2:N) = [C(X1,X2)/STD(X1)/STD(X2|X1),..,C(X1,XN)/STD(X1)/STD(XN|XN-1,..,X1)]
!   R(2,2) = 1
!   R(2,3:N) =[C(X2,X3|X1)/STD(X2|X1)/STD(X3|X2,X1),..,C(X2,XN|X1)/STD(X2|X1)/STD(XN|XN-1,..,X1)]
!                  ....
!                  etc.
!               3b) The lower triangular part of R contains the
!               normalized conditional standard deviations (which is
!               used in the reqression approximation C1C2), i.e.,
!               R(2:N,1)  = [STD(X2|X1) STD(X3|X1),....,STD(XN|X1) ]/STD(X1)
!               R(3:N,2)  = [STD(X3|X1,X2),....,STD(XN|X1,X2) ]/STD(X2|X1)
!               .....
!               etc.
!     Cm     = Conditional mean given Xc
!     A,B    = lower and upper integration limits length Nt+Nd
!     INFIN  = INTEGER, array of integration limits flags:  length Nt+Nd   (in)
!             if INFIN(I) < 0, Ith limits are (-infinity, infinity);
!             if INFIN(I) = 0, Ith limits are (-infinity, B(I)];
!             if INFIN(I) = 1, Ith limits are [A(I), infinity);
!             if INFIN(I) = 2, Ith limits are [A(I), B(I)].
!    INDEX1  = permutation index vector, i.e., giving the indices to the
!              variables original place.
!    INFIS   = Number of redundant variables of Xt
!    INFISD  = Number of redundant variables of Xd
!    NDIM    = Number of relevant dimensions to integrate. This is the
!             same as the rank of the submatrix of Cov([Xt,Xd]) minus
!             the INFIS variables of Xt and INFISD variables of Xd.
!    Y       = working array
!    CDI     = Cholesky diagonal elements which contains either
!               CDI(J) = STD(Xj | X1,...,Xj-1,Xc) if Xj is stochastic given
!                        X1,...Xj, Xc
!              or
!               CDI(J) = COV(Xj,Xk | X1,..,Xk-1,Xc  )/STD(Xk | X1,..,Xk-1,Xc)
!               if Xj is determinstically determined given X1,..,Xk,Xc
!               for some k<j.
!
!     Subroutine to sort integration limits and determine Cholesky
!     factor.
!
!     Note: COVSRT works only on the upper triangular part of R
      LOGICAL,                          INTENT(in)    :: BCVSRT
      INTEGER,                          INTENT(in)    :: Nt,Nd
      DOUBLE PRECISION, DIMENSION(:,:), INTENT(inout) :: R
      DOUBLE PRECISION, DIMENSION(:  ), INTENT(inout) :: Cm,A,B
      INTEGER,          DIMENSION(:  ), INTENT(inout) :: INFI
      INTEGER,          DIMENSION(:  ), INTENT(inout) :: INDEX1
      INTEGER,                          INTENT(out) :: INFIS,INFISD,NDIM
      DOUBLE PRECISION, DIMENSION(:  ), INTENT(out)   :: Y, CDI
      double precision :: covErr
!     Local variables
      INTEGER :: N,N1,I, J, K, L, JMIN, Ndleft
      INTEGER ::  K1, K0, Nullity,INFJ,FINA,FINB
      DOUBLE PRECISION :: SUMSQ, AJ, BJ, TMP,D, E,EPSL
      DOUBLE PRECISION :: AA, Ca, Pa, APJ, PRBJ,RMAX
      DOUBLE PRECISION :: CVDIAG, AMIN, BMIN, PRBMIN
      DOUBLE PRECISION :: LTOL,TOL, ZERO,xCut
      LOGICAL          :: isOK = .TRUE.
      LOGICAL          :: isXd = .FALSE.
      LOGICAL          :: isXt = .FALSE.
      LOGICAL          :: chkLim
!      PARAMETER ( SQTWPI = 2.506628274631001D0, )
      PARAMETER (ZERO = 0.D0, TOL = 1D-16)

      xCut = mInfinity
!     xCut = MIN(ABS(xCutOff),8.0D0)
      INFIS  = 0
      INFISD = 0
      Ndim   = 0
      Ndleft = Nd
      N      = Nt + Nd

      LTOL = mSmall
      TMP  = ZERO
      DO I = 1, N
         IF (R(I,I).GT.TMP) TMP = R(I,I)
      ENDDO
      EPSL = tmp*MAX(mCovEps,N*mSmall) !tmp
      !IF (N < 10) EPSL = MIN(1D-10,EPSL)
      !LTOL = EPSL
!      EPSL = MAX(EPS2,LTOL)
      IF (TMP.GT.EPSL) THEN
         DO I = 1, N
            IF ((INFI(I)  <  0).OR.(R(I,I)<=LTOL)) THEN
               IF (INDEX1(I)<=Nt)  THEN
                  INFIS = INFIS+1
               ELSEIF (R(I,I)<=LTOL) THEN
                  INFISD = INFISD+1
               ENDIF
            ENDIF
         END DO
      ELSE
         LTOL       = EPSL
         INFIS      = Nt
         INFISD     = Nd
         R(1:N,1:N) = ZERO
      ENDIF

      covErr = 20.d0*LTOL

      N1  = N-INFIS-INFISD
      CDI(N1+1:N) = gZERO
      !PRINT *,'COVSRT'
      !CALL PRINTCOF(N,A,B,INFI,R,INDEX1)

!     Move any redundant variables of Xd to innermost positions.
      LP3: DO I = N, N-INFISD+1, -1
         isXt = (INDEX1(I)<=Nt)
         IF ( (R(I,I) > LTOL) .OR. (isXt)) THEN
            DO J = 1,I-1
               isXd = (INDEX1(J)>Nt)
               IF ( (R(J,J) <= LTOL) .AND.isXd) THEN
                  CALL RCSWAP(J, I, N, N, R,INDEX1,Cm, A, B, INFI)
                  !GO TO 10
                  CYCLE LP3
               ENDIF
            END DO
         ENDIF
! 10
      END DO LP3
!
!     Move any doubly infinite limits or any redundant of Xt to the next
!     innermost positions.
!
      LP4: DO I = N-INFISD, N1+1, -1
         isXd = (INDEX1(I)>Nt)
         IF ( ((INFI(I) > -1).AND.(R(I,I) > LTOL))
     &        .OR. isXd) THEN
            DO J = 1,I-1
               isXt = (INDEX1(J)<=Nt)
               IF ( (INFI(J)  <  0 .OR. (R(J,J)<= LTOL))
     &              .AND. (isXt)) THEN
                  CALL RCSWAP( J, I, N,N, R,INDEX1,Cm, A, B, INFI)
                  !GO TO 15
                  CYCLE LP4
               ENDIF
            END DO
         ENDIF
!15
      END DO LP4

!      CALL mexprintf('Before sorting')
!      CALL PRINTCOF(N,A,B,INFI,R,INDEX1)
!      CALL PRINTVEC(CDI,'CDI')
!      CALL PRINTVEC(Cm,'Cm')

      IF ( N1 <= 0 ) GOTO 200
!
!     Sort remaining limits and determine Cholesky factor.
!
      Y(1:N1) = gZERO
      K       = 1
      Ndleft  = Nd - INFISD
      Nullity = 0
      DO  WHILE (K .LE. N1)

!     IF (Ndim.EQ.3) EPSL = MAX(EPS2,1D-10)
!     Determine the integration limits for variable with minimum
!     expected probability and interchange that variable with Kth.

         K0     = K - Nullity
         PRBMIN = gTWO
         JMIN   = K
         CVDIAG = ZERO
         RMAX   = ZERO
         IF ((Ndleft>0) .OR. (NDIM < Nd+mNIT)) THEN
            DO J = K, N1
               isXd = (INDEX1(J)>Nt)
               isOK = ((NDIM <= Nd+mNIT).OR.isXd)
               IF ( R(J,J) <= K0*K0*EPSL .OR. (.NOT. isOK)) THEN
                  RMAX = max(RMAX,ABS(R(J,J)))
               ELSE
                  TMP = ZERO    ! =  conditional mean of Y(I) given Y(1:I-1)
                  DO I = 1, K0 - 1
                     TMP = TMP + R(I,J)*Y(I)
                  END DO
                  SUMSQ = SQRT( R(J,J))

                  IF (INFI(J)>-1) THEN
                                ! May have infinite int. limits if Nd>0
                     IF (INFI(J).NE.0) THEN
                        AJ = ( A(J) - TMP )/SUMSQ
                     ENDIF
                     IF (INFI(J).NE.1) THEN
                        BJ = ( B(J) - TMP )/SUMSQ
                     ENDIF
                  ENDIF
                  IF (isXd) THEN
                     AA = (Cm(J)+TMP)/SUMSQ ! inflection point
                     CALL EXLMS(AA,AJ,BJ,INFI(J),D,E,Ca,Pa)
                     PRBJ = E - D
                  ELSE
                                !CALL MVNLMS( AJ, BJ, INFI(J), D, E )
                     CALL MVNLIMITS(AJ,BJ,INFI(J),APJ,PRBJ)
                  ENDIF
                                !IF ( EMIN + D .GE. E + DMIN ) THEN
                  IF ( PRBJ  <  PRBMIN ) THEN
                     JMIN = J
                     AMIN = AJ
                     BMIN = BJ
                     PRBMIN = MAX(PRBJ,ZERO)
                     CVDIAG = SUMSQ
                  ENDIF
               ENDIF
            END DO
         END IF
!
!     Compute Ith column of Cholesky factor.
!     Compute expected value for Ith integration variable (without
!     considering the jacobian) and
!     scale Ith covariance matrix row and limits.
!
! 40
         IF ( CVDIAG.GT.TOL) THEN
            isXd = (INDEX1(JMIN)>Nt)
            IF (isXd) THEN
               Ndleft = Ndleft - 1
            ELSEIF (BCVSRT.EQV..FALSE..AND.(PRBMIN+LTOL>=gONE)) THEN
!BCVSRT.EQ.
               J = 1
               AJ = R(J,JMIN)*xCut
               BJ = - (R(J,JMIN))*xCut
               if (INFI(J)>=0) then
                  if (INFI(J).ne.0) then
                     AJ = -(R(J,JMIN))*MAX(A(J),-xCut)
                  endif
                  if (INFI(J).ne.1) then
                     BJ = - (R(J,JMIN))*MIN(B(J),xCut)
                  endif
               endif
               if (R(J,JMIN)<gZERO) THEN
                  CALL SWAP(AJ,BJ)
               endif
               INFJ = INFI(JMIN)
               AJ   = A(JMIN)+AJ
               BJ   = B(JMIN)+BJ

               D = gZERO
               DO J = 2, K0-1
                  D = D + ABS(R(J,JMIN))
               END DO

               AJ = (AJ + D*xCut)/CVDIAG
               BJ = (BJ - D*xCut)/CVDIAG
               CALL ADJLIMITS(AJ,BJ,INFJ)
               IF (INFJ < 0) THEN
                  !variable is redundnant
                  !                     CALL mexPrintf('AdjLim'//CHAR(10))
                  IF ( JMIN < N1 ) THEN
                   CALL RCSWAP( JMIN,N1,N1,N,R,INDEX1,Cm,A,B,INFI)
                     ! move conditional standarddeviations
                     R(JMIN,1:K0-1) = R(N1,1:K0-1)

                     Y(JMIN) = Y(N1)
                  ENDIF
                  R(1:N1,N1)     = gZERO
                  R(N1,1:N1)     = gZERO
                  Y(N1)   = gZERO
                  INFIS = INFIS+1
                  N1    = N1-1
                  GOTO 100
               END IF
            ENDIF
            NDIM = NDIM + 1     !Number of relevant dimensions to integrate

            IF ( K < JMIN ) THEN
               CALL RCSWAP( K, JMIN, N1,N, R,INDEX1,Cm, A, B, INFI)
               ! SWAP conditional standarddeviations
               DO J = 1,K0-1  !MIN(K0, K-1)
                  CALL SWAP(R(K,J),R(JMIN,J))
               END DO
            END IF

            R(K0,K) = CVDIAG
            CDI(K)  = CVDIAG     ! Store the diagonal element
            DO I = K0+1,K
               R(I,K) = gZERO;
               R(K,I) = gZERO
            END DO

            K1 = K
            I  = K1 + 1
            DO WHILE (I <= N1)
               TMP = ZERO
               DO J = 1, K0 - 1
                  !tmp = tmp + L(i,j).*L(k1,j)
                  TMP = TMP + R(J,I)*R(J,K1)
               END DO
                  ! Cov(Xk,Xi|X1,X2,...Xk-1)/STD(Xk|X1,X2,...Xk-1)
               R(K0,I)  = (R(K1,I) - TMP) /CVDIAG
                  ! Var(Xi|X1,X2,...Xk)
               R(I,I) = R(I,I) - R(K0,I) * R(K0,I)

               IF (R(I,I).GT.LTOL) THEN
                  R(I,K0) = SQRT(R(I,I)) ! STD(Xi|X1,X2,...Xk)
               ELSE   !!IF (R(I,I) .LE. LTOL) THEN !TOL
                                !CALL mexprintf('Singular')
                  isXd = (index1(I)>Nt)
                  if (isXd) then
                     Ndleft = Ndleft - 1
                  ELSEIF (BCVSRT.EQV..FALSE.) THEN
!     BCVSRT.EQ.
                     J = 1
                     AJ = R(J,I)*xCut
                     BJ = - (R(J,I))*xCut
                     if (INFI(J)>=0) then
                        if (INFI(J).ne.0) then
                           AJ = -(R(J,I))*MAX(A(J),-xCut)
                        endif
                        if (INFI(J).ne.1) then
                           BJ = - (R(J,I))*MIN(B(J),xCut)
                        endif
                     endif
                     if (R(J,I)<gZERO) THEN
                        CALL SWAP(AJ,BJ)
                     endif
                     INFJ = INFI(I)
                     AJ   = A(I)+AJ
                     BJ   = B(I)+BJ

                     D = gZERO
                     DO J = 2, K0
                        D = D + ABS(R(J,I))
                     END DO

                     AJ = (AJ + D*xCut)-mXcutOff
                     BJ = (BJ - D*xCut)+mXcutOff
                     !call printvar(Aj,'Aj')
                     !call printvar(Bj,'Bj')
                     CALL ADJLIMITS(AJ,BJ,INFJ)
                     IF (INFJ < 0) THEN
                                !variable is redundnant
                        !CALL mexPrintf('AdjLim'//CHAR(10))
                        IF ( I < N1 ) THEN
                           CALL RCSWAP( I,N1,N1,N,R,INDEX1,Cm,A,B,INFI)
                                ! move conditional standarddeviations
                           R(I,1:K0-1) = R(N1,1:K0-1)

                           Y(I) = Y(N1)
                        ENDIF
                        R(1:N1,N1)     = gZERO
                        R(N1,1:N1)     = gZERO
                        Y(N1)   = gZERO
                        INFIS = INFIS+1
                        N1    = N1-1

                        !CALL mexprintf('covsrt updated N1')
                        !call printvar(INFIS,' Infis')
                        GOTO 75
                     END IF
                  END IF
                  IF (mNIT>100) THEN
                     R(I,K0) = gZERO
                  ELSE
                     R(I,K0) = MAX(SQRT(MAX(R(I,I), gZERO)),LTOL)
                  ENDIF
                  Nullity = Nullity + 1
                  K  = K + 1
                  IF (K  <  I) THEN
                     CALL RCSWAP( K, I, N1,N,R,INDEX1,Cm, A, B, INFI)
                     ! SWAP conditional standarddeviations
                     DO J = 1, K0
                        CALL SWAP(R(K,J),R(I,J))
                     END DO
                  ENDIF
                  chkLim = .FALSE. !((.not.isXd).AND.(BCVSRT.EQ..FALSE.))
                  L = INFIS
                  CALL  RCSCALE(chkLim,K,K0,N1,N,K1,INFIS,CDI,Cm,
     &                 R,A,B,INFI,INDEX1)
                  if (L.ne.INFIS) THEN
                     K = K - 1
                     I = I - 1
                  ENDIF
               END IF
               I = I + 1
 75            CONTINUE
            END DO
            INFJ = INFI(K1)

            IF (K1 .EQ.1) THEN
               FINA = 0
               FINB = 0
               IF (INFJ.GE.0) THEN
                  IF  (INFJ.NE.0) FINA = 1
                  IF  (INFJ.NE.1) FINB = 1
               ENDIF
               CALL C1C2(K1+1,N1,K0,A,B,INFI, Y, R,
     &              AMIN, BMIN, FINA,FINB)
               INFJ = 2*FINA+FINB-1
               CALL MVNLIMITS(AMIN,BMIN,INFJ,APJ,PRBMIN)
            ENDIF

            Y(K0) = gettmean(AMIN,BMIN,INFJ,PRBMIN)


            R( K0, K1 ) = R( K0, K1 ) / CVDIAG
            DO J = 1, K0 - 1
               ! conditional covariances
               R( J, K1 ) = R( J, K1 ) / CVDIAG
               ! conditional standard dev.s used in regression eq.
               R( K1, J ) = R( K1, J ) / CVDIAG
            END DO

            A( K1 ) = A( K1 )/CVDIAG
            B( K1 ) = B( K1 )/CVDIAG

            K  = K  + 1
100         CONTINUE
         ELSE
            covErr = RMAX
            R(K:N1,K:N1) = gZERO
            I = K
            DO WHILE (I <= N1)
!  Scale  covariance matrix rows and limits
!  If the conditional covariance matrix diagonal entry is zero,
!  permute limits and/or rows, if necessary.
               chkLim = ((index1(I)<=Nt).AND.(BCVSRT.EQV..FALSE.))
               L = INFIS
               CALL RCSCALE(chkLim,I,K0-1,N1,N,K1,INFIS,CDI,Cm,
     &              R,A,B,INFI,INDEX1)
               if (L.EQ.INFIS) I = I + 1
            END DO
            Nullity = N1 - K0 + 1
            GOTO 200  !RETURN
         END IF
      END DO
 200  CONTINUE
      IF (Ndim .GT. 0) THEN  ! N1<K
         ! K1 = index to the last stochastic varible to integrate
         ! If last stoch. variable is Xt: reduce dimension of integral by 1
         IF (ALL(INDEX1(K1:N1).LE.Nt)) Ndim = Ndim-1
      ENDIF
!      CALL mexprintf('After sorting')

!      CALL PRINTCOF(N,A,B,INFI,R,INDEX1)
!      CALL printvar(A(1),'A1')
!      CALL printvar(B(1),'B1')
!      CALL printvar(INFIS,'INFIS')
!      CALL PRINTVEC(CDI,'CDI')
!      CALL PRINTVEC(Y,'Y')
!      CALL PRINTVEC(AA1,'AA1')
!      CALL PRINTVEC(BB1,'BB1')
!      CALL PRINTVAR(NDIM,TXT='NDIM')
!      CALL PRINTVAR(NIT,TXT='NIT')
!      DEALLOCATE(AA1)
!      DEALLOCATE(BB1)
      RETURN
      END SUBROUTINE COVSRT

      SUBROUTINE COVSRT1(BCVSRT, Nt,Nd,R,Cm,A,B,INFI,INDEX1,
     &     INFIS,INFISD, NDIM, Y, CDI )
      USE FIMOD
      USE SWAPMOD
!      USE GLOBALCONST
!      USE GLOBALDATA, ONLY : EPS2,NIT,xCutOff,Nc1c2
      IMPLICIT NONE
!COVSRT1  sort integration limits and determine Cholesky factor.
!
!     Nt, Nd = size info about Xt and Xd variables.
!     R      = Covariance/Cholesky factored matrix for [Xt,Xd,Xc] (in)
!              On input:
!              Note: Only upper triangular part is needed/used.
!               1a) the first upper triangular the Nt + Nd times Nt + Nd
!                   block contains COV([Xt,Xd]|Xc)
!                   (conditional covariance matrix for Xt and Xd given Xc)
!               2a) The upper triangular part of the Nt+Nd+Nc times Nc
!                   last block contains the cholesky matrix for Xc, i.e.,
!
!              On output:
!               1b) part 2a) mentioned above is unchanged, only necessary
!                   permutations according to INDEX1 is done.
!               2b) part 1a) mentioned above is changed to a special
!                   form of cholesky matrix: (N = Nt+Nd-INFIS-INFISD)
!                  R(1,1) = 1
!                  R(1,2:N) = [COV(X1,X2)/STD(X1),....COV(X1,XN)/STD(X1)]
!                  R(2,2) = 1
!                  R(2,3:N) = [COV(X2,X3)/STD(X2|X1),....COV(X2,XN)/STD(X2|X1)]
!                  ....
!                  etc.
!               3b) The lower triangular part of R contains the
!               conditional standard deviations, i.e.,
!               R(2:N,1)  = [STD(X2|X1) STD(X3|X1),....,STD(XN|X1) ]
!               R(3:N,2)  = [STD(X3|X1,X2),....,STD(XN|X1,X2) ]
!               .....
!               etc.
!
!     Cm     = Conditional mean given Xc
!     A,B    = lower and upper integration limits length Nt+Nd
!     INFIN  = INTEGER, array of integration limits flags:  length Nt+Nd   (in)
!             if INFIN(I) < 0, Ith limits are (-infinity, infinity);
!             if INFIN(I) = 0, Ith limits are (-infinity, B(I)];
!             if INFIN(I) = 1, Ith limits are [A(I), infinity);
!             if INFIN(I) = 2, Ith limits are [A(I), B(I)].
!    INDEX1  = permutation index vector
!    INFIS   = Number of redundant variables of Xt
!    INFISD  = Number of redundant variables of Xd
!    NDIM    = Number of relevant dimensions to integrate. This is the
!             same as the rank of the submatrix of Cov([Xt,Xd]) minus
!             the INFIS variables of Xt and INFISD variables of Xd.
!    Y       = working array
!    CDI     = Cholesky diagonal elements which contains either
!               CDI(J) = STD(Xj| X1,,,Xj-1,Xc) if Xj is stochastic given
!                X1,...Xj, Xc
!              or
!               CDI(J) = COV(Xj,Xk|X1,..,Xk-1,Xc  )/STD(Xk| X1,,,Xk-1,Xc)
!               if Xj is determinstically determined given X1,..,Xk,Xc
!               for some k<j.
!
!     Subroutine to sort integration limits and determine Cholesky
!     factor.
!
!     Note: COVSRT1 works only on the upper triangular part of R
      LOGICAL,                          INTENT(in)    :: BCVSRT
      INTEGER,                          INTENT(in)    :: Nt,Nd
      DOUBLE PRECISION, DIMENSION(:,:), INTENT(inout) :: R
      DOUBLE PRECISION, DIMENSION(:  ), INTENT(inout) :: Cm,A,B
      INTEGER,          DIMENSION(:  ), INTENT(inout) :: INFI
      INTEGER,          DIMENSION(:  ), INTENT(inout) :: INDEX1
      INTEGER,                          INTENT(out) :: INFIS,INFISD,NDIM
      DOUBLE PRECISION, DIMENSION(:  ), INTENT(out)   :: Y, CDI
!     Local variables
      INTEGER :: N,N1,I, J, K, L, JMIN,Ndleft
      INTEGER ::  K1, K0, Nullity, INFJ, FINA, FINB
      DOUBLE PRECISION :: SUMSQ, AJ, BJ, TMP,D, E,EPSL
      DOUBLE PRECISION :: AA, Ca, Pa, APJ, PRBJ, RMAX, xCut
      DOUBLE PRECISION :: CVDIAG, AMIN, BMIN, PRBMIN
      DOUBLE PRECISION :: LTOL,TOL, ZERO,EIGHT
!      INTEGER, PARAMETER :: NMAX = 1500
!      DOUBLE PRECISION, DIMENSION(NMAX) :: AP,BP
!      INTEGER, DIMENSION(NMAX) :: INFP
!      INTEGER :: Nabp
      LOGICAL          :: isOK = .TRUE.
      LOGICAL          :: isXd = .FALSE.
      LOGICAL          :: isXt = .FALSE.
      LOGICAL          :: chkLim
!      PARAMETER ( SQTWPI = 2.506628274631001D0)
      PARAMETER (ZERO = 0.D0,EIGHT = 8.D0, TOL=1.0D-16)

      xCut = MIN(ABS(mXcutOff),EIGHT)
      INFIS  = 0
      INFISD = 0
      Ndim   = 0
      Ndleft = Nd
      N      = Nt+Nd

!     IF (N < 10) EPSL = MIN(1D-10,EPS2)
      LTOL = TOL
      TMP  = ZERO
      DO I = 1, N
         IF (R(I,I).GT.TMP) TMP = R(I,I)
      ENDDO

      EPSL = MAX(mCovEps,N*TMP*mSmall)
      IF (TMP.GT.EPSL) THEN
         DO I = 1, N
            IF ((INFI(I)  <  0).OR.(R(I,I).LE.LTOL)) THEN
               IF (INDEX1(I).LE.Nt)  THEN
                  INFIS = INFIS+1
               ELSEIF (R(I,I).LE.LTOL) THEN
                  INFISD = INFISD+1
               ENDIF
            ENDIF
         END DO
      ELSE
         !CALL PRINTCOF(N,A,B,INFI,R,INDEX1)
         !CALL PRINTVEC(CDI)
         !CALL PRINTVEC(Cm)
         LTOL   = EPSL
         INFIS  = Nt
         INFISD = Nd
         R(1:N,1:N) = ZERO
      ENDIF
      N1  = N-INFIS-INFISD
      CDI(N1+1:N) = ZERO
!     PRINT *,'COVSRT'
!     CALL PRINTCOF(N,A,B,INFI,R,INDEX1)

!     Move any redundant variables of Xd to innermost positions.
      LP1:  DO I = N, N-INFISD+1, -1
         isXt = (INDEX1(I).LE.Nt)
         IF ( R(I,I) .GT. LTOL .OR. isXt) THEN
            DO J = 1,I-1
               isXd = (INDEX1(J).GT.Nt)
               IF ( R(J,J) .LE. LTOL .AND. isXd) THEN
                  CALL RCSWAP( J, I, N,N, R,INDEX1,Cm, A, B, INFI)
                  !GO TO 10
                  CYCLE LP1
               ENDIF
            END DO
         ENDIF
! 10
      END DO LP1
!
!     Move any doubly infinite limits or any redundant of Xt to the next
!     innermost positions.
!
      LP2: DO I = N-INFISD, N1+1, -1
         isXd = (INDEX1(I).GT.Nt)
         IF ( ((INFI(I) .GE. 0).AND. (R(I,I).GT. LTOL) )
     &        .OR. isXd) THEN
            DO J = 1,I-1
               isXt = (INDEX1(J).LE.Nt)
               IF ( (INFI(J)  <  0 .OR. (R(J,J).LE. LTOL))
     &              .AND. isXt) THEN
                  CALL RCSWAP( J, I, N,N, R,INDEX1,Cm, A, B, INFI)
                  !GO TO 15
                  CYCLE LP2
               ENDIF
            END DO
         ENDIF
 !15
      END DO LP2

      IF ( N1 .LE. 0 ) RETURN
!      CALL mexprintf('Before sorting')
!      CALL PRINTCOF(N,A,B,INFI,R,INDEX1)


!     Sort remaining limits and determine Cholesky factor.
      Y(1:N1) = ZERO
      K  = 1
!      N1  = N-INFIS-INFISD
      Ndleft = Nd - INFISD
      Nullity = 0

!      Nabp  = 0
!      AP(1:N1) = ZERO
!      BP(1:N1) = zero
      DO  WHILE (K .LE. N1)

!     Determine the integration limits for variable with minimum
!     expected probability and interchange that variable with Kth.

         K0     = K-Nullity
         PRBMIN = 2.d0
         JMIN   = K
         CVDIAG = ZERO
         RMAX   = ZERO
         IF (Ndleft.GT.0 .OR. NDIM < Nd+mNIT) THEN
            DO J = K,N1
               isXd = (INDEX1(J).GT.Nt)
               isOK = ((NDIM <= Nd+mNIT).OR.isXd)
               IF ( R(J,J) .LE. K0*K0*EPSL.OR. (.NOT. isOK)) THEN
                  RMAX = max(RMAX,R(J,J))
               ELSE
                  TMP = Y(J) ! =  the conditional mean of Y(J) given Y(1:J-1)
                  SUMSQ = SQRT( R(J,J))

                  IF (INFI(J) < 0) GO TO 30 ! May have infinite int. limits if Nd>0
                  IF (INFI(J).NE.0) THEN
                     AJ = ( A(J) - TMP )/SUMSQ
                  ENDIF
                  IF (INFI(J).NE.1) THEN
                     BJ = ( B(J) - TMP )/SUMSQ
                  ENDIF
 30               IF (INDEX1(J).GT.Nt) THEN
                     AA = (Cm(J)+TMP)/SUMSQ ! inflection point
                     CALL EXLMS(AA,AJ,BJ,INFI(J),D,E,Ca,Pa)
                     PRBJ = E-D
                  ELSE
                                !CALL MVNLMS( AJ, BJ, INFI(J), D, E )
                     CALL MVNLIMITS(AJ,BJ,INFI(J),APJ,PRBJ)
                  ENDIF
                  IF ( PRBJ  <  PRBMIN ) THEN
                     JMIN = J
                     AMIN = AJ
                     BMIN = BJ
                     PRBMIN = MAX(PRBJ,ZERO)
                     CVDIAG = SUMSQ
                  ENDIF
               ENDIF
            END DO
         END IF
!
!     Compute Ith column of Cholesky factor.
!     Compute expected value for Ith integration variable (without
!     considering the jacobian) and
!     scale Ith covariance matrix row and limits.
!
!40
         IF ( CVDIAG.GT.TOL) THEN
            IF (INDEX1(JMIN).GT.Nt) THEN
               Ndleft = Ndleft-1
            ELSE
               IF (BCVSRT.EQV..FALSE..AND.(PRBMIN+LTOL.GE.gONE)) THEN
!BCVSRT.EQ.
                  I = 1
                  AJ = R(I,JMIN)*xCut
                  BJ = - (R(I,JMIN))*xCut
                  if (INFI(1)>=0) then
                     if (INFI(1).ne.0) then
                        AJ = -(R(I,JMIN))*MAX(A(I),-xCut)
                     endif
                     if (INFI(1).ne.1) then
                        BJ = - (R(I,JMIN))*MIN(B(I),xCut)
                     endif
                  endif
                  if (R(I,JMIN)<gZERO) THEN
                     CALL SWAP(AJ,BJ)
                  endif
                  INFJ = INFI(JMIN)
                  AJ   = A(JMIN)+AJ
                  BJ   = B(JMIN)+BJ

                  D = gZERO
                  DO I = 2, K0-1
                     D = D + ABS(R(I,JMIN))
                  END DO



                  AJ = (AJ + D*xCut)/CVDIAG
                  BJ = (BJ - D*xCut)/CVDIAG
                  CALL ADJLIMITS(AJ,BJ,INFJ)
                  IF (INFJ < 0) THEN
                     !variable is redundnant
                     !CALL mexPrintf('AdjLim'//CHAR(10))
                   IF ( JMIN < N1 ) THEN
                   CALL RCSWAP( JMIN, N1, N1,N, R,INDEX1,Cm, A, B, INFI)
                                ! SWAP conditional standarddeviations
                      DO I = 1,K0-1
                         CALL SWAP(R(JMIN,I),R(N1,I))
                      END DO
                      CALL SWAP(Y(N1),Y(JMIN))
                   ENDIF
                   INFIS = INFIS+1
                   N1    = N1-1
                   GOTO 100
                  END IF
                ENDIF
            ENDIF
            NDIM = NDIM + 1     !Number of relevant dimensions to integrate

            IF ( K < JMIN ) THEN

               CALL RCSWAP( K, JMIN, N1,N, R,INDEX1,Cm, A, B, INFI)
               ! SWAP conditional standarddeviations
               DO J=1,K0-1
                  CALL SWAP(R(K,J),R(JMIN,J))
               END DO
               CALL SWAP(Y(K),Y(JMIN))
            END IF


            R(K0,K:N1) = R(K0,K:N1)/CVDIAG
            R(K0,K) = CVDIAG
            CDI(K)  = CVDIAG     ! Store the diagonal element
            DO I = K0+1,K
               R(I,K) = ZERO
               R(K,I) = ZERO
            END DO

            K1  = K
            !IF (K .EQ. N1) GOTO 200

!  Cov(Xi,Xj|Xk,Xk+1,..,Xn)=
!             Cov(Xi,Xj|Xk+1,..,Xn) -
 !             Cov(Xi,Xk|Xk+1,..Xn)*Cov(Xj,Xk|Xk+1,..Xn)
            I = K1 +1
            DO WHILE (I <= N1)
                ! Var(Xj | Xk,Xk+1,...,Xn)
               R(I,I)  =  R(I,I) - R(K0,I)*R(K0,I)
               IF (R(I,I).GT.LTOL) THEN
                  R(I,K0)     = SQRT(R(I,I)) ! STD(Xi|X1,X2,...Xk)
                  R(I,I+1:N1) = R(I,I+1:N1) - R(K0,I+1:N1)*R(K0,I)
               ELSE
                  R(I,K0) = MAX(SQRT(MAX(R(I,I), gZERO)),LTOL)
                  Nullity = Nullity + 1
                  K  = K + 1
                  IF (K  <  I) THEN
                     CALL RCSWAP( K, I, N1,N,R,INDEX1,Cm, A, B, INFI)
                     ! SWAP conditional standarddeviations
                     DO J=1,K0
                        CALL SWAP(R(K,J),R(I,J))
                     END DO
                     CALL SWAP(Y(K),Y(I))
                  ENDIF
                  isXd = (INDEX1(K).GT.Nt)
                  IF (isXd) Ndleft = Ndleft-1
                  chkLim = ((.not.isXd).AND.(BCVSRT.EQV..FALSE.))
                  L = INFIS
                  CALL  RCSCALE(chkLim,K,K0,N1,N,K1,INFIS,CDI,Cm,
     &                 R,A,B,INFI,INDEX1,Y)
                  IF (L.NE.INFIS) I = I - 1
               END IF
               I = I +1
            END DO
            INFJ = INFI(K1)
            IF (K0 == 1) THEN
               FINA = 0
               FINB = 0
               IF (INFJ.GE.0) THEN
                  IF  (INFJ.NE.0) FINA = 1
                  IF  (INFJ.NE.1) FINB = 1
               ENDIF
               CALL C1C2(K1+1,N1,K0,A,B,INFI, Y, R,
     &              AMIN, BMIN, FINA,FINB)
               INFJ = 2*FINA+FINB-1
               CALL MVNLIMITS(AMIN,BMIN,INFJ,APJ,PRBMIN)
            ENDIF
            Y(K0) = GETTMEAN(AMIN,BMIN,INFJ,PRBMIN)

                     ! conditional mean (expectation)
                                ! E(Y(K+1:N)|Y(1),Y(2),...,Y(K))
            Y(K+1:N1) = Y(K+1:N1)+Y(K0)*R(K0,K+1:N1)
            R(K0,K1)  = R(K0,K1)/CVDIAG ! conditional covariances
            DO J = 1, K0 - 1
               R(J,K1) = R(J,K1)/CVDIAG ! conditional covariances
               R(K1,J) = R(K1,J)/CVDIAG ! conditional standard dev.s used in regression eq.
            END DO
            A(K1) = A(K1)/CVDIAG
            B(K1) = B(K1)/CVDIAG

            K  = K  + 1
 100        CONTINUE
         ELSE
            R(K:N1,K:N1) = gZERO
!            CALL PRINTCOF(N,A,B,INFI,R,INDEX1)
            I = K
            DO WHILE (I <= N1)
!  Scale  covariance matrix rows and limits
!  If the conditional covariance matrix diagonal entry is zero,
!  permute limits and/or rows, if necessary.
               chkLim = ((index1(I)<=Nt).AND.(BCVSRT.EQV..FALSE.))
               L = INFIS
               CALL RCSCALE(chkLim,I,K0-1,N1,N,K1,INFIS,CDI,Cm,
     &              R,A,B,INFI,INDEX1)
               if (L.EQ.INFIS) I = I + 1
            END DO
            Nullity = N1 - K0 + 1
            GOTO 200  !RETURN
         END IF
      END DO

 200  CONTINUE
      IF (Ndim .GT. 0) THEN     ! N1<K
         ! K1 = index to the last stochastic varible to integrate
         IF (ALL(INDEX1(K1:N1).LE.Nt)) Ndim = Ndim - 1
      ENDIF
!      CALL mexprintf('After sorting')
!      CALL PRINTCOF(N,A,B,INFI,R,INDEX1)
!      CALL PRINTVEC(CDI)
!      CALL PRINTVAR(NDIM,TXT='NDIM')
      RETURN
      END SUBROUTINE COVSRT1

      SUBROUTINE RCSWAP( P, Q, N,Ntd, C,IND,Cm, A, B, INFIN )
      USE SWAPMOD
      IMPLICIT NONE
! RCSWAP  Swaps rows and columns P and Q in situ, with P <= Q.
!
!
!   CALL  RCSWAP( P, Q, N, Ntd, C,IND A, B, INFIN, Cm)
!
!    P, Q  = row/column number to swap P<=Q<=N
!    N     = Number of significant variables of [Xt,Xd]
!    Ntd   = length(Xt)+length(Xd)
!    C     = upper triangular cholesky factor.Cov([Xt,Xd,Xc]) size Ntdc X Ntdc
!    IND   = permutation index vector. (index to variables original place).
!    Cm    = conditional mean
!    A,B   = lower and upper integration limit, respectively.
!    INFIN = if INFIN(I) < 0, Ith limits are (-infinity, infinity);
!            if INFIN(I) = 0, Ith limits are (-infinity, B(I)];
!            if INFIN(I) = 1, Ith limits are [A(I), infinity);
!            if INFIN(I) = 2, Ith limits are [A(I), B(I)].
!
! NOTE: RCSWAP works only on the upper triangular part of C
      DOUBLE PRECISION, DIMENSION(:,:),INTENT(inout) :: C
      INTEGER, DIMENSION(:),INTENT(inout) :: IND
      INTEGER, DIMENSION(:),          OPTIONAL,INTENT(inout) :: INFIN
      DOUBLE PRECISION, DIMENSION(:), OPTIONAL,INTENT(inout) :: A,B,Cm
      INTEGER,INTENT(in) :: P, Q, N, Ntd
! local variable
      INTEGER :: J, Ntdc
      LOGICAL :: isXc
      IF (PRESENT(Cm))    CALL SWAP( Cm(P), Cm(Q) )
      IF (PRESENT(A))     CALL SWAP( A(P), A(Q) )
      IF (PRESENT(B))     CALL SWAP( B(P), B(Q) )
      IF (PRESENT(INFIN)) CALL SWAP(INFIN(P),INFIN(Q))

      CALL SWAP(IND(P),IND(Q))

      CALL SWAP( C(P,P), C(Q,Q) )
      DO J = 1, P-1
         CALL SWAP( C(J,P), C(J,Q) )
      END DO
      DO J = P+1, Q-1
         CALL SWAP( C(P,J), C(J,Q) )
      END DO
      DO J = Q+1, N
         CALL SWAP( C(P,J), C(Q,J) )
      END DO
      Ntdc = SIZE(C,DIM=1)
      isXc = (N < Ntdc)
      IF (isXc) THEN
         !Swap row P and Q of Xc variables
         DO J = Ntd+1, Ntdc
            CALL SWAP( C(P,J), C(Q,J) )
         END DO
      ENDIF
      RETURN
      END SUBROUTINE RCSWAP
      end module RINDMOD