1 #include "cs_mex.h"
2 /* cs_lusol: solve A*x=b using a sparse LU factorization */
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])3 void mexFunction
4 (
5 int nargout,
6 mxArray *pargout [ ],
7 int nargin,
8 const mxArray *pargin [ ]
9 )
10 {
11 double tol ;
12 CS_INT order ;
13 if (nargout > 1 || nargin < 2 || nargin > 4)
14 {
15 mexErrMsgTxt ("Usage: x = cs_lusol(A,b,order,tol)") ;
16 }
17 order = (nargin < 3) ? 2 : mxGetScalar (pargin [2]) ;
18 order = CS_MAX (order, 0) ;
19 order = CS_MIN (order, 3) ;
20 if (nargin == 2)
21 {
22 tol = 1 ; /* normal partial pivoting */
23 }
24 else if (nargin == 3)
25 {
26 tol = (order == 1) ? 0.001 : 1 ; /* tol = 0.001 for amd(A+A') */
27 }
28 else
29 {
30 tol = mxGetScalar (pargin [3]) ;
31 }
32 if (mxIsComplex (pargin [0]) || mxIsComplex (pargin [1]))
33 {
34 #ifndef NCOMPLEX
35 cs_cl *A, Amatrix ;
36 cs_complex_t *x ;
37 A = cs_cl_mex_get_sparse (&Amatrix, 1, pargin [0]) ; /* get A */
38 x = cs_cl_mex_get_double (A->n, pargin [1]) ; /* x = b */
39 if (!cs_cl_lusol (order, A, x, tol)) /* x = A\x */
40 {
41 mexErrMsgTxt ("failed (singular or out of memory)") ;
42 }
43 cs_cl_free (A->x) ; /* complex copy no longer needed */
44 pargout [0] = cs_cl_mex_put_double (A->n, x) ; /* return x */
45 #else
46 mexErrMsgTxt ("complex matrices not supported") ;
47 #endif
48 }
49 else
50 {
51 cs_dl *A, Amatrix ;
52 double *x, *b ;
53 A = cs_dl_mex_get_sparse (&Amatrix, 1, 1, pargin [0]) ; /* get A */
54 b = cs_dl_mex_get_double (A->n, pargin [1]) ; /* get b */
55 x = cs_dl_mex_put_double (A->n, b, &(pargout [0])) ; /* x = b */
56 if (!cs_dl_lusol (order, A, x, tol)) /* x = A\x */
57 {
58 mexErrMsgTxt ("failed (singular or out of memory)") ;
59 }
60 }
61 }
62