C     function L = myluinc(A,milu,no_steps)
C
C     Simplified incomplete LU factorisation method, which preserves the
C     sparsity pattern of A in the factorisation. milu is a flag which
C     if set will implement modified ilu, in which discarded fill-in
C     entries are subtracted from the diagonal to preserve the row sum.
C
C     Written by George Beckett - MATLAB 6.0 Fortran MEX
C     Thursday 31st May 2001.
C
C
      SUBROUTINE mexFunction(nlhs, plhs, nrhs, prhs)
      INTEGER plhs(*), prhs(*)
      INTEGER nlhs, nrhs
      INTEGER norows, nocols, nnzs
      INTEGER milu, steps
      INTEGER pSrci, pSrcj, pSrcr, ptemp
      INTEGER pDsti, pDstj, pDstr
C
C     Validate input and output arguments.
C
      IF((nrhs .EQ. 0) .OR. (nrhs .GT. 3)) THEN
         CALL mexErrMsgTxt('Invalid no. inputs.')
      ELSEIF(nlhs .GT. 1) THEN
         CALL mexErrMsgTxt('Invalid no. output arguments')
      ENDIF
C
C     Check that A is sparse, real, and has correct dimensions.
C
      IF(mxIsComplex(prhs(1)) .EQ. 1)
     +     CALL mexErrMsgTxt('Input A must be a real matrix')
C
      IF(mxIsSparse(prhs(1)) .EQ. 0)
     +     CALL mexErrMsgTxt('Input A must be a sparse matrix')
C
C     Read in data from input matrix
C
      norows = mxGetM(prhs(1))
      nocols = mxGetN(prhs(1))
      nnzs = mxGetNzmax(prhs(1))
      pSrci = mxGetIr(prhs(1))
      pSrcj = mxGetJc(prhs(1))
      pSrcr = mxGetPr(prhs(1))
C
      IF(norows .GT. nocols)
     +     CALL mexErrMsgTxt('In input A, rows>columns.')
C
C     Check to see if the user has specified an optional second
C     argument, i.e. milu. Default value is zero.
C
      milu = 0
C
      IF(nrhs .GT. 1) THEN
C     Check that the input is a scalar
         IF(mxGetM(prhs(2))*mxGetN(prhs(2)) .EQ. 1) THEN
            ptemp = mxGetPr(prhs(2))
            CALL getscalar(milu, %val(ptemp))
         ELSE
            CALL mexErrMsgTxt('Input milu must be a scalar.')
         ENDIF
      ENDIF
C
C     Check to see if the user has specified an optional third
C     argument. Default is norows-1
C
      steps = norows - 1
C
      IF(nrhs .EQ. 3) THEN
C     Check that the input is a scalar
         IF(mxGetM(prhs(3))*mxGetN(prhs(3)) .EQ. 1) THEN
            ptemp = mxGetPr(prhs(3))
            CALL getscalar(steps, %val(ptemp))
         ELSE
            CALL mexErrMsgTxt('Input no_steps must be a scalar.')
         ENDIF
      ENDIF
C
C     Create sparse output matrix, setup pointers, and copy input
C     matrix into it.
C
      plhs(1) = mxCreateSparse(norows, nocols, nnzs, 0)
      pDsti = mxGetIr(plhs(1))
      pDstj = mxGetJc(plhs(1))
      pDstr = mxGetPr(plhs(1))
C
      CALL copymtx(nocols, nnzs, %val(pSrci), %val(pSrcj), %val(pSrcr),
     +     %val(pDsti), %val(pDstj), %val(pDstr))
C
C     Now we compute the ilu factorisation of the matrix.
C
      CALL crtilu(milu, steps, norows, nocols, nnzs, %val(pDsti), 
     +     %val(pDstj), %val(pDstr))
C
      RETURN
      END
C
C
C     Subroutine to retrieve scalar from an input mxArray
C
C
      SUBROUTINE getscalar(dest, pSrc)
      REAL*8 temp
      REAL*8 pSrc(1)
      INTEGER dest
C
      temp = pSrc(1)
      dest = IFIX(temp)
C
      RETURN
      END
C
C
C     Subroutine to copy one sparse matrix into another.
C
C
      SUBROUTINE copymtx(nocols, nnzs, pAi, pAj, pAr, pBi, pBj, pBr)
      INTEGER nocols, nnzs, kk
      INTEGER pAi(nnzs), pAj(nocols+1), pBi(nnzs), pBj(nocols+1)
      REAL*8 pAr(nnzs), pBr(nnzs)
C
      DO 10 kk=1,nocols+1
         pBj(kk) = pAj(kk)
 10   CONTINUE
C
      DO 20 KK=1,nnzs
         pBi(kk) = pAi(kk)
         pBr(kk) = pAr(kk)
 20   CONTINUE
C
      RETURN
      END
C
C
C     Subroutine to create ilu decomposition of input matrix,
C     which is overwritten in the process. This one
C     takes a bit of following!!! 
C
C
      SUBROUTINE crtilu(milu, steps, norows, nocols, nnzs,
     +     pAi, pAj, pAr)
      INTEGER milu, steps, norows, nocols, nnzs
      INTEGER pAi(0:nnzs-1), pAj(0:nocols)
      REAL*8 pAr(0:nnzs-1)
C     
      INTEGER prow, crow, ccol
      INTEGER diaidx, subidx, rjidx, ijidx, iiidx, kk
      REAL*8 d, e
C
C     We step along the sub-diagonal elements one at a time.
C
      DO 10 prow = 0,steps-1
C     Search for appropriate diagonal entry in sparsity pattern.
         diaidx = -1

         DO 20 kk=pAj(prow),pAj(prow+1)-1
            IF(pAi(kk) .EQ. prow) THEN
               diaidx = kk
               GOTO 21
            ENDIF
 20      CONTINUE

C     Check that we found a diagonal entry.
 21      IF(diaidx .EQ. -1)
     +       CALL mexErrMsgTxt('Zero on diagonal has occured.')
C         
         d = 1.0/pAr(diaidx)
C
C     Now step down each sub-diagonal entry in column `prow'.
C
         DO 30 subidx = diaidx+1,pAj(prow+1)-1
            crow = pAi(subidx)
            e = pAr(subidx)*d
            pAr(subidx) = e

            DO 40 ccol = prow+1, nocols-1
C
C     Check that (prow,ccol) is nonzero.
C
               rjidx = -1
               
               DO 50 kk=pAj(ccol),pAj(ccol+1)-1
                  IF(pAi(kk) .EQ. prow) THEN
                     rjidx = kk
                     GOTO 51
                  ENDIF
 50            CONTINUE
C
C     If (prow,ccol) not in sparsity pattern then abort the current
C     loop.
C
 51            IF(rjidx .EQ. -1) GOTO 40
C               
C     Check that (crow,ccol) is nonzero.
C
               ijidx = -1
C
               DO 60 kk = rjidx+1,pAj(ccol+1)-1
                  IF(pAi(kk) .EQ. crow) THEN
                     ijidx = kk
                     GOTO 61
                  ENDIF
 60            CONTINUE
C
C     If (crow,ccol) not in sparsity pattern then either abort the
C     current loop or if milu is set add contribution to the appropriate
C     diagonal (crow,crow). */
C
 61            IF(ijidx .EQ. -1) THEN
                  IF(milu .EQ. 0) GOTO 40
C
                  iiidx = -1
C
                  DO 70 kk = pAj(crow),pAj(crow+1)-1
                     IF(pAi(kk) .EQ. crow) THEN
                        iiidx = kk
                        GOTO 71
                     ENDIF
 70               CONTINUE
C
 71               IF(iiidx .EQ. -1)
     +                 CALL mexErrMsgTxt('No diagonal entry for milu.')
C                  
                  pAr(iiidx) = pAr(iiidx) - e*pAr(rjidx)
               ELSE
C     Perform the elimination
                  pAr(ijidx) = pAr(ijidx) - e*pAr(rjidx)
               ENDIF
 40         CONTINUE
 30      CONTINUE
 10   CONTINUE
C
      RETURN
      END