1 #include "cs_mex.h"
2 /* cs_add: sparse matrix addition */
3
4 #ifndef NCOMPLEX
get_complex(const mxArray * a)5 static cs_complex_t get_complex (const mxArray *a)
6 {
7 cs_complex_t s = mxGetScalar (a) ;
8 if (mxIsComplex (a))
9 {
10 double *z = mxGetPi (a) ;
11 s += I * z [0] ;
12 }
13 return (s) ;
14 }
15 #endif
16
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])17 void mexFunction
18 (
19 int nargout,
20 mxArray *pargout [ ],
21 int nargin,
22 const mxArray *pargin [ ]
23 )
24 {
25 if (nargout > 1 || nargin < 2 || nargin > 4)
26 {
27 mexErrMsgTxt ("Usage: C = cs_add(A,B,alpha,beta)") ;
28 }
29 if (mxIsComplex (pargin [0]) || mxIsComplex (pargin [1])
30 || (nargin > 2 && mxIsComplex (pargin [2]))
31 || (nargin > 3 && mxIsComplex (pargin [3])))
32 {
33 #ifndef NCOMPLEX
34 cs_complex_t alpha, beta ;
35 cs_cl Amatrix, Bmatrix, *A, *B, *C, *D ;
36 A = cs_cl_mex_get_sparse (&Amatrix, 0, pargin [0]) ; /* get A */
37 B = cs_cl_mex_get_sparse (&Bmatrix, 0, pargin [1]) ; /* get B */
38 alpha = (nargin < 3) ? 1 : get_complex (pargin [2]) ; /* get alpha */
39 beta = (nargin < 4) ? 1 : get_complex (pargin [3]) ; /* get beta */
40 C = cs_cl_add (A,B,alpha,beta) ; /* C = alpha*A + beta *B */
41 cs_cl_dropzeros (C) ; /* drop zeros */
42 D = cs_cl_transpose (C, 1) ; /* sort result via double transpose */
43 cs_cl_spfree (C) ;
44 C = cs_cl_transpose (D, 1) ;
45 cs_cl_spfree (D) ;
46 pargout [0] = cs_cl_mex_put_sparse (&C) ; /* return C */
47 #else
48 mexErrMsgTxt ("complex matrices not supported") ;
49 #endif
50 }
51 else
52 {
53 double alpha, beta ;
54 cs_dl Amatrix, Bmatrix, *A, *B, *C, *D ;
55 A = cs_dl_mex_get_sparse (&Amatrix, 0, 1, pargin [0]) ; /* get A */
56 B = cs_dl_mex_get_sparse (&Bmatrix, 0, 1, pargin [1]) ; /* get B */
57 alpha = (nargin < 3) ? 1 : mxGetScalar (pargin [2]) ; /* get alpha */
58 beta = (nargin < 4) ? 1 : mxGetScalar (pargin [3]) ; /* get beta */
59 C = cs_dl_add (A,B,alpha,beta) ; /* C = alpha*A + beta *B */
60 cs_dl_dropzeros (C) ; /* drop zeros */
61 D = cs_dl_transpose (C, 1) ; /* sort result via double transpose */
62 cs_dl_spfree (C) ;
63 C = cs_dl_transpose (D, 1) ;
64 cs_dl_spfree (D) ;
65 pargout [0] = cs_dl_mex_put_sparse (&C) ; /* return C */
66 }
67 }
68