2
/* cs_lsolve: x=L\b. L must be sparse and lower triangular. b must be a
3
* full or sparse vector. x is full or sparse, depending on b.
5
* Time taken is O(flop count), which may be less than n if b is sparse,
6
* depending on L and b.
8
* This function works with MATLAB 7.2, but is not perfectly compatible with
9
* the requirements of a MATLAB mexFunction when b is sparse. X is returned
10
* as an unsorted sparse vector. Also, this mexFunction temporarily modifies
11
* its input, L, by modifying L->p (in the cs_dfs function) and then restoring
12
* it. This could be corrected by creating a copy of L->p
13
* (see cs_dmperm_mex.c), but this would take O(n) time, destroying the
14
* O(flop count) time complexity of this function.
16
* Note that b cannot be sparse complex. This function does not support
17
* sparse complex L and b because the sparse x=L\b only accesses part of the
18
* matrix L. Converting L from a MATLAB complex matrix to a CXSparse complex
19
* matrix requires all of L to be accessed, defeating the purpose of this
22
* L can be sparse complex, but in that case b must be full real or complex,
31
const mxArray *pargin [ ]
34
CS_INT top, nz, p, *xi, n ;
35
if (nargout > 1 || nargin != 2)
37
mexErrMsgTxt ("Usage: x = cs_lsolve(L,b)") ;
39
if (mxIsSparse (pargin [1]))
41
cs_dl Lmatrix, Bmatrix, *L, *B, *X ;
43
if (mxIsComplex (pargin [0]) || mxIsComplex (pargin [1]))
45
mexErrMsgTxt ("sparse complex case not supported") ;
47
L = cs_dl_mex_get_sparse (&Lmatrix, 1, 1, pargin [0]) ;/* get L */
49
B = cs_dl_mex_get_sparse (&Bmatrix, 0, 1, pargin [1]) ;/* get sparse b*/
50
cs_mex_check (0, n, 1, 0, 1, 1, pargin [1]) ;
51
xi = cs_dl_malloc (2*n, sizeof (CS_INT)) ; /* get workspace */
52
x = cs_dl_malloc (n, sizeof (double)) ;
53
top = cs_dl_spsolve (L, B, 0, xi, x, NULL, 1) ; /* x = L\b */
54
X = cs_dl_spalloc (n, 1, n-top, 1, 0) ; /* create sparse x*/
57
for (p = top ; p < n ; p++)
60
X->x [nz++] = x [xi [p]] ;
63
pargout [0] = cs_dl_mex_put_sparse (&X) ;
67
else if (mxIsComplex (pargin [0]) || mxIsComplex (pargin [1]))
70
cs_cl Lmatrix, Bmatrix, *L, *B, *X ;
72
L = cs_cl_mex_get_sparse (&Lmatrix, 1, pargin [0]) ; /* get L */
74
x = cs_cl_mex_get_double (n, pargin [1]) ; /* x = b */
75
cs_cl_lsolve (L, x) ; /* x = L\x */
77
pargout [0] = cs_cl_mex_put_double (n, x) ; /* return x */
79
mexErrMsgTxt ("complex matrices not supported") ;
86
L = cs_dl_mex_get_sparse (&Lmatrix, 1, 1, pargin [0]) ; /* get L */
88
b = cs_dl_mex_get_double (n, pargin [1]) ; /* get b */
89
x = cs_dl_mex_put_double (n, b, &(pargout [0])) ; /* x = b */
90
cs_dl_lsolve (L, x) ; /* x = L\x */