1 #include "cs_mex.h"
2 /* z = cs_gaxpy (A,x,y) computes z = A*x+y */
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 if (nargout > 1 || nargin != 3)
12 {
13 mexErrMsgTxt ("Usage: z = cs_gaxpy(A,x,y)") ;
14 }
15 if (mxIsComplex (pargin [0]) || mxIsComplex (pargin [1])
16 || mxIsComplex (pargin [2]))
17 {
18 #ifndef NCOMPLEX
19 cs_cl Amatrix, *A ;
20 cs_complex_t *x, *z ;
21 A = cs_cl_mex_get_sparse (&Amatrix, 0, pargin [0]) ;/* get A */
22 x = cs_cl_mex_get_double (A->n, pargin [1]) ; /* get x */
23 z = cs_cl_mex_get_double (A->m, pargin [2]) ; /* z = y */
24 cs_cl_gaxpy (A, x, z) ; /* z = z + A*x */
25 cs_free (x) ;
26 cs_free (A->x) ;
27 pargout [0] = cs_cl_mex_put_double (A->m, z) ; /* return z */
28 #else
29 mexErrMsgTxt ("complex matrices not supported") ;
30 #endif
31 }
32 else
33 {
34 cs_dl Amatrix, *A ;
35 double *x, *y, *z ;
36 A = cs_dl_mex_get_sparse (&Amatrix, 0, 1, pargin [0]) ;/* get A */
37 x = cs_dl_mex_get_double (A->n, pargin [1]) ; /* get x */
38 y = cs_dl_mex_get_double (A->m, pargin [2]) ; /* get y */
39 z = cs_dl_mex_put_double (A->m, y, &(pargout [0])) ; /* z = y */
40 cs_dl_gaxpy (A, x, z) ; /* z = z + A*x */
41 }
42 }
43