1 #include "cs_mex.h"
2 /* cs_permute: permute a sparse matrix */
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 CS_INT ignore, *P, *Q, *Pinv, m, n ;
12 if (nargout > 1 || nargin != 3)
13 {
14 mexErrMsgTxt ("Usage: C = cs_permute(A,p,q)") ;
15 }
16 m = mxGetM (pargin [0]) ;
17 n = mxGetN (pargin [0]) ;
18 P = cs_dl_mex_get_int (m, pargin [1], &ignore, 1) ; /* get P */
19 Q = cs_dl_mex_get_int (n, pargin [2], &ignore, 1) ; /* get Q */
20 Pinv = cs_pinv (P, m) ; /* P = Pinv' */
21 if (mxIsComplex (pargin [0]))
22 {
23 #ifndef NCOMPLEX
24 cs_cl Amatrix, *A, *C, *D ;
25 A = cs_cl_mex_get_sparse (&Amatrix, 0, pargin [0]) ; /* get A */
26 C = cs_cl_permute (A, Pinv, Q, 1) ; /* C = A(p,q) */
27 cs_cl_free (A->x) ;
28 D = cs_cl_transpose (C, 1) ; /* sort C via double transpose */
29 cs_cl_spfree (C) ;
30 C = cs_cl_transpose (D, 1) ;
31 cs_cl_spfree (D) ;
32 pargout [0] = cs_cl_mex_put_sparse (&C) ; /* return C */
33 #else
34 mexErrMsgTxt ("complex matrices not supported") ;
35 #endif
36 }
37 else
38 {
39 cs_dl Amatrix, *A, *C, *D ;
40 A = cs_dl_mex_get_sparse (&Amatrix, 0, 1, pargin [0]) ; /* get A */
41 C = cs_dl_permute (A, Pinv, Q, 1) ; /* C = A(p,q) */
42 D = cs_dl_transpose (C, 1) ; /* sort C via double transpose */
43 cs_dl_spfree (C) ;
44 C = cs_dl_transpose (D, 1) ;
45 cs_dl_spfree (D) ;
46 pargout [0] = cs_dl_mex_put_sparse (&C) ; /* return C */
47 }
48 cs_free (Pinv) ;
49 cs_free (P) ;
50 cs_free (Q) ;
51 }
52