2
/* cs_lu: sparse LU factorization, with optional fill-reducing ordering */
8
const mxArray *pargin [ ]
11
CS_INT n, order, *p, *Q ;
13
if (nargout > 4 || nargin > 3 || nargin < 1)
15
mexErrMsgTxt ("Usage: [L,U,p,q] = cs_lu (A,tol)") ;
17
if (nargin == 2) /* determine tol and ordering */
19
tol = mxGetScalar (pargin [1]) ;
20
order = (nargout == 4) ? 1 : 0 ; /* amd (A+A'), or natural */
25
order = (nargout == 4) ? 2 : 0 ; /* amd(S'*S) w/dense rows or I */
27
if (mxIsComplex (pargin [0]))
32
cs_cl Amatrix, *L, *U, *A, *D ;
33
A = cs_cl_mex_get_sparse (&Amatrix, 1, pargin [0]) ; /* get A */
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) ;
42
N->L = cs_cl_transpose (D, 1) ;
44
cs_cl_dropzeros (N->U) ; /* drop zeros from U and sort it */
45
D = cs_cl_transpose (N->U, 1) ;
47
N->U = cs_cl_transpose (D, 1) ;
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 */
54
if (nargout == 4) pargout [3] = cs_dl_mex_put_int (S->q, n, 1, 0) ;
58
mexErrMsgTxt ("complex matrices not supported") ;
65
cs_dl Amatrix, *L, *U, *A, *D ;
66
A = cs_dl_mex_get_sparse (&Amatrix, 1, 1, pargin [0]) ; /* get A */
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) ;
74
N->L = cs_dl_transpose (D, 1) ;
76
cs_dl_dropzeros (N->U) ; /* drop zeros from U and sort it */
77
D = cs_dl_transpose (N->U, 1) ;
79
N->U = cs_dl_transpose (D, 1) ;
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 */
86
if (nargout == 4) pargout [3] = cs_dl_mex_put_int (S->q, n, 1, 0) ;