1 // =============================================================================
2 // === spqr_qmult mexFunction ==================================================
3 // =============================================================================
4 
5 #include "spqr_mx.hpp"
6 
7 //  Multiply Q times X, where Q is stored as a struct, in Householder form.
8 //
9 //  method = 0: Y = Q'*X    default
10 //  method = 1: Y = Q*X
11 //  method = 2: Y = X*Q'
12 //  method = 3: Y = X*Q
13 //
14 //  Usage:
15 //
16 //  Y = spqr_qmult (Q,X,method) ;
17 //
18 //  where Q is the struct from [Q,R,E] = spqr (A,opts) with
19 //  opts.Q = 'Householder'
20 //
21 //  Q.H: Householder vectors (m-by-nh).  In each column, the nonzero entry with
22 //      the smallest row index must be equal to 1.0.
23 //  Q.Tau: Householder coefficients (1-by-nh).
24 //  Q.P: inverse row permutation.  P [i] = k if row i of X is row k of H and Y.
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     mxArray *Hmatlab, *Tau, *P ;
35     Long *HPinv, *Yp, *Yi ;
36     double *Hx, *Xx, *Tx, *Px, dummy ;
37     Long m, n, k, nh, nb, p, i, method, mh, gotP, X_is_sparse, is_complex, hnz,
38         tnz, xnz, inuse, count ;
39     cholmod_sparse *Ysparse, *H, Hmatrix, *Xsparse, Xsmatrix ;
40     cholmod_dense *Ydense, *Xdense, Xdmatrix, *HTau, HTau_matrix ;
41     cholmod_common Common, *cc ;
42 
43     // -------------------------------------------------------------------------
44     // start CHOLMOD and set parameters
45     // -------------------------------------------------------------------------
46 
47     cc = &Common ;
48     cholmod_l_start (cc) ;
49     spqr_mx_config (SPUMONI, cc) ;
50 
51     // -------------------------------------------------------------------------
52     // check inputs
53     // -------------------------------------------------------------------------
54 
55     // nargin can be 2 or 3
56     // nargout can be 0 or 1
57 
58     if (nargout > 1)
59     {
60         mexErrMsgIdAndTxt ("MATLAB:maxlhs", "Too many output arguments") ;
61     }
62     if (nargin < 2)
63     {
64         mexErrMsgIdAndTxt ("MATLAB:minrhs", "Not enough input arguments") ;
65     }
66     if (nargin > 3)
67     {
68         mexErrMsgIdAndTxt ("MATLAB:maxrhs", "Too many input arguments") ;
69     }
70 
71     if (!mxIsStruct (pargin [0]))
72     {
73         mexErrMsgIdAndTxt ("QR:invalidInput", "invalid Q (must be a struct)") ;
74     }
75 
76     // -------------------------------------------------------------------------
77     // get H, Tau, and P from the Q struct
78     // -------------------------------------------------------------------------
79 
80     i = mxGetFieldNumber (pargin [0], "H") ;
81     if (i < 0)
82     {
83         mexErrMsgIdAndTxt ("QR:invalidInput", "invalid Q struct") ;
84     }
85     Hmatlab = mxGetFieldByNumber (pargin [0], 0, i) ;
86     nh = mxGetN (Hmatlab) ;
87     if (!mxIsSparse (Hmatlab))
88     {
89         mexErrMsgIdAndTxt ("QR:invalidInput", "H must be sparse") ;
90     }
91     i = mxGetFieldNumber (pargin [0], "Tau") ;
92     if (i < 0)
93     {
94         mexErrMsgIdAndTxt ("QR:invalidInput", "invalid Q struct") ;
95     }
96     Tau = mxGetFieldByNumber (pargin [0], 0, i) ;
97     if (nh != mxGetNumberOfElements (Tau))
98     {
99         mexErrMsgIdAndTxt ("QR:invalidInput",
100             "H and Tau must have the same number of columns") ;
101     }
102 
103     is_complex = mxIsComplex (Tau) || mxIsComplex (Hmatlab) ||
104         mxIsComplex (pargin [1]) ;
105 
106     // -------------------------------------------------------------------------
107     // get the Householder vectors
108     // -------------------------------------------------------------------------
109 
110     H = spqr_mx_get_sparse (Hmatlab, &Hmatrix, &dummy) ;
111     mh = H->nrow ;
112     Hx = spqr_mx_merge_if_complex (Hmatlab, is_complex, &hnz, cc) ;
113     if (is_complex)
114     {
115         // H has been converted from real or zomplex to complex
116         H->x = Hx ;
117         H->z = NULL ;
118         H->xtype = CHOLMOD_COMPLEX ;
119     }
120 
121     // -------------------------------------------------------------------------
122     // get Tau
123     // -------------------------------------------------------------------------
124 
125     HTau = spqr_mx_get_dense (Tau, &HTau_matrix, &dummy) ;
126     Tx = spqr_mx_merge_if_complex (Tau, is_complex, &tnz, cc) ;
127     if (is_complex)
128     {
129         // HTau has been converted from real or zomplex to complex
130         HTau->x = Tx ;
131         HTau->z = NULL ;
132         HTau->xtype = CHOLMOD_COMPLEX ;
133     }
134 
135     // -------------------------------------------------------------------------
136     // get method
137     // -------------------------------------------------------------------------
138 
139     if (nargin < 3)
140     {
141         method = 0 ;
142     }
143     else
144     {
145         method = (Long) mxGetScalar (pargin [2]) ;
146         if (method < 0 || method > 3)
147         {
148             mexErrMsgIdAndTxt ("QR:invalidInput", "invalid method") ;
149         }
150     }
151 
152     // -------------------------------------------------------------------------
153     // get X
154     // -------------------------------------------------------------------------
155 
156     m = mxGetM (pargin [1]) ;
157     n = mxGetN (pargin [1]) ;
158     X_is_sparse = mxIsSparse (pargin [1]) ;
159     Xsparse = NULL ;
160     if (X_is_sparse)
161     {
162         Xsparse = spqr_mx_get_sparse (pargin [1], &Xsmatrix, &dummy) ;
163     }
164     else
165     {
166         Xdense = spqr_mx_get_dense (pargin [1], &Xdmatrix, &dummy) ;
167     }
168     Xx = spqr_mx_merge_if_complex (pargin [1], is_complex, &xnz, cc) ;
169     if (is_complex)
170     {
171         // X has been converted from real or zomplex to complex
172         if (X_is_sparse)
173         {
174             Xsparse->x = Xx ;
175             Xsparse->xtype = CHOLMOD_COMPLEX ;
176         }
177         else
178         {
179             Xdense->x = Xx ;
180             Xdense->xtype = CHOLMOD_COMPLEX ;
181         }
182     }
183 
184     if (method == 0 || method == 1)
185     {
186         if (mh != m)
187         {
188             mexErrMsgIdAndTxt ("QR:invalidInput",
189                 "H and X must have same number of rows") ;
190         }
191     }
192     else
193     {
194         if (mh != n)
195         {
196             mexErrMsgIdAndTxt ("QR:invalidInput",
197                 "# of cols of X must equal # of rows of H") ;
198         }
199     }
200 
201     // -------------------------------------------------------------------------
202     // get P
203     // -------------------------------------------------------------------------
204 
205     i = mxGetFieldNumber (pargin [0], "P") ;
206     gotP = (i >= 0) ;
207     HPinv = NULL ;
208 
209     if (gotP)
210     {
211         // get P from the H struct
212         P = mxGetFieldByNumber (pargin [0], 0, i) ;
213         if (mxGetNumberOfElements (P) != mh)
214         {
215             mexErrMsgIdAndTxt ("QR:invalidInput",
216                 "P must be a vector of length equal to # rows of H") ;
217         }
218         HPinv = (Long *) cholmod_l_malloc (mh, sizeof (Long), cc) ;
219         Px = mxGetPr (P) ;
220         for (i = 0 ; i < mh ; i++)
221         {
222             HPinv [i] = (Long) (Px [i] - 1) ;
223             if (HPinv [i] < 0 || HPinv [i] >= mh)
224             {
225                 mexErrMsgIdAndTxt ("QR:invalidInput", "invalid permutation") ;
226             }
227         }
228     }
229 
230     // -------------------------------------------------------------------------
231     // Y = Q'*X, Q*X, X*Q or X*Q'
232     // -------------------------------------------------------------------------
233 
234     if (is_complex)
235     {
236         if (X_is_sparse)
237         {
238             Ysparse = SuiteSparseQR_qmult <Complex> (method, H,
239                 HTau, HPinv, Xsparse, cc) ;
240             pargout [0] = spqr_mx_put_sparse (&Ysparse, cc) ;
241         }
242         else
243         {
244             Ydense = SuiteSparseQR_qmult <Complex> (method, H,
245                 HTau, HPinv, Xdense, cc) ;
246             pargout [0] = spqr_mx_put_dense (&Ydense, cc) ;
247         }
248     }
249     else
250     {
251         if (X_is_sparse)
252         {
253             Ysparse = SuiteSparseQR_qmult <double> (method, H,
254                 HTau, HPinv, Xsparse, cc) ;
255             pargout [0] = spqr_mx_put_sparse (&Ysparse, cc) ;
256         }
257         else
258         {
259             Ydense = SuiteSparseQR_qmult <double> (method, H,
260                 HTau, HPinv, Xdense, cc) ;
261             pargout [0] = spqr_mx_put_dense (&Ydense, cc) ;
262         }
263     }
264 
265     // -------------------------------------------------------------------------
266     // free workspace
267     // -------------------------------------------------------------------------
268 
269     cholmod_l_free (mh, sizeof (Long), HPinv, cc) ;
270 
271     if (is_complex)
272     {
273         // free the merged copies of the real parts of the H and Tau matrices
274         cholmod_l_free (hnz, sizeof (Complex), Hx, cc) ;
275         cholmod_l_free (tnz, sizeof (Complex), Tx, cc) ;
276         cholmod_l_free (xnz, sizeof (Complex), Xx, cc) ;
277     }
278     cholmod_l_finish (cc) ;
279 
280 #if 0
281     // malloc count for testing only ...
282     spqr_mx_get_usage (pargout [0], 1, &inuse, &count, cc) ;
283     if (inuse != cc->memory_inuse || count != cc->malloc_count)
284     {
285         mexErrMsgIdAndTxt ("QR:internalError", "memory leak!") ;
286     }
287 #endif
288 }
289