1 #include "cs_mex.h"
2 /* cs_usolve: x=U\b. U must be sparse and upper triangular. b must be a
3 * full or sparse vector. x is full or sparse, depending on b.
4 *
5 * Time taken is O(flop count), which may be less than n if b is sparse,
6 * depending on U and b.
7 *
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, U, by modifying U->p (in the cs_dfs function) and then restoring
12 * it. This could be corrected by creating a copy of U->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.
15 *
16 * Note that b cannot be sparse complex. This function does not support
17 * sparse complex U and b because the sparse x=U\b only accesses part of the
18 * matrix U. Converting U from a MATLAB complex matrix to a CXSparse complex
19 * matrix requires all of U to be accessed, defeating the purpose of this
20 * function.
21 *
22 * U can be sparse complex, but in that case b must be full real or complex,
23 * not sparse.
24 */
25
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])26 void mexFunction
27 (
28 int nargout,
29 mxArray *pargout [ ],
30 int nargin,
31 const mxArray *pargin [ ]
32 )
33 {
34 CS_INT top, nz, p, *xi, n ;
35 if (nargout > 1 || nargin != 2)
36 {
37 mexErrMsgTxt ("Usage: x = cs_usolve(U,b)") ;
38 }
39 if (mxIsSparse (pargin [1]))
40 {
41 cs_dl Umatrix, Bmatrix, *U, *B, *X ;
42 double *x ;
43 if (mxIsComplex (pargin [0]) || mxIsComplex (pargin [1]))
44 {
45 mexErrMsgTxt ("sparse complex case not supported") ;
46 }
47 U = cs_dl_mex_get_sparse (&Umatrix, 1, 1, pargin [0]) ;/* get U */
48 n = U->n ;
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 (U, B, 0, xi, x, NULL, 0) ; /* x = U\b */
54 X = cs_dl_spalloc (n, 1, n-top, 1, 0) ; /* create sparse x*/
55 X->p [0] = 0 ;
56 nz = 0 ;
57 for (p = top ; p < n ; p++)
58 {
59 X->i [nz] = xi [p] ;
60 X->x [nz++] = x [xi [p]] ;
61 }
62 X->p [1] = nz ;
63 pargout [0] = cs_dl_mex_put_sparse (&X) ;
64 cs_free (x) ;
65 cs_free (xi) ;
66 }
67 else if (mxIsComplex (pargin [0]) || mxIsComplex (pargin [1]))
68 {
69 #ifndef NCOMPLEX
70 cs_cl Umatrix, *U ;
71 cs_complex_t *x ;
72 U = cs_cl_mex_get_sparse (&Umatrix, 1, pargin [0]) ; /* get U */
73 n = U->n ;
74 x = cs_cl_mex_get_double (n, pargin [1]) ; /* x = b */
75 cs_cl_usolve (U, x) ; /* x = U\x */
76 cs_free (U->x) ;
77 pargout [0] = cs_cl_mex_put_double (n, x) ; /* return x */
78 #else
79 mexErrMsgTxt ("complex matrices not supported") ;
80 #endif
81 }
82 else
83 {
84 cs_dl Umatrix, *U ;
85 double *x, *b ;
86 U = cs_dl_mex_get_sparse (&Umatrix, 1, 1, pargin [0]) ; /* get U */
87 n = U->n ;
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_usolve (U, x) ; /* x = U\x */
91 }
92 }
93