/* function [L] = myluinc(A,milu,no_steps)

   Simplified incomplete LU factorisation method, which preserves the
   sparsity pattern of A in the factorisation. milu is a flag which if
   set will implement modified ilu, in which discarded fill-in entries
   are subtracted from the diagonal to preserve the row sum.

   Written by George Beckett
   Thursday 31st May 2001.
*/

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], \
		 int nrhs, const mxArray *prhs[]){
  int no_rows, no_cols, no_nonzeros, milu=0, no_steps;
  int *pAi=NULL, *pAj=NULL, *pSrci=NULL, *pSrcj=NULL;
  int pivot_row, current_row, current_col;
  int diag_index, sub_index, rj_index, ij_index, ii_index;
  int k;

  double *pAr=NULL, *pSrcr=NULL;
  double d, e;

  /* Validate no. inputs and outputs. */
  if(!nrhs || (nrhs>3))
    mexErrMsgTxt("Invalid no. inputs.");
  else if(nlhs>1)
    mexErrMsgTxt("Invalid no. outputs.");
 
  /* Check that A is sparse, real, and has correct dimensions. */
  if(mxIsComplex(prhs[0]))
    mexErrMsgTxt("Input A must be a real matrix.");

  if(!mxIsSparse(prhs[0]))
    mexErrMsgTxt("Input A must be a sparse matrix.");

  no_rows = mxGetM(prhs[0]);
  no_cols = mxGetN(prhs[0]);
  no_nonzeros = mxGetNzmax(prhs[0]);
  pSrci = mxGetIr(prhs[0]);
  pSrcj = mxGetJc(prhs[0]);
  pSrcr = mxGetPr(prhs[0]);

  if(no_rows>no_cols)
    mexErrMsgTxt("In input A, no. rows is greater than no. columns.");

  /* Check if the user has specified an optional second argument,
     i.e. milu. Default value is 0. */
  if(nrhs>1){
    if(!mxIsEmpty(prhs[1])){
      if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) ||
	 (mxGetN(prhs[1])*mxGetM(prhs[1]) != 1))
	mexErrMsgTxt("Input milu must be scalar.");
      
      milu = (int) mxGetScalar(prhs[1]);
    }
  } 


  /* Check to see if the user has specified an optional third
     argument. Default is norows-1. */

  no_steps = no_rows-1;

  if(nrhs==3){
    if(!mxIsEmpty(prhs[2])){
      if(!mxIsDouble(prhs[2]) || mxIsComplex(prhs[2]) ||
	 (mxGetN(prhs[2])*mxGetM(prhs[2]) != 1))
	mexErrMsgTxt("Input no_steps must be scalar.");
      
      no_steps = (int) mxGetScalar(prhs[2]);
    }
  } 
  
  /* Create sparse output matrix, setup pointers, and copy input
     matrix into output. */
  plhs[0] = mxCreateSparse(no_rows, no_cols, no_nonzeros, mxREAL);
  pAi = mxGetIr(plhs[0]);
  pAj = mxGetJc(plhs[0]);
  pAr = mxGetPr(plhs[0]);

  for(k=0; k<no_cols+1; k++)
    pAj[k] = pSrcj[k];

  for(k=0; k<no_nonzeros; k++){
    pAi[k] = pSrci[k];
    pAr[k] = pSrcr[k];
  }
      
  /* Now, we compute the ilu factorisation of the matrix. This one
     takes a bit of following!!! */

  /* We step along the sub-diagonal elements one at a time. */
  for(pivot_row=0; pivot_row<no_steps; pivot_row++){
    /* Search for appropriate diagonal entry in sparsity pattern. */
    diag_index=-1;

    for(k=pAj[pivot_row]; k<pAj[pivot_row+1]; k++){
      if(pAi[k]==pivot_row){
	diag_index=k;
	break;
      }
    }

    /* Check that we found a diagonal entry. */
    if(diag_index==-1)
      mexErrMsgTxt("Zero on diagonal has occured.");

    d = 1.0/pAr[diag_index];

    /* Now step down each sub-diagonal entry in column `pivot_row'. */
    for(sub_index=diag_index+1; sub_index<pAj[pivot_row+1];
	sub_index++){
      current_row = pAi[sub_index];
      e = pAr[sub_index]*d;
      pAr[sub_index] = e;

      for(current_col=pivot_row+1; current_col<no_cols;
	  current_col++){

	/* Check that (pivot_row,current_col) is nonzero. */
	rj_index=-1;
	
	for(k=pAj[current_col]; k<pAj[current_col+1]; k++)
	  if(pAi[k]==pivot_row){
	    rj_index=k;
	    break;
	  }

	/* If (pivot_row,current_col) not in sparsity pattern then
	   abort the current loop. */
	if(rj_index==-1) continue;

	/* Check that (current_row,current_col) is nonzero. */
	ij_index=-1;

	for(k=rj_index+1; k<pAj[current_col+1]; k++)
	  if(pAi[k]==current_row){
	    ij_index=k;
	    break;
	  }

	/* If (pivot_row,current_col) not in sparsity pattern then
	   either abort the current loop or if milu is set add
	   contribution to the appropriate diagonal
	   (current_row,current_row). */
	if(ij_index==-1){
	  if(!milu)
	    continue;
	  else{
	    /* Firstly find the appropriate entry. */
	    ii_index = -1;
	    
	    for(k=pAj[current_row]; k<pAj[current_row+1]; k++)
	      if(pAi[k]==current_row){
		ii_index = k;
		break;
	      }
	    
	    if(ii_index == -1)
	      mexErrMsgTxt("No diagonal entry for milu.");

	    pAr[ii_index] = pAr[ii_index] - e*pAr[rj_index];
	  }
	}
	else
	  /* Now perform the elimination. */
	  pAr[ij_index] = pAr[ij_index] - e*pAr[rj_index];
      }
    }
   
    /* mlfDisp(mlfFull(plhs[0])); */
  }

  return;
}