|
5164
|
1 /* ========================================================================== */ |
|
|
2 /* === UMFPACK mexFunction ================================================== */ |
|
|
3 /* ========================================================================== */ |
|
|
4 |
|
|
5 /* -------------------------------------------------------------------------- */ |
|
|
6 /* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis. CISE Dept, */ |
|
|
7 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */ |
|
|
8 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */ |
|
|
9 /* -------------------------------------------------------------------------- */ |
|
|
10 |
|
|
11 /* |
|
|
12 MATLAB interface for umfpack. |
|
|
13 |
|
|
14 Factor or solve a sparse linear system, returning either the solution |
|
|
15 x to Ax=b or A'x'=b', or the factorization LU=P(R\A)Q or LU=PAQ. A must be |
|
|
16 sparse, with nonzero dimensions, but it may be complex, singular, and/or |
|
|
17 rectangular. b must be a dense n-by-1 vector (real or complex). |
|
|
18 L is unit lower triangular, U is upper triangular, and R is diagonal. |
|
|
19 P and Q are permutation matrices (permutations of an identity matrix). |
|
|
20 |
|
|
21 The matrix A is scaled, by default. Each row i is divided by r (i), where |
|
|
22 r (i) is the sum of the absolute values of the entries in that row. The |
|
|
23 scaled matrix has an infinity norm of 1. The scale factors r (i) are |
|
|
24 returned in a diagonal sparse matrix. If the factorization is: |
|
|
25 |
|
|
26 [L, U, P, Q, R] = umfpack (A) ; |
|
|
27 |
|
|
28 then the factorization is |
|
|
29 |
|
|
30 L*U = P * (R \ A) * Q |
|
|
31 |
|
|
32 This is safer than returning a matrix R such that L*U = P*R*A*Q, because |
|
|
33 it avoids the division by small entries. If r(i) is subnormal, multiplying |
|
|
34 by 1/r(i) would result in an IEEE Infinity, but dividing by r(i) is safe. |
|
|
35 |
|
|
36 The factorization |
|
|
37 |
|
|
38 [L, U, P, Q] = umfpack (A) ; |
|
|
39 |
|
|
40 returns LU factors such that L*U = P*A*Q, with no scaling. |
|
|
41 |
|
|
42 See umfpack.m, umfpack_details.m, and umfpack.h for details. |
|
|
43 |
|
|
44 Note that this mexFunction accesses only the user-callable UMFPACK routines. |
|
|
45 Thus, is also provides another example of how user C code can access |
|
|
46 UMFPACK. |
|
|
47 |
|
|
48 If NO_TRANSPOSE_FORWARD_SLASH is not defined at compile time, then the |
|
|
49 forward slash (/) operator acts almost like x = b/A in MATLAB 6.1. It is |
|
|
50 solved by factorizing the array transpose, and then x = (A.'\b.').' is |
|
|
51 solved. This is the default behavior (for historical reasons), since |
|
|
52 factorizing A can behave perform much differently than factorizing its |
|
|
53 transpose. |
|
|
54 |
|
|
55 If NO_TRANSPOSE_FORWARD_SLASH is defined at compile time, then the forward |
|
|
56 slash operator does not act like x=b/A in MATLAB 6.1. It is solved by |
|
|
57 factorizing A, and then solving via the transposed L and U matrices. |
|
|
58 The solution is still x = (A.'\b.').', except that A is factorized instead |
|
|
59 of A.'. |
|
|
60 |
|
|
61 Modified for v4.3.1, Jan 10, 2005: default has been changed to |
|
|
62 NO_TRANSPOSE_FORWARD_SLASH, to test iterative refinement for b/A. |
|
|
63 v4.4: added method for computing the determinant. |
|
|
64 */ |
|
|
65 #define NO_TRANSPOSE_FORWARD_SLASH /* default has changed for v4.3.1 */ |
|
|
66 |
|
|
67 #include "umfpack.h" |
|
|
68 #include "mex.h" |
|
|
69 #include "matrix.h" |
|
|
70 #include <string.h> |
|
|
71 #include <math.h> |
|
|
72 #include <float.h> |
|
|
73 |
|
|
74 #define MIN(a,b) (((a) < (b)) ? (a) : (b)) |
|
|
75 #define MAX(a,b) (((a) > (b)) ? (a) : (b)) |
|
|
76 #define STRING_MATCH(s1,s2) (strcmp ((s1), (s2)) == 0) |
|
|
77 #ifndef TRUE |
|
|
78 #define TRUE (1) |
|
|
79 #endif |
|
|
80 #ifndef FALSE |
|
|
81 #define FALSE (0) |
|
|
82 #endif |
|
|
83 |
|
|
84 /* ========================================================================== */ |
|
|
85 /* === error ================================================================ */ |
|
|
86 /* ========================================================================== */ |
|
|
87 |
|
|
88 /* Return an error message */ |
|
|
89 |
|
|
90 static void error |
|
|
91 ( |
|
|
92 char *s, |
|
|
93 int A_is_complex, |
|
|
94 int nargout, |
|
|
95 mxArray *pargout [ ], |
|
|
96 double Control [UMFPACK_CONTROL], |
|
|
97 double Info [UMFPACK_INFO], |
|
|
98 int status, |
|
|
99 int do_info |
|
|
100 ) |
|
|
101 { |
|
|
102 int i ; |
|
|
103 double *OutInfo ; |
|
|
104 if (A_is_complex) |
|
|
105 { |
|
|
106 umfpack_zi_report_status (Control, status) ; |
|
|
107 umfpack_zi_report_info (Control, Info) ; |
|
|
108 } |
|
|
109 else |
|
|
110 { |
|
|
111 umfpack_di_report_status (Control, status) ; |
|
|
112 umfpack_di_report_info (Control, Info) ; |
|
|
113 } |
|
|
114 if (do_info > 0) |
|
|
115 { |
|
|
116 /* return Info */ |
|
|
117 pargout [do_info] = mxCreateDoubleMatrix (1, UMFPACK_INFO, mxREAL) ; |
|
|
118 OutInfo = mxGetPr (pargout [do_info]) ; |
|
|
119 for (i = 0 ; i < UMFPACK_INFO ; i++) |
|
|
120 { |
|
|
121 OutInfo [i] = Info [i] ; |
|
|
122 } |
|
|
123 } |
|
|
124 mexErrMsgTxt (s) ; |
|
|
125 } |
|
|
126 |
|
|
127 |
|
|
128 /* ========================================================================== */ |
|
|
129 /* === UMFPACK ============================================================== */ |
|
|
130 /* ========================================================================== */ |
|
|
131 |
|
|
132 void mexFunction |
|
|
133 ( |
|
|
134 int nargout, /* number of outputs */ |
|
|
135 mxArray *pargout [ ], /* output arguments */ |
|
|
136 int nargin, /* number of inputs */ |
|
|
137 const mxArray *pargin [ ] /* input arguments */ |
|
|
138 ) |
|
|
139 { |
|
|
140 |
|
|
141 /* ---------------------------------------------------------------------- */ |
|
|
142 /* local variables */ |
|
|
143 /* ---------------------------------------------------------------------- */ |
|
|
144 |
|
|
145 double Info [UMFPACK_INFO], Control [UMFPACK_CONTROL], dx, dz, dexp ; |
|
|
146 double *Lx, *Lz, *Ux, *Uz, *Ax, *Az, *Bx, *Bz, *Xx, *Xz, *User_Control, |
|
|
147 *p, *q, *OutInfo, *p1, *p2, *p3, *p4, *Ltx, *Ltz, *Rs, *Px, *Qx ; |
|
|
148 void *Symbolic, *Numeric ; |
|
|
149 int *Lp, *Li, *Up, *Ui, *Ap, *Ai, *P, *Q, do_solve, lnz, unz, nn, i, |
|
|
150 transpose, size, do_info, do_numeric, *Front_npivcol, op, k, *Rp, *Ri, |
|
|
151 *Front_parent, *Chain_start, *Chain_maxrows, *Chain_maxcols, nz, status, |
|
|
152 nfronts, nchains, *Ltp, *Ltj, *Qinit, print_level, status2, no_scale, |
|
|
153 *Front_1strow, *Front_leftmostdesc, n_row, n_col, n_inner, sys, |
|
|
154 ignore1, ignore2, ignore3, A_is_complex, B_is_complex, X_is_complex, |
|
|
155 *Pp, *Pi, *Qp, *Qi, do_recip, do_det ; |
|
|
156 mxArray *Amatrix, *Bmatrix, *User_Control_matrix, *User_Qinit ; |
|
|
157 char *operator, *operation ; |
|
|
158 mxComplexity Atype, Xtype ; |
|
|
159 char warning [200] ; |
|
|
160 |
|
|
161 #ifndef NO_TRANSPOSE_FORWARD_SLASH |
|
|
162 int *Cp, *Ci ; |
|
|
163 double *Cx, *Cz ; |
|
|
164 #endif |
|
|
165 |
|
|
166 /* ---------------------------------------------------------------------- */ |
|
|
167 /* get inputs A, b, and the operation to perform */ |
|
|
168 /* ---------------------------------------------------------------------- */ |
|
|
169 |
|
|
170 User_Control_matrix = (mxArray *) NULL ; |
|
|
171 User_Qinit = (mxArray *) NULL ; |
|
|
172 |
|
|
173 do_info = 0 ; |
|
|
174 do_solve = FALSE ; |
|
|
175 do_numeric = TRUE ; |
|
|
176 transpose = FALSE ; |
|
|
177 no_scale = FALSE ; |
|
|
178 do_det = FALSE ; |
|
|
179 |
|
|
180 /* find the operator */ |
|
|
181 op = 0 ; |
|
|
182 for (i = 0 ; i < nargin ; i++) |
|
|
183 { |
|
|
184 if (mxIsChar (pargin [i])) |
|
|
185 { |
|
|
186 op = i ; |
|
|
187 break ; |
|
|
188 } |
|
|
189 } |
|
|
190 |
|
|
191 if (op > 0) |
|
|
192 { |
|
|
193 operator = mxArrayToString (pargin [op]) ; |
|
|
194 |
|
|
195 if (STRING_MATCH (operator, "\\")) |
|
|
196 { |
|
|
197 |
|
|
198 /* -------------------------------------------------------------- */ |
|
|
199 /* matrix left divide, x = A\b */ |
|
|
200 /* -------------------------------------------------------------- */ |
|
|
201 |
|
|
202 /* |
|
|
203 [x, Info] = umfpack (A, '\', b) ; |
|
|
204 [x, Info] = umfpack (A, '\', b, Control) ; |
|
|
205 [x, Info] = umfpack (A, Qinit, '\', b, Control) ; |
|
|
206 [x, Info] = umfpack (A, Qinit, '\', b) ; |
|
|
207 */ |
|
|
208 |
|
|
209 operation = "x = A\\b" ; |
|
|
210 do_solve = TRUE ; |
|
|
211 Amatrix = (mxArray *) pargin [0] ; |
|
|
212 Bmatrix = (mxArray *) pargin [op+1] ; |
|
|
213 |
|
|
214 if (nargout == 2) |
|
|
215 { |
|
|
216 do_info = 1 ; |
|
|
217 } |
|
|
218 if (op == 2) |
|
|
219 { |
|
|
220 User_Qinit = (mxArray *) pargin [1] ; |
|
|
221 } |
|
|
222 if ((op == 1 && nargin == 4) || (op == 2 && nargin == 5)) |
|
|
223 { |
|
|
224 User_Control_matrix = (mxArray *) pargin [nargin-1] ; |
|
|
225 } |
|
|
226 if (nargin < 3 || nargin > 5 || nargout > 2) |
|
|
227 { |
|
|
228 mexErrMsgTxt ("wrong number of arguments") ; |
|
|
229 } |
|
|
230 |
|
|
231 } |
|
|
232 else if (STRING_MATCH (operator, "/")) |
|
|
233 { |
|
|
234 |
|
|
235 /* -------------------------------------------------------------- */ |
|
|
236 /* matrix right divide, x = b/A */ |
|
|
237 /* -------------------------------------------------------------- */ |
|
|
238 |
|
|
239 /* |
|
|
240 [x, Info] = umfpack (b, '/', A) ; |
|
|
241 [x, Info] = umfpack (b, '/', A, Control) ; |
|
|
242 [x, Info] = umfpack (b, '/', A, Qinit) ; |
|
|
243 [x, Info] = umfpack (b, '/', A, Qinit, Control) ; |
|
|
244 */ |
|
|
245 |
|
|
246 operation = "x = b/A" ; |
|
|
247 do_solve = TRUE ; |
|
|
248 transpose = TRUE ; |
|
|
249 Amatrix = (mxArray *) pargin [2] ; |
|
|
250 Bmatrix = (mxArray *) pargin [0] ; |
|
|
251 |
|
|
252 if (nargout == 2) |
|
|
253 { |
|
|
254 do_info = 1 ; |
|
|
255 } |
|
|
256 if (nargin == 5) |
|
|
257 { |
|
|
258 User_Qinit = (mxArray *) pargin [3] ; |
|
|
259 User_Control_matrix = (mxArray *) pargin [4] ; |
|
|
260 } |
|
|
261 else if (nargin == 4) |
|
|
262 { |
|
|
263 /* Control is k-by-1 where k > 1, Qinit is 1-by-n */ |
|
|
264 if (mxGetM (pargin [3]) == 1) |
|
|
265 { |
|
|
266 User_Qinit = (mxArray *) pargin [3] ; |
|
|
267 } |
|
|
268 else |
|
|
269 { |
|
|
270 User_Control_matrix = (mxArray *) pargin [3] ; |
|
|
271 } |
|
|
272 } |
|
|
273 else if (nargin < 3 || nargin > 5 || nargout > 2) |
|
|
274 { |
|
|
275 mexErrMsgTxt ("wrong number of arguments") ; |
|
|
276 } |
|
|
277 |
|
|
278 } |
|
|
279 else if (STRING_MATCH (operator, "symbolic")) |
|
|
280 { |
|
|
281 |
|
|
282 /* -------------------------------------------------------------- */ |
|
|
283 /* symbolic factorization only */ |
|
|
284 /* -------------------------------------------------------------- */ |
|
|
285 |
|
|
286 /* |
|
|
287 [P Q Fr Ch Info] = umfpack (A, 'symbolic') ; |
|
|
288 [P Q Fr Ch Info] = umfpack (A, 'symbolic', Control) ; |
|
|
289 [P Q Fr Ch Info] = umfpack (A, Qinit, 'symbolic') ; |
|
|
290 [P Q Fr Ch Info] = umfpack (A, Qinit, 'symbolic', Control) ; |
|
|
291 */ |
|
|
292 |
|
|
293 operation = "symbolic factorization" ; |
|
|
294 do_numeric = FALSE ; |
|
|
295 Amatrix = (mxArray *) pargin [0] ; |
|
|
296 |
|
|
297 if (nargout == 5) |
|
|
298 { |
|
|
299 do_info = 4 ; |
|
|
300 } |
|
|
301 if (op == 2) |
|
|
302 { |
|
|
303 User_Qinit = (mxArray *) pargin [1] ; |
|
|
304 } |
|
|
305 if ((op == 1 && nargin == 3) || (op == 2 && nargin == 4)) |
|
|
306 { |
|
|
307 User_Control_matrix = (mxArray *) pargin [nargin-1] ; |
|
|
308 } |
|
|
309 if (nargin < 2 || nargin > 4 || nargout > 5 || nargout < 4) |
|
|
310 { |
|
|
311 mexErrMsgTxt ("wrong number of arguments") ; |
|
|
312 } |
|
|
313 |
|
|
314 } |
|
|
315 else if (STRING_MATCH (operator, "det")) |
|
|
316 { |
|
|
317 |
|
|
318 /* -------------------------------------------------------------- */ |
|
|
319 /* compute the determinant */ |
|
|
320 /* -------------------------------------------------------------- */ |
|
|
321 |
|
|
322 /* |
|
|
323 * [det] = umfpack (A, 'det') ; |
|
|
324 * [dmantissa dexp] = umfpack (A, 'det') ; |
|
|
325 */ |
|
|
326 |
|
|
327 operation = "determinant" ; |
|
|
328 do_det = TRUE ; |
|
|
329 Amatrix = (mxArray *) pargin [0] ; |
|
|
330 if (nargin > 2 || nargout > 2) |
|
|
331 { |
|
|
332 mexErrMsgTxt ("wrong number of arguments") ; |
|
|
333 } |
|
|
334 |
|
|
335 } |
|
|
336 else |
|
|
337 { |
|
|
338 mexErrMsgTxt ("operator must be '/', '\\', or 'symbolic'") ; |
|
|
339 } |
|
|
340 mxFree (operator) ; |
|
|
341 |
|
|
342 } |
|
|
343 else if (nargin > 0) |
|
|
344 { |
|
|
345 |
|
|
346 /* ------------------------------------------------------------------ */ |
|
|
347 /* LU factorization */ |
|
|
348 /* ------------------------------------------------------------------ */ |
|
|
349 |
|
|
350 /* |
|
|
351 with scaling: |
|
|
352 [L, U, P, Q, R, Info] = umfpack (A) ; |
|
|
353 [L, U, P, Q, R, Info] = umfpack (A, Qinit) ; |
|
|
354 |
|
|
355 scaling determined by Control settings: |
|
|
356 [L, U, P, Q, R, Info] = umfpack (A, Control) ; |
|
|
357 [L, U, P, Q, R, Info] = umfpack (A, Qinit, Control) ; |
|
|
358 |
|
|
359 with no scaling: |
|
|
360 [L, U, P, Q] = umfpack (A) ; |
|
|
361 [L, U, P, Q] = umfpack (A, Control) ; |
|
|
362 [L, U, P, Q] = umfpack (A, Qinit) ; |
|
|
363 [L, U, P, Q] = umfpack (A, Qinit, Control) ; |
|
|
364 */ |
|
|
365 |
|
|
366 operation = "numeric factorization" ; |
|
|
367 Amatrix = (mxArray *) pargin [0] ; |
|
|
368 |
|
|
369 no_scale = nargout <= 4 ; |
|
|
370 |
|
|
371 if (nargout == 6) |
|
|
372 { |
|
|
373 do_info = 5 ; |
|
|
374 } |
|
|
375 if (nargin == 3) |
|
|
376 { |
|
|
377 User_Qinit = (mxArray *) pargin [1] ; |
|
|
378 User_Control_matrix = (mxArray *) pargin [2] ; |
|
|
379 } |
|
|
380 else if (nargin == 2) |
|
|
381 { |
|
|
382 /* Control is k-by-1 where k > 1, Qinit is 1-by-n */ |
|
|
383 if (mxGetM (pargin [1]) == 1) |
|
|
384 { |
|
|
385 User_Qinit = (mxArray *) pargin [1] ; |
|
|
386 } |
|
|
387 else |
|
|
388 { |
|
|
389 User_Control_matrix = (mxArray *) pargin [1] ; |
|
|
390 } |
|
|
391 } |
|
|
392 else if (nargin > 3 || nargout > 6 || nargout < 4) |
|
|
393 { |
|
|
394 mexErrMsgTxt ("wrong number of arguments") ; |
|
|
395 } |
|
|
396 |
|
|
397 } |
|
|
398 else |
|
|
399 { |
|
|
400 |
|
|
401 /* ------------------------------------------------------------------ */ |
|
|
402 /* return default control settings */ |
|
|
403 /* ------------------------------------------------------------------ */ |
|
|
404 |
|
|
405 /* |
|
|
406 Control = umfpack ; |
|
|
407 umfpack ; |
|
|
408 */ |
|
|
409 |
|
|
410 if (nargout > 1) |
|
|
411 { |
|
|
412 mexErrMsgTxt ("wrong number of arguments") ; |
|
|
413 } |
|
|
414 |
|
|
415 pargout [0] = mxCreateDoubleMatrix (UMFPACK_CONTROL, 1, mxREAL) ; |
|
|
416 User_Control = mxGetPr (pargout [0]) ; |
|
|
417 umfpack_di_defaults (User_Control) ; |
|
|
418 |
|
|
419 return ; |
|
|
420 } |
|
|
421 |
|
|
422 /* ---------------------------------------------------------------------- */ |
|
|
423 /* check inputs */ |
|
|
424 /* ---------------------------------------------------------------------- */ |
|
|
425 |
|
|
426 if (mxGetNumberOfDimensions (Amatrix) != 2) |
|
|
427 { |
|
|
428 mexErrMsgTxt ("input matrix A must be 2-dimensional") ; |
|
|
429 } |
|
|
430 n_row = mxGetM (Amatrix) ; |
|
|
431 n_col = mxGetN (Amatrix) ; |
|
|
432 nn = MAX (n_row, n_col) ; |
|
|
433 n_inner = MIN (n_row, n_col) ; |
|
|
434 if (do_solve && n_row != n_col) |
|
|
435 { |
|
|
436 mexErrMsgTxt ("input matrix A must square for '\\' or '/'") ; |
|
|
437 } |
|
|
438 if (!mxIsSparse (Amatrix)) |
|
|
439 { |
|
|
440 mexErrMsgTxt ("input matrix A must be sparse") ; |
|
|
441 } |
|
|
442 if (n_row == 0 || n_col == 0) |
|
|
443 { |
|
|
444 mexErrMsgTxt ("input matrix A cannot have zero rows or zero columns") ; |
|
|
445 } |
|
|
446 |
|
|
447 /* The real/complex status of A determines which version to use, */ |
|
|
448 /* (umfpack_di_* or umfpack_zi_*). */ |
|
|
449 A_is_complex = mxIsComplex (Amatrix) ; |
|
|
450 Atype = A_is_complex ? mxCOMPLEX : mxREAL ; |
|
|
451 Ap = mxGetJc (Amatrix) ; |
|
|
452 Ai = mxGetIr (Amatrix) ; |
|
|
453 Ax = mxGetPr (Amatrix) ; |
|
|
454 Az = mxGetPi (Amatrix) ; |
|
|
455 |
|
|
456 if (do_solve) |
|
|
457 { |
|
|
458 |
|
|
459 if (n_row != n_col) |
|
|
460 { |
|
|
461 mexErrMsgTxt ("A must be square for \\ or /") ; |
|
|
462 } |
|
|
463 if (transpose) |
|
|
464 { |
|
|
465 if (mxGetM (Bmatrix) != 1 || mxGetN (Bmatrix) != nn) |
|
|
466 { |
|
|
467 mexErrMsgTxt ("b has the wrong dimensions") ; |
|
|
468 } |
|
|
469 } |
|
|
470 else |
|
|
471 { |
|
|
472 if (mxGetM (Bmatrix) != nn || mxGetN (Bmatrix) != 1) |
|
|
473 { |
|
|
474 mexErrMsgTxt ("b has the wrong dimensions") ; |
|
|
475 } |
|
|
476 } |
|
|
477 if (mxGetNumberOfDimensions (Bmatrix) != 2) |
|
|
478 { |
|
|
479 mexErrMsgTxt ("input matrix b must be 2-dimensional") ; |
|
|
480 } |
|
|
481 if (mxIsSparse (Bmatrix)) |
|
|
482 { |
|
|
483 mexErrMsgTxt ("input matrix b cannot be sparse") ; |
|
|
484 } |
|
|
485 if (mxGetClassID (Bmatrix) != mxDOUBLE_CLASS) |
|
|
486 { |
|
|
487 mexErrMsgTxt ("input matrix b must double precision matrix") ; |
|
|
488 } |
|
|
489 |
|
|
490 B_is_complex = mxIsComplex (Bmatrix) ; |
|
|
491 Bx = mxGetPr (Bmatrix) ; |
|
|
492 Bz = mxGetPi (Bmatrix) ; |
|
|
493 |
|
|
494 X_is_complex = A_is_complex || B_is_complex ; |
|
|
495 Xtype = X_is_complex ? mxCOMPLEX : mxREAL ; |
|
|
496 } |
|
|
497 |
|
|
498 /* ---------------------------------------------------------------------- */ |
|
|
499 /* set the Control parameters */ |
|
|
500 /* ---------------------------------------------------------------------- */ |
|
|
501 |
|
|
502 if (A_is_complex) |
|
|
503 { |
|
|
504 umfpack_zi_defaults (Control) ; |
|
|
505 } |
|
|
506 else |
|
|
507 { |
|
|
508 umfpack_di_defaults (Control) ; |
|
|
509 } |
|
|
510 if (User_Control_matrix) |
|
|
511 { |
|
|
512 if (mxGetClassID (User_Control_matrix) != mxDOUBLE_CLASS || |
|
|
513 mxIsSparse (User_Control_matrix)) |
|
|
514 { |
|
|
515 mexErrMsgTxt ("Control must be a dense real matrix") ; |
|
|
516 } |
|
|
517 size = UMFPACK_CONTROL ; |
|
|
518 size = MIN (size, mxGetNumberOfElements (User_Control_matrix)) ; |
|
|
519 User_Control = mxGetPr (User_Control_matrix) ; |
|
|
520 for (i = 0 ; i < size ; i++) |
|
|
521 { |
|
|
522 Control [i] = User_Control [i] ; |
|
|
523 } |
|
|
524 } |
|
|
525 |
|
|
526 if (no_scale) |
|
|
527 { |
|
|
528 /* turn off scaling for [L, U, P, Q] = umfpack (A) ; |
|
|
529 * ignoring the input value of Control (24) for the usage |
|
|
530 * [L, U, P, Q] = umfpack (A, Control) ; */ |
|
|
531 Control [UMFPACK_SCALE] = UMFPACK_SCALE_NONE ; |
|
|
532 } |
|
|
533 |
|
|
534 if (mxIsNaN (Control [UMFPACK_PRL])) |
|
|
535 { |
|
|
536 print_level = UMFPACK_DEFAULT_PRL ; |
|
|
537 } |
|
|
538 else |
|
|
539 { |
|
|
540 print_level = (int) Control [UMFPACK_PRL] ; |
|
|
541 } |
|
|
542 |
|
|
543 Control [UMFPACK_PRL] = print_level ; |
|
|
544 |
|
|
545 /* ---------------------------------------------------------------------- */ |
|
|
546 /* get Qinit, if present */ |
|
|
547 /* ---------------------------------------------------------------------- */ |
|
|
548 |
|
|
549 if (User_Qinit) |
|
|
550 { |
|
|
551 if (mxGetM (User_Qinit) != 1 || mxGetN (User_Qinit) != n_col) |
|
|
552 { |
|
|
553 mexErrMsgTxt ("Qinit must be 1-by-n_col") ; |
|
|
554 } |
|
|
555 if (mxGetNumberOfDimensions (User_Qinit) != 2) |
|
|
556 { |
|
|
557 mexErrMsgTxt ("input Qinit must be 2-dimensional") ; |
|
|
558 } |
|
|
559 if (mxIsComplex (User_Qinit)) |
|
|
560 { |
|
|
561 mexErrMsgTxt ("input Qinit must not be complex") ; |
|
|
562 } |
|
|
563 if (mxGetClassID (User_Qinit) != mxDOUBLE_CLASS) |
|
|
564 { |
|
|
565 mexErrMsgTxt ("input Qinit must be a double matrix") ; |
|
|
566 } |
|
|
567 if (mxIsSparse (User_Qinit)) |
|
|
568 { |
|
|
569 mexErrMsgTxt ("input Qinit must be dense") ; |
|
|
570 } |
|
|
571 Qinit = (int *) mxMalloc (n_col * sizeof (int)) ; |
|
|
572 p = mxGetPr (User_Qinit) ; |
|
|
573 for (k = 0 ; k < n_col ; k++) |
|
|
574 { |
|
|
575 /* convert from 1-based to 0-based indexing */ |
|
|
576 Qinit [k] = ((int) (p [k])) - 1 ; |
|
|
577 } |
|
|
578 |
|
|
579 } |
|
|
580 else |
|
|
581 { |
|
|
582 /* umfpack_*_qsymbolic will call colamd to get Qinit. This is the */ |
|
|
583 /* same as calling umfpack_*_symbolic with Qinit set to NULL*/ |
|
|
584 Qinit = (int *) NULL ; |
|
|
585 } |
|
|
586 |
|
|
587 /* ---------------------------------------------------------------------- */ |
|
|
588 /* report the inputs A and Qinit */ |
|
|
589 /* ---------------------------------------------------------------------- */ |
|
|
590 |
|
|
591 if (print_level >= 2) |
|
|
592 { |
|
|
593 /* print the operation */ |
|
|
594 mexPrintf ("\numfpack: %s\n", operation) ; |
|
|
595 } |
|
|
596 |
|
|
597 if (A_is_complex) |
|
|
598 { |
|
|
599 umfpack_zi_report_control (Control) ; |
|
|
600 if (print_level >= 3) mexPrintf ("\nA: ") ; |
|
|
601 (void) umfpack_zi_report_matrix (n_row, n_col, Ap, Ai, Ax, Az, |
|
|
602 1, Control) ; |
|
|
603 if (Qinit) |
|
|
604 { |
|
|
605 if (print_level >= 3) mexPrintf ("\nQinit: ") ; |
|
|
606 (void) umfpack_zi_report_perm (n_col, Qinit, Control) ; |
|
|
607 } |
|
|
608 } |
|
|
609 else |
|
|
610 { |
|
|
611 umfpack_di_report_control (Control) ; |
|
|
612 if (print_level >= 3) mexPrintf ("\nA: ") ; |
|
|
613 (void) umfpack_di_report_matrix (n_row, n_col, Ap, Ai, Ax, |
|
|
614 1, Control) ; |
|
|
615 if (Qinit) |
|
|
616 { |
|
|
617 if (print_level >= 3) mexPrintf ("\nQinit: ") ; |
|
|
618 (void) umfpack_di_report_perm (n_col, Qinit, Control) ; |
|
|
619 } |
|
|
620 } |
|
|
621 |
|
|
622 #ifndef NO_TRANSPOSE_FORWARD_SLASH |
|
|
623 /* ---------------------------------------------------------------------- */ |
|
|
624 /* create the array transpose for x = b/A */ |
|
|
625 /* ---------------------------------------------------------------------- */ |
|
|
626 |
|
|
627 if (transpose) |
|
|
628 { |
|
|
629 /* note that in this case A will be square (nn = n_row = n_col) */ |
|
|
630 /* x = (A.'\b.').' will be computed */ |
|
|
631 |
|
|
632 /* make sure Ci and Cx exist, avoid malloc of zero-sized arrays. */ |
|
|
633 nz = MAX (Ap [nn], 1) ; |
|
|
634 |
|
|
635 Cp = (int *) mxMalloc ((nn+1) * sizeof (int)) ; |
|
|
636 Ci = (int *) mxMalloc (nz * sizeof (int)) ; |
|
|
637 Cx = (double *) mxMalloc (nz * sizeof (double)) ; |
|
|
638 if (A_is_complex) |
|
|
639 { |
|
|
640 Cz = (double *) mxMalloc (nz * sizeof (double)) ; |
|
|
641 status = umfpack_zi_transpose (nn, nn, Ap, Ai, Ax, Az, |
|
|
642 (int *) NULL, (int *) NULL, Cp, Ci, Cx, Cz, FALSE) ; |
|
|
643 } |
|
|
644 else |
|
|
645 { |
|
|
646 status = umfpack_di_transpose (nn, nn, Ap, Ai, Ax, |
|
|
647 (int *) NULL, (int *) NULL, Cp, Ci, Cx) ; |
|
|
648 } |
|
|
649 |
|
|
650 if (status != UMFPACK_OK) |
|
|
651 { |
|
|
652 error ("transpose of A failed", A_is_complex, nargout, pargout, |
|
|
653 Control, Info, status, do_info); |
|
|
654 return ; |
|
|
655 } |
|
|
656 |
|
|
657 /* modify pointers so that C will be factorized and solved, not A */ |
|
|
658 Ap = Cp ; |
|
|
659 Ai = Ci ; |
|
|
660 Ax = Cx ; |
|
|
661 if (A_is_complex) |
|
|
662 { |
|
|
663 Az = Cz ; |
|
|
664 } |
|
|
665 } |
|
|
666 #endif |
|
|
667 |
|
|
668 /* ---------------------------------------------------------------------- */ |
|
|
669 /* perform the symbolic factorization */ |
|
|
670 /* ---------------------------------------------------------------------- */ |
|
|
671 |
|
|
672 if (A_is_complex) |
|
|
673 { |
|
|
674 status = umfpack_zi_qsymbolic (n_row, n_col, Ap, Ai, Ax, Az, |
|
|
675 Qinit, &Symbolic, Control, Info) ; |
|
|
676 } |
|
|
677 else |
|
|
678 { |
|
|
679 status = umfpack_di_qsymbolic (n_row, n_col, Ap, Ai, Ax, |
|
|
680 Qinit, &Symbolic, Control, Info) ; |
|
|
681 } |
|
|
682 |
|
|
683 if (Qinit) |
|
|
684 { |
|
|
685 mxFree (Qinit) ; |
|
|
686 } |
|
|
687 |
|
|
688 if (status < 0) |
|
|
689 { |
|
|
690 error ("symbolic factorization failed", A_is_complex, nargout, pargout, |
|
|
691 Control, Info, status, do_info) ; |
|
|
692 return ; |
|
|
693 } |
|
|
694 |
|
|
695 /* ---------------------------------------------------------------------- */ |
|
|
696 /* report the Symbolic object */ |
|
|
697 /* ---------------------------------------------------------------------- */ |
|
|
698 |
|
|
699 if (A_is_complex) |
|
|
700 { |
|
|
701 (void) umfpack_zi_report_symbolic (Symbolic, Control) ; |
|
|
702 } |
|
|
703 else |
|
|
704 { |
|
|
705 (void) umfpack_di_report_symbolic (Symbolic, Control) ; |
|
|
706 } |
|
|
707 |
|
|
708 /* ---------------------------------------------------------------------- */ |
|
|
709 /* perform numeric factorization, or just return symbolic factorization */ |
|
|
710 /* ---------------------------------------------------------------------- */ |
|
|
711 |
|
|
712 if (do_numeric) |
|
|
713 { |
|
|
714 |
|
|
715 /* ------------------------------------------------------------------ */ |
|
|
716 /* perform the numeric factorization */ |
|
|
717 /* ------------------------------------------------------------------ */ |
|
|
718 |
|
|
719 if (A_is_complex) |
|
|
720 { |
|
|
721 status = umfpack_zi_numeric (Ap, Ai, Ax, Az, Symbolic, &Numeric, |
|
|
722 Control, Info) ; |
|
|
723 } |
|
|
724 else |
|
|
725 { |
|
|
726 status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, |
|
|
727 Control, Info) ; |
|
|
728 } |
|
|
729 |
|
|
730 /* ------------------------------------------------------------------ */ |
|
|
731 /* free the symbolic factorization */ |
|
|
732 /* ------------------------------------------------------------------ */ |
|
|
733 |
|
|
734 if (A_is_complex) |
|
|
735 { |
|
|
736 umfpack_zi_free_symbolic (&Symbolic) ; |
|
|
737 } |
|
|
738 else |
|
|
739 { |
|
|
740 umfpack_di_free_symbolic (&Symbolic) ; |
|
|
741 } |
|
|
742 |
|
|
743 /* ------------------------------------------------------------------ */ |
|
|
744 /* report the Numeric object */ |
|
|
745 /* ------------------------------------------------------------------ */ |
|
|
746 |
|
|
747 if (status < 0) |
|
|
748 { |
|
|
749 error ("numeric factorization failed", A_is_complex, nargout, |
|
|
750 pargout, Control, Info, status, do_info); |
|
|
751 return ; |
|
|
752 } |
|
|
753 |
|
|
754 if (A_is_complex) |
|
|
755 { |
|
|
756 (void) umfpack_zi_report_numeric (Numeric, Control) ; |
|
|
757 } |
|
|
758 else |
|
|
759 { |
|
|
760 (void) umfpack_di_report_numeric (Numeric, Control) ; |
|
|
761 } |
|
|
762 |
|
|
763 /* ------------------------------------------------------------------ */ |
|
|
764 /* return the solution, determinant, or the factorization */ |
|
|
765 /* ------------------------------------------------------------------ */ |
|
|
766 |
|
|
767 if (do_solve) |
|
|
768 { |
|
|
769 /* -------------------------------------------------------------- */ |
|
|
770 /* solve Ax=b or A'x'=b', and return just the solution x */ |
|
|
771 /* -------------------------------------------------------------- */ |
|
|
772 |
|
|
773 #ifndef NO_TRANSPOSE_FORWARD_SLASH |
|
|
774 if (transpose) |
|
|
775 { |
|
|
776 /* A.'x.'=b.' gives the same x=b/A as solving A'x'=b' */ |
|
|
777 /* since C=A.' was factorized, solve with sys = UMFPACK_A */ |
|
|
778 /* since x and b are vectors, x.' and b.' are implicit */ |
|
|
779 pargout [0] = mxCreateDoubleMatrix (1, nn, Xtype) ; |
|
|
780 } |
|
|
781 else |
|
|
782 { |
|
|
783 pargout [0] = mxCreateDoubleMatrix (nn, 1, Xtype) ; |
|
|
784 } |
|
|
785 sys = UMFPACK_A ; |
|
|
786 #else |
|
|
787 if (transpose) |
|
|
788 { |
|
|
789 /* If A is real, A'x=b is the same as A.'x=b. */ |
|
|
790 /* x and b are vectors, so x and b are the same as x' and b'. */ |
|
|
791 /* If A is complex, then A.'x.'=b.' gives the same solution x */ |
|
|
792 /* as the complex conjugate transpose. If we used the A'x=b */ |
|
|
793 /* option in umfpack_*_solve, we would have to form b' on */ |
|
|
794 /* input and x' on output (negating the imaginary part). */ |
|
|
795 /* We can save this work by just using the A.'x=b option in */ |
|
|
796 /* umfpack_*_solve. Then, forming x.' and b.' is implicit, */ |
|
|
797 /* since x and b are just vectors anyway. */ |
|
|
798 /* In both cases, the system to solve is A.'x=b */ |
|
|
799 pargout [0] = mxCreateDoubleMatrix (1, nn, Xtype) ; |
|
|
800 sys = UMFPACK_Aat ; |
|
|
801 } |
|
|
802 else |
|
|
803 { |
|
|
804 pargout [0] = mxCreateDoubleMatrix (nn, 1, Xtype) ; |
|
|
805 sys = UMFPACK_A ; |
|
|
806 } |
|
|
807 #endif |
|
|
808 |
|
|
809 /* -------------------------------------------------------------- */ |
|
|
810 /* print the right-hand-side, B */ |
|
|
811 /* -------------------------------------------------------------- */ |
|
|
812 |
|
|
813 if (print_level >= 3) mexPrintf ("\nright-hand side, b: ") ; |
|
|
814 if (B_is_complex) |
|
|
815 { |
|
|
816 (void) umfpack_zi_report_vector (nn, Bx, Bz, Control) ; |
|
|
817 } |
|
|
818 else |
|
|
819 { |
|
|
820 (void) umfpack_di_report_vector (nn, Bx, Control) ; |
|
|
821 } |
|
|
822 |
|
|
823 /* -------------------------------------------------------------- */ |
|
|
824 /* solve the system */ |
|
|
825 /* -------------------------------------------------------------- */ |
|
|
826 |
|
|
827 Xx = mxGetPr (pargout [0]) ; |
|
|
828 Xz = mxGetPi (pargout [0]) ; |
|
|
829 status2 = UMFPACK_OK ; |
|
|
830 |
|
|
831 if (A_is_complex) |
|
|
832 { |
|
|
833 if (!B_is_complex) |
|
|
834 { |
|
|
835 /* umfpack_zi_solve expects a complex B */ |
|
|
836 Bz = (double *) mxCalloc (nn, sizeof (double)) ; |
|
|
837 } |
|
|
838 status = umfpack_zi_solve (sys, Ap, Ai, Ax, Az, Xx, Xz, Bx, Bz, |
|
|
839 Numeric, Control, Info) ; |
|
|
840 if (!B_is_complex) |
|
|
841 { |
|
|
842 mxFree (Bz) ; |
|
|
843 } |
|
|
844 } |
|
|
845 else |
|
|
846 { |
|
|
847 if (B_is_complex) |
|
|
848 { |
|
|
849 /* Ax=b when b is complex and A is sparse can be split */ |
|
|
850 /* into two systems, A*xr=br and A*xi=bi, where r denotes */ |
|
|
851 /* the real part and i the imaginary part of x and b. */ |
|
|
852 status2 = umfpack_di_solve (sys, Ap, Ai, Ax, Xz, Bz, |
|
|
853 Numeric, Control, Info) ; |
|
|
854 } |
|
|
855 status = umfpack_di_solve (sys, Ap, Ai, Ax, Xx, Bx, |
|
|
856 Numeric, Control, Info) ; |
|
|
857 } |
|
|
858 |
|
|
859 #ifndef NO_TRANSPOSE_FORWARD_SLASH |
|
|
860 /* -------------------------------------------------------------- */ |
|
|
861 /* free the transposed matrix C */ |
|
|
862 /* -------------------------------------------------------------- */ |
|
|
863 |
|
|
864 if (transpose) |
|
|
865 { |
|
|
866 mxFree (Cp) ; |
|
|
867 mxFree (Ci) ; |
|
|
868 mxFree (Cx) ; |
|
|
869 if (A_is_complex) |
|
|
870 { |
|
|
871 mxFree (Cz) ; |
|
|
872 } |
|
|
873 } |
|
|
874 #endif |
|
|
875 |
|
|
876 /* -------------------------------------------------------------- */ |
|
|
877 /* free the Numeric object */ |
|
|
878 /* -------------------------------------------------------------- */ |
|
|
879 |
|
|
880 if (A_is_complex) |
|
|
881 { |
|
|
882 umfpack_zi_free_numeric (&Numeric) ; |
|
|
883 } |
|
|
884 else |
|
|
885 { |
|
|
886 umfpack_di_free_numeric (&Numeric) ; |
|
|
887 } |
|
|
888 |
|
|
889 /* -------------------------------------------------------------- */ |
|
|
890 /* check error status */ |
|
|
891 /* -------------------------------------------------------------- */ |
|
|
892 |
|
|
893 if (status < 0 || status2 < 0) |
|
|
894 { |
|
|
895 mxDestroyArray (pargout [0]) ; |
|
|
896 error ("solve failed", A_is_complex, nargout, pargout, Control, |
|
|
897 Info, status, do_info) ; |
|
|
898 return ; |
|
|
899 } |
|
|
900 |
|
|
901 /* -------------------------------------------------------------- */ |
|
|
902 /* print the solution, X */ |
|
|
903 /* -------------------------------------------------------------- */ |
|
|
904 |
|
|
905 if (print_level >= 3) mexPrintf ("\nsolution, x: ") ; |
|
|
906 if (X_is_complex) |
|
|
907 { |
|
|
908 (void) umfpack_zi_report_vector (nn, Xx, Xz, Control) ; |
|
|
909 } |
|
|
910 else |
|
|
911 { |
|
|
912 (void) umfpack_di_report_vector (nn, Xx, Control) ; |
|
|
913 } |
|
|
914 |
|
|
915 /* -------------------------------------------------------------- */ |
|
|
916 /* warn about singular or near-singular matrices */ |
|
|
917 /* -------------------------------------------------------------- */ |
|
|
918 |
|
|
919 /* no warning is given if Control (1) is zero */ |
|
|
920 |
|
|
921 if (Control [UMFPACK_PRL] >= 1) |
|
|
922 { |
|
|
923 if (status == UMFPACK_WARNING_singular_matrix) |
|
|
924 { |
|
|
925 sprintf (warning, "matrix is singular\n" |
|
|
926 "Try increasing Control (%d) and Control (%d).\n" |
|
|
927 "(Suppress this warning with Control (%d) = 0.)\n", |
|
|
928 1+UMFPACK_PIVOT_TOLERANCE, |
|
|
929 1+UMFPACK_SYM_PIVOT_TOLERANCE, |
|
|
930 1+UMFPACK_PRL) ; |
|
|
931 mexWarnMsgTxt (warning) ; |
|
|
932 } |
|
|
933 else if (Info [UMFPACK_RCOND] < DBL_EPSILON) |
|
|
934 { |
|
|
935 sprintf (warning, "matrix is nearly singular, rcond = %g\n" |
|
|
936 "Try increasing Control (%d) and Control (%d).\n" |
|
|
937 "(Suppress this warning with Control (%d) = 0.)\n", |
|
|
938 Info [UMFPACK_RCOND], |
|
|
939 1+UMFPACK_PIVOT_TOLERANCE, |
|
|
940 1+UMFPACK_SYM_PIVOT_TOLERANCE, |
|
|
941 1+UMFPACK_PRL) ; |
|
|
942 mexWarnMsgTxt (warning) ; |
|
|
943 } |
|
|
944 } |
|
|
945 |
|
|
946 } |
|
|
947 else if (do_det) |
|
|
948 { |
|
|
949 |
|
|
950 /* -------------------------------------------------------------- */ |
|
|
951 /* get the determinant */ |
|
|
952 /* -------------------------------------------------------------- */ |
|
|
953 |
|
|
954 if (nargout == 2) |
|
|
955 { |
|
|
956 /* [det dexp] = umfpack (A, 'det') ; |
|
|
957 * return determinant in the form det * 10^dexp */ |
|
|
958 p = &dexp ; |
|
|
959 } |
|
|
960 else |
|
|
961 { |
|
|
962 /* [det] = umfpack (A, 'det') ; |
|
|
963 * return determinant as a single scalar (overflow or |
|
|
964 * underflow is much more likely) */ |
|
|
965 p = (double *) NULL ; |
|
|
966 } |
|
|
967 if (A_is_complex) |
|
|
968 { |
|
|
969 status = umfpack_zi_get_determinant (&dx, &dz, p, |
|
|
970 Numeric, Info) ; |
|
|
971 umfpack_zi_free_numeric (&Numeric) ; |
|
|
972 } |
|
|
973 else |
|
|
974 { |
|
|
975 status = umfpack_di_get_determinant (&dx, p, |
|
|
976 Numeric, Info) ; |
|
|
977 umfpack_di_free_numeric (&Numeric) ; |
|
|
978 dz = 0 ; |
|
|
979 } |
|
|
980 if (status < 0) |
|
|
981 { |
|
|
982 error ("extracting LU factors failed", A_is_complex, nargout, |
|
|
983 pargout, Control, Info, status, do_info) ; |
|
|
984 } |
|
|
985 if (A_is_complex) |
|
|
986 { |
|
|
987 pargout [0] = mxCreateDoubleMatrix (1, 1, mxCOMPLEX) ; |
|
|
988 p = mxGetPr (pargout [0]) ; |
|
|
989 *p = dx ; |
|
|
990 p = mxGetPi (pargout [0]) ; |
|
|
991 *p = dz ; |
|
|
992 } |
|
|
993 else |
|
|
994 { |
|
|
995 pargout [0] = mxCreateDoubleMatrix (1, 1, mxREAL) ; |
|
|
996 p = mxGetPr (pargout [0]) ; |
|
|
997 *p = dx ; |
|
|
998 } |
|
|
999 if (nargout == 2) |
|
|
1000 { |
|
|
1001 pargout [1] = mxCreateDoubleMatrix (1, 1, mxREAL) ; |
|
|
1002 p = mxGetPr (pargout [1]) ; |
|
|
1003 *p = dexp ; |
|
|
1004 } |
|
|
1005 |
|
|
1006 } |
|
|
1007 else |
|
|
1008 { |
|
|
1009 |
|
|
1010 /* -------------------------------------------------------------- */ |
|
|
1011 /* get L, U, P, Q, and r */ |
|
|
1012 /* -------------------------------------------------------------- */ |
|
|
1013 |
|
|
1014 if (A_is_complex) |
|
|
1015 { |
|
|
1016 status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1, &ignore2, |
|
|
1017 &ignore3, Numeric) ; |
|
|
1018 } |
|
|
1019 else |
|
|
1020 { |
|
|
1021 status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2, |
|
|
1022 &ignore3, Numeric) ; |
|
|
1023 } |
|
|
1024 |
|
|
1025 if (status < 0) |
|
|
1026 { |
|
|
1027 if (A_is_complex) |
|
|
1028 { |
|
|
1029 umfpack_zi_free_numeric (&Numeric) ; |
|
|
1030 } |
|
|
1031 else |
|
|
1032 { |
|
|
1033 umfpack_di_free_numeric (&Numeric) ; |
|
|
1034 } |
|
|
1035 error ("extracting LU factors failed", A_is_complex, nargout, |
|
|
1036 pargout, Control, Info, status, do_info) ; |
|
|
1037 return ; |
|
|
1038 } |
|
|
1039 |
|
|
1040 /* avoid malloc of zero-sized arrays */ |
|
|
1041 lnz = MAX (lnz, 1) ; |
|
|
1042 unz = MAX (unz, 1) ; |
|
|
1043 |
|
|
1044 /* get temporary space, for the *** ROW *** form of L */ |
|
|
1045 Ltp = (int *) mxMalloc ((n_row+1) * sizeof (int)) ; |
|
|
1046 Ltj = (int *) mxMalloc (lnz * sizeof (int)) ; |
|
|
1047 Ltx = (double *) mxMalloc (lnz * sizeof (double)) ; |
|
|
1048 if (A_is_complex) |
|
|
1049 { |
|
|
1050 Ltz = (double *) mxMalloc (lnz * sizeof (double)) ; |
|
|
1051 } |
|
|
1052 else |
|
|
1053 { |
|
|
1054 Ltz = (double *) NULL ; |
|
|
1055 } |
|
|
1056 |
|
|
1057 /* create permanent copy of the output matrix U */ |
|
|
1058 pargout [1] = mxCreateSparse (n_inner, n_col, unz, Atype) ; |
|
|
1059 Up = mxGetJc (pargout [1]) ; |
|
|
1060 Ui = mxGetIr (pargout [1]) ; |
|
|
1061 Ux = mxGetPr (pargout [1]) ; |
|
|
1062 Uz = mxGetPi (pargout [1]) ; |
|
|
1063 |
|
|
1064 /* temporary space for the integer permutation vectors */ |
|
|
1065 P = (int *) mxMalloc (n_row * sizeof (int)) ; |
|
|
1066 Q = (int *) mxMalloc (n_col * sizeof (int)) ; |
|
|
1067 |
|
|
1068 /* get scale factors, if requested */ |
|
|
1069 status2 = UMFPACK_OK ; |
|
|
1070 if (!no_scale) |
|
|
1071 { |
|
|
1072 /* create a diagonal sparse matrix for the scale factors */ |
|
|
1073 pargout [4] = mxCreateSparse (n_row, n_row, n_row, mxREAL) ; |
|
|
1074 Rp = mxGetJc (pargout [4]) ; |
|
|
1075 Ri = mxGetIr (pargout [4]) ; |
|
|
1076 for (i = 0 ; i < n_row ; i++) |
|
|
1077 { |
|
|
1078 Rp [i] = i ; |
|
|
1079 Ri [i] = i ; |
|
|
1080 } |
|
|
1081 Rp [n_row] = n_row ; |
|
|
1082 Rs = mxGetPr (pargout [4]) ; |
|
|
1083 } |
|
|
1084 else |
|
|
1085 { |
|
|
1086 Rs = (double *) NULL ; |
|
|
1087 } |
|
|
1088 |
|
|
1089 /* get Lt, U, P, Q, and Rs from the numeric object */ |
|
|
1090 if (A_is_complex) |
|
|
1091 { |
|
|
1092 status = umfpack_zi_get_numeric (Ltp, Ltj, Ltx, Ltz, Up, Ui, Ux, |
|
|
1093 Uz, P, Q, (double *) NULL, (double *) NULL, |
|
|
1094 &do_recip, Rs, Numeric) ; |
|
|
1095 umfpack_zi_free_numeric (&Numeric) ; |
|
|
1096 } |
|
|
1097 else |
|
|
1098 { |
|
|
1099 status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Ui, |
|
|
1100 Ux, P, Q, (double *) NULL, |
|
|
1101 &do_recip, Rs, Numeric) ; |
|
|
1102 umfpack_di_free_numeric (&Numeric) ; |
|
|
1103 } |
|
|
1104 |
|
|
1105 /* for the mexFunction, -DNRECIPROCAL must be set, |
|
|
1106 * so do_recip must be FALSE */ |
|
|
1107 |
|
|
1108 if (status < 0 || status2 < 0 || do_recip) |
|
|
1109 { |
|
|
1110 mxFree (Ltp) ; |
|
|
1111 mxFree (Ltj) ; |
|
|
1112 mxFree (Ltx) ; |
|
|
1113 if (Ltz) mxFree (Ltz) ; |
|
|
1114 mxFree (P) ; |
|
|
1115 mxFree (Q) ; |
|
|
1116 mxDestroyArray (pargout [1]) ; |
|
|
1117 error ("extracting LU factors failed", A_is_complex, nargout, |
|
|
1118 pargout, Control, Info, status, do_info) ; |
|
|
1119 return ; |
|
|
1120 } |
|
|
1121 |
|
|
1122 /* create sparse permutation matrix for P */ |
|
|
1123 pargout [2] = mxCreateSparse (n_row, n_row, n_row, mxREAL) ; |
|
|
1124 Pp = mxGetJc (pargout [2]) ; |
|
|
1125 Pi = mxGetIr (pargout [2]) ; |
|
|
1126 Px = mxGetPr (pargout [2]) ; |
|
|
1127 for (k = 0 ; k < n_row ; k++) |
|
|
1128 { |
|
|
1129 Pp [k] = k ; |
|
|
1130 Px [k] = 1 ; |
|
|
1131 Pi [P [k]] = k ; |
|
|
1132 } |
|
|
1133 Pp [n_row] = n_row ; |
|
|
1134 |
|
|
1135 /* create sparse permutation matrix for Q */ |
|
|
1136 pargout [3] = mxCreateSparse (n_col, n_col, n_col, mxREAL) ; |
|
|
1137 Qp = mxGetJc (pargout [3]) ; |
|
|
1138 Qi = mxGetIr (pargout [3]) ; |
|
|
1139 Qx = mxGetPr (pargout [3]) ; |
|
|
1140 for (k = 0 ; k < n_col ; k++) |
|
|
1141 { |
|
|
1142 Qp [k] = k ; |
|
|
1143 Qx [k] = 1 ; |
|
|
1144 Qi [k] = Q [k] ; |
|
|
1145 } |
|
|
1146 Qp [n_col] = n_col ; |
|
|
1147 |
|
|
1148 /* permanent copy of L */ |
|
|
1149 pargout [0] = mxCreateSparse (n_row, n_inner, lnz, Atype) ; |
|
|
1150 Lp = mxGetJc (pargout [0]) ; |
|
|
1151 Li = mxGetIr (pargout [0]) ; |
|
|
1152 Lx = mxGetPr (pargout [0]) ; |
|
|
1153 Lz = mxGetPi (pargout [0]) ; |
|
|
1154 |
|
|
1155 /* convert L from row form to column form */ |
|
|
1156 if (A_is_complex) |
|
|
1157 { |
|
|
1158 /* non-conjugate array transpose */ |
|
|
1159 status = umfpack_zi_transpose (n_inner, n_row, Ltp, Ltj, Ltx, |
|
|
1160 Ltz, (int *) NULL, (int *) NULL, Lp, Li, Lx, Lz, FALSE) ; |
|
|
1161 } |
|
|
1162 else |
|
|
1163 { |
|
|
1164 status = umfpack_di_transpose (n_inner, n_row, Ltp, Ltj, Ltx, |
|
|
1165 (int *) NULL, (int *) NULL, Lp, Li, Lx) ; |
|
|
1166 } |
|
|
1167 |
|
|
1168 mxFree (Ltp) ; |
|
|
1169 mxFree (Ltj) ; |
|
|
1170 mxFree (Ltx) ; |
|
|
1171 if (Ltz) mxFree (Ltz) ; |
|
|
1172 |
|
|
1173 if (status < 0) |
|
|
1174 { |
|
|
1175 mxFree (P) ; |
|
|
1176 mxFree (Q) ; |
|
|
1177 mxDestroyArray (pargout [0]) ; |
|
|
1178 mxDestroyArray (pargout [1]) ; |
|
|
1179 mxDestroyArray (pargout [2]) ; |
|
|
1180 mxDestroyArray (pargout [3]) ; |
|
|
1181 error ("constructing L failed", A_is_complex, nargout, pargout, |
|
|
1182 Control, Info, status, do_info) ; |
|
|
1183 return ; |
|
|
1184 } |
|
|
1185 |
|
|
1186 /* -------------------------------------------------------------- */ |
|
|
1187 /* print L, U, P, and Q */ |
|
|
1188 /* -------------------------------------------------------------- */ |
|
|
1189 |
|
|
1190 if (A_is_complex) |
|
|
1191 { |
|
|
1192 if (print_level >= 3) mexPrintf ("\nL: ") ; |
|
|
1193 (void) umfpack_zi_report_matrix (n_row, n_inner, Lp, Li, |
|
|
1194 Lx, Lz, 1, Control) ; |
|
|
1195 if (print_level >= 3) mexPrintf ("\nU: ") ; |
|
|
1196 (void) umfpack_zi_report_matrix (n_inner, n_col, Up, Ui, |
|
|
1197 Ux, Uz, 1, Control) ; |
|
|
1198 if (print_level >= 3) mexPrintf ("\nP: ") ; |
|
|
1199 (void) umfpack_zi_report_perm (n_row, P, Control) ; |
|
|
1200 if (print_level >= 3) mexPrintf ("\nQ: ") ; |
|
|
1201 (void) umfpack_zi_report_perm (n_col, Q, Control) ; |
|
|
1202 } |
|
|
1203 else |
|
|
1204 { |
|
|
1205 if (print_level >= 3) mexPrintf ("\nL: ") ; |
|
|
1206 (void) umfpack_di_report_matrix (n_row, n_inner, Lp, Li, |
|
|
1207 Lx, 1, Control) ; |
|
|
1208 if (print_level >= 3) mexPrintf ("\nU: ") ; |
|
|
1209 (void) umfpack_di_report_matrix (n_inner, n_col, Up, Ui, |
|
|
1210 Ux, 1, Control) ; |
|
|
1211 if (print_level >= 3) mexPrintf ("\nP: ") ; |
|
|
1212 (void) umfpack_di_report_perm (n_row, P, Control) ; |
|
|
1213 if (print_level >= 3) mexPrintf ("\nQ: ") ; |
|
|
1214 (void) umfpack_di_report_perm (n_col, Q, Control) ; |
|
|
1215 } |
|
|
1216 |
|
|
1217 mxFree (P) ; |
|
|
1218 mxFree (Q) ; |
|
|
1219 |
|
|
1220 } |
|
|
1221 |
|
|
1222 } |
|
|
1223 else |
|
|
1224 { |
|
|
1225 |
|
|
1226 /* ------------------------------------------------------------------ */ |
|
|
1227 /* return the symbolic factorization */ |
|
|
1228 /* ------------------------------------------------------------------ */ |
|
|
1229 |
|
|
1230 Q = (int *) mxMalloc (n_col * sizeof (int)) ; |
|
|
1231 P = (int *) mxMalloc (n_row * sizeof (int)) ; |
|
|
1232 Front_npivcol = (int *) mxMalloc ((nn+1) * sizeof (int)) ; |
|
|
1233 Front_parent = (int *) mxMalloc ((nn+1) * sizeof (int)) ; |
|
|
1234 Front_1strow = (int *) mxMalloc ((nn+1) * sizeof (int)) ; |
|
|
1235 Front_leftmostdesc = (int *) mxMalloc ((nn+1) * sizeof (int)) ; |
|
|
1236 Chain_start = (int *) mxMalloc ((nn+1) * sizeof (int)) ; |
|
|
1237 Chain_maxrows = (int *) mxMalloc ((nn+1) * sizeof (int)) ; |
|
|
1238 Chain_maxcols = (int *) mxMalloc ((nn+1) * sizeof (int)) ; |
|
|
1239 |
|
|
1240 if (A_is_complex) |
|
|
1241 { |
|
|
1242 status = umfpack_zi_get_symbolic (&ignore1, &ignore2, &ignore3, |
|
|
1243 &nz, &nfronts, &nchains, P, Q, Front_npivcol, |
|
|
1244 Front_parent, Front_1strow, Front_leftmostdesc, |
|
|
1245 Chain_start, Chain_maxrows, Chain_maxcols, Symbolic) ; |
|
|
1246 umfpack_zi_free_symbolic (&Symbolic) ; |
|
|
1247 } |
|
|
1248 else |
|
|
1249 { |
|
|
1250 status = umfpack_di_get_symbolic (&ignore1, &ignore2, &ignore3, |
|
|
1251 &nz, &nfronts, &nchains, P, Q, Front_npivcol, |
|
|
1252 Front_parent, Front_1strow, Front_leftmostdesc, |
|
|
1253 Chain_start, Chain_maxrows, Chain_maxcols, Symbolic) ; |
|
|
1254 umfpack_di_free_symbolic (&Symbolic) ; |
|
|
1255 } |
|
|
1256 |
|
|
1257 if (status < 0) |
|
|
1258 { |
|
|
1259 mxFree (P) ; |
|
|
1260 mxFree (Q) ; |
|
|
1261 mxFree (Front_npivcol) ; |
|
|
1262 mxFree (Front_parent) ; |
|
|
1263 mxFree (Front_1strow) ; |
|
|
1264 mxFree (Front_leftmostdesc) ; |
|
|
1265 mxFree (Chain_start) ; |
|
|
1266 mxFree (Chain_maxrows) ; |
|
|
1267 mxFree (Chain_maxcols) ; |
|
|
1268 error ("extracting symbolic factors failed", A_is_complex, nargout, |
|
|
1269 pargout, Control, Info, status, do_info) ; |
|
|
1270 return ; |
|
|
1271 } |
|
|
1272 |
|
|
1273 /* create sparse permutation matrix for P */ |
|
|
1274 pargout [0] = mxCreateSparse (n_row, n_row, n_row, mxREAL) ; |
|
|
1275 Pp = mxGetJc (pargout [0]) ; |
|
|
1276 Pi = mxGetIr (pargout [0]) ; |
|
|
1277 Px = mxGetPr (pargout [0]) ; |
|
|
1278 for (k = 0 ; k < n_row ; k++) |
|
|
1279 { |
|
|
1280 Pp [k] = k ; |
|
|
1281 Px [k] = 1 ; |
|
|
1282 Pi [P [k]] = k ; |
|
|
1283 } |
|
|
1284 Pp [n_row] = n_row ; |
|
|
1285 |
|
|
1286 /* create sparse permutation matrix for Q */ |
|
|
1287 pargout [1] = mxCreateSparse (n_col, n_col, n_col, mxREAL) ; |
|
|
1288 Qp = mxGetJc (pargout [1]) ; |
|
|
1289 Qi = mxGetIr (pargout [1]) ; |
|
|
1290 Qx = mxGetPr (pargout [1]) ; |
|
|
1291 for (k = 0 ; k < n_col ; k++) |
|
|
1292 { |
|
|
1293 Qp [k] = k ; |
|
|
1294 Qx [k] = 1 ; |
|
|
1295 Qi [k] = Q [k] ; |
|
|
1296 } |
|
|
1297 Qp [n_col] = n_col ; |
|
|
1298 |
|
|
1299 /* create Fr */ |
|
|
1300 pargout [2] = mxCreateDoubleMatrix (nfronts+1, 4, mxREAL) ; |
|
|
1301 |
|
|
1302 p1 = mxGetPr (pargout [2]) ; |
|
|
1303 p2 = p1 + nfronts + 1 ; |
|
|
1304 p3 = p2 + nfronts + 1 ; |
|
|
1305 p4 = p3 + nfronts + 1 ; |
|
|
1306 for (i = 0 ; i <= nfronts ; i++) |
|
|
1307 { |
|
|
1308 /* convert parent, 1strow, and leftmostdesc to 1-based */ |
|
|
1309 p1 [i] = (double) (Front_npivcol [i]) ; |
|
|
1310 p2 [i] = (double) (Front_parent [i] + 1) ; |
|
|
1311 p3 [i] = (double) (Front_1strow [i] + 1) ; |
|
|
1312 p4 [i] = (double) (Front_leftmostdesc [i] + 1) ; |
|
|
1313 } |
|
|
1314 |
|
|
1315 /* create Ch */ |
|
|
1316 pargout [3] = mxCreateDoubleMatrix (nchains+1, 3, mxREAL) ; |
|
|
1317 p1 = mxGetPr (pargout [3]) ; |
|
|
1318 p2 = p1 + nchains + 1 ; |
|
|
1319 p3 = p2 + nchains + 1 ; |
|
|
1320 for (i = 0 ; i < nchains ; i++) |
|
|
1321 { |
|
|
1322 p1 [i] = (double) (Chain_start [i] + 1) ; /* convert to 1-based */ |
|
|
1323 p2 [i] = (double) (Chain_maxrows [i]) ; |
|
|
1324 p3 [i] = (double) (Chain_maxcols [i]) ; |
|
|
1325 } |
|
|
1326 p1 [nchains] = Chain_start [nchains] + 1 ; |
|
|
1327 p2 [nchains] = 0 ; |
|
|
1328 p3 [nchains] = 0 ; |
|
|
1329 |
|
|
1330 mxFree (P) ; |
|
|
1331 mxFree (Q) ; |
|
|
1332 mxFree (Front_npivcol) ; |
|
|
1333 mxFree (Front_parent) ; |
|
|
1334 mxFree (Front_1strow) ; |
|
|
1335 mxFree (Front_leftmostdesc) ; |
|
|
1336 mxFree (Chain_start) ; |
|
|
1337 mxFree (Chain_maxrows) ; |
|
|
1338 mxFree (Chain_maxcols) ; |
|
|
1339 } |
|
|
1340 |
|
|
1341 /* ---------------------------------------------------------------------- */ |
|
|
1342 /* report Info */ |
|
|
1343 /* ---------------------------------------------------------------------- */ |
|
|
1344 |
|
|
1345 if (A_is_complex) |
|
|
1346 { |
|
|
1347 umfpack_zi_report_info (Control, Info) ; |
|
|
1348 } |
|
|
1349 else |
|
|
1350 { |
|
|
1351 umfpack_di_report_info (Control, Info) ; |
|
|
1352 } |
|
|
1353 |
|
|
1354 if (do_info > 0) |
|
|
1355 { |
|
|
1356 /* return Info */ |
|
|
1357 pargout [do_info] = mxCreateDoubleMatrix (1, UMFPACK_INFO, mxREAL) ; |
|
|
1358 OutInfo = mxGetPr (pargout [do_info]) ; |
|
|
1359 for (i = 0 ; i < UMFPACK_INFO ; i++) |
|
|
1360 { |
|
|
1361 OutInfo [i] = Info [i] ; |
|
|
1362 } |
|
|
1363 } |
|
|
1364 } |
