#include "fintrf.h"
!C=======================================================================
!C	Version 2011-11-14 19:15:14
!C	Lin Renwen 
!C 	linrenwen@gmail.com
!C 	SUFE SHANGHAI
!C=======================================================================

! Purpose:  (see Peter Boettcher also)
! 
! Multiple-page A and Multiple-page B
! for index = 1:end
!     A(:,:, index) \ B(:,:, index);
! end
! 
! Multiple-page A and Single-page B
! for index = 1:end
!     A(:,:, index) \ B
! end
! 
! where index may be multi-dimensional;

! USAGE
! [C,D] = ndml(A, B)

! Hardware and Software
! AMD 
! Windows 7 32 bit 
! Matlab R2010a
! Visual Studio 2005
! Intel Fortran 11.0

! SOME HINT FOR CODING
!    For ND array 
!    1 prhs->xmp->x->xfp1D->xfp3D->xPagefp2D
!    2 plhs<-ymp<-yfp1D<-yfp3D<-yPagefp2D
!
!    1 numX  
!    2 dimsX -> numDims
!    3 dimsX -> mX nX -> strideX (numX) -> numPageX


! HOW TO COMPILE
!     mex  -g ndmm.F MatlabAPImx.obj MatlabAPImex.obj MatlabAPImat.obj 
!     D:\matlab2010\extern\lib\win32\microsoft\libmwlapack.lib 
!     D:\matlab2010\extern\lib\win32\microsoft\libmwblas.lib


! CALLS
!     libmwlapack.lib and James Tursa's Matlab API toolbox
!     We can use DGESV  in libmwlapack;
!     where INFO > 0, INPLIES THERE IS A PROBLEM OF SIGULARITY.

!C=======================================================================

!C gateway subroutine
!
!
!#ifndef mwSize
!#define mwSize integer*4
!#endif
!
!#ifndef mwPointer
!#define mwPointer integer*4
!#endif
!
!#ifndef mwIndex
!#define mwIndex integer*4
!#endif

!C      
!C=======================================================================   
      subroutine MEXFUNCTION(nlhs, plhs, nrhs, prhs) 
      implicit none
!-ARG MEXFUN
      mwPointer plhs(*), prhs(*)
      integer*4 nlhs, nrhs

!C     Function declarations:
      mwPointer mxGetDimensions
      mwSize  mxGetNumberOfElements
      mwSize mxGetNumberOfDimensions
      mwPointer mxCreateNumericArray
      mwPointer mxGetPr
      integer mxIsDouble
      integer*4 mxClassIDFromClassName
!      mwPointer mxIsInt8 
      mwPointer mxGetM, mxGetN
      mwPointer mxGetString    
      
      integer*4 classid
      integer*4 complexflag

!-LOC C/Y A B alpha beta TRANSA,TRANSB
      mwPointer Amp, Bmp, Cmp !ymp
      real(8), target, allocatable :: Cfp3D(:,:,:)
      real(8), target, allocatable :: Afp3D(:,:,:)
      real(8), target, allocatable :: Bfp3D(:,:,:)
      real*8, pointer :: ap2D(:,:),bp2D(:,:),cp2D(:,:)
      
!    1 numX  
!    2 dimsX -> numDims
!    3 dimsX -> mX nX -> strideX (numX) -> numPageX
      
      mwsize numC
      mwsize, allocatable :: dimsC(:)
      mwsize numDimsC ! dimsX -> numDims
      mwsize mC, nC, strideC, numPagesC 
      !dimsX -> mX nX -> strideX (numX) -> numPageX
      
      mwsize numA
      mwsize, allocatable :: dimsA(:)
      mwsize numDimsA ! dimsX -> numDims
      mwsize mA, nA, strideA, numPagesA 
      !dimsX -> mX nX -> strideX (numX) -> numPageX
      
      mwsize numB
      mwsize, allocatable :: dimsB(:)
      mwsize numDimsB ! dimsX -> numDims
      mwsize mB, nB, strideB, numPagesB 
      !dimsX -> mX nX -> strideX (numX) -> numPageX
  
! LOC TMP
      mwsize ii, jj, kk, tmp1, tmp2, status
! .. Parameters ..
      real(8) ONE,ZERO
      PARAMETER (ONE=1.0, ZERO=0.0)
      double precision one1, zero1

!C PART I

!C=======================================================================
! PEREPARE
!    1 numX
!    2 dimsX -> numDimsX
!    3 dimsX -> mX nX -> strideX (numX) -> numPagesX      

      numA = mxGetNumberOfElements(prhs(1))
      numDimsA = mxGetNumberOfDimensions(prhs(1))
      allocate(dimsA(numDimsA))
!      call mxCopyPtrToInteger4(mxGetDimensions(prhs(2)), dimsA,numDimsA)
      call mxCopyPtrToPtrArray(mxGetDimensions(prhs(1)), dimsA,numDimsA)
      mA = dimsA(1)
      nA = dimsA(2)
      strideA = mA * nA
      numPagesA = numA / strideA
      
      numB = mxGetNumberOfElements(prhs(2)) !
      numDimsB = mxGetNumberOfDimensions(prhs(2))
      allocate(dimsB(numDimsB))
      call mxCopyPtrToPtrArray(mxGetDimensions(prhs(2)), dimsB,numDimsB)
      mB = dimsB(1)
      nB = dimsB(2)
      strideB = mB * nB
      numPagesB = numB / strideB
      
      numC = mA * nB * numPagesA !
      numDimsC = numDimsA
      allocate(dimsC(numDimsC))
      dimsC = dimsA 
      dimsC(2) = dimsB(2) !
      mC = dimsC(1)
      nC = dimsC(2)
      strideC = mC * nC
      numPagesC = numC / strideC
      
! CHECK INPUT
      if( nrhs /= 2 .or. mxIsDouble(prhs(1)) == 0) then 
          call mexErrMsgTxt("Need 2 double 3-infD input") 
      endif 
      
      if( numPagesA /= numPagesB .or. nA /= mB) then 
          call mexErrMsgTxt("1. numPagesA /= numPagesB or 2. nA /= mB") 
      endif

!C BODY
!C=======================================================================
 ! parepare fun

      classid = mxClassIDFromClassName('double')
      complexflag = 0
      
      plhs(1) = 
     & mxCreateNumericArray(numDimsC, dimsC, classid, complexflag)
      Cmp = mxGetpr(plhs(1))
      Amp = mxGetpr(prhs(1))
      Bmp = mxGetpr(prhs(2))

      allocate(Afp3D(mA,nA,numPagesA))
      allocate(Bfp3D(mB,nB,numPagesB))
      allocate(Cfp3D(mC,nC,numPagesC))
      
      call mxCopyPtrToPtrArray(Amp, Afp3D, numA)
      call mxCopyPtrToPtrArray(Bmp, Bfp3D, numB)
      call mxCopyPtrToPtrArray(Cmp, Cfp3D, numC)

!       do kk = 1,numPagesC
!                  Cfp3D(1:mC,1:nC,kk) = matmul(Afp3D(1:mA,1:nA,kk),
!     &       Bfp3D(1:mB,1:nB,kk));
!       enddo

      do kk = 1,numPagesC
          ap2D => Afp3D(:,:,kk)
          bp2D => Bfp3D(:,:,kk)
          cp2D => Cfp3D(:,:,kk)
          cp2D = matmul(ap2D,bp2D);
!      call dgemm('N','N',mA,nB,nA,1.0,ap2D,mA,bp2D,nA,0.0,cp2D,mA)
!       dgemm('N', 'N', mA, nB, nA, 1.0, ap2D, mA, bp2D, nA, 0.0, cp2D, mA)
!          call dgemm(ch1, ch2, m, n, p, one, %val(A), m, %val(B), p, zero, %val(C), m)
      enddo

      call mxCopyReal8ToPtr(Cfp3D, Cmp, numC)

      DEallocate(Afp3D)
      DEallocate(Bfp3D)
      DEallocate(Cfp3D)
      
      deallocate(dimsA)
      deallocate(dimsB)
      deallocate(dimsC)

      return 
      end subroutine mexFunction

!C APPENDIX      
!C=======================================================================

 
!C=======================================================================
