Mercurial > octave-nkf

comparison liboctave/UMFPACK/UMFPACK/MATLAB/umfpackmex.c @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children
comparison
equal deleted inserted replaced
5163:9f3299378193 5164:57077d0ddc8e
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 }