lusubs.c

ZEUS is a family of Eulerian (grid based) Magneto-Hydrodynamic codes (MHD) for use in astrophysics,

#include "mex.h"

void mexFunction( int nlhs,       mxArray *plhs[],
                  int nrhs, const mxArray *prhs[]
                 )
                
{   /* Forward/back solve function for LU factorisations with permutations
     *       
     *  Inputs:
     *           
     *           L,U - n x n sparse lower and upper factors
     *           P,Q - n x n sparse perumtation matrices
     *           b   - n x 1 dense vector
     *           
     *  Outputs:
     *   
     *           y = Q*( U\L\ (P*b) )
     *
     *  Darren Engwirda 2006 with thanks to Tim Davis       
     */  

    double *s, *b, *x, *y;
    int *jc, *ir, n, i, k;
    
    /* Check I/O number */
    if (nlhs!=1) {
        mexErrMsgTxt("Wrong number of outputs");
    }
    if (nrhs!=5) {
        mexErrMsgTxt("Wrong number of inputs");
    }
    
    /*** Brief error checking is only done on "b" ***/
    
    n = mxGetM(prhs[4]);
    
    if ((n!=mxGetM(prhs[0])) || (mxGetN(prhs[4])!=1)) {
        mexErrMsgTxt("Wrong input dimensions");
    }

    
    /* Alloacte output */
    plhs[0] = mxCreateDoubleMatrix(n, 1, mxREAL);
    
    /* Pointers to I/O */
    b = mxGetPr(prhs[4]);       /* Rhs vector    */
    y = mxGetPr(plhs[0]);       /* Output vector */
    
    /* Working variable */ 
    x = mxMalloc(n * sizeof(double));
    
    
    /* Make x = P*b where P is the permutation matrix */
    
    ir = mxGetIr(prhs[2]);      /* Permutation information from P */
    
    for (i=0; i<n; i++) {       
        x[ir[i]] = b[i];        
    }

    
    /* Solve L*x = x */
    
    s = mxGetPr(prhs[0]); ir = mxGetIr(prhs[0]); jc = mxGetJc(prhs[0]);
    
    for (i=0; i<n; i++) {
    
        int stop = jc[i+1];
        double temp = x[i] / s[jc[i]];
        
	    x[i] = temp;

	    for (k=jc[i]+1; k<stop; k++) {

	        x[ir[k]] -= s[k] * temp;
	    }
    }
    
    
    /* Solve U*x = x */
    
    s = mxGetPr(prhs[1]); ir = mxGetIr(prhs[1]); jc = mxGetJc(prhs[1]);
    
    for (i=n-1; i>=0; i--) {
	
        int stop = jc[i+1]-1;
        double temp = x[i] / s[stop];
    
	    x[i] = temp;

	    for (k=jc[i]; k<stop; k++)	{
	    
		    x[ir[k]] -= s[k] * temp;
		}
    }
    
    
    /* Make y = Q*x where Q is the permutation matrix */
    
    ir = mxGetIr(prhs[3]);      /* Permutation information from Q */
    
    for (i=0; i<n; i++) {       
        y[ir[i]] = x[i];        
    }
    
    
    /* Free allocated memory */
    mxFree(x);
    
    /* End */
}