1 #include "cs_mex.h"
2 /* cs_qr: sparse QR factorization */
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 m, n, order, *p ;
12 if (nargout > 5 || nargin != 1)
13 {
14 mexErrMsgTxt ("Usage: [V,beta,p,R,q] = cs_qr(A)") ;
15 }
16 order = (nargout == 5) ? 3 : 0 ; /* determine ordering */
17 m = mxGetM (pargin [0]) ;
18 n = mxGetN (pargin [0]) ;
19 if (m < n) mexErrMsgTxt ("A must have # rows >= # columns") ;
20 if (mxIsComplex (pargin [0]))
21 {
22 #ifndef NCOMPLEX
23 cs_cls *S ;
24 cs_cln *N ;
25 cs_cl Amatrix, *A, *D ;
26 A = cs_cl_mex_get_sparse (&Amatrix, 0, pargin [0]) ; /* get A */
27 S = cs_cl_sqr (order, A, 1) ; /* symbolic QR ordering & analysis*/
28 N = cs_cl_qr (A, S) ; /* numeric QR factorization */
29 cs_free (A->x) ;
30 if (!N) mexErrMsgTxt ("qr failed") ;
31 cs_cl_dropzeros (N->L) ; /* drop zeros from V and sort */
32 D = cs_cl_transpose (N->L, 1) ;
33 cs_cl_spfree (N->L) ;
34 N->L = cs_cl_transpose (D, 1) ;
35 cs_cl_spfree (D) ;
36 cs_cl_dropzeros (N->U) ; /* drop zeros from R and sort */
37 D = cs_cl_transpose (N->U, 1) ;
38 cs_cl_spfree (N->U) ;
39 N->U = cs_cl_transpose (D, 1) ;
40 cs_cl_spfree (D) ;
41 m = N->L->m ; /* m may be larger now */
42 p = cs_cl_pinv (S->pinv, m) ; /* p = pinv' */
43 pargout [0] = cs_cl_mex_put_sparse (&(N->L)) ; /* return V */
44 cs_dl_mex_put_double (n, N->B, &(pargout [1])) ; /* return beta */
45 pargout [2] = cs_dl_mex_put_int (p, m, 1, 1) ; /* return p */
46 pargout [3] = cs_cl_mex_put_sparse (&(N->U)) ; /* return R */
47 pargout [4] = cs_dl_mex_put_int (S->q, n, 1, 0) ; /* return q */
48 cs_cl_nfree (N) ;
49 cs_cl_sfree (S) ;
50 #else
51 mexErrMsgTxt ("complex matrices not supported") ;
52 #endif
53 }
54 else
55 {
56 cs_dls *S ;
57 cs_dln *N ;
58 cs_dl Amatrix, *A, *D ;
59 A = cs_dl_mex_get_sparse (&Amatrix, 0, 1, pargin [0]) ; /* get A */
60 S = cs_dl_sqr (order, A, 1) ; /* symbolic QR ordering & analysis*/
61 N = cs_dl_qr (A, S) ; /* numeric QR factorization */
62 if (!N) mexErrMsgTxt ("qr failed") ;
63 cs_dl_dropzeros (N->L) ; /* drop zeros from V and sort */
64 D = cs_dl_transpose (N->L, 1) ;
65 cs_dl_spfree (N->L) ;
66 N->L = cs_dl_transpose (D, 1) ;
67 cs_dl_spfree (D) ;
68 cs_dl_dropzeros (N->U) ; /* drop zeros from R and sort */
69 D = cs_dl_transpose (N->U, 1) ;
70 cs_dl_spfree (N->U) ;
71 N->U = cs_dl_transpose (D, 1) ;
72 cs_dl_spfree (D) ;
73 m = N->L->m ; /* m may be larger now */
74 p = cs_dl_pinv (S->pinv, m) ; /* p = pinv' */
75 pargout [0] = cs_dl_mex_put_sparse (&(N->L)) ; /* return V */
76 cs_dl_mex_put_double (n, N->B, &(pargout [1])) ; /* return beta */
77 pargout [2] = cs_dl_mex_put_int (p, m, 1, 1) ; /* return p */
78 pargout [3] = cs_dl_mex_put_sparse (&(N->U)) ; /* return R */
79 pargout [4] = cs_dl_mex_put_int (S->q, n, 1, 0) ; /* return q */
80 cs_dl_nfree (N) ;
81 cs_dl_sfree (S) ;
82 }
83 }
84