1 #include "cs_mex.h"
2 /* cs_lu: sparse LU factorization, with optional fill-reducing ordering */
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 CS_INT n, order, *p ;
12 double tol ;
13 if (nargout > 4 || nargin > 3 || nargin < 1)
14 {
15 mexErrMsgTxt ("Usage: [L,U,p,q] = cs_lu (A,tol)") ;
16 }
17 if (nargin == 2) /* determine tol and ordering */
18 {
19 tol = mxGetScalar (pargin [1]) ;
20 order = (nargout == 4) ? 1 : 0 ; /* amd (A+A'), or natural */
21 }
22 else
23 {
24 tol = 1 ;
25 order = (nargout == 4) ? 2 : 0 ; /* amd(S'*S) w/dense rows or I */
26 }
27 if (mxIsComplex (pargin [0]))
28 {
29 #ifndef NCOMPLEX
30 cs_cls *S ;
31 cs_cln *N ;
32 cs_cl Amatrix, *A, *D ;
33 A = cs_cl_mex_get_sparse (&Amatrix, 1, pargin [0]) ; /* get A */
34 n = A->n ;
35 S = cs_cl_sqr (order, A, 0) ; /* symbolic ordering, no QR bound */
36 N = cs_cl_lu (A, S, tol) ; /* numeric factorization */
37 if (!N) mexErrMsgTxt ("cs_lu failed (singular, or out of memory)") ;
38 cs_cl_free (A->x) ; /* complex copy no longer needed */
39 cs_cl_dropzeros (N->L) ; /* drop zeros from L and sort it */
40 D = cs_cl_transpose (N->L, 1) ;
41 cs_cl_spfree (N->L) ;
42 N->L = cs_cl_transpose (D, 1) ;
43 cs_cl_spfree (D) ;
44 cs_cl_dropzeros (N->U) ; /* drop zeros from U and sort it */
45 D = cs_cl_transpose (N->U, 1) ;
46 cs_cl_spfree (N->U) ;
47 N->U = cs_cl_transpose (D, 1) ;
48 cs_cl_spfree (D) ;
49 p = cs_cl_pinv (N->pinv, n) ; /* p=pinv' */
50 pargout [0] = cs_cl_mex_put_sparse (&(N->L)) ; /* return L */
51 pargout [1] = cs_cl_mex_put_sparse (&(N->U)) ; /* return U */
52 pargout [2] = cs_dl_mex_put_int (p, n, 1, 1) ; /* return p */
53 /* return Q */
54 if (nargout == 4) pargout [3] = cs_dl_mex_put_int (S->q, n, 1, 0) ;
55 cs_cl_nfree (N) ;
56 cs_cl_sfree (S) ;
57 #else
58 mexErrMsgTxt ("complex matrices not supported") ;
59 #endif
60 }
61 else
62 {
63 cs_dls *S ;
64 cs_dln *N ;
65 cs_dl Amatrix, *A, *D ;
66 A = cs_dl_mex_get_sparse (&Amatrix, 1, 1, pargin [0]) ; /* get A */
67 n = A->n ;
68 S = cs_dl_sqr (order, A, 0) ; /* symbolic ordering, no QR bound */
69 N = cs_dl_lu (A, S, tol) ; /* numeric factorization */
70 if (!N) mexErrMsgTxt ("cs_lu failed (singular, or out of memory)") ;
71 cs_dl_dropzeros (N->L) ; /* drop zeros from L and sort it */
72 D = cs_dl_transpose (N->L, 1) ;
73 cs_dl_spfree (N->L) ;
74 N->L = cs_dl_transpose (D, 1) ;
75 cs_dl_spfree (D) ;
76 cs_dl_dropzeros (N->U) ; /* drop zeros from U and sort it */
77 D = cs_dl_transpose (N->U, 1) ;
78 cs_dl_spfree (N->U) ;
79 N->U = cs_dl_transpose (D, 1) ;
80 cs_dl_spfree (D) ;
81 p = cs_dl_pinv (N->pinv, n) ; /* p=pinv' */
82 pargout [0] = cs_dl_mex_put_sparse (&(N->L)) ; /* return L */
83 pargout [1] = cs_dl_mex_put_sparse (&(N->U)) ; /* return U */
84 pargout [2] = cs_dl_mex_put_int (p, n, 1, 1) ; /* return p */
85 /* return Q */
86 if (nargout == 4) pargout [3] = cs_dl_mex_put_int (S->q, n, 1, 0) ;
87 cs_dl_nfree (N) ;
88 cs_dl_sfree (S) ;
89 }
90 }
91