1 /* ========================================================================== */
2 /* === UMFPACK mexFunction ================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* Copyright (c) 2005-2012 by Timothy A. Davis, http://www.suitesparse.com. */
7 /* All Rights Reserved. See ../Doc/License.txt for License. */
8 /* -------------------------------------------------------------------------- */
9
10 /*
11 MATLAB interface for umfpack.
12
13 Factor or solve a sparse linear system, returning either the solution
14 x to Ax=b or A'x'=b', or the factorization LU=P(R\A)Q or LU=PAQ. A must be
15 sparse, with nonzero dimensions, but it may be complex, singular, and/or
16 rectangular. b must be a dense n-by-1 vector (real or complex).
17 L is unit lower triangular, U is upper triangular, and R is diagonal.
18 P and Q are permutation matrices (permutations of an identity matrix).
19
20 The matrix A is scaled, by default. Each row i is divided by r (i), where
21 r (i) is the sum of the absolute values of the entries in that row. The
22 scaled matrix has an infinity norm of 1. The scale factors r (i) are
23 returned in a diagonal sparse matrix. If the factorization is:
24
25 [L, U, P, Q, R] = umfpack (A) ;
26
27 then the factorization is
28
29 L*U = P * (R \ A) * Q
30
31 This is safer than returning a matrix R such that L*U = P*R*A*Q, because
32 it avoids the division by small entries. If r(i) is subnormal, multiplying
33 by 1/r(i) would result in an IEEE Infinity, but dividing by r(i) is safe.
34
35 The factorization
36
37 [L, U, P, Q] = umfpack (A) ;
38
39 returns LU factors such that L*U = P*A*Q, with no scaling.
40
41 See umfpack.m, umfpack_details.m, and umfpack.h for details.
42
43 Note that this mexFunction accesses only the user-callable UMFPACK routines.
44 Thus, is also provides another example of how user C code can access
45 UMFPACK.
46
47 Unlike MATLAB, x=b/A is solved by factorizing A, and then solving via the
48 transposed L and U matrices. The solution is still x = (A.'\b.').', except
49 that A is factorized instead of A.'.
50
51 v5.1: port to 64-bit MATLAB
52 */
53
54 #include "umfpack.h"
55 #include "mex.h"
56 #include "matrix.h"
57 #include <string.h>
58 #include <math.h>
59 #include <float.h>
60
61 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
62 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
63 #define MATCH(s1,s2) (strcmp ((s1), (s2)) == 0)
64 #ifndef TRUE
65 #define TRUE (1)
66 #endif
67 #ifndef FALSE
68 #define FALSE (0)
69 #endif
70 #define Long SuiteSparse_long
71
72 /* ========================================================================== */
73 /* === error ================================================================ */
74 /* ========================================================================== */
75
76 /* Return an error message */
77
error(char * s,Long A_is_complex,int nargout,mxArray * pargout[],double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO],Long status)78 static void error
79 (
80 char *s,
81 Long A_is_complex,
82 int nargout,
83 mxArray *pargout [ ],
84 double Control [UMFPACK_CONTROL],
85 double Info [UMFPACK_INFO],
86 Long status
87 )
88 {
89 Long i ;
90 double *Out_Info ;
91 if (A_is_complex)
92 {
93 umfpack_zl_report_status (Control, status) ;
94 umfpack_zl_report_info (Control, Info) ;
95 }
96 else
97 {
98 umfpack_dl_report_status (Control, status) ;
99 umfpack_dl_report_info (Control, Info) ;
100 }
101
102 mexErrMsgTxt (s) ;
103 }
104
105
106 /* ========================================================================== */
107 /* === get_option =========================================================== */
108 /* ========================================================================== */
109
110 /* get a single string or numeric option from the MATLAB options struct */
111
get_option(const mxArray * mxopts,const char * field,double * x,Long * x_present,char ** s)112 static int get_option
113 (
114 /* inputs: */
115 const mxArray *mxopts, /* the MATLAB struct */
116 const char *field, /* the field to get from the MATLAB struct */
117
118 /* outputs: */
119 double *x, /* double value of the field, if present */
120 Long *x_present, /* true if double x is present */
121 char **s /* char value of the field, if present; */
122 /* must be mxFree'd by caller when done */
123 )
124 {
125 Long f ;
126 mxArray *p ;
127
128 /* find the field number */
129 if (mxopts == NULL || mxIsEmpty (mxopts) || !mxIsStruct (mxopts))
130 {
131 /* mxopts is not present, or [ ], or not a struct */
132 f = -1 ;
133 }
134 else
135 {
136 /* f will be -1 if the field is not present */
137 f = mxGetFieldNumber (mxopts, field) ;
138 }
139
140 /* get the field, or NULL if not present */
141 if (f == -1)
142 {
143 p = NULL ;
144 }
145 else
146 {
147 p = mxGetFieldByNumber (mxopts, 0, f) ;
148 }
149
150 *x_present = FALSE ;
151 if (s != NULL)
152 {
153 *s = NULL ;
154 }
155
156 if (p == NULL)
157 {
158 /* option not present */
159 return (TRUE) ;
160 }
161 if (mxIsNumeric (p))
162 {
163 /* option is numeric */
164 if (x == NULL)
165 {
166 mexPrintf ("opts.%s field must be a string\n", field) ;
167 mexErrMsgIdAndTxt ("UMFPACK:invalidInput", "invalid option") ;
168 }
169 *x = mxGetScalar (p) ;
170 *x_present = TRUE ;
171 }
172 else if (mxIsChar (p))
173 {
174 /* option is a MATLAB string; convert it to a C-style string */
175 if (s == NULL)
176 {
177 mexPrintf ("opts.%s field must be a numeric value\n", field) ;
178 mexErrMsgIdAndTxt ("UMFPACK:invalidInput", "invalid option") ;
179 }
180 *s = mxArrayToString (p) ;
181 }
182 return (TRUE) ;
183 }
184
185
186 /* ========================================================================== */
187 /* === get_all_options ====================================================== */
188 /* ========================================================================== */
189
190 /* get all the options from the MATLAB struct.
191
192 opts.prl >= 0, default 1 (errors only)
193 opts.strategy 'auto', 'unsymmetric', 'symmetric', default auto
194 opts.ordering 'amd' AMD for A+A', COLAMD for A'A
195 'default' use CHOLMOD (AMD then METIS; take best fount)
196 'metis' use METIS
197 'none' no fill-reducing ordering
198 'given' use Qinit (this is default if Qinit present)
199 'best' try AMD/COLAMD, METIS, and NESDIS; take best
200 opts.tol default 0.1
201 opts.symtol default 0.001
202 opts.scale row scaling: 'none', 'sum', 'max'
203 opts.irstep max # of steps of iterative refinement, default 2
204 opts.singletons 'yes','no' default 'yes'
205 */
206
get_all_options(const mxArray * mxopts,double * Control)207 int get_all_options
208 (
209 const mxArray *mxopts,
210 double *Control
211 )
212 {
213 double x ;
214 char *s ;
215 Long x_present, i, info_details ;
216
217 /* ---------------------------------------------------------------------- */
218 /* prl: an integer, default 1 */
219 /* ---------------------------------------------------------------------- */
220
221 get_option (mxopts, "prl", &x, &x_present, NULL) ;
222 Control [UMFPACK_PRL] = x_present ? ((Long) x) : UMFPACK_DEFAULT_PRL ;
223 if (mxIsNaN (Control [UMFPACK_PRL]))
224 {
225 Control [UMFPACK_PRL] = UMFPACK_DEFAULT_PRL ;
226 }
227
228 /* ---------------------------------------------------------------------- */
229 /* strategy: a string */
230 /* ---------------------------------------------------------------------- */
231
232 get_option (mxopts, "strategy", NULL, &x_present, &s) ;
233 if (s != NULL)
234 {
235 if (MATCH (s, "auto"))
236 {
237 Control [UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO ;
238 }
239 else if (MATCH (s, "unsymmetric"))
240 {
241 Control [UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC ;
242 }
243 else if (MATCH (s, "symmetric"))
244 {
245 Control [UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC ;
246 }
247 else if (MATCH (s, "default"))
248 {
249 Control [UMFPACK_STRATEGY] = UMFPACK_DEFAULT_STRATEGY ;
250 }
251 else
252 {
253 mexErrMsgIdAndTxt ("UMFPACK:invalidInput","invalid strategy") ;
254 }
255 mxFree (s) ;
256 }
257
258 /* ---------------------------------------------------------------------- */
259 /* ordering: a string */
260 /* ---------------------------------------------------------------------- */
261
262 get_option (mxopts, "ordering", NULL, &x_present, &s) ;
263 if (s != NULL)
264 {
265 if (MATCH (s, "amd") || MATCH (s, "colamd") || MATCH (s,"bestamd"))
266 {
267 Control [UMFPACK_ORDERING] = UMFPACK_ORDERING_AMD ;
268 }
269 #ifndef NCHOLMOD
270 else if (MATCH (s, "fixed") || MATCH (s, "none") || MATCH (s,"natural"))
271 {
272 Control [UMFPACK_ORDERING] = UMFPACK_ORDERING_NONE ;
273 }
274 else if (MATCH (s, "metis"))
275 {
276 Control [UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS ;
277 }
278 else if (MATCH (s, "cholmod"))
279 {
280 Control [UMFPACK_ORDERING] = UMFPACK_ORDERING_CHOLMOD ;
281 }
282 else if (MATCH (s, "best"))
283 {
284 Control [UMFPACK_ORDERING] = UMFPACK_ORDERING_BEST ;
285 }
286 #endif
287 else if (MATCH (s, "given"))
288 {
289 Control [UMFPACK_ORDERING] = UMFPACK_ORDERING_GIVEN ;
290 }
291 else if (MATCH (s, "default"))
292 {
293 Control [UMFPACK_ORDERING] = UMFPACK_DEFAULT_ORDERING ;
294 }
295 else
296 {
297 mexErrMsgIdAndTxt ("UMFPACK:invalidInput","invalid ordering") ;
298 }
299 mxFree (s) ;
300 }
301
302 /* ---------------------------------------------------------------------- */
303 /* scale: a string */
304 /* ---------------------------------------------------------------------- */
305
306 get_option (mxopts, "scale", NULL, &x_present, &s) ;
307 if (s != NULL)
308 {
309 if (MATCH (s, "none"))
310 {
311 Control [UMFPACK_SCALE] = UMFPACK_SCALE_NONE ;
312 }
313 else if (MATCH (s, "sum"))
314 {
315 Control [UMFPACK_SCALE] = UMFPACK_SCALE_SUM ;
316 }
317 else if (MATCH (s, "max"))
318 {
319 Control [UMFPACK_SCALE] = UMFPACK_SCALE_MAX ;
320 }
321 else if (MATCH (s, "default"))
322 {
323 Control [UMFPACK_SCALE] = UMFPACK_DEFAULT_SCALE ;
324 }
325 else
326 {
327 mexErrMsgIdAndTxt ("UMFPACK:invalidInput","invalid scale") ;
328 }
329 mxFree (s) ;
330 }
331
332 /* ---------------------------------------------------------------------- */
333 /* tol: a double */
334 /* ---------------------------------------------------------------------- */
335
336 get_option (mxopts, "tol", &x, &x_present, NULL) ;
337 Control [UMFPACK_PIVOT_TOLERANCE]
338 = x_present ? x : UMFPACK_DEFAULT_PIVOT_TOLERANCE ;
339
340 /* ---------------------------------------------------------------------- */
341 /* symtol: a double */
342 /* ---------------------------------------------------------------------- */
343
344 get_option (mxopts, "symtol", &x, &x_present, NULL) ;
345 Control [UMFPACK_SYM_PIVOT_TOLERANCE]
346 = x_present ? x : UMFPACK_DEFAULT_SYM_PIVOT_TOLERANCE ;
347
348 /* ---------------------------------------------------------------------- */
349 /* irstep: an integer */
350 /* ---------------------------------------------------------------------- */
351
352 get_option (mxopts, "irstep", &x, &x_present, NULL) ;
353 Control [UMFPACK_IRSTEP] = x_present ? x : UMFPACK_DEFAULT_IRSTEP ;
354
355 /* ---------------------------------------------------------------------- */
356 /* singletons: a string */
357 /* ---------------------------------------------------------------------- */
358
359 get_option (mxopts, "singletons", NULL, &x_present, &s) ;
360 if (s != NULL)
361 {
362 if (MATCH (s, "enable"))
363 {
364 Control [UMFPACK_SINGLETONS] = TRUE ;
365 }
366 else if (MATCH (s, "disable"))
367 {
368 Control [UMFPACK_SINGLETONS] = FALSE ;
369 }
370 else if (MATCH (s, "default"))
371 {
372 Control [UMFPACK_SINGLETONS] = UMFPACK_DEFAULT_SINGLETONS ;
373 }
374 mxFree (s) ;
375 }
376
377 /* ---------------------------------------------------------------------- */
378 /* details: an int */
379 /* ---------------------------------------------------------------------- */
380
381 get_option (mxopts, "details", &x, &x_present, NULL) ;
382 info_details = x_present ? x : 0 ;
383 return (info_details) ;
384 }
385
386
387 /* ========================================================================== */
388 /* === umfpack_mx_info_details ============================================== */
389 /* ========================================================================== */
390
391 /* Return detailed info struct; useful for UMFPACK development only */
392
393 #define XFIELD(x) mxSetFieldByNumber (info, 0, k++, mxCreateDoubleScalar (x))
394 #define SFIELD(s) mxSetFieldByNumber (info, 0, k++, mxCreateString (s))
395 #define YESNO(x) ((x) ? "yes" : "no")
396
umfpack_mx_info_details(double * Control,double * Info)397 mxArray *umfpack_mx_info_details /* return a struct with info statistics */
398 (
399 double *Control,
400 double *Info
401 )
402 {
403 Long k = 0 ;
404 mxArray *info ;
405 Long sizeof_unit = (Long) Info [UMFPACK_SIZE_OF_UNIT] ;
406
407 const char *info_struct [ ] =
408 {
409 "control_prl",
410 "control_dense_row",
411 "control_dense_col",
412 "control_tol",
413 "control_block_size",
414 "control_strategy",
415 "control_alloc_init",
416 "control_irstep",
417 "umfpack_compiled_with_BLAS",
418 "control_ordering",
419 "control_singletons",
420 "control_fixQ",
421 "control_amd_dense",
422 "control_symtol",
423 "control_scale",
424 "control_front_alloc",
425 "control_droptol",
426 "control_aggressive",
427
428 "status",
429
430 "nrow",
431 "ncol",
432 "anz",
433
434 "sizeof_unit",
435 "sizeof_int",
436 "sizeof_long",
437 "sizeof_pointer",
438 "sizeof_entry",
439 "number_of_dense_rows",
440 "number_of_empty_rows",
441 "number_of_dense_cols",
442 "number_of_empty_cols",
443
444 "number_of_memory_defragmentations_during_symbolic_analysis",
445 "memory_usage_for_symbolic_analysis_in_bytes",
446 "size_of_symbolic_factorization_in_bytes",
447 "symbolic_time",
448 "symbolic_walltime",
449 "strategy_used",
450 "ordering_used",
451 "Qfixed",
452 "diagonol_pivots_preferred",
453
454 "number_of_column_singletons",
455 "number_of_row_singletons",
456
457 /* only computed if symmetric strategy used */
458 "symmetric_strategy_S_size",
459 "symmetric_strategy_S_symmetric",
460 "symmetric_strategy_pattern_symmetry",
461 "symmetric_strategy_nnz_A_plus_AT",
462 "symmetric_strategy_nnz_diag",
463 "symmetric_strategy_lunz",
464 "symmetric_strategy_flops",
465 "symmetric_strategy_ndense",
466 "symmetric_strategy_dmax",
467
468 "estimated_size_of_numeric_factorization_in_bytes",
469 "estimated_peak_memory_in_bytes",
470 "estimated_number_of_floating_point_operations",
471 "estimated_number_of_entries_in_L",
472 "estimated_number_of_entries_in_U",
473 "estimated_variable_init_in_bytes",
474 "estimated_variable_peak_in_bytes",
475 "estimated_variable_final_in_bytes",
476 "estimated_number_of_entries_in_largest_frontal_matrix",
477 "estimated_largest_frontal_matrix_row_dimension",
478 "estimated_largest_frontal_matrix_col_dimension",
479
480 "size_of_numeric_factorization_in_bytes",
481 "total_memory_usage_in_bytes",
482 "number_of_floating_point_operations",
483 "number_of_entries_in_L",
484 "number_of_entries_in_U",
485 "variable_init_in_bytes",
486 "variable_peak_in_bytes",
487 "variable_final_in_bytes",
488 "number_of_entries_in_largest_frontal_matrix",
489 "largest_frontal_matrix_row_dimension",
490 "largest_frontal_matrix_col_dimension",
491
492 "number_of_memory_defragmentations_during_numeric_factorization",
493 "number_of_memory_reallocations_during_numeric_factorization",
494 "number_of_costly_memory_reallocations_during_numeric_factorization",
495 "number_of_integers_in_compressed_pattern_of_L_and_U",
496 "number_of_entries_in_LU_data_structure",
497 "numeric_time",
498 "nnz_diagonal_of_U",
499 "rcond_estimate",
500 "scaling_used",
501 "min_abs_row_sum_of_A",
502 "max_abs_row_sum_of_A",
503 "min_abs_diagonal_of_U",
504 "max_abs_diagonal_of_U",
505 "alloc_init_used",
506 "number_of_forced_updates",
507 "numeric_walltime",
508 "symmetric_strategy_number_off_diagonal_pivots",
509
510 "number_of_entries_in_L_including_dropped_entries",
511 "number_of_entries_in_U_including_dropped_entries",
512 "number_of_small_entries_dropped_from_L_and_U",
513
514 "number_of_iterative_refinement_steps_taken",
515 "number_of_iterative_refinement_steps_attempted",
516 "omega1",
517 "omega2",
518 "solve_flops",
519 "solve_time",
520 "solve_walltime"
521 } ;
522
523 info = mxCreateStructMatrix (1, 1, 100, info_struct) ;
524
525 XFIELD (Control [UMFPACK_PRL]) ;
526 XFIELD (Control [UMFPACK_DENSE_ROW]) ;
527 XFIELD (Control [UMFPACK_DENSE_COL]) ;
528 XFIELD (Control [UMFPACK_PIVOT_TOLERANCE]) ;
529 XFIELD (Control [UMFPACK_BLOCK_SIZE]) ;
530
531 switch ((int) Control [UMFPACK_STRATEGY])
532 {
533 case UMFPACK_STRATEGY_UNSYMMETRIC: SFIELD ("unsymmetric") ; break ;
534 case UMFPACK_STRATEGY_SYMMETRIC: SFIELD ("symmetric") ; break ;
535 default:
536 case UMFPACK_DEFAULT_STRATEGY: SFIELD ("auto") ; break ;
537 }
538
539 XFIELD (Control [UMFPACK_ALLOC_INIT]) ;
540 XFIELD (Control [UMFPACK_IRSTEP]) ;
541 SFIELD (YESNO (Control [UMFPACK_COMPILED_WITH_BLAS])) ;
542
543 switch ((int) Control [UMFPACK_ORDERING])
544 {
545 case UMFPACK_ORDERING_NONE: SFIELD ("none") ; break ;
546 case UMFPACK_ORDERING_AMD: SFIELD ("amd") ; break ;
547 case UMFPACK_ORDERING_METIS: SFIELD ("metis") ; break ;
548 default:
549 case UMFPACK_ORDERING_CHOLMOD: SFIELD ("cholmod") ; break ;
550 case UMFPACK_ORDERING_BEST: SFIELD ("best") ; break ;
551 case UMFPACK_ORDERING_GIVEN: SFIELD ("given") ; break ;
552 }
553
554 SFIELD (YESNO (Control [UMFPACK_SINGLETONS])) ;
555
556 if (Control [UMFPACK_FIXQ] > 0)
557 {
558 SFIELD ("forced true") ;
559 }
560 else if (Control [UMFPACK_FIXQ] < 0)
561 {
562 SFIELD ("forced false") ;
563 }
564 else
565 {
566 SFIELD ("auto") ;
567 }
568
569 XFIELD (Control [UMFPACK_AMD_DENSE]) ;
570 XFIELD (Control [UMFPACK_SYM_PIVOT_TOLERANCE]) ;
571
572 switch ((int) Control [UMFPACK_SCALE])
573 {
574 case UMFPACK_SCALE_NONE: SFIELD ("none") ; break ;
575 case UMFPACK_SCALE_MAX: SFIELD ("max") ; break ;
576 default:
577 case UMFPACK_SCALE_SUM: SFIELD ("sum") ; break ;
578 }
579
580 XFIELD (Control [UMFPACK_FRONT_ALLOC_INIT]) ;
581 XFIELD (Control [UMFPACK_DROPTOL]) ;
582 SFIELD (YESNO (Control [UMFPACK_AGGRESSIVE])) ;
583
584 switch ((int) Info [UMFPACK_STATUS])
585 {
586 case UMFPACK_OK:
587 SFIELD ("ok") ; break ;
588 case UMFPACK_WARNING_singular_matrix:
589 SFIELD ("singular matrix") ; break ;
590 case UMFPACK_WARNING_determinant_underflow:
591 SFIELD ("determinant underflow") ; break ;
592 case UMFPACK_WARNING_determinant_overflow:
593 SFIELD ("determinant overflow") ; break ;
594 case UMFPACK_ERROR_out_of_memory:
595 SFIELD ("out of memory") ; break ;
596 case UMFPACK_ERROR_invalid_Numeric_object:
597 SFIELD ("invalid numeric LU object") ; break ;
598 case UMFPACK_ERROR_invalid_Symbolic_object:
599 SFIELD ("invalid symbolic LU object") ; break ;
600 case UMFPACK_ERROR_argument_missing:
601 SFIELD ("argument missing") ; break ;
602 case UMFPACK_ERROR_n_nonpositive:
603 SFIELD ("n < 0") ; break ;
604 case UMFPACK_ERROR_invalid_matrix:
605 SFIELD ("invalid matrix") ; break ;
606 case UMFPACK_ERROR_different_pattern:
607 SFIELD ("pattern changed") ; break ;
608 case UMFPACK_ERROR_invalid_system:
609 SFIELD ("invalid system") ; break ;
610 case UMFPACK_ERROR_invalid_permutation:
611 SFIELD ("invalid permutation") ; break ;
612 case UMFPACK_ERROR_internal_error:
613 SFIELD ("internal error; contact DrTimothyAldenDavis@gmail.com") ;
614 break ;
615 case UMFPACK_ERROR_file_IO:
616 SFIELD ("file I/O error") ; break ;
617 case UMFPACK_ERROR_ordering_failed:
618 SFIELD ("ordering failed") ; break ;
619 default:
620 if (Info [UMFPACK_STATUS] < 0)
621 {
622 SFIELD ("unknown error") ;
623 }
624 else
625 {
626 SFIELD ("unknown warning") ;
627 }
628 break ;
629 }
630
631 XFIELD (Info [UMFPACK_NROW]) ;
632 XFIELD (Info [UMFPACK_NCOL]) ;
633 XFIELD (Info [UMFPACK_NZ]) ;
634
635 XFIELD (Info [UMFPACK_SIZE_OF_UNIT]) ;
636 XFIELD (Info [UMFPACK_SIZE_OF_INT]) ;
637 XFIELD (Info [UMFPACK_SIZE_OF_LONG]) ;
638 XFIELD (Info [UMFPACK_SIZE_OF_POINTER]) ;
639 XFIELD (Info [UMFPACK_SIZE_OF_ENTRY]) ;
640 XFIELD (Info [UMFPACK_NDENSE_ROW]) ;
641 XFIELD (Info [UMFPACK_NEMPTY_ROW]) ;
642 XFIELD (Info [UMFPACK_NDENSE_COL]) ;
643 XFIELD (Info [UMFPACK_NEMPTY_COL]) ;
644
645 XFIELD (Info [UMFPACK_SYMBOLIC_DEFRAG]) ;
646 XFIELD (Info [UMFPACK_SYMBOLIC_PEAK_MEMORY] * sizeof_unit) ;
647 XFIELD (Info [UMFPACK_SYMBOLIC_SIZE] * sizeof_unit) ;
648 XFIELD (Info [UMFPACK_SYMBOLIC_TIME]) ;
649 XFIELD (Info [UMFPACK_SYMBOLIC_WALLTIME]) ;
650
651 switch ((int) Info [UMFPACK_STRATEGY_USED])
652 {
653 default:
654 case UMFPACK_STRATEGY_UNSYMMETRIC: SFIELD ("unsymmetric") ; break ;
655 case UMFPACK_STRATEGY_SYMMETRIC: SFIELD ("symmetric") ; break ;
656 }
657
658 switch ((int) Info [UMFPACK_ORDERING_USED])
659 {
660 case UMFPACK_ORDERING_AMD: SFIELD ("amd") ; break ;
661 case UMFPACK_ORDERING_METIS: SFIELD ("metis") ; break ;
662 case UMFPACK_ORDERING_GIVEN: SFIELD ("given") ; break ;
663 default: SFIELD ("none") ; break ;
664 }
665
666 SFIELD (YESNO (Info [UMFPACK_QFIXED])) ;
667 SFIELD (YESNO (Info [UMFPACK_DIAG_PREFERRED])) ;
668
669 XFIELD (Info [UMFPACK_COL_SINGLETONS]) ;
670 XFIELD (Info [UMFPACK_ROW_SINGLETONS]) ;
671
672 /* only computed if symmetric ordering is used */
673 XFIELD (Info [UMFPACK_N2]) ;
674 SFIELD (YESNO (Info [UMFPACK_S_SYMMETRIC])) ;
675 XFIELD (Info [UMFPACK_PATTERN_SYMMETRY]) ;
676 XFIELD (Info [UMFPACK_NZ_A_PLUS_AT]) ;
677 XFIELD (Info [UMFPACK_NZDIAG]) ;
678 XFIELD (Info [UMFPACK_SYMMETRIC_LUNZ]) ;
679 XFIELD (Info [UMFPACK_SYMMETRIC_FLOPS]) ;
680 XFIELD (Info [UMFPACK_SYMMETRIC_NDENSE]) ;
681 XFIELD (Info [UMFPACK_SYMMETRIC_DMAX]) ;
682
683 XFIELD (Info [UMFPACK_NUMERIC_SIZE_ESTIMATE] * sizeof_unit) ;
684 XFIELD (Info [UMFPACK_PEAK_MEMORY_ESTIMATE] * sizeof_unit) ;
685 XFIELD (Info [UMFPACK_FLOPS_ESTIMATE]) ;
686 XFIELD (Info [UMFPACK_LNZ_ESTIMATE]) ;
687 XFIELD (Info [UMFPACK_UNZ_ESTIMATE]) ;
688 XFIELD (Info [UMFPACK_VARIABLE_INIT_ESTIMATE] * sizeof_unit) ;
689 XFIELD (Info [UMFPACK_VARIABLE_PEAK_ESTIMATE] * sizeof_unit) ;
690 XFIELD (Info [UMFPACK_VARIABLE_FINAL_ESTIMATE] * sizeof_unit) ;
691 XFIELD (Info [UMFPACK_MAX_FRONT_SIZE_ESTIMATE]) ;
692 XFIELD (Info [UMFPACK_MAX_FRONT_NROWS_ESTIMATE]) ;
693 XFIELD (Info [UMFPACK_MAX_FRONT_NCOLS_ESTIMATE]) ;
694
695 XFIELD (Info [UMFPACK_NUMERIC_SIZE] * sizeof_unit) ;
696 XFIELD (Info [UMFPACK_PEAK_MEMORY] * sizeof_unit) ;
697 XFIELD (Info [UMFPACK_FLOPS]) ;
698 XFIELD (Info [UMFPACK_LNZ]) ;
699 XFIELD (Info [UMFPACK_UNZ]) ;
700 XFIELD (Info [UMFPACK_VARIABLE_INIT] * sizeof_unit) ;
701 XFIELD (Info [UMFPACK_VARIABLE_PEAK] * sizeof_unit) ;
702 XFIELD (Info [UMFPACK_VARIABLE_FINAL] * sizeof_unit) ;
703 XFIELD (Info [UMFPACK_MAX_FRONT_SIZE]) ;
704 XFIELD (Info [UMFPACK_MAX_FRONT_NROWS]) ;
705 XFIELD (Info [UMFPACK_MAX_FRONT_NCOLS]) ;
706
707 XFIELD (Info [UMFPACK_NUMERIC_DEFRAG]) ;
708 XFIELD (Info [UMFPACK_NUMERIC_REALLOC]) ;
709 XFIELD (Info [UMFPACK_NUMERIC_COSTLY_REALLOC]) ;
710 XFIELD (Info [UMFPACK_COMPRESSED_PATTERN]) ;
711 XFIELD (Info [UMFPACK_LU_ENTRIES]) ;
712 XFIELD (Info [UMFPACK_NUMERIC_TIME]) ;
713 XFIELD (Info [UMFPACK_UDIAG_NZ]) ;
714 XFIELD (Info [UMFPACK_RCOND]) ;
715
716 switch ((int) Info [UMFPACK_WAS_SCALED])
717 {
718 case UMFPACK_SCALE_NONE: SFIELD ("none") ; break ;
719 case UMFPACK_SCALE_MAX: SFIELD ("max") ; break ;
720 default:
721 case UMFPACK_SCALE_SUM: SFIELD ("sum") ; break ;
722 }
723
724 XFIELD (Info [UMFPACK_RSMIN]) ;
725 XFIELD (Info [UMFPACK_RSMAX]) ;
726 XFIELD (Info [UMFPACK_UMIN]) ;
727 XFIELD (Info [UMFPACK_UMAX]) ;
728 XFIELD (Info [UMFPACK_ALLOC_INIT_USED]) ;
729 XFIELD (Info [UMFPACK_FORCED_UPDATES]) ;
730 XFIELD (Info [UMFPACK_NUMERIC_WALLTIME]) ;
731 XFIELD (Info [UMFPACK_NOFF_DIAG]) ;
732
733 XFIELD (Info [UMFPACK_ALL_LNZ]) ;
734 XFIELD (Info [UMFPACK_ALL_UNZ]) ;
735 XFIELD (Info [UMFPACK_NZDROPPED]) ;
736
737 XFIELD (Info [UMFPACK_IR_TAKEN]) ;
738 XFIELD (Info [UMFPACK_IR_ATTEMPTED]) ;
739 XFIELD (Info [UMFPACK_OMEGA1]) ;
740 XFIELD (Info [UMFPACK_OMEGA2]) ;
741 XFIELD (Info [UMFPACK_SOLVE_FLOPS]) ;
742 XFIELD (Info [UMFPACK_SOLVE_TIME]) ;
743 XFIELD (Info [UMFPACK_SOLVE_WALLTIME]) ;
744
745 return (info) ;
746 }
747
748
749 /* ========================================================================== */
750 /* === umfpack_mx_info_user ================================================= */
751 /* ========================================================================== */
752
753 /* Return user-friendly info struct */
754
umfpack_mx_info_user(double * Control,double * Info,Long do_solve)755 mxArray *umfpack_mx_info_user /* return a struct with info statistics */
756 (
757 double *Control,
758 double *Info,
759 Long do_solve
760 )
761 {
762 Long k = 0 ;
763 mxArray *info ;
764 Long sizeof_unit = (Long) Info [UMFPACK_SIZE_OF_UNIT] ;
765
766 const char *info_struct [ ] =
767 {
768 "analysis_time",
769 "strategy_used",
770 "ordering_used",
771 "memory_usage_in_bytes",
772 "factorization_flop_count",
773 "nnz_in_L_plus_U",
774 "rcond_estimate",
775 "factorization_time",
776 /* if solve */
777 "iterative_refinement_steps",
778 "solve_flop_count",
779 "solve_time"
780 } ;
781
782 info = mxCreateStructMatrix (1, 1, do_solve ? 11 : 8, info_struct) ;
783
784 XFIELD (Info [UMFPACK_SYMBOLIC_WALLTIME]) ;
785
786 switch ((int) Info [UMFPACK_STRATEGY_USED])
787 {
788 default:
789 case UMFPACK_STRATEGY_UNSYMMETRIC: SFIELD ("unsymmetric") ; break ;
790 case UMFPACK_STRATEGY_SYMMETRIC: SFIELD ("symmetric") ; break ;
791 }
792
793 switch ((int) Info [UMFPACK_ORDERING_USED])
794 {
795 case UMFPACK_ORDERING_AMD: SFIELD ("amd") ; break ;
796 case UMFPACK_ORDERING_METIS: SFIELD ("metis") ; break ;
797 case UMFPACK_ORDERING_GIVEN: SFIELD ("given") ; break ;
798 default: SFIELD ("none") ; break ;
799 }
800
801 XFIELD (Info [UMFPACK_PEAK_MEMORY] * sizeof_unit) ;
802 XFIELD (Info [UMFPACK_FLOPS]) ;
803 XFIELD (Info [UMFPACK_LNZ] + Info [UMFPACK_UNZ] - Info [UMFPACK_UDIAG_NZ]) ;
804 XFIELD (Info [UMFPACK_RCOND]) ;
805 XFIELD (Info [UMFPACK_NUMERIC_WALLTIME]) ;
806
807 if (do_solve)
808 {
809 XFIELD (Info [UMFPACK_IR_TAKEN]) ;
810 XFIELD (Info [UMFPACK_SOLVE_FLOPS]) ;
811 XFIELD (Info [UMFPACK_SOLVE_WALLTIME]) ;
812 }
813
814 return (info) ;
815 }
816
817
818 /* ========================================================================== */
819 /* === umfpack_mx_defaults ================================================== */
820 /* ========================================================================== */
821
822 /* Return a struct with default Control settings (except opts.details). */
823
umfpack_mx_defaults(void)824 mxArray *umfpack_mx_defaults ( void )
825 {
826 mxArray *opts ;
827 const char *opts_struct [ ] =
828 {
829 "prl",
830 "strategy",
831 "ordering",
832 "tol",
833 "symtol",
834 "scale",
835 "irstep",
836 "singletons"
837 } ;
838 opts = mxCreateStructMatrix (1, 1, 8, opts_struct) ;
839 mxSetFieldByNumber (opts, 0, 0,
840 mxCreateDoubleScalar (UMFPACK_DEFAULT_PRL)) ;
841 mxSetFieldByNumber (opts, 0, 1, mxCreateString ("auto")) ;
842 mxSetFieldByNumber (opts, 0, 2, mxCreateString ("default")) ;
843 mxSetFieldByNumber (opts, 0, 3,
844 mxCreateDoubleScalar (UMFPACK_DEFAULT_PIVOT_TOLERANCE)) ;
845 mxSetFieldByNumber (opts, 0, 4,
846 mxCreateDoubleScalar (UMFPACK_DEFAULT_SYM_PIVOT_TOLERANCE)) ;
847 mxSetFieldByNumber (opts, 0, 5, mxCreateString ("sum")) ;
848 mxSetFieldByNumber (opts, 0, 6,
849 mxCreateDoubleScalar (UMFPACK_DEFAULT_IRSTEP)) ;
850 mxSetFieldByNumber (opts, 0, 7, mxCreateString ("enable")) ;
851 return (opts) ;
852 }
853
854
855 /* ========================================================================== */
856 /* === UMFPACK ============================================================== */
857 /* ========================================================================== */
858
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])859 void mexFunction
860 (
861 int nargout, /* number of outputs */
862 mxArray *pargout [ ], /* output arguments */
863 int nargin, /* number of inputs */
864 const mxArray *pargin [ ] /* input arguments */
865 )
866 {
867
868 /* ---------------------------------------------------------------------- */
869 /* local variables */
870 /* ---------------------------------------------------------------------- */
871
872 double Info [UMFPACK_INFO], Control [UMFPACK_CONTROL], dx, dz, dexp ;
873 double *Lx, *Lz, *Ux, *Uz, *Ax, *Az, *Bx, *Bz, *Xx, *Xz, *User_Control,
874 *p, *q, *Out_Info, *p1, *p2, *p3, *p4, *Ltx, *Ltz, *Rs, *Px, *Qx ;
875 void *Symbolic, *Numeric ;
876 Long *Lp, *Li, *Up, *Ui, *Ap, *Ai, *P, *Q, do_solve, lnz, unz, nn, i,
877 transpose, size, do_info, do_numeric, *Front_npivcol, op, k, *Rp, *Ri,
878 *Front_parent, *Chain_start, *Chain_maxrows, *Chain_maxcols, nz, status,
879 nfronts, nchains, *Ltp, *Ltj, *Qinit, print_level, status2, no_scale,
880 *Front_1strow, *Front_leftmostdesc, n_row, n_col, n_inner, sys,
881 ignore1, ignore2, ignore3, A_is_complex, B_is_complex, X_is_complex,
882 *Pp, *Pi, *Qp, *Qi, do_recip, do_det ;
883 mxArray *Amatrix, *Bmatrix, *User_Control_struct, *User_Qinit ;
884 char *operator, *operation ;
885 mxComplexity Atype, Xtype ;
886 char warning [200] ;
887 int info_details ;
888
889 /* ---------------------------------------------------------------------- */
890 /* get inputs A, b, and the operation to perform */
891 /* ---------------------------------------------------------------------- */
892
893 if (nargin > 1 && mxIsStruct (pargin [nargin-1]))
894 {
895 User_Control_struct = (mxArray *) (pargin [nargin-1]) ;
896 }
897 else
898 {
899 User_Control_struct = (mxArray *) NULL ;
900 }
901
902 User_Qinit = (mxArray *) NULL ;
903
904 do_info = 0 ;
905 do_solve = FALSE ;
906 do_numeric = TRUE ;
907 transpose = FALSE ;
908 no_scale = FALSE ;
909 do_det = FALSE ;
910
911 /* find the operator */
912 op = 0 ;
913 for (i = 0 ; i < nargin ; i++)
914 {
915 if (mxIsChar (pargin [i]))
916 {
917 op = i ;
918 break ;
919 }
920 }
921
922 if (op > 0)
923 {
924 operator = mxArrayToString (pargin [op]) ;
925
926 if (MATCH (operator, "\\"))
927 {
928
929 /* -------------------------------------------------------------- */
930 /* matrix left divide, x = A\b */
931 /* -------------------------------------------------------------- */
932
933 /*
934 [x, Info] = umfpack (A, '\', b) ;
935 [x, Info] = umfpack (A, '\', b, Control) ;
936 [x, Info] = umfpack (A, Qinit, '\', b) ;
937 [x, Info] = umfpack (A, Qinit, '\', b, Control) ;
938 */
939
940 operation = "x = A\\b" ;
941 do_solve = TRUE ;
942 Amatrix = (mxArray *) pargin [0] ;
943 Bmatrix = (mxArray *) pargin [op+1] ;
944
945 if (nargout == 2)
946 {
947 do_info = 1 ;
948 }
949 if (op == 2)
950 {
951 User_Qinit = (mxArray *) pargin [1] ;
952 }
953 if (nargin < 3 || nargin > 5 || nargout > 2)
954 {
955 mexErrMsgTxt ("wrong number of arguments") ;
956 }
957
958 }
959 else if (MATCH (operator, "/"))
960 {
961
962 /* -------------------------------------------------------------- */
963 /* matrix right divide, x = b/A */
964 /* -------------------------------------------------------------- */
965
966 /*
967 [x, Info] = umfpack (b, '/', A) ;
968 [x, Info] = umfpack (b, '/', A, Control) ;
969 [x, Info] = umfpack (b, '/', A, Qinit) ;
970 [x, Info] = umfpack (b, '/', A, Qinit, Control) ;
971 */
972
973 operation = "x = b/A" ;
974 do_solve = TRUE ;
975 transpose = TRUE ;
976 Amatrix = (mxArray *) pargin [2] ;
977 Bmatrix = (mxArray *) pargin [0] ;
978
979 if (nargout == 2)
980 {
981 do_info = 1 ;
982 }
983 if (nargin >= 4 && mxIsDouble (pargin [3]))
984 {
985 User_Qinit = (mxArray *) pargin [3] ;
986 }
987 if (nargin < 3 || nargin > 5 || nargout > 2)
988 {
989 mexErrMsgTxt ("wrong number of arguments") ;
990 }
991
992 }
993 else if (MATCH (operator, "symbolic"))
994 {
995
996 /* -------------------------------------------------------------- */
997 /* symbolic factorization only */
998 /* -------------------------------------------------------------- */
999
1000 /*
1001 [P Q Fr Ch Info] = umfpack (A, 'symbolic') ;
1002 [P Q Fr Ch Info] = umfpack (A, 'symbolic', Control) ;
1003 [P Q Fr Ch Info] = umfpack (A, Qinit, 'symbolic') ;
1004 [P Q Fr Ch Info] = umfpack (A, Qinit, 'symbolic', Control) ;
1005 */
1006
1007 operation = "symbolic factorization" ;
1008 do_numeric = FALSE ;
1009 Amatrix = (mxArray *) pargin [0] ;
1010
1011 if (nargout == 5)
1012 {
1013 do_info = 4 ;
1014 }
1015 if (op == 2)
1016 {
1017 User_Qinit = (mxArray *) pargin [1] ;
1018 }
1019 if (nargin < 2 || nargin > 4 || nargout > 5 || nargout < 4)
1020 {
1021 mexErrMsgTxt ("wrong number of arguments") ;
1022 }
1023
1024 }
1025 else if (MATCH (operator, "det"))
1026 {
1027
1028 /* -------------------------------------------------------------- */
1029 /* compute the determinant */
1030 /* -------------------------------------------------------------- */
1031
1032 /*
1033 * [det] = umfpack (A, 'det') ;
1034 * [dmantissa dexp] = umfpack (A, 'det') ;
1035 */
1036
1037 operation = "determinant" ;
1038 do_det = TRUE ;
1039 Amatrix = (mxArray *) pargin [0] ;
1040 if (nargin > 2 || nargout > 2)
1041 {
1042 mexErrMsgTxt ("wrong number of arguments") ;
1043 }
1044
1045 }
1046 else
1047 {
1048 mexErrMsgTxt ("operator must be '/', '\\', or 'symbolic'") ;
1049 }
1050 mxFree (operator) ;
1051
1052 }
1053 else if (nargin > 0)
1054 {
1055
1056 /* ------------------------------------------------------------------ */
1057 /* LU factorization */
1058 /* ------------------------------------------------------------------ */
1059
1060 /*
1061 with scaling:
1062 [L, U, P, Q, R, Info] = umfpack (A) ;
1063 [L, U, P, Q, R, Info] = umfpack (A, Qinit) ;
1064
1065 scaling determined by Control settings:
1066 [L, U, P, Q, R, Info] = umfpack (A, Control) ;
1067 [L, U, P, Q, R, Info] = umfpack (A, Qinit, Control) ;
1068
1069 with no scaling:
1070 [L, U, P, Q] = umfpack (A) ;
1071 [L, U, P, Q] = umfpack (A, Control) ;
1072 [L, U, P, Q] = umfpack (A, Qinit) ;
1073 [L, U, P, Q] = umfpack (A, Qinit, Control) ;
1074 */
1075
1076 operation = "numeric factorization" ;
1077 Amatrix = (mxArray *) pargin [0] ;
1078
1079 no_scale = (nargout <= 4) ;
1080
1081 if (nargout == 6)
1082 {
1083 do_info = 5 ;
1084 }
1085 if (nargin >= 2 && mxIsDouble (pargin [1]))
1086 {
1087 User_Qinit = (mxArray *) pargin [1] ;
1088 }
1089 if (nargin > 3 || nargout > 6 || nargout < 4)
1090 {
1091 mexErrMsgTxt ("wrong number of arguments") ;
1092 }
1093
1094 }
1095 else
1096 {
1097
1098 /* ------------------------------------------------------------------ */
1099 /* return default control settings */
1100 /* ------------------------------------------------------------------ */
1101
1102 /*
1103 Control = umfpack ;
1104 umfpack ;
1105 */
1106
1107 if (nargout > 1)
1108 {
1109 mexErrMsgTxt ("wrong number of arguments") ;
1110 }
1111
1112 /* return default opts struct */
1113 pargout [0] = umfpack_mx_defaults ( ) ;
1114 return ;
1115 }
1116
1117 /* ---------------------------------------------------------------------- */
1118 /* check inputs */
1119 /* ---------------------------------------------------------------------- */
1120
1121 if (mxGetNumberOfDimensions (Amatrix) != 2)
1122 {
1123 mexErrMsgTxt ("input matrix A must be 2-dimensional") ;
1124 }
1125 n_row = mxGetM (Amatrix) ;
1126 n_col = mxGetN (Amatrix) ;
1127 nn = MAX (n_row, n_col) ;
1128 n_inner = MIN (n_row, n_col) ;
1129 if (do_solve && n_row != n_col)
1130 {
1131 mexErrMsgTxt ("input matrix A must square for '\\' or '/'") ;
1132 }
1133 if (!mxIsSparse (Amatrix))
1134 {
1135 mexErrMsgTxt ("input matrix A must be sparse") ;
1136 }
1137 if (n_row == 0 || n_col == 0)
1138 {
1139 mexErrMsgTxt ("input matrix A cannot have zero rows or zero columns") ;
1140 }
1141
1142 /* The real/complex status of A determines which version to use, */
1143 /* (umfpack_dl_* or umfpack_zl_*). */
1144 A_is_complex = mxIsComplex (Amatrix) ;
1145 Atype = A_is_complex ? mxCOMPLEX : mxREAL ;
1146 Ap = (Long *) mxGetJc (Amatrix) ;
1147 Ai = (Long *) mxGetIr (Amatrix) ;
1148 Ax = mxGetPr (Amatrix) ;
1149 Az = mxGetPi (Amatrix) ;
1150
1151 if (do_solve)
1152 {
1153
1154 if (n_row != n_col)
1155 {
1156 mexErrMsgTxt ("A must be square for \\ or /") ;
1157 }
1158 if (transpose)
1159 {
1160 if (mxGetM (Bmatrix) != 1 || mxGetN (Bmatrix) != nn)
1161 {
1162 mexErrMsgTxt ("b has the wrong dimensions") ;
1163 }
1164 }
1165 else
1166 {
1167 if (mxGetM (Bmatrix) != nn || mxGetN (Bmatrix) != 1)
1168 {
1169 mexErrMsgTxt ("b has the wrong dimensions") ;
1170 }
1171 }
1172 if (mxGetNumberOfDimensions (Bmatrix) != 2)
1173 {
1174 mexErrMsgTxt ("input matrix b must be 2-dimensional") ;
1175 }
1176 if (mxIsSparse (Bmatrix))
1177 {
1178 mexErrMsgTxt ("input matrix b cannot be sparse") ;
1179 }
1180 if (mxGetClassID (Bmatrix) != mxDOUBLE_CLASS)
1181 {
1182 mexErrMsgTxt ("input matrix b must double precision matrix") ;
1183 }
1184
1185 B_is_complex = mxIsComplex (Bmatrix) ;
1186 Bx = mxGetPr (Bmatrix) ;
1187 Bz = mxGetPi (Bmatrix) ;
1188
1189 X_is_complex = A_is_complex || B_is_complex ;
1190 Xtype = X_is_complex ? mxCOMPLEX : mxREAL ;
1191 }
1192
1193 /* ---------------------------------------------------------------------- */
1194 /* set the Control parameters */
1195 /* ---------------------------------------------------------------------- */
1196
1197 if (A_is_complex)
1198 {
1199 umfpack_zl_defaults (Control) ;
1200 }
1201 else
1202 {
1203 umfpack_dl_defaults (Control) ;
1204 }
1205
1206 info_details = 0 ;
1207 if (User_Control_struct != NULL)
1208 {
1209 info_details = get_all_options (User_Control_struct, Control) ;
1210 }
1211
1212 if (no_scale)
1213 {
1214 /* turn off scaling for [L, U, P, Q] = umfpack (A) ;
1215 * ignoring the input value of Control (24) for the usage
1216 * [L, U, P, Q] = umfpack (A, Control) ; */
1217 Control [UMFPACK_SCALE] = UMFPACK_SCALE_NONE ;
1218 }
1219
1220 print_level = (Long) Control [UMFPACK_PRL] ;
1221
1222 /* ---------------------------------------------------------------------- */
1223 /* get Qinit, if present */
1224 /* ---------------------------------------------------------------------- */
1225
1226 if (User_Qinit)
1227 {
1228 if (mxGetM (User_Qinit) != 1 || mxGetN (User_Qinit) != n_col)
1229 {
1230 mexErrMsgTxt ("Qinit must be 1-by-n_col") ;
1231 }
1232 if (mxGetNumberOfDimensions (User_Qinit) != 2)
1233 {
1234 mexErrMsgTxt ("input Qinit must be 2-dimensional") ;
1235 }
1236 if (mxIsComplex (User_Qinit))
1237 {
1238 mexErrMsgTxt ("input Qinit must not be complex") ;
1239 }
1240 if (mxGetClassID (User_Qinit) != mxDOUBLE_CLASS)
1241 {
1242 mexErrMsgTxt ("input Qinit must be a double matrix") ;
1243 }
1244 if (mxIsSparse (User_Qinit))
1245 {
1246 mexErrMsgTxt ("input Qinit must be dense") ;
1247 }
1248 Qinit = (Long *) mxMalloc (n_col * sizeof (Long)) ;
1249 p = mxGetPr (User_Qinit) ;
1250 for (k = 0 ; k < n_col ; k++)
1251 {
1252 /* convert from 1-based to 0-based */
1253 Qinit [k] = ((Long) (p [k])) - 1 ;
1254 }
1255 Control [UMFPACK_ORDERING] = UMFPACK_ORDERING_GIVEN ;
1256 }
1257 else
1258 {
1259 /* umfpack_*_qsymbolic will call colamd to get Qinit. This is the */
1260 /* same as calling umfpack_*_symbolic with Qinit set to NULL*/
1261 Qinit = (Long *) NULL ;
1262 }
1263
1264 /* ---------------------------------------------------------------------- */
1265 /* report the inputs A and Qinit */
1266 /* ---------------------------------------------------------------------- */
1267
1268 if (print_level >= 2)
1269 {
1270 /* print the operation */
1271 mexPrintf ("\numfpack: %s\n", operation) ;
1272 }
1273
1274 if (A_is_complex)
1275 {
1276 umfpack_zl_report_control (Control) ;
1277 if (print_level >= 3) mexPrintf ("\nA: ") ;
1278 (void) umfpack_zl_report_matrix (n_row, n_col, Ap, Ai, Ax, Az,
1279 1, Control) ;
1280 if (Qinit)
1281 {
1282 if (print_level >= 3) mexPrintf ("\nQinit: ") ;
1283 (void) umfpack_zl_report_perm (n_col, Qinit, Control) ;
1284 }
1285 }
1286 else
1287 {
1288 umfpack_dl_report_control (Control) ;
1289 if (print_level >= 3) mexPrintf ("\nA: ") ;
1290 (void) umfpack_dl_report_matrix (n_row, n_col, Ap, Ai, Ax,
1291 1, Control) ;
1292 if (Qinit)
1293 {
1294 if (print_level >= 3) mexPrintf ("\nQinit: ") ;
1295 (void) umfpack_dl_report_perm (n_col, Qinit, Control) ;
1296 }
1297 }
1298
1299 /* ---------------------------------------------------------------------- */
1300 /* perform the symbolic factorization */
1301 /* ---------------------------------------------------------------------- */
1302
1303 if (A_is_complex)
1304 {
1305 status = umfpack_zl_qsymbolic (n_row, n_col, Ap, Ai, Ax, Az,
1306 Qinit, &Symbolic, Control, Info) ;
1307 }
1308 else
1309 {
1310 status = umfpack_dl_qsymbolic (n_row, n_col, Ap, Ai, Ax,
1311 Qinit, &Symbolic, Control, Info) ;
1312 }
1313
1314 if (Qinit)
1315 {
1316 mxFree (Qinit) ;
1317 }
1318
1319 if (status < 0)
1320 {
1321 error ("symbolic factorization failed", A_is_complex, nargout, pargout,
1322 Control, Info, status) ;
1323 return ;
1324 }
1325
1326 /* ---------------------------------------------------------------------- */
1327 /* report the Symbolic object */
1328 /* ---------------------------------------------------------------------- */
1329
1330 if (A_is_complex)
1331 {
1332 (void) umfpack_zl_report_symbolic (Symbolic, Control) ;
1333 }
1334 else
1335 {
1336 (void) umfpack_dl_report_symbolic (Symbolic, Control) ;
1337 }
1338
1339 /* ---------------------------------------------------------------------- */
1340 /* perform numeric factorization, or just return symbolic factorization */
1341 /* ---------------------------------------------------------------------- */
1342
1343 if (do_numeric)
1344 {
1345
1346 /* ------------------------------------------------------------------ */
1347 /* perform the numeric factorization */
1348 /* ------------------------------------------------------------------ */
1349
1350 if (A_is_complex)
1351 {
1352 status = umfpack_zl_numeric (Ap, Ai, Ax, Az, Symbolic, &Numeric,
1353 Control, Info) ;
1354 }
1355 else
1356 {
1357 status = umfpack_dl_numeric (Ap, Ai, Ax, Symbolic, &Numeric,
1358 Control, Info) ;
1359 }
1360
1361 /* ------------------------------------------------------------------ */
1362 /* free the symbolic factorization */
1363 /* ------------------------------------------------------------------ */
1364
1365 if (A_is_complex)
1366 {
1367 umfpack_zl_free_symbolic (&Symbolic) ;
1368 }
1369 else
1370 {
1371 umfpack_dl_free_symbolic (&Symbolic) ;
1372 }
1373
1374 /* ------------------------------------------------------------------ */
1375 /* report the Numeric object */
1376 /* ------------------------------------------------------------------ */
1377
1378 if (status < 0)
1379 {
1380 error ("numeric factorization failed", A_is_complex, nargout,
1381 pargout, Control, Info, status);
1382 return ;
1383 }
1384
1385 if (A_is_complex)
1386 {
1387 (void) umfpack_zl_report_numeric (Numeric, Control) ;
1388 }
1389 else
1390 {
1391 (void) umfpack_dl_report_numeric (Numeric, Control) ;
1392 }
1393
1394 /* ------------------------------------------------------------------ */
1395 /* return the solution, determinant, or the factorization */
1396 /* ------------------------------------------------------------------ */
1397
1398 if (do_solve)
1399 {
1400 /* -------------------------------------------------------------- */
1401 /* solve Ax=b or A'x'=b', and return just the solution x */
1402 /* -------------------------------------------------------------- */
1403
1404 if (transpose)
1405 {
1406 /* If A is real, A'x=b is the same as A.'x=b. */
1407 /* x and b are vectors, so x and b are the same as x' and b'. */
1408 /* If A is complex, then A.'x.'=b.' gives the same solution x */
1409 /* as the complex conjugate transpose. If we used the A'x=b */
1410 /* option in umfpack_*_solve, we would have to form b' on */
1411 /* input and x' on output (negating the imaginary part). */
1412 /* We can save this work by just using the A.'x=b option in */
1413 /* umfpack_*_solve. Then, forming x.' and b.' is implicit, */
1414 /* since x and b are just vectors anyway. */
1415 /* In both cases, the system to solve is A.'x=b */
1416 pargout [0] = mxCreateDoubleMatrix (1, nn, Xtype) ;
1417 sys = UMFPACK_Aat ;
1418 }
1419 else
1420 {
1421 pargout [0] = mxCreateDoubleMatrix (nn, 1, Xtype) ;
1422 sys = UMFPACK_A ;
1423 }
1424
1425 /* -------------------------------------------------------------- */
1426 /* print the right-hand-side, B */
1427 /* -------------------------------------------------------------- */
1428
1429 if (print_level >= 3) mexPrintf ("\nright-hand side, b: ") ;
1430 if (B_is_complex)
1431 {
1432 (void) umfpack_zl_report_vector (nn, Bx, Bz, Control) ;
1433 }
1434 else
1435 {
1436 (void) umfpack_dl_report_vector (nn, Bx, Control) ;
1437 }
1438
1439 /* -------------------------------------------------------------- */
1440 /* solve the system */
1441 /* -------------------------------------------------------------- */
1442
1443 Xx = mxGetPr (pargout [0]) ;
1444 Xz = mxGetPi (pargout [0]) ;
1445 status2 = UMFPACK_OK ;
1446
1447 if (A_is_complex)
1448 {
1449 if (!B_is_complex)
1450 {
1451 /* umfpack_zl_solve expects a complex B */
1452 Bz = (double *) mxCalloc (nn, sizeof (double)) ;
1453 }
1454 status = umfpack_zl_solve (sys, Ap, Ai, Ax, Az, Xx, Xz, Bx, Bz,
1455 Numeric, Control, Info) ;
1456 if (!B_is_complex)
1457 {
1458 mxFree (Bz) ;
1459 }
1460 }
1461 else
1462 {
1463 if (B_is_complex)
1464 {
1465 /* Ax=b when b is complex and A is sparse can be split */
1466 /* into two systems, A*xr=br and A*xi=bi, where r denotes */
1467 /* the real part and i the imaginary part of x and b. */
1468 status2 = umfpack_dl_solve (sys, Ap, Ai, Ax, Xz, Bz,
1469 Numeric, Control, Info) ;
1470 }
1471 status = umfpack_dl_solve (sys, Ap, Ai, Ax, Xx, Bx,
1472 Numeric, Control, Info) ;
1473 }
1474
1475 /* -------------------------------------------------------------- */
1476 /* free the Numeric object */
1477 /* -------------------------------------------------------------- */
1478
1479 if (A_is_complex)
1480 {
1481 umfpack_zl_free_numeric (&Numeric) ;
1482 }
1483 else
1484 {
1485 umfpack_dl_free_numeric (&Numeric) ;
1486 }
1487
1488 /* -------------------------------------------------------------- */
1489 /* check error status */
1490 /* -------------------------------------------------------------- */
1491
1492 if (status < 0 || status2 < 0)
1493 {
1494 mxDestroyArray (pargout [0]) ;
1495 error ("solve failed", A_is_complex, nargout, pargout, Control,
1496 Info, status) ;
1497 return ;
1498 }
1499
1500 /* -------------------------------------------------------------- */
1501 /* print the solution, X */
1502 /* -------------------------------------------------------------- */
1503
1504 if (print_level >= 3) mexPrintf ("\nsolution, x: ") ;
1505 if (X_is_complex)
1506 {
1507 (void) umfpack_zl_report_vector (nn, Xx, Xz, Control) ;
1508 }
1509 else
1510 {
1511 (void) umfpack_dl_report_vector (nn, Xx, Control) ;
1512 }
1513
1514 /* -------------------------------------------------------------- */
1515 /* warn about singular or near-singular matrices */
1516 /* -------------------------------------------------------------- */
1517
1518 /* no warning is given if Control (1) is zero */
1519
1520 if (print_level >= 1)
1521 {
1522 if (status == UMFPACK_WARNING_singular_matrix)
1523 {
1524 mexWarnMsgTxt (
1525 "matrix is singular\n"
1526 "Try increasing opts.tol and opts.symtol.\n"
1527 "Suppress this warning with opts.prl=0\n") ;
1528 }
1529 else if (Info [UMFPACK_RCOND] < DBL_EPSILON)
1530 {
1531 sprintf (warning, "matrix is nearly singular, rcond = %g\n"
1532 "Try increasing opts.tol and opts.symtol.\n"
1533 "Suppress this warning with opts.prl=0\n",
1534 Info [UMFPACK_RCOND]) ;
1535 mexWarnMsgTxt (warning) ;
1536 }
1537 }
1538
1539 }
1540 else if (do_det)
1541 {
1542
1543 /* -------------------------------------------------------------- */
1544 /* get the determinant */
1545 /* -------------------------------------------------------------- */
1546
1547 if (nargout == 2)
1548 {
1549 /* [det dexp] = umfpack (A, 'det') ;
1550 * return determinant in the form det * 10^dexp */
1551 p = &dexp ;
1552 }
1553 else
1554 {
1555 /* [det] = umfpack (A, 'det') ;
1556 * return determinant as a single scalar (overflow or
1557 * underflow is much more likely) */
1558 p = (double *) NULL ;
1559 }
1560 if (A_is_complex)
1561 {
1562 status = umfpack_zl_get_determinant (&dx, &dz, p,
1563 Numeric, Info) ;
1564 umfpack_zl_free_numeric (&Numeric) ;
1565 }
1566 else
1567 {
1568 status = umfpack_dl_get_determinant (&dx, p,
1569 Numeric, Info) ;
1570 umfpack_dl_free_numeric (&Numeric) ;
1571 dz = 0 ;
1572 }
1573 if (status < 0)
1574 {
1575 error ("extracting LU factors failed", A_is_complex, nargout,
1576 pargout, Control, Info, status) ;
1577 }
1578 if (A_is_complex)
1579 {
1580 pargout [0] = mxCreateDoubleMatrix (1, 1, mxCOMPLEX) ;
1581 p = mxGetPr (pargout [0]) ;
1582 *p = dx ;
1583 p = mxGetPi (pargout [0]) ;
1584 *p = dz ;
1585 }
1586 else
1587 {
1588 pargout [0] = mxCreateDoubleMatrix (1, 1, mxREAL) ;
1589 p = mxGetPr (pargout [0]) ;
1590 *p = dx ;
1591 }
1592 if (nargout == 2)
1593 {
1594 pargout [1] = mxCreateDoubleMatrix (1, 1, mxREAL) ;
1595 p = mxGetPr (pargout [1]) ;
1596 *p = dexp ;
1597 }
1598
1599 }
1600 else
1601 {
1602
1603 /* -------------------------------------------------------------- */
1604 /* get L, U, P, Q, and r */
1605 /* -------------------------------------------------------------- */
1606
1607 if (A_is_complex)
1608 {
1609 status = umfpack_zl_get_lunz (&lnz, &unz, &ignore1, &ignore2,
1610 &ignore3, Numeric) ;
1611 }
1612 else
1613 {
1614 status = umfpack_dl_get_lunz (&lnz, &unz, &ignore1, &ignore2,
1615 &ignore3, Numeric) ;
1616 }
1617
1618 if (status < 0)
1619 {
1620 if (A_is_complex)
1621 {
1622 umfpack_zl_free_numeric (&Numeric) ;
1623 }
1624 else
1625 {
1626 umfpack_dl_free_numeric (&Numeric) ;
1627 }
1628 error ("extracting LU factors failed", A_is_complex, nargout,
1629 pargout, Control, Info, status) ;
1630 return ;
1631 }
1632
1633 /* avoid malloc of zero-sized arrays */
1634 lnz = MAX (lnz, 1) ;
1635 unz = MAX (unz, 1) ;
1636
1637 /* get temporary space, for the *** ROW *** form of L */
1638 Ltp = (Long *) mxMalloc ((n_row+1) * sizeof (Long)) ;
1639 Ltj = (Long *) mxMalloc (lnz * sizeof (Long)) ;
1640 Ltx = (double *) mxMalloc (lnz * sizeof (double)) ;
1641 if (A_is_complex)
1642 {
1643 Ltz = (double *) mxMalloc (lnz * sizeof (double)) ;
1644 }
1645 else
1646 {
1647 Ltz = (double *) NULL ;
1648 }
1649
1650 /* create permanent copy of the output matrix U */
1651 pargout [1] = mxCreateSparse (n_inner, n_col, unz, Atype) ;
1652 Up = (Long *) mxGetJc (pargout [1]) ;
1653 Ui = (Long *) mxGetIr (pargout [1]) ;
1654 Ux = mxGetPr (pargout [1]) ;
1655 Uz = mxGetPi (pargout [1]) ;
1656
1657 /* temporary space for the integer permutation vectors */
1658 P = (Long *) mxMalloc (n_row * sizeof (Long)) ;
1659 Q = (Long *) mxMalloc (n_col * sizeof (Long)) ;
1660
1661 /* get scale factors, if requested */
1662 status2 = UMFPACK_OK ;
1663 if (!no_scale)
1664 {
1665 /* create a diagonal sparse matrix for the scale factors */
1666 pargout [4] = mxCreateSparse (n_row, n_row, n_row, mxREAL) ;
1667 Rp = (Long *) mxGetJc (pargout [4]) ;
1668 Ri = (Long *) mxGetIr (pargout [4]) ;
1669 for (i = 0 ; i < n_row ; i++)
1670 {
1671 Rp [i] = i ;
1672 Ri [i] = i ;
1673 }
1674 Rp [n_row] = n_row ;
1675 Rs = mxGetPr (pargout [4]) ;
1676 }
1677 else
1678 {
1679 Rs = (double *) NULL ;
1680 }
1681
1682 /* get Lt, U, P, Q, and Rs from the numeric object */
1683 if (A_is_complex)
1684 {
1685 status = umfpack_zl_get_numeric (Ltp, Ltj, Ltx, Ltz, Up, Ui, Ux,
1686 Uz, P, Q, (double *) NULL, (double *) NULL,
1687 &do_recip, Rs, Numeric) ;
1688 umfpack_zl_free_numeric (&Numeric) ;
1689 }
1690 else
1691 {
1692 status = umfpack_dl_get_numeric (Ltp, Ltj, Ltx, Up, Ui,
1693 Ux, P, Q, (double *) NULL,
1694 &do_recip, Rs, Numeric) ;
1695 umfpack_dl_free_numeric (&Numeric) ;
1696 }
1697
1698 /* for the mexFunction, -DNRECIPROCAL must be set,
1699 * so do_recip must be FALSE */
1700
1701 if (status < 0 || status2 < 0 || do_recip)
1702 {
1703 mxFree (Ltp) ;
1704 mxFree (Ltj) ;
1705 mxFree (Ltx) ;
1706 if (Ltz) mxFree (Ltz) ;
1707 mxFree (P) ;
1708 mxFree (Q) ;
1709 mxDestroyArray (pargout [1]) ;
1710 error ("extracting LU factors failed", A_is_complex, nargout,
1711 pargout, Control, Info, status) ;
1712 return ;
1713 }
1714
1715 /* create sparse permutation matrix for P */
1716 pargout [2] = mxCreateSparse (n_row, n_row, n_row, mxREAL) ;
1717 Pp = (Long *) mxGetJc (pargout [2]) ;
1718 Pi = (Long *) mxGetIr (pargout [2]) ;
1719 Px = mxGetPr (pargout [2]) ;
1720 for (k = 0 ; k < n_row ; k++)
1721 {
1722 Pp [k] = k ;
1723 Px [k] = 1 ;
1724 Pi [P [k]] = k ;
1725 }
1726 Pp [n_row] = n_row ;
1727
1728 /* create sparse permutation matrix for Q */
1729 pargout [3] = mxCreateSparse (n_col, n_col, n_col, mxREAL) ;
1730 Qp = (Long *) mxGetJc (pargout [3]) ;
1731 Qi = (Long *) mxGetIr (pargout [3]) ;
1732 Qx = mxGetPr (pargout [3]) ;
1733 for (k = 0 ; k < n_col ; k++)
1734 {
1735 Qp [k] = k ;
1736 Qx [k] = 1 ;
1737 Qi [k] = Q [k] ;
1738 }
1739 Qp [n_col] = n_col ;
1740
1741 /* permanent copy of L */
1742 pargout [0] = mxCreateSparse (n_row, n_inner, lnz, Atype) ;
1743 Lp = (Long *) mxGetJc (pargout [0]) ;
1744 Li = (Long *) mxGetIr (pargout [0]) ;
1745 Lx = mxGetPr (pargout [0]) ;
1746 Lz = mxGetPi (pargout [0]) ;
1747
1748 /* convert L from row form to column form */
1749 if (A_is_complex)
1750 {
1751 /* non-conjugate array transpose */
1752 status = umfpack_zl_transpose (n_inner, n_row, Ltp, Ltj, Ltx,
1753 Ltz, (Long *) NULL, (Long *) NULL, Lp, Li, Lx, Lz,
1754 FALSE) ;
1755 }
1756 else
1757 {
1758 status = umfpack_dl_transpose (n_inner, n_row, Ltp, Ltj, Ltx,
1759 (Long *) NULL, (Long *) NULL, Lp, Li, Lx) ;
1760 }
1761
1762 mxFree (Ltp) ;
1763 mxFree (Ltj) ;
1764 mxFree (Ltx) ;
1765 if (Ltz) mxFree (Ltz) ;
1766
1767 if (status < 0)
1768 {
1769 mxFree (P) ;
1770 mxFree (Q) ;
1771 mxDestroyArray (pargout [0]) ;
1772 mxDestroyArray (pargout [1]) ;
1773 mxDestroyArray (pargout [2]) ;
1774 mxDestroyArray (pargout [3]) ;
1775 error ("constructing L failed", A_is_complex, nargout, pargout,
1776 Control, Info, status) ;
1777 return ;
1778 }
1779
1780 /* -------------------------------------------------------------- */
1781 /* print L, U, P, and Q */
1782 /* -------------------------------------------------------------- */
1783
1784 if (A_is_complex)
1785 {
1786 if (print_level >= 3) mexPrintf ("\nL: ") ;
1787 (void) umfpack_zl_report_matrix (n_row, n_inner, Lp, Li,
1788 Lx, Lz, 1, Control) ;
1789 if (print_level >= 3) mexPrintf ("\nU: ") ;
1790 (void) umfpack_zl_report_matrix (n_inner, n_col, Up, Ui,
1791 Ux, Uz, 1, Control) ;
1792 if (print_level >= 3) mexPrintf ("\nP: ") ;
1793 (void) umfpack_zl_report_perm (n_row, P, Control) ;
1794 if (print_level >= 3) mexPrintf ("\nQ: ") ;
1795 (void) umfpack_zl_report_perm (n_col, Q, Control) ;
1796 }
1797 else
1798 {
1799 if (print_level >= 3) mexPrintf ("\nL: ") ;
1800 (void) umfpack_dl_report_matrix (n_row, n_inner, Lp, Li,
1801 Lx, 1, Control) ;
1802 if (print_level >= 3) mexPrintf ("\nU: ") ;
1803 (void) umfpack_dl_report_matrix (n_inner, n_col, Up, Ui,
1804 Ux, 1, Control) ;
1805 if (print_level >= 3) mexPrintf ("\nP: ") ;
1806 (void) umfpack_dl_report_perm (n_row, P, Control) ;
1807 if (print_level >= 3) mexPrintf ("\nQ: ") ;
1808 (void) umfpack_dl_report_perm (n_col, Q, Control) ;
1809 }
1810
1811 mxFree (P) ;
1812 mxFree (Q) ;
1813
1814 }
1815
1816 }
1817 else
1818 {
1819
1820 /* ------------------------------------------------------------------ */
1821 /* return the symbolic factorization */
1822 /* ------------------------------------------------------------------ */
1823
1824 Q = (Long *) mxMalloc (n_col * sizeof (Long)) ;
1825 P = (Long *) mxMalloc (n_row * sizeof (Long)) ;
1826 Front_npivcol = (Long *) mxMalloc ((nn+1) * sizeof (Long)) ;
1827 Front_parent = (Long *) mxMalloc ((nn+1) * sizeof (Long)) ;
1828 Front_1strow = (Long *) mxMalloc ((nn+1) * sizeof (Long)) ;
1829 Front_leftmostdesc = (Long *) mxMalloc ((nn+1) * sizeof (Long)) ;
1830 Chain_start = (Long *) mxMalloc ((nn+1) * sizeof (Long)) ;
1831 Chain_maxrows = (Long *) mxMalloc ((nn+1) * sizeof (Long)) ;
1832 Chain_maxcols = (Long *) mxMalloc ((nn+1) * sizeof (Long)) ;
1833
1834 if (A_is_complex)
1835 {
1836 status = umfpack_zl_get_symbolic (&ignore1, &ignore2, &ignore3,
1837 &nz, &nfronts, &nchains, P, Q, Front_npivcol,
1838 Front_parent, Front_1strow, Front_leftmostdesc,
1839 Chain_start, Chain_maxrows, Chain_maxcols, Symbolic) ;
1840 umfpack_zl_free_symbolic (&Symbolic) ;
1841 }
1842 else
1843 {
1844 status = umfpack_dl_get_symbolic (&ignore1, &ignore2, &ignore3,
1845 &nz, &nfronts, &nchains, P, Q, Front_npivcol,
1846 Front_parent, Front_1strow, Front_leftmostdesc,
1847 Chain_start, Chain_maxrows, Chain_maxcols, Symbolic) ;
1848 umfpack_dl_free_symbolic (&Symbolic) ;
1849 }
1850
1851 if (status < 0)
1852 {
1853 mxFree (P) ;
1854 mxFree (Q) ;
1855 mxFree (Front_npivcol) ;
1856 mxFree (Front_parent) ;
1857 mxFree (Front_1strow) ;
1858 mxFree (Front_leftmostdesc) ;
1859 mxFree (Chain_start) ;
1860 mxFree (Chain_maxrows) ;
1861 mxFree (Chain_maxcols) ;
1862 error ("extracting symbolic factors failed", A_is_complex, nargout,
1863 pargout, Control, Info, status) ;
1864 return ;
1865 }
1866
1867 /* create sparse permutation matrix for P */
1868 pargout [0] = mxCreateSparse (n_row, n_row, n_row, mxREAL) ;
1869 Pp = (Long *) mxGetJc (pargout [0]) ;
1870 Pi = (Long *) mxGetIr (pargout [0]) ;
1871 Px = mxGetPr (pargout [0]) ;
1872 for (k = 0 ; k < n_row ; k++)
1873 {
1874 Pp [k] = k ;
1875 Px [k] = 1 ;
1876 Pi [P [k]] = k ;
1877 }
1878 Pp [n_row] = n_row ;
1879
1880 /* create sparse permutation matrix for Q */
1881 pargout [1] = mxCreateSparse (n_col, n_col, n_col, mxREAL) ;
1882 Qp = (Long *) mxGetJc (pargout [1]) ;
1883 Qi = (Long *) mxGetIr (pargout [1]) ;
1884 Qx = mxGetPr (pargout [1]) ;
1885 for (k = 0 ; k < n_col ; k++)
1886 {
1887 Qp [k] = k ;
1888 Qx [k] = 1 ;
1889 Qi [k] = Q [k] ;
1890 }
1891 Qp [n_col] = n_col ;
1892
1893 /* create Fr */
1894 pargout [2] = mxCreateDoubleMatrix (nfronts+1, 4, mxREAL) ;
1895
1896 p1 = mxGetPr (pargout [2]) ;
1897 p2 = p1 + nfronts + 1 ;
1898 p3 = p2 + nfronts + 1 ;
1899 p4 = p3 + nfronts + 1 ;
1900 for (i = 0 ; i <= nfronts ; i++)
1901 {
1902 /* convert parent, 1strow, and leftmostdesc to 1-based */
1903 p1 [i] = (double) (Front_npivcol [i]) ;
1904 p2 [i] = (double) (Front_parent [i] + 1) ;
1905 p3 [i] = (double) (Front_1strow [i] + 1) ;
1906 p4 [i] = (double) (Front_leftmostdesc [i] + 1) ;
1907 }
1908
1909 /* create Ch */
1910 pargout [3] = mxCreateDoubleMatrix (nchains+1, 3, mxREAL) ;
1911 p1 = mxGetPr (pargout [3]) ;
1912 p2 = p1 + nchains + 1 ;
1913 p3 = p2 + nchains + 1 ;
1914 for (i = 0 ; i < nchains ; i++)
1915 {
1916 p1 [i] = (double) (Chain_start [i] + 1) ; /* convert to 1-based */
1917 p2 [i] = (double) (Chain_maxrows [i]) ;
1918 p3 [i] = (double) (Chain_maxcols [i]) ;
1919 }
1920 p1 [nchains] = Chain_start [nchains] + 1 ;
1921 p2 [nchains] = 0 ;
1922 p3 [nchains] = 0 ;
1923
1924 mxFree (P) ;
1925 mxFree (Q) ;
1926 mxFree (Front_npivcol) ;
1927 mxFree (Front_parent) ;
1928 mxFree (Front_1strow) ;
1929 mxFree (Front_leftmostdesc) ;
1930 mxFree (Chain_start) ;
1931 mxFree (Chain_maxrows) ;
1932 mxFree (Chain_maxcols) ;
1933 }
1934
1935 /* ---------------------------------------------------------------------- */
1936 /* report Info */
1937 /* ---------------------------------------------------------------------- */
1938
1939 if (A_is_complex)
1940 {
1941 umfpack_zl_report_info (Control, Info) ;
1942 }
1943 else
1944 {
1945 umfpack_dl_report_info (Control, Info) ;
1946 }
1947
1948 if (do_info > 0)
1949 {
1950 /* return Info */
1951 if (info_details > 0)
1952 {
1953 pargout [do_info] = umfpack_mx_info_details (Control, Info) ;
1954 }
1955 else
1956 {
1957 pargout [do_info] = umfpack_mx_info_user (Control, Info, do_solve) ;
1958 }
1959 }
1960 }
1961