1 /* ========================================================================== */
2 /* === MATLAB/cholmod_matlab ================================================ */
3 /* ========================================================================== */
4
5 /* Utility routines for the CHOLMOD MATLAB mexFunctions.
6 *
7 * If CHOLMOD runs out of memory, MATLAB will terminate the mexFunction
8 * immediately since it uses mxMalloc (see sputil_config, below). Likewise,
9 * if mxCreate* or mxMalloc (as called in this file) fails, MATLAB will also
10 * terminate the mexFunction. When this occurs, MATLAB frees all allocated
11 * memory, so we don't have to worry about memory leaks. If this were not the
12 * case, the routines in this file would suffer from memory leaks whenever an
13 * error occurred.
14 */
15
16 #include "cholmod_matlab.h"
17
18 #ifndef INT64_T
19 #define INT64_T long long
20 #endif
21
22 /* This file pointer is used for the mread and mwrite mexFunctions. It must
23 * be a global variable, because the file pointer is not passed to the
24 * sputil_error_handler function when an error occurs. */
25 FILE *sputil_file = NULL ;
26
27 /* ========================================================================== */
28 /* === sputil_config ======================================================== */
29 /* ========================================================================== */
30
31 /* Define function pointers and other parameters for a mexFunction */
32
sputil_config(Long spumoni,cholmod_common * cm)33 void sputil_config (Long spumoni, cholmod_common *cm)
34 {
35 /* cholmod_l_solve must return a real or zomplex X for MATLAB */
36 cm->prefer_zomplex = TRUE ;
37
38 /* printing and error handling */
39 if (spumoni == 0)
40 {
41 /* do not print anything from within CHOLMOD */
42 cm->print = -1 ;
43 SuiteSparse_config.printf_func = NULL ;
44 }
45 else
46 {
47 /* spumoni = 1: print warning and error messages. cholmod_l_print_*
48 * routines will print a one-line summary of each object printed.
49 * spumoni = 2: also print a short summary of each object.
50 */
51 cm->print = spumoni + 2 ;
52 /* SuiteSparse_config.printf_func = mexPrintf : */
53 }
54
55 /* error handler */
56 cm->error_handler = sputil_error_handler ;
57
58 /* Turn off METIS memory guard. It is not needed, because mxMalloc will
59 * safely terminate the mexFunction and free any workspace without killing
60 * all of MATLAB. This assumes cholmod_make was used to compile CHOLMOD
61 * for MATLAB. */
62 cm->metis_memory = 0.0 ;
63 }
64
65
66 /* ========================================================================== */
67 /* === sputil_error_handler ================================================= */
68 /* ========================================================================== */
69
sputil_error_handler(int status,const char * file,int line,const char * message)70 void sputil_error_handler (int status, const char *file, int line,
71 const char *message)
72 {
73 if (status < CHOLMOD_OK)
74 {
75 /*
76 mexPrintf ("ERROR: file %s line %d, status %d\n", file, line, status) ;
77 */
78 if (sputil_file != NULL)
79 {
80 fclose (sputil_file) ;
81 sputil_file = NULL ;
82 }
83 mexErrMsgTxt (message) ;
84 }
85 /*
86 else
87 {
88 mexPrintf ("Warning: file %s line %d, status %d\n", file, line, status);
89 }
90 */
91 }
92
93
94 /* ========================================================================== */
95 /* === sputil_get_sparse ==================================================== */
96 /* ========================================================================== */
97
98 /* Create a shallow CHOLMOD copy of a MATLAB sparse matrix. No memory is
99 * allocated. The resulting matrix A must not be modified.
100 */
101
sputil_get_sparse(const mxArray * Amatlab,cholmod_sparse * A,double * dummy,Long stype)102 cholmod_sparse *sputil_get_sparse
103 (
104 const mxArray *Amatlab, /* MATLAB version of the matrix */
105 cholmod_sparse *A, /* CHOLMOD version of the matrix */
106 double *dummy, /* a pointer to a valid scalar double */
107 Long stype /* -1: lower, 0: unsymmetric, 1: upper */
108 )
109 {
110 Long *Ap ;
111 A->nrow = mxGetM (Amatlab) ;
112 A->ncol = mxGetN (Amatlab) ;
113 A->p = (Long *) mxGetJc (Amatlab) ;
114 A->i = (Long *) mxGetIr (Amatlab) ;
115 Ap = A->p ;
116 A->nzmax = Ap [A->ncol] ;
117 A->packed = TRUE ;
118 A->sorted = TRUE ;
119 A->nz = NULL ;
120 A->itype = CHOLMOD_LONG ; /* was CHOLMOD_INT in v1.6 and earlier */
121 A->dtype = CHOLMOD_DOUBLE ;
122 A->stype = stype ;
123
124 #ifndef MATLAB6p1_OR_EARLIER
125
126 if (mxIsLogical (Amatlab))
127 {
128 A->x = NULL ;
129 A->z = NULL ;
130 A->xtype = CHOLMOD_PATTERN ;
131 }
132 else if (mxIsEmpty (Amatlab))
133 {
134 /* this is not dereferenced, but the existence (non-NULL) of these
135 * pointers is checked in CHOLMOD */
136 A->x = dummy ;
137 A->z = dummy ;
138 A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
139 }
140 else if (mxIsDouble (Amatlab))
141 {
142 A->x = mxGetPr (Amatlab) ;
143 A->z = mxGetPi (Amatlab) ;
144 A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
145 }
146 else
147 {
148 /* only logical and complex/real double matrices supported */
149 sputil_error (ERROR_INVALID_TYPE, 0) ;
150 }
151
152 #else
153
154 if (mxIsEmpty (Amatlab))
155 {
156 /* this is not dereferenced, but the existence (non-NULL) of these
157 * pointers is checked in CHOLMOD */
158 A->x = dummy ;
159 A->z = dummy ;
160 A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
161 }
162 else
163 {
164 /* in MATLAB 6.1, the matrix is sparse, so it must be double */
165 A->x = mxGetPr (Amatlab) ;
166 A->z = mxGetPi (Amatlab) ;
167 A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
168 }
169
170 #endif
171
172 return (A) ;
173 }
174
175
176 /* ========================================================================== */
177 /* === sputil_get_dense ===================================================== */
178 /* ========================================================================== */
179
180 /* Create a shallow CHOLMOD copy of a MATLAB dense matrix. No memory is
181 * allocated. Only double (real and zomplex) matrices are supported. The
182 * resulting matrix B must not be modified.
183 */
184
sputil_get_dense(const mxArray * Amatlab,cholmod_dense * A,double * dummy)185 cholmod_dense *sputil_get_dense
186 (
187 const mxArray *Amatlab, /* MATLAB version of the matrix */
188 cholmod_dense *A, /* CHOLMOD version of the matrix */
189 double *dummy /* a pointer to a valid scalar double */
190 )
191 {
192 A->nrow = mxGetM (Amatlab) ;
193 A->ncol = mxGetN (Amatlab) ;
194 A->d = A->nrow ;
195 A->nzmax = A->nrow * A->ncol ;
196 A->dtype = CHOLMOD_DOUBLE ;
197
198 if (mxIsEmpty (Amatlab))
199 {
200 A->x = dummy ;
201 A->z = dummy ;
202 }
203 else if (mxIsDouble (Amatlab))
204 {
205 A->x = mxGetPr (Amatlab) ;
206 A->z = mxGetPi (Amatlab) ;
207 }
208 else
209 {
210 /* only full double matrices supported by sputil_get_dense */
211 sputil_error (ERROR_INVALID_TYPE, 0) ;
212 }
213 A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
214
215 return (A) ;
216 }
217
218
219 /* ========================================================================== */
220 /* === sputil_get_sparse_pattern ============================================ */
221 /* ========================================================================== */
222
223 /* Create a CHOLMOD_PATTERN sparse matrix for a MATLAB matrix, depending on the
224 * type:
225 *
226 * (1) MATLAB full real double: duplicate CHOLMOD_REAL sparse matrix.
227 * (2) MATLAB full complex double: duplicate CHOLMOD_ZOMPLEX sparse matrix.
228 * (3) MATLAB full logical: duplicate CHOLMOD_PATTERN sparse matrix.
229 * (4) MATLAB sparse real double: shallow CHOLMOD_REAL copy.
230 * (5) MATLAB sparse complex double: shallow CHOLMOD_ZOMPLEX copy.
231 * (6) MATLAB sparse logical: shallow CHOLMOD_PATTERN copy.
232 *
233 * A shallow copy or duplicate is returned; the shallow copy must not be freed.
234 * For a shallow copy, the return value A is the same as Ashallow. For a
235 * complete duplicate, A and Ashallow will differ.
236 */
237
sputil_get_sparse_pattern(const mxArray * Amatlab,cholmod_sparse * Ashallow,double * dummy,cholmod_common * cm)238 cholmod_sparse *sputil_get_sparse_pattern
239 (
240 const mxArray *Amatlab, /* MATLAB version of the matrix */
241 cholmod_sparse *Ashallow, /* shallow CHOLMOD version of the matrix */
242 double *dummy, /* a pointer to a valid scalar double */
243 cholmod_common *cm
244 )
245 {
246 cholmod_sparse *A = NULL ;
247
248 if (!mxIsSparse (Amatlab))
249 {
250
251 /* ------------------------------------------------------------------ */
252 /* A = sparse (X) where X is full */
253 /* ------------------------------------------------------------------ */
254
255 if (mxIsDouble (Amatlab))
256 {
257
258 /* -------------------------------------------------------------- */
259 /* convert full double X into sparse matrix A (pattern only) */
260 /* -------------------------------------------------------------- */
261
262 cholmod_dense Xmatrix, *X ;
263 X = sputil_get_dense (Amatlab, &Xmatrix, dummy) ;
264 A = cholmod_l_dense_to_sparse (X, FALSE, cm) ;
265
266 }
267
268 #ifndef MATLAB6p1_OR_EARLIER
269
270 else if (mxIsLogical (Amatlab))
271 {
272
273 /* -------------------------------------------------------------- */
274 /* convert full logical MATLAB matrix into CHOLMOD_PATTERN */
275 /* -------------------------------------------------------------- */
276
277 /* (this is copied and modified from t_cholmod_dense.c) */
278
279 char *x ;
280 Long *Ap, *Ai ;
281 Long nrow, ncol, i, j, nz, nzmax, p ;
282
283 /* -------------------------------------------------------------- */
284 /* count the number of nonzeros in the result */
285 /* -------------------------------------------------------------- */
286
287 nrow = mxGetM (Amatlab) ;
288 ncol = mxGetN (Amatlab) ;
289 x = (char *) mxGetData (Amatlab) ;
290 nzmax = nrow * ncol ;
291 for (nz = 0, j = 0 ; j < nzmax ; j++)
292 {
293 if (x [j])
294 {
295 nz++ ;
296 }
297 }
298
299 /* -------------------------------------------------------------- */
300 /* allocate the result A */
301 /* -------------------------------------------------------------- */
302
303 A = cholmod_l_allocate_sparse (nrow, ncol, nz, TRUE, TRUE, 0,
304 CHOLMOD_PATTERN, cm) ;
305
306 if (cm->status < CHOLMOD_OK)
307 {
308 return (NULL) ; /* out of memory */
309 }
310 Ap = A->p ;
311 Ai = A->i ;
312
313 /* -------------------------------------------------------------- */
314 /* copy the full logical matrix into the sparse matrix A */
315 /* -------------------------------------------------------------- */
316
317 p = 0 ;
318 for (j = 0 ; j < ncol ; j++)
319 {
320 Ap [j] = p ;
321 for (i = 0 ; i < nrow ; i++)
322 {
323 if (x [i+j*nrow])
324 {
325 Ai [p++] = i ;
326 }
327 }
328 }
329 /* ASSERT (p == nz) ; */
330 Ap [ncol] = nz ;
331 }
332
333 #endif
334
335 else
336 {
337 /* only double and logical matrices supported */
338 sputil_error (ERROR_INVALID_TYPE, 0) ;
339 }
340
341 }
342 else
343 {
344
345 /* ------------------------------------------------------------------ */
346 /* create a shallow copy of sparse matrix A (default stype is zero) */
347 /* ------------------------------------------------------------------ */
348
349 A = sputil_get_sparse (Amatlab, Ashallow, dummy, 0) ;
350 A->x = NULL ;
351 A->z = NULL ;
352 A->xtype = CHOLMOD_PATTERN ;
353 }
354
355 return (A) ;
356 }
357
358
359 /* ========================================================================== */
360 /* === sputil_put_sparse ==================================================== */
361 /* ========================================================================== */
362
363 /* Creates a true MATLAB version of a CHOLMOD sparse matrix. The CHOLMOD sparse
364 * matrix is destroyed. Both real and zomplex matrices are supported.
365 */
366
sputil_put_sparse(cholmod_sparse ** Ahandle,cholmod_common * cm)367 mxArray *sputil_put_sparse
368 (
369 cholmod_sparse **Ahandle, /* CHOLMOD version of the matrix */
370 cholmod_common *cm
371 )
372 {
373 mxArray *Amatlab ;
374 cholmod_sparse *A ;
375 A = *Ahandle ;
376 if (!A || A->x == NULL) mexErrMsgTxt ("invalid sparse matrix") ;
377 Amatlab = mxCreateSparse (0, 0, 0,
378 (A->xtype != CHOLMOD_REAL) ? mxCOMPLEX: mxREAL) ;
379 mxSetM (Amatlab, A->nrow) ;
380 mxSetN (Amatlab, A->ncol) ;
381 mxSetNzmax (Amatlab, A->nzmax) ;
382 MXFREE (mxGetJc (Amatlab)) ;
383 MXFREE (mxGetIr (Amatlab)) ;
384 MXFREE (mxGetPr (Amatlab)) ;
385 mxSetJc (Amatlab, A->p) ;
386 mxSetIr (Amatlab, A->i) ;
387 mxSetPr (Amatlab, A->x) ;
388 mexMakeMemoryPersistent (A->p) ;
389 mexMakeMemoryPersistent (A->i) ;
390 mexMakeMemoryPersistent (A->x) ;
391 if (A->xtype != CHOLMOD_REAL)
392 {
393 if (A->z == NULL) mexErrMsgTxt ("invalid complex sparse matrix") ;
394 MXFREE (mxGetPi (Amatlab)) ;
395 mxSetPi (Amatlab, A->z) ;
396 mexMakeMemoryPersistent (A->z) ;
397 }
398 A->p = NULL ;
399 A->i = NULL ;
400 A->x = NULL ;
401 A->z = NULL ;
402 cholmod_l_free_sparse (Ahandle, cm) ;
403 return (Amatlab) ;
404 }
405
406
407 /* ========================================================================== */
408 /* === sputil_put_dense ===================================================== */
409 /* ========================================================================== */
410
411 /* Creates a true MATLAB version of a CHOLMOD dense matrix. The CHOLMOD dense
412 * matrix is destroyed. Both real and zomplex matrices are supported.
413 */
414
sputil_put_dense(cholmod_dense ** Ahandle,cholmod_common * cm)415 mxArray *sputil_put_dense
416 (
417 cholmod_dense **Ahandle, /* CHOLMOD version of the matrix */
418 cholmod_common *cm
419 )
420 {
421 mxArray *Amatlab ;
422 cholmod_dense *A ;
423 A = *Ahandle ;
424 if (!A || A->x == NULL) mexErrMsgTxt ("invalid dense matrix") ;
425 Amatlab = mxCreateDoubleMatrix (0, 0,
426 (A->xtype != CHOLMOD_REAL) ? mxCOMPLEX: mxREAL) ;
427 mxSetM (Amatlab, A->nrow) ;
428 mxSetN (Amatlab, A->ncol) ;
429 MXFREE (mxGetPr (Amatlab)) ;
430 mxSetPr (Amatlab, A->x) ;
431 mexMakeMemoryPersistent (A->x) ;
432 if (A->xtype != CHOLMOD_REAL)
433 {
434 if (A->z == NULL) mexErrMsgTxt ("invalid complex dense matrix") ;
435 MXFREE (mxGetPi (Amatlab)) ;
436 mxSetPi (Amatlab, A->z) ;
437 mexMakeMemoryPersistent (A->z) ;
438 }
439 A->x = NULL ;
440 A->z = NULL ;
441 cholmod_l_free_dense (Ahandle, cm) ;
442 return (Amatlab) ;
443 }
444
445
446 /* ========================================================================== */
447 /* === sputil_put_int ======================================================= */
448 /* ========================================================================== */
449
450 /* Convert a Long vector into a double mxArray */
451
sputil_put_int(Long * P,Long n,Long one_based)452 mxArray *sputil_put_int
453 (
454 Long *P, /* vector to convert */
455 Long n, /* length of P */
456 Long one_based /* 1 if convert from 0-based to 1-based, 0 otherwise */
457 )
458 {
459 double *p ;
460 mxArray *Q ;
461 Long i ;
462 Q = mxCreateDoubleMatrix (1, n, mxREAL) ;
463 p = mxGetPr (Q) ;
464 for (i = 0 ; i < n ; i++)
465 {
466 p [i] = (double) (P [i] + one_based) ;
467 }
468 return (Q) ;
469 }
470
471
472 /* ========================================================================== */
473 /* === sputil_error ========================================================= */
474 /* ========================================================================== */
475
476 /* An integer is out of range, or other error has occurred. */
477
sputil_error(Long error,Long is_index)478 void sputil_error
479 (
480 Long error, /* kind of error */
481 Long is_index /* TRUE if a matrix index, FALSE if a matrix dimension */
482 )
483 {
484 if (error == ERROR_TOO_SMALL)
485 {
486 mexErrMsgTxt (is_index ?
487 "sparse: index into matrix must be positive" :
488 "sparse: sparse matrix sizes must be non-negative integers") ;
489 }
490 else if (error == ERROR_HUGE)
491 {
492 mexErrMsgTxt (is_index ?
493 "sparse: index into matrix is too large" :
494 "sparse: sparse matrix size is too large") ;
495 }
496 else if (error == ERROR_NOT_INTEGER)
497 {
498 mexErrMsgTxt (is_index ?
499 "sparse: index into matrix must be an integer" :
500 "sparse: sparse matrix size must be an integer") ;
501 }
502 else if (error == ERROR_TOO_LARGE)
503 {
504 mexErrMsgTxt ("sparse: index exceeds matrix dimensions") ;
505 }
506 else if (error == ERROR_USAGE)
507 {
508 mexErrMsgTxt (
509 "Usage:\n"
510 "A = sparse (S)\n"
511 "A = sparse (i,j,s,m,n,nzmax)\n"
512 "A = sparse (i,j,s,m,n)\n"
513 "A = sparse (i,j,s)\n"
514 "A = sparse (m,n)\n") ;
515 }
516 else if (error == ERROR_LENGTH)
517 {
518 mexErrMsgTxt ("sparse: vectors must be the same lengths") ;
519 }
520 else if (error == ERROR_INVALID_TYPE)
521 {
522 mexErrMsgTxt ("matrix class not supported") ;
523 }
524 }
525
526
527 /* ========================================================================== */
528 /* === sputil_double_to_int ================================================= */
529 /* ========================================================================== */
530
531 /* convert a double into an integer */
532
sputil_double_to_int(double x,Long is_index,Long n)533 Long sputil_double_to_int /* returns integer value of x */
534 (
535 double x, /* double value to convert */
536 Long is_index, /* TRUE if a matrix index, FALSE if a matrix dimension */
537 Long n /* if a matrix index, x cannot exceed this dimension,
538 * except that -1 is treated as infinity */
539 )
540 {
541 Long i ;
542 if (x > INT_MAX)
543 {
544 /* x is way too big for an integer */
545 sputil_error (ERROR_HUGE, is_index) ;
546 }
547 else if (x < 0)
548 {
549 /* x must be non-negative */
550 sputil_error (ERROR_TOO_SMALL, is_index) ;
551 }
552 i = (Long) x ;
553 if (x != (double) i)
554 {
555 /* x must be an integer */
556 sputil_error (ERROR_NOT_INTEGER, is_index) ;
557 }
558 if (is_index)
559 {
560 if (i < 1)
561 {
562 sputil_error (ERROR_TOO_SMALL, is_index) ;
563 }
564 else if (i > n && n != EMPTY)
565 {
566 sputil_error (ERROR_TOO_LARGE, is_index) ;
567 }
568 }
569 return (i) ;
570 }
571
572
573 /* ========================================================================== */
574 /* === sputil_nelements ===================================================== */
575 /* ========================================================================== */
576
577 /* return the number of elements in an mxArray. Trigger an error on integer
578 * overflow (in case the argument is sparse) */
579
sputil_nelements(const mxArray * arg)580 Long sputil_nelements (const mxArray *arg)
581 {
582 double size ;
583 const Long *dims ;
584 Long k, ndims ;
585 ndims = mxGetNumberOfDimensions (arg) ;
586 dims = (Long *) mxGetDimensions (arg) ;
587 size = 1 ;
588 for (k = 0 ; k < ndims ; k++)
589 {
590 size *= dims [k] ;
591 }
592 return (sputil_double_to_int (size, FALSE, 0)) ;
593 }
594
595
596 /* ========================================================================== */
597 /* === sputil_get_double ==================================================== */
598 /* ========================================================================== */
599
sputil_get_double(const mxArray * arg)600 double sputil_get_double (const mxArray *arg)
601 {
602 if (sputil_nelements (arg) < 1)
603 {
604 /* [] is not a scalar, but its value is zero so that
605 * sparse ([],[],[]) is a 0-by-0 matrix */
606 return (0) ;
607 }
608 return (mxGetScalar (arg)) ;
609 }
610
611
612 /* ========================================================================== */
613 /* === sputil_get_integer =================================================== */
614 /* ========================================================================== */
615
616 /* return an argument as a non-negative integer scalar, or -1 if error */
617
sputil_get_integer(const mxArray * arg,Long is_index,Long n)618 Long sputil_get_integer
619 (
620 const mxArray *arg, /* MATLAB argument to convert */
621 Long is_index, /* TRUE if an index, FALSE if a matrix dimension */
622 Long n /* maximum value, if an index */
623 )
624 {
625 double x = sputil_get_double (arg) ;
626 if (mxIsInf (x) || mxIsNaN (x))
627 {
628 /* arg is Inf or NaN, return -1 */
629 return (EMPTY) ;
630 }
631 return (sputil_double_to_int (x, is_index, n)) ;
632 }
633
634
635 /* ========================================================================== */
636 /* === sputil_trim ========================================================== */
637 /* ========================================================================== */
638
639 /* Remove columns k to n-1 from a sparse matrix S, leaving columns 0 to k-1.
640 * S must be packed (there can be no S->nz array). This condition is not
641 * checked, since only packed matrices are passed to this routine. */
642
sputil_trim(cholmod_sparse * S,Long k,cholmod_common * cm)643 void sputil_trim
644 (
645 cholmod_sparse *S,
646 Long k,
647 cholmod_common *cm
648 )
649 {
650 Long *Sp ;
651 Long ncol ;
652 size_t n1, nznew ;
653
654 if (S == NULL)
655 {
656 return ;
657 }
658
659 ncol = S->ncol ;
660 if (k < 0 || k >= ncol)
661 {
662 /* do not modify S */
663 return ;
664 }
665
666 /* reduce S->p in size. This cannot fail. */
667 n1 = ncol + 1 ;
668 S->p = cholmod_l_realloc (k+1, sizeof (Long), S->p, &n1, cm) ;
669
670 /* get the new number of entries in S */
671 Sp = S->p ;
672 nznew = Sp [k] ;
673
674 /* reduce S->i, S->x, and S->z (if present) to size nznew */
675 cholmod_l_reallocate_sparse (nznew, S, cm) ;
676
677 /* S now has only k columns */
678 S->ncol = k ;
679 }
680
681
682 /* ========================================================================== */
683 /* === sputil_extract_zeros ================================================= */
684 /* ========================================================================== */
685
686 /* Create a sparse binary (real double) matrix Z that contains the pattern
687 * of explicit zeros in the sparse real/zomplex double matrix A. */
688
sputil_extract_zeros(cholmod_sparse * A,cholmod_common * cm)689 cholmod_sparse *sputil_extract_zeros
690 (
691 cholmod_sparse *A,
692 cholmod_common *cm
693 )
694 {
695 Long *Ap, *Ai, *Zp, *Zi ;
696 double *Ax, *Az, *Zx ;
697 Long j, p, nzeros = 0, is_complex, pz, nrow, ncol ;
698 cholmod_sparse *Z ;
699
700 if (A == NULL || A->xtype == CHOLMOD_PATTERN || A->xtype == CHOLMOD_COMPLEX)
701 {
702 /* only sparse real/zomplex double matrices supported */
703 sputil_error (ERROR_INVALID_TYPE, 0) ;
704 }
705
706 Ap = A->p ;
707 Ai = A->i ;
708 Ax = A->x ;
709 Az = A->z ;
710 ncol = A->ncol ;
711 nrow = A->nrow ;
712 is_complex = (A->xtype == CHOLMOD_ZOMPLEX) ;
713
714 /* count the number of zeros in a sparse matrix A */
715 for (j = 0 ; j < ncol ; j++)
716 {
717 for (p = Ap [j] ; p < Ap [j+1] ; p++)
718 {
719 if (CHOLMOD_IS_ZERO (Ax [p]) &&
720 ((is_complex) ? CHOLMOD_IS_ZERO (Az [p]) : TRUE))
721 {
722 nzeros++ ;
723 }
724 }
725 }
726
727 /* allocate the Z matrix with space for all the zero entries */
728 Z = cholmod_l_spzeros (nrow, ncol, nzeros, CHOLMOD_REAL, cm) ;
729
730 /* extract the zeros from A and store them in Z as binary values */
731 if (nzeros > 0)
732 {
733 Zp = Z->p ;
734 Zi = Z->i ;
735 Zx = Z->x ;
736 pz = 0 ;
737 for (j = 0 ; j < ncol ; j++)
738 {
739 Zp [j] = pz ;
740 for (p = Ap [j] ; p < Ap [j+1] ; p++)
741 {
742 if (CHOLMOD_IS_ZERO (Ax [p]) &&
743 ((is_complex) ? CHOLMOD_IS_ZERO (Az [p]) : TRUE))
744 {
745 Zi [pz] = Ai [p] ;
746 Zx [pz] = 1 ;
747 pz++ ;
748 }
749 }
750 }
751 Zp [ncol] = pz ;
752 }
753
754 return (Z) ;
755 }
756
757
758 /* ========================================================================== */
759 /* === sputil_drop_zeros ==================================================== */
760 /* ========================================================================== */
761
762 /* Drop zeros from a packed CHOLMOD sparse matrix (zomplex or real). This is
763 * very similar to CHOLMOD/MatrixOps/cholmod_drop, except that this routine has
764 * no tolerance parameter and it can handle zomplex matrices. NaN's are left
765 * in the matrix. If this is used on the sparse matrix version of the factor
766 * L, then the update/downdate methods cannot be applied to L (ldlupdate).
767 * Returns the number of entries dropped.
768 */
769
sputil_drop_zeros(cholmod_sparse * S)770 Long sputil_drop_zeros
771 (
772 cholmod_sparse *S
773 )
774 {
775 double sik, zik ;
776 Long *Sp, *Si ;
777 double *Sx, *Sz ;
778 Long pdest, k, ncol, p, pend, nz ;
779
780 if (S == NULL)
781 {
782 return (0) ;
783 }
784
785 Sp = S->p ;
786 Si = S->i ;
787 Sx = S->x ;
788 Sz = S->z ;
789 pdest = 0 ;
790 ncol = S->ncol ;
791 nz = Sp [ncol] ;
792
793 if (S->xtype == CHOLMOD_ZOMPLEX)
794 {
795 for (k = 0 ; k < ncol ; k++)
796 {
797 p = Sp [k] ;
798 pend = Sp [k+1] ;
799 Sp [k] = pdest ;
800 for ( ; p < pend ; p++)
801 {
802 sik = Sx [p] ;
803 zik = Sz [p] ;
804 if (CHOLMOD_IS_NONZERO (sik) || CHOLMOD_IS_NONZERO (zik))
805 {
806 if (p != pdest)
807 {
808 Si [pdest] = Si [p] ;
809 Sx [pdest] = sik ;
810 Sz [pdest] = zik ;
811 }
812 pdest++ ;
813 }
814 }
815 }
816 }
817 else
818 {
819 for (k = 0 ; k < ncol ; k++)
820 {
821 p = Sp [k] ;
822 pend = Sp [k+1] ;
823 Sp [k] = pdest ;
824 for ( ; p < pend ; p++)
825 {
826 sik = Sx [p] ;
827 if (CHOLMOD_IS_NONZERO (sik))
828 {
829 if (p != pdest)
830 {
831 Si [pdest] = Si [p] ;
832 Sx [pdest] = sik ;
833 }
834 pdest++ ;
835 }
836 }
837 }
838 }
839 Sp [ncol] = pdest ;
840 return (nz - pdest) ;
841 }
842
843
844 /* ========================================================================== */
845 /* === sputil_copy_ij ======================================================= */
846 /* ========================================================================== */
847
848 /* copy i or j arguments into an Long vector. For small integer types, i and
849 * and j can be returned with negative entries; this error condition is caught
850 * later, in cholmod_triplet_to_sparse.
851 *
852 * Future work: if the mxClassID matches the default integer (INT32 for 32-bit
853 * MATLAB and INT64 for 64-bit), then it would save memory to patch in the
854 * vector with a pointer copy, rather than making a copy of the whole vector.
855 * This would require that the 1-based i and j vectors be converted on the fly
856 * to 0-based vectors in cholmod_triplet_to_sparse.
857 */
858
sputil_copy_ij(Long is_scalar,Long scalar,void * vector,mxClassID category,Long nz,Long n,Long * I)859 Long sputil_copy_ij /* returns the dimension, n */
860 (
861 Long is_scalar, /* TRUE if argument is a scalar, FALSE otherwise */
862 Long scalar, /* scalar value of the argument */
863 void *vector, /* vector value of the argument */
864 mxClassID category, /* type of vector */
865 Long nz, /* length of output vector I */
866 Long n, /* maximum dimension, EMPTY if not yet known */
867 Long *I /* vector of length nz to copy into */
868 )
869 {
870 Long i, k, ok, ok2, ok3, n2 ;
871
872 if (is_scalar)
873 {
874 n2 = scalar ;
875 if (n == EMPTY)
876 {
877 n = scalar ;
878 }
879 i = scalar - 1 ;
880 for (k = 0 ; k < nz ; k++)
881 {
882 I [k] = i ;
883 }
884 }
885 else
886 {
887 /* copy double input into Long vector (convert to 0-based) */
888 ok = TRUE ;
889 ok2 = TRUE ;
890 ok3 = TRUE ;
891 n2 = 0 ;
892
893 switch (category)
894 {
895
896 #ifndef MATLAB6p1_OR_EARLIER
897
898 /* MATLAB 6.1 or earlier do not have mxLOGICAL_CLASS */
899 case mxLOGICAL_CLASS:
900
901 for (k = 0 ; k < nz ; k++)
902 {
903 i = (Long) (((mxLogical *) vector) [k]) ;
904 I [k] = i - 1 ;
905 n2 = MAX (n2, i) ;
906 }
907 break ;
908
909 #endif
910
911 case mxCHAR_CLASS:
912
913 for (k = 0 ; k < nz ; k++)
914 {
915 i = (Long) (((mxChar *) vector) [k]) ;
916 I [k] = i - 1 ;
917 n2 = MAX (n2, i) ;
918 }
919 break ;
920
921 case mxINT8_CLASS:
922
923 for (k = 0 ; k < nz ; k++)
924 {
925 i = (Long) (((INT8_T *) vector) [k]) ;
926 I [k] = i - 1 ;
927 n2 = MAX (n2, i) ;
928 }
929 break ;
930
931 case mxUINT8_CLASS:
932
933 for (k = 0 ; k < nz ; k++)
934 {
935 i = (Long) (((UINT8_T *) vector) [k]) ;
936 I [k] = i - 1 ;
937 n2 = MAX (n2, i) ;
938 }
939 break ;
940
941 case mxINT16_CLASS:
942
943 for (k = 0 ; k < nz ; k++)
944 {
945 i = (Long) (((INT16_T *) vector) [k]) ;
946 I [k] = i - 1 ;
947 n2 = MAX (n2, i) ;
948 }
949 break ;
950
951 case mxUINT16_CLASS:
952
953 for (k = 0 ; k < nz ; k++)
954 {
955 i = (Long) (((UINT16_T *) vector) [k]) ;
956 I [k] = i - 1 ;
957 n2 = MAX (n2, i) ;
958 }
959 break ;
960
961 case mxINT32_CLASS:
962
963 for (k = 0 ; k < nz ; k++)
964 {
965 i = (Long) (((INT32_T *) vector) [k]) ;
966 I [k] = i - 1 ;
967 n2 = MAX (n2, i) ;
968 }
969 break ;
970
971 case mxUINT32_CLASS:
972
973 for (k = 0 ; ok3 && k < nz ; k++)
974 {
975 double y = (((UINT32_T *) vector) [k]) ;
976 i = (Long) y ;
977 ok3 = (y < Long_max) ;
978 I [k] = i - 1 ;
979 n2 = MAX (n2, i) ;
980 }
981 break ;
982
983 case mxINT64_CLASS:
984
985 for (k = 0 ; ok2 && ok3 && k < nz ; k++)
986 {
987 INT64_T y = ((INT64_T *) vector) [k] ;
988 i = (Long) y ;
989 ok2 = (y > 0) ;
990 ok3 = (y < Long_max) ;
991 I [k] = i - 1 ;
992 n2 = MAX (n2, i) ;
993 }
994 break ;
995
996 case mxUINT64_CLASS:
997
998 for (k = 0 ; ok2 && ok3 && k < nz ; k++)
999 {
1000 unsigned INT64_T y = ((unsigned INT64_T *) vector) [k] ;
1001 i = (Long) y ;
1002 ok2 = (y > 0) ;
1003 ok3 = (y < Long_max) ;
1004 I [k] = i - 1 ;
1005 n2 = MAX (n2, i) ;
1006 }
1007 break ;
1008
1009 case mxSINGLE_CLASS:
1010
1011 for (k = 0 ; ok && ok2 && ok3 && k < nz ; k++)
1012 {
1013 float y = ((float *) vector) [k] ;
1014 i = (Long) y ;
1015 ok = (y == (float) i) ;
1016 ok2 = (y > 0) ;
1017 ok3 = (y < Long_max) ;
1018 I [k] = i - 1 ;
1019 n2 = MAX (n2, i) ;
1020 }
1021 break ;
1022
1023 case mxDOUBLE_CLASS:
1024
1025 for (k = 0 ; ok && ok2 && ok3 && k < nz ; k++)
1026 {
1027 double y = ((double *) vector) [k] ;
1028 i = (Long) y ;
1029 ok = (y == (double) i) ;
1030 ok2 = (y > 0) ;
1031 ok3 = (y < Long_max) ;
1032 I [k] = i - 1 ;
1033 n2 = MAX (n2, i) ;
1034 }
1035
1036 break ;
1037
1038 default:
1039
1040 sputil_error (ERROR_INVALID_TYPE, FALSE) ;
1041 break ;
1042 }
1043
1044 if (!ok)
1045 {
1046 sputil_error (ERROR_NOT_INTEGER, TRUE) ;
1047 }
1048
1049 if (!ok2)
1050 {
1051 sputil_error (ERROR_TOO_SMALL, TRUE) ;
1052 }
1053
1054 if (!ok3)
1055 {
1056 sputil_error (ERROR_HUGE, TRUE) ;
1057 }
1058
1059 }
1060 return ((n == EMPTY) ? n2 : n) ;
1061 }
1062
1063
1064 /* ========================================================================== */
1065 /* === sputil_dense_to_sparse =============================================== */
1066 /* ========================================================================== */
1067
1068 /* Convert a dense matrix of any numeric type into a
1069 * sparse double or sparse logical matrix.
1070 */
1071
1072 #define COUNT_NZ \
1073 { \
1074 for (j = 0 ; j < ncol ; j++) \
1075 { \
1076 for (i = 0 ; i < nrow ; i++) \
1077 { \
1078 xij = X [i + j*nrow] ; \
1079 if (CHOLMOD_IS_NONZERO (xij)) \
1080 { \
1081 nz++ ; \
1082 } \
1083 } \
1084 } \
1085 }
1086
1087 #define COPY_DENSE_TO_SPARSE(stype) \
1088 { \
1089 stype *Sx ; \
1090 Sp = (Long *) mxGetJc (S) ; \
1091 Si = (Long *) mxGetIr (S) ; \
1092 Sx = (stype *) mxGetData (S) ; \
1093 nz = 0 ; \
1094 for (j = 0 ; j < ncol ; j++) \
1095 { \
1096 Sp [j] = nz ; \
1097 for (i = 0 ; i < nrow ; i++) \
1098 { \
1099 xij = X [i + j*nrow] ; \
1100 if (CHOLMOD_IS_NONZERO (xij)) \
1101 { \
1102 Si [nz] = i ; \
1103 Sx [nz] = (stype) xij ; \
1104 nz++ ; \
1105 } \
1106 } \
1107 } \
1108 Sp [ncol] = nz ; \
1109 }
1110
1111 #define DENSE_TO_SPARSE(type) \
1112 { \
1113 type *X, xij ; \
1114 X = (type *) mxGetData (arg) ; \
1115 COUNT_NZ ; \
1116 S = mxCreateSparse (nrow, ncol, nz, mxREAL) ; \
1117 COPY_DENSE_TO_SPARSE (double) ; \
1118 }
1119
sputil_dense_to_sparse(const mxArray * arg)1120 mxArray *sputil_dense_to_sparse (const mxArray *arg)
1121 {
1122 mxArray *S = NULL ;
1123 Long *Sp, *Si ;
1124 Long nrow, ncol, nz, i, j ;
1125
1126 nrow = mxGetM (arg) ;
1127 ncol = mxGetN (arg) ;
1128 nz = 0 ;
1129
1130 if (mxIsComplex (arg))
1131 {
1132
1133 /* ------------------------------------------------------------------ */
1134 /* convert a complex dense matrix into a complex sparse matrix */
1135 /* ------------------------------------------------------------------ */
1136
1137 double xij, zij ;
1138 double *X, *Z, *Sx, *Sz ;
1139
1140 if (mxGetClassID (arg) != mxDOUBLE_CLASS)
1141 {
1142 /* A complex matrix can have any class (int8, int16, single, etc),
1143 * but this function only supports complex double. This condition
1144 * is not checked in the caller. */
1145 sputil_error (ERROR_INVALID_TYPE, FALSE) ;
1146 }
1147
1148 X = mxGetPr (arg) ;
1149 Z = mxGetPi (arg) ;
1150 for (j = 0 ; j < ncol ; j++)
1151 {
1152 for (i = 0 ; i < nrow ; i++)
1153 {
1154 xij = X [i + j*nrow] ;
1155 zij = Z [i + j*nrow] ;
1156 if (CHOLMOD_IS_NONZERO (xij) || CHOLMOD_IS_NONZERO (zij))
1157 {
1158 nz++ ;
1159 }
1160 }
1161 }
1162 S = mxCreateSparse (nrow, ncol, nz, mxCOMPLEX) ;
1163 Sp = (Long *) mxGetJc (S) ;
1164 Si = (Long *) mxGetIr (S) ;
1165 Sx = mxGetPr (S) ;
1166 Sz = mxGetPi (S) ;
1167 nz = 0 ;
1168 for (j = 0 ; j < ncol ; j++)
1169 {
1170 Sp [j] = nz ;
1171 for (i = 0 ; i < nrow ; i++)
1172 {
1173 xij = X [i + j*nrow] ;
1174 zij = Z [i + j*nrow] ;
1175 if (CHOLMOD_IS_NONZERO (xij) || CHOLMOD_IS_NONZERO (zij))
1176 {
1177 Si [nz] = i ;
1178 Sx [nz] = xij ;
1179 Sz [nz] = zij ;
1180 nz++ ;
1181 }
1182 }
1183 }
1184 Sp [ncol] = nz ;
1185
1186 }
1187 else
1188 {
1189
1190 /* ------------------------------------------------------------------ */
1191 /* convert real matrix (any class) to sparse double or logical */
1192 /* ------------------------------------------------------------------ */
1193
1194 switch (mxGetClassID (arg))
1195 {
1196
1197 #ifndef MATLAB6p1_OR_EARLIER
1198
1199 /* MATLAB 6.1 or earlier do not have mxLOGICAL_CLASS */
1200 case mxLOGICAL_CLASS:
1201 {
1202 mxLogical *X, xij ;
1203 X = (mxLogical *) mxGetData (arg) ;
1204 COUNT_NZ ;
1205 S = mxCreateSparseLogicalMatrix (nrow, ncol, nz) ;
1206 COPY_DENSE_TO_SPARSE (mxLogical) ;
1207 }
1208 break ;
1209 #endif
1210
1211 case mxCHAR_CLASS:
1212
1213 DENSE_TO_SPARSE (mxChar) ;
1214 break ;
1215
1216 case mxINT8_CLASS:
1217
1218 DENSE_TO_SPARSE (char) ;
1219 break ;
1220
1221 case mxUINT8_CLASS:
1222
1223 DENSE_TO_SPARSE (unsigned char) ;
1224 break ;
1225
1226 case mxINT16_CLASS:
1227
1228 DENSE_TO_SPARSE (short) ;
1229 break ;
1230
1231 case mxUINT16_CLASS:
1232
1233 DENSE_TO_SPARSE (unsigned short) ;
1234 break ;
1235
1236 case mxINT32_CLASS:
1237
1238 DENSE_TO_SPARSE (INT32_T) ;
1239 break ;
1240
1241 case mxUINT32_CLASS:
1242
1243 DENSE_TO_SPARSE (unsigned INT32_T) ;
1244 break ;
1245
1246 case mxINT64_CLASS:
1247
1248 DENSE_TO_SPARSE (INT64_T) ;
1249 break ;
1250
1251 case mxUINT64_CLASS:
1252
1253 DENSE_TO_SPARSE (unsigned INT64_T) ;
1254 break ;
1255
1256 case mxSINGLE_CLASS:
1257
1258 DENSE_TO_SPARSE (float) ;
1259 break ;
1260
1261 case mxDOUBLE_CLASS:
1262
1263 DENSE_TO_SPARSE (double) ;
1264 break ;
1265
1266 default:
1267
1268 sputil_error (ERROR_INVALID_TYPE, FALSE) ;
1269 break ;
1270 }
1271 }
1272
1273 return (S) ;
1274 }
1275
1276
1277 /* ========================================================================== */
1278 /* === sputil_triplet_to_sparse ============================================= */
1279 /* ========================================================================== */
1280
1281 /* Convert a triplet form into a sparse matrix. If complex, s must be double.
1282 * If real, s can be of any class.
1283 */
1284
sputil_triplet_to_sparse(Long nrow,Long ncol,Long nz,Long nzmax,Long i_is_scalar,Long i,void * i_vector,mxClassID i_class,Long j_is_scalar,Long j,void * j_vector,mxClassID j_class,Long s_is_scalar,double x,double z,void * x_vector,double * z_vector,mxClassID s_class,Long s_complex,cholmod_common * cm)1285 cholmod_sparse *sputil_triplet_to_sparse
1286 (
1287 Long nrow, Long ncol, Long nz, Long nzmax,
1288 Long i_is_scalar, Long i, void *i_vector, mxClassID i_class,
1289 Long j_is_scalar, Long j, void *j_vector, mxClassID j_class,
1290 Long s_is_scalar, double x, double z, void *x_vector, double *z_vector,
1291 mxClassID s_class, Long s_complex,
1292 cholmod_common *cm
1293 )
1294 {
1295 double dummy = 0 ;
1296 cholmod_triplet *T ;
1297 cholmod_sparse *S ;
1298 double *Tx, *Tz ;
1299 Long *Ti, *Tj ;
1300 Long k, x_patch ;
1301
1302 /* ---------------------------------------------------------------------- */
1303 /* allocate the triplet form */
1304 /* ---------------------------------------------------------------------- */
1305
1306 /* Note that nrow and ncol may be EMPTY; this is not an error condition.
1307 * Allocate the numerical part of T only if s is a scalar. */
1308 x_patch = (!s_is_scalar && (s_class == mxDOUBLE_CLASS || s_complex)) ;
1309
1310 T = cholmod_l_allocate_triplet (MAX (0,nrow), MAX (0,ncol), nz, 0,
1311 x_patch ? CHOLMOD_PATTERN :
1312 (s_complex ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL), cm) ;
1313 Ti = T->i ;
1314 Tj = T->j ;
1315 Tx = T->x ;
1316 Tz = T->z ;
1317
1318 /* ---------------------------------------------------------------------- */
1319 /* fill the triplet form */
1320 /* ---------------------------------------------------------------------- */
1321
1322 if (s_is_scalar)
1323 {
1324
1325 /* ------------------------------------------------------------------ */
1326 /* fill T->x and T->z with a scalar value */
1327 /* ------------------------------------------------------------------ */
1328
1329 for (k = 0 ; k < nz ; k++)
1330 {
1331 Tx [k] = x ;
1332 }
1333 if (s_complex)
1334 {
1335 for (k = 0 ; k < nz ; k++)
1336 {
1337 Tz [k] = z ;
1338 }
1339 }
1340
1341 }
1342 else
1343 {
1344
1345 /* ------------------------------------------------------------------ */
1346 /* copy x/z_vector into T->x and T->z, and convert to double */
1347 /* ------------------------------------------------------------------ */
1348
1349 if (s_complex)
1350 {
1351
1352 /* Patch in s as the numerical values of the triplet matrix.
1353 * Note that T->x and T->z must not be free'd when done. */
1354 T->x = (x_vector == NULL) ? &dummy : x_vector ;
1355 T->z = (z_vector == NULL) ? &dummy : z_vector ;
1356 T->xtype = CHOLMOD_ZOMPLEX ;
1357
1358 }
1359 else switch (s_class)
1360 {
1361
1362 #ifndef MATLAB6p1_OR_EARLIER
1363
1364 /* MATLAB 6.1 or earlier do not have mxLOGICAL_CLASS */
1365 case mxLOGICAL_CLASS:
1366
1367 for (k = 0 ; k < nz ; k++)
1368 {
1369 Tx [k] = (double) (((mxLogical *) x_vector) [k]) ;
1370 }
1371 break ;
1372
1373 #endif
1374
1375 case mxCHAR_CLASS:
1376
1377 for (k = 0 ; k < nz ; k++)
1378 {
1379 Tx [k] = (double) (((mxChar *) x_vector) [k]) ;
1380 }
1381 break ;
1382
1383 case mxINT8_CLASS:
1384
1385 for (k = 0 ; k < nz ; k++)
1386 {
1387 Tx [k] = (double) (((INT8_T *) x_vector) [k]) ;
1388 }
1389 break ;
1390
1391 case mxUINT8_CLASS:
1392
1393 for (k = 0 ; k < nz ; k++)
1394 {
1395 Tx [k] = (double) (((UINT8_T *) x_vector) [k]) ;
1396 }
1397 break ;
1398
1399 case mxINT16_CLASS:
1400
1401 for (k = 0 ; k < nz ; k++)
1402 {
1403 Tx [k] = (double) (((INT16_T *) x_vector) [k]) ;
1404 }
1405 break ;
1406
1407 case mxUINT16_CLASS:
1408
1409 for (k = 0 ; k < nz ; k++)
1410 {
1411 Tx [k] = (double) (((UINT16_T *) x_vector) [k]) ;
1412 }
1413 break ;
1414
1415 case mxINT32_CLASS:
1416
1417 for (k = 0 ; k < nz ; k++)
1418 {
1419 Tx [k] = (double) (((INT32_T *) x_vector) [k]) ;
1420 }
1421 break ;
1422
1423 case mxUINT32_CLASS:
1424
1425 for (k = 0 ; k < nz ; k++)
1426 {
1427 Tx [k] = (double) (((UINT32_T *) x_vector) [k]) ;
1428 }
1429 break ;
1430
1431 case mxINT64_CLASS:
1432
1433 for (k = 0 ; k < nz ; k++)
1434 {
1435 Tx [k] = (double) (((INT64_T *) x_vector) [k]) ;
1436 }
1437 break ;
1438
1439 case mxUINT64_CLASS:
1440
1441 for (k = 0 ; k < nz ; k++)
1442 {
1443 Tx [k] = (double) (((unsigned INT64_T *) x_vector) [k]) ;
1444 }
1445 break ;
1446
1447 case mxSINGLE_CLASS:
1448
1449 for (k = 0 ; k < nz ; k++)
1450 {
1451 Tx [k] = (double) (((float *) x_vector) [k]) ;
1452 }
1453 break ;
1454
1455 case mxDOUBLE_CLASS:
1456
1457 /* Patch in s as the numerical values of the triplet matrix.
1458 * Note that T->x must not be free'd when done. */
1459 T->x = (x_vector == NULL) ? &dummy : x_vector ;
1460 T->xtype = CHOLMOD_REAL ;
1461 break ;
1462
1463 default:
1464
1465 sputil_error (ERROR_INVALID_TYPE, FALSE) ;
1466 break ;
1467 }
1468 }
1469
1470 /* copy i in to the integer vector T->i */
1471 nrow = sputil_copy_ij (i_is_scalar, i, i_vector, i_class, nz, nrow, Ti) ;
1472
1473 /* copy j in to the integer vector T->j */
1474 ncol = sputil_copy_ij (j_is_scalar, j, j_vector, j_class, nz, ncol, Tj) ;
1475
1476 /* nrow and ncol are known */
1477 T->nrow = nrow ;
1478 T->ncol = ncol ;
1479 T->nnz = nz ;
1480
1481 /* ---------------------------------------------------------------------- */
1482 /* convert triplet to sparse matrix */
1483 /* ---------------------------------------------------------------------- */
1484
1485 /* If the triplet matrix T is invalid, or if CHOLMOD runs out of memory,
1486 * then S is NULL. */
1487 S = cholmod_l_triplet_to_sparse (T, nzmax, cm) ;
1488
1489 /* ---------------------------------------------------------------------- */
1490 /* free workspace */
1491 /* ---------------------------------------------------------------------- */
1492
1493 /* do not free T->x or T->z if it points to input x_vector */
1494 if (x_patch)
1495 {
1496 T->x = NULL ;
1497 T->z = NULL ;
1498 T->xtype = CHOLMOD_PATTERN ;
1499 }
1500 cholmod_l_free_triplet (&T, cm) ;
1501 return (S) ;
1502 }
1503
1504
1505 /* ========================================================================== */
1506 /* === sputil_copy_sparse =================================================== */
1507 /* ========================================================================== */
1508
1509 /* copy a sparse matrix, S = sparse(A), dropping any zero entries and ensuring
1510 * the nzmax(S) == nnz(S). Explicit zero entries in A "cannot" occur, in
1511 * the current version of MATLAB ... but a user mexFunction might generate a
1512 * matrix with explicit zeros. This function ensures S=sparse(A) drops those
1513 * explicit zeros. */
1514
sputil_copy_sparse(const mxArray * A)1515 mxArray *sputil_copy_sparse (const mxArray *A)
1516 {
1517 double aij, zij ;
1518 mxArray *S ;
1519 double *Ax, *Az, *Sx, *Sz ;
1520 Long *Ap, *Ai, *Sp, *Si ;
1521 Long anz, snz, p, j, nrow, ncol, pend ;
1522
1523 #ifndef MATLAB6p1_OR_EARLIER
1524
1525 /* MATLAB 6.1 or earlier : all sparse matrices are OK */
1526 if (! (mxGetClassID (A) == mxLOGICAL_CLASS ||
1527 mxGetClassID (A) == mxDOUBLE_CLASS))
1528 {
1529 /* Only sparse logical and real/complex double matrices supported.
1530 * This condition is not checked in the caller. */
1531 sputil_error (ERROR_INVALID_TYPE, 0) ;
1532 }
1533
1534 #endif
1535
1536 nrow = mxGetM (A) ;
1537 ncol = mxGetN (A) ;
1538 Ap = (Long *) mxGetJc (A) ;
1539 Ai = (Long *) mxGetIr (A) ;
1540 anz = Ap [ncol] ;
1541
1542 snz = 0 ;
1543
1544 #ifndef MATLAB6p1_OR_EARLIER
1545
1546 /* MATLAB 6.1 or earlier do not have mxLOGICAL_CLASS */
1547 if (mxIsLogical (A))
1548 {
1549
1550 /* ------------------------------------------------------------------ */
1551 /* copy a sparse logical matrix */
1552 /* ------------------------------------------------------------------ */
1553
1554 /* count the number of nonzeros in A */
1555 mxLogical *Al, *Sl ;
1556 Al = mxGetLogicals (A) ;
1557 for (p = 0 ; p < anz ; p++)
1558 {
1559 if (Al [p])
1560 {
1561 snz++ ;
1562 }
1563 }
1564
1565 /* allocate S */
1566 S = mxCreateSparseLogicalMatrix (nrow, ncol, snz) ;
1567 Sp = (Long *) mxGetJc (S) ;
1568 Si = (Long *) mxGetIr (S) ;
1569 Sl = mxGetLogicals (S) ;
1570
1571 /* copy A into S, dropping zero entries */
1572 snz = 0 ;
1573 for (j = 0 ; j < ncol ; j++)
1574 {
1575 Sp [j] = snz ;
1576 pend = Ap [j+1] ;
1577 for (p = Ap [j] ; p < pend ; p++)
1578 {
1579 if (Al [p])
1580 {
1581 Si [snz] = Ai [p] ;
1582 Sl [snz] = 1 ;
1583 snz++ ;
1584 }
1585 }
1586 }
1587
1588 }
1589 else
1590
1591 #endif
1592
1593 if (mxIsComplex (A))
1594 {
1595
1596 /* ------------------------------------------------------------------ */
1597 /* copy a sparse complex double matrix */
1598 /* ------------------------------------------------------------------ */
1599
1600 /* count the number of nonzeros in A */
1601 Ax = mxGetPr (A) ;
1602 Az = mxGetPi (A) ;
1603 for (p = 0 ; p < anz ; p++)
1604 {
1605 aij = Ax [p] ;
1606 zij = Az [p] ;
1607 if (CHOLMOD_IS_NONZERO (aij) || CHOLMOD_IS_NONZERO (zij))
1608 {
1609 snz++ ;
1610 }
1611 }
1612
1613 /* allocate S */
1614 S = mxCreateSparse (nrow, ncol, snz, mxCOMPLEX) ;
1615 Sp = (Long *) mxGetJc (S) ;
1616 Si = (Long *) mxGetIr (S) ;
1617 Sx = mxGetPr (S) ;
1618 Sz = mxGetPi (S) ;
1619
1620 /* copy A into S, dropping zero entries */
1621 snz = 0 ;
1622 for (j = 0 ; j < ncol ; j++)
1623 {
1624 Sp [j] = snz ;
1625 pend = Ap [j+1] ;
1626 for (p = Ap [j] ; p < pend ; p++)
1627 {
1628 aij = Ax [p] ;
1629 zij = Az [p] ;
1630 if (CHOLMOD_IS_NONZERO (aij) || CHOLMOD_IS_NONZERO (zij))
1631 {
1632 Si [snz] = Ai [p] ;
1633 Sx [snz] = aij ;
1634 Sz [snz] = zij ;
1635 snz++ ;
1636 }
1637 }
1638 }
1639
1640 }
1641 else
1642 {
1643
1644 /* ------------------------------------------------------------------ */
1645 /* copy a sparse real double matrix */
1646 /* ------------------------------------------------------------------ */
1647
1648 /* count the number of nonzeros in A */
1649 Ax = mxGetPr (A) ;
1650 for (p = 0 ; p < anz ; p++)
1651 {
1652 aij = Ax [p] ;
1653 if (CHOLMOD_IS_NONZERO (aij))
1654 {
1655 snz++ ;
1656 }
1657 }
1658
1659 /* allocate S */
1660 S = mxCreateSparse (nrow, ncol, snz, mxREAL) ;
1661 Sp = (Long *) mxGetJc (S) ;
1662 Si = (Long *) mxGetIr (S) ;
1663 Sx = mxGetPr (S) ;
1664
1665 /* copy A into S, dropping zero entries */
1666 snz = 0 ;
1667 for (j = 0 ; j < ncol ; j++)
1668 {
1669 Sp [j] = snz ;
1670 pend = Ap [j+1] ;
1671 for (p = Ap [j] ; p < pend ; p++)
1672 {
1673 aij = Ax [p] ;
1674 if (CHOLMOD_IS_NONZERO (aij))
1675 {
1676 Si [snz] = Ai [p] ;
1677 Sx [snz] = aij ;
1678 snz++ ;
1679 }
1680 }
1681 }
1682 }
1683
1684 Sp [ncol] = snz ;
1685 return (S) ;
1686 }
1687
1688
1689 /* ========================================================================== */
1690 /* === sputil_sparse_to_dense =============================================== */
1691 /* ========================================================================== */
1692
1693 /* convert a sparse double or logical array to a dense double array */
1694
sputil_sparse_to_dense(const mxArray * S)1695 mxArray *sputil_sparse_to_dense (const mxArray *S)
1696 {
1697 mxArray *X ;
1698 double *Sx, *Sz, *Xx, *Xz ;
1699 Long *Sp, *Si ;
1700 Long nrow, ncol, i, j, p, pend, j2 ;
1701
1702 #ifndef MATLAB6p1_OR_EARLIER
1703
1704 /* MATLAB 6.1 or earlier : all sparse matrices are OK */
1705 if (! (mxGetClassID (S) == mxLOGICAL_CLASS ||
1706 mxGetClassID (S) == mxDOUBLE_CLASS))
1707 {
1708 /* only sparse logical and real/complex double matrices supported */
1709 sputil_error (ERROR_INVALID_TYPE, 0) ;
1710 }
1711
1712 #endif
1713
1714 nrow = mxGetM (S) ;
1715 ncol = mxGetN (S) ;
1716 Sp = (Long *) mxGetJc (S) ;
1717 Si = (Long *) mxGetIr (S) ;
1718
1719 #ifndef MATLAB6p1_OR_EARLIER
1720
1721 /* MATLAB 6.1 or earlier do not have mxLOGICAL_CLASS */
1722 if (mxIsLogical (S))
1723 {
1724 /* logical */
1725 mxLogical *Sl ;
1726 Sl = (mxLogical *) mxGetData (S) ;
1727 X = mxCreateDoubleMatrix (nrow, ncol, mxREAL) ;
1728 Xx = mxGetPr (X) ;
1729 for (j = 0 ; j < ncol ; j++)
1730 {
1731 pend = Sp [j+1] ;
1732 j2 = j*nrow ;
1733 for (p = Sp [j] ; p < pend ; p++)
1734 {
1735 Xx [Si [p] + j2] = (double) (Sl [p]) ;
1736 }
1737 }
1738 }
1739 else
1740
1741 #endif
1742
1743 if (mxIsComplex (S))
1744 {
1745 /* complex */
1746 Sx = mxGetPr (S) ;
1747 Sz = mxGetPi (S) ;
1748 X = mxCreateDoubleMatrix (nrow, ncol, mxCOMPLEX) ;
1749 Xx = mxGetPr (X) ;
1750 Xz = mxGetPi (X) ;
1751 for (j = 0 ; j < ncol ; j++)
1752 {
1753 pend = Sp [j+1] ;
1754 j2 = j*nrow ;
1755 for (p = Sp [j] ; p < pend ; p++)
1756 {
1757 i = Si [p] ;
1758 Xx [i + j2] = Sx [p] ;
1759 Xz [i + j2] = Sz [p] ;
1760 }
1761 }
1762 }
1763 else
1764 {
1765 /* real */
1766 Sx = mxGetPr (S) ;
1767 X = mxCreateDoubleMatrix (nrow, ncol, mxREAL) ;
1768 Xx = mxGetPr (X) ;
1769 for (j = 0 ; j < ncol ; j++)
1770 {
1771 pend = Sp [j+1] ;
1772 j2 = j*nrow ;
1773 for (p = Sp [j] ; p < pend ; p++)
1774 {
1775 Xx [Si [p] + j2] = Sx [p] ;
1776 }
1777 }
1778 }
1779
1780 return (X) ;
1781 }
1782
1783
1784 /* ========================================================================== */
1785 /* === sputil_check_ijvector ================================================ */
1786 /* ========================================================================== */
1787
1788 /* Check a sparse i or j input argument */
1789
sputil_check_ijvector(const mxArray * arg)1790 void sputil_check_ijvector (const mxArray *arg)
1791 {
1792 if (mxIsComplex (arg))
1793 {
1794 /* i and j cannot be complex */
1795 sputil_error (ERROR_NOT_INTEGER, TRUE) ;
1796 }
1797 if (mxIsSparse (arg))
1798 {
1799 /* the i and j arguments for sparse(i,j,s,...) can be sparse, but if so
1800 * they must have no zero entries. */
1801 double mn, m, nz ;
1802 Long *p, n ;
1803 m = (double) mxGetM (arg) ;
1804 n = mxGetN (arg) ;
1805 mn = m*n ;
1806 p = (Long *) mxGetJc (arg) ;
1807 nz = p [n] ;
1808 if (mn != nz)
1809 {
1810 /* i or j contains at least one zero, which is invalid */
1811 sputil_error (ERROR_TOO_SMALL, TRUE) ;
1812 }
1813 }
1814 }
1815
1816
1817 /* ========================================================================== */
1818 /* === sputil_sparse ======================================================== */
1819 /* ========================================================================== */
1820
1821 /* Implements the sparse2 mexFunction */
1822
sputil_sparse(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])1823 void sputil_sparse
1824 (
1825 int nargout,
1826 mxArray *pargout [ ],
1827 int nargin,
1828 const mxArray *pargin [ ]
1829 )
1830 {
1831 double x, z ;
1832 double *z_vector ;
1833 void *i_vector, *j_vector, *x_vector ;
1834 mxArray *s_array ;
1835 cholmod_sparse *S, *Z ;
1836 cholmod_common Common, *cm ;
1837 Long nrow, ncol, k, nz, i_is_scalar, j_is_scalar, s_is_sparse,
1838 s_is_scalar, ilen, jlen, slen, nzmax, i, j, s_complex, ndropped ;
1839 mxClassID i_class, j_class, s_class ;
1840
1841 /* ---------------------------------------------------------------------- */
1842 /* start CHOLMOD and set defaults */
1843 /* ---------------------------------------------------------------------- */
1844
1845 cm = &Common ;
1846 cholmod_l_start (cm) ;
1847 sputil_config (SPUMONI, cm) ;
1848
1849 /* ---------------------------------------------------------------------- */
1850 /* get inputs */
1851 /* ---------------------------------------------------------------------- */
1852
1853 if (nargout > 2 || nargin > 6 || nargin == 4 || nargin == 0)
1854 {
1855 sputil_error (ERROR_USAGE, FALSE) ;
1856 }
1857
1858 /* ---------------------------------------------------------------------- */
1859 /* convert inputs into a sparse matrix S */
1860 /* ---------------------------------------------------------------------- */
1861
1862 S = NULL ;
1863 Z = NULL ;
1864
1865 if (nargin == 1)
1866 {
1867
1868 /* ------------------------------------------------------------------ */
1869 /* S = sparse (A) where A is sparse or full */
1870 /* ------------------------------------------------------------------ */
1871
1872 nrow = mxGetM (pargin [0]) ;
1873 ncol = mxGetN (pargin [0]) ;
1874
1875 if (mxIsSparse (pargin [0]))
1876 {
1877
1878 /* -------------------------------------------------------------- */
1879 /* S = sparse (A) where A is sparse (double, complex, or logical) */
1880 /* -------------------------------------------------------------- */
1881
1882 pargout [0] = sputil_copy_sparse (pargin [0]) ;
1883
1884 }
1885 else
1886 {
1887
1888 /* -------------------------------------------------------------- */
1889 /* S = sparse (A) where A is full (real or complex) */
1890 /* -------------------------------------------------------------- */
1891
1892 /* A can be of any numeric type (mxLogical, int8, ..., double),
1893 * except that if A is complex, it must also be double. */
1894 pargout [0] = sputil_dense_to_sparse (pargin [0]) ;
1895 }
1896
1897 }
1898 else if (nargin == 2)
1899 {
1900
1901 /* ------------------------------------------------------------------ */
1902 /* S = sparse (m,n) */
1903 /* ------------------------------------------------------------------ */
1904
1905 Long *Sp ;
1906 nrow = sputil_get_integer (pargin [0], FALSE, 0) ;
1907 ncol = sputil_get_integer (pargin [1], FALSE, 0) ;
1908 pargout [0] = mxCreateSparse (nrow, ncol, 1, mxREAL) ;
1909 Sp = (Long *) mxGetJc (pargout [0]) ;
1910 Sp [0] = 0 ;
1911
1912 }
1913 else
1914 {
1915
1916 /* ------------------------------------------------------------------ */
1917 /* S = sparse (i,j,s), sparse (i,j,s,m,n) or sparse (i,j,s,m,n,nzmax) */
1918 /* ------------------------------------------------------------------ */
1919
1920 /* i, j, and s can be of any numeric type */
1921
1922 /* ensure i and j are valid (i and j cannot be complex) */
1923 sputil_check_ijvector (pargin [0]) ;
1924 sputil_check_ijvector (pargin [1]) ;
1925
1926 /* convert s from sparse to dense (double), if necessary */
1927 s_is_sparse = mxIsSparse (pargin [2]) ;
1928 if (s_is_sparse)
1929 {
1930 /* s must be double (real/complex) or logical */
1931 s_array = sputil_sparse_to_dense (pargin [2]) ;
1932 }
1933 else
1934 {
1935 s_array = (mxArray *) pargin [2] ;
1936 }
1937
1938 /* s is now full. It can be any class, except if complex it must also
1939 * be double */
1940 s_class = mxGetClassID (s_array) ;
1941 s_complex = mxIsComplex (s_array) ;
1942 if (s_complex && s_class != mxDOUBLE_CLASS)
1943 {
1944 /* for complex case, only double class is supported */
1945 sputil_error (ERROR_INVALID_TYPE, 0) ;
1946 }
1947
1948 /* get sizes of inputs */
1949 ilen = sputil_nelements (pargin [0]) ;
1950 jlen = sputil_nelements (pargin [1]) ;
1951 slen = sputil_nelements (s_array) ;
1952
1953 /* if i, j, s are scalars, they "float" to sizes of non-scalar args */
1954 i_is_scalar = (ilen == 1) ;
1955 j_is_scalar = (jlen == 1) ;
1956 s_is_scalar = (slen == 1) ;
1957
1958 /* find the length */
1959 if (!i_is_scalar)
1960 {
1961 /* if i is not a scalar, let it determine the length */
1962 nz = ilen ;
1963 }
1964 else if (!j_is_scalar)
1965 {
1966 /* otherwise, if j is not a scalar, let it determine the length */
1967 nz = jlen ;
1968 }
1969 else
1970 {
1971 /* finally, i and j are both scalars, so let s determine length */
1972 nz = slen ;
1973 }
1974
1975 /* make sure the sizes are compatible */
1976 if (!((i_is_scalar || ilen == nz) &&
1977 (j_is_scalar || jlen == nz) &&
1978 (s_is_scalar || slen == nz)))
1979 {
1980 sputil_error (ERROR_LENGTH, FALSE) ;
1981 }
1982
1983 if (nargin > 4)
1984 {
1985 nrow = sputil_get_integer (pargin [3], FALSE, 0) ;
1986 ncol = sputil_get_integer (pargin [4], FALSE, 0) ;
1987 }
1988 else
1989 {
1990 /* nrow and ncol will be discovered by scanning i and j */
1991 nrow = EMPTY ;
1992 ncol = EMPTY ;
1993 }
1994
1995 if (nargin > 5)
1996 {
1997 nzmax = sputil_get_integer (pargin [5], FALSE, 0) ;
1998 nzmax = MAX (nzmax, nz) ;
1999 }
2000 else
2001 {
2002 nzmax = nz ;
2003 }
2004
2005 /* ------------------------------------------------------------------ */
2006 /* convert triplet form to sparse form */
2007 /* ------------------------------------------------------------------ */
2008
2009 i = i_is_scalar ? sputil_get_integer (pargin [0], TRUE, nrow) : 0 ;
2010 i_vector = mxGetData (pargin [0]) ;
2011 i_class = mxGetClassID (pargin [0]) ;
2012
2013 j = j_is_scalar ? sputil_get_integer (pargin [1], TRUE, ncol) : 0 ;
2014 j_vector = mxGetData (pargin [1]) ;
2015 j_class = mxGetClassID (pargin [1]) ;
2016
2017 x_vector = mxGetData (s_array) ;
2018 z_vector = mxGetPi (s_array) ;
2019 x = sputil_get_double (s_array) ;
2020 z = (s_complex && z_vector != NULL) ? (z_vector [0]) : 0 ;
2021
2022 S = sputil_triplet_to_sparse (nrow, ncol, nz, nzmax,
2023 i_is_scalar, i, i_vector, i_class,
2024 j_is_scalar, j, j_vector, j_class,
2025 s_is_scalar, x, z, x_vector, z_vector,
2026 s_class, s_complex,
2027 cm) ;
2028
2029 /* set nzmax(S) to nnz(S), unless nzmax is specified on input */
2030 if (nargin <= 5 && S != NULL)
2031 {
2032 cholmod_l_reallocate_sparse (cholmod_l_nnz (S, cm), S, cm) ;
2033 }
2034
2035 if (nargout > 1)
2036 {
2037 /* return a binary pattern of the explicit zero entries, for the
2038 * [S Z] = sparse(i,j,x, ...) form. */
2039 Z = sputil_extract_zeros (S, cm) ;
2040 }
2041
2042 /* drop explicit zeros from S */
2043 ndropped = sputil_drop_zeros (S) ;
2044
2045 /* if entries dropped, set nzmax(S) to nnz(S), unless nzmax specified */
2046 if (ndropped > 0 && nargin <= 5 && S != NULL)
2047 {
2048 cholmod_l_reallocate_sparse (cholmod_l_nnz (S, cm), S, cm) ;
2049 }
2050
2051 if (s_is_sparse)
2052 {
2053 mxDestroyArray (s_array) ;
2054 }
2055 }
2056
2057 /* ---------------------------------------------------------------------- */
2058 /* convert S into a MATLAB sparse matrix */
2059 /* ---------------------------------------------------------------------- */
2060
2061 k = 0 ;
2062 if (S != NULL)
2063 {
2064
2065 #ifndef MATLAB6p1_OR_EARLIER
2066
2067 /* MATLAB 6.1 or earlier do not have mxLOGICAL_CLASS */
2068 if (mxIsLogical (pargin [2]))
2069 {
2070 /* copy S into a MATLAB sparse logical matrix */
2071 mxLogical *s_logical ;
2072 pargout [0] = mxCreateSparseLogicalMatrix (0, 0, 0) ;
2073 s_logical = cholmod_l_malloc (S->nzmax, sizeof (mxLogical), cm) ;
2074 for (k = 0 ; k < (Long) (S->nzmax) ; k++)
2075 {
2076 s_logical [k] = 1 ;
2077 }
2078 MXFREE (mxGetData (pargout [0])) ;
2079 mxSetData (pargout [0], s_logical) ;
2080 mexMakeMemoryPersistent (s_logical) ;
2081 k++ ;
2082 }
2083 else
2084
2085 #endif
2086 if (S->x == NULL) mexErrMsgTxt ("invalid sparse matrix") ;
2087
2088 if (mxIsComplex (pargin [2]))
2089 {
2090 /* copy S into a MATLAB sparse complex double matrix */
2091 if (S->z == NULL) mexErrMsgTxt ("invalid complex sparse matrix") ;
2092 pargout [0] = mxCreateSparse (0, 0, 0, mxCOMPLEX) ;
2093 MXFREE (mxGetPr (pargout [0])) ;
2094 MXFREE (mxGetPi (pargout [0])) ;
2095 mxSetPr (pargout [0], S->x) ;
2096 mxSetPi (pargout [0], S->z) ;
2097 mexMakeMemoryPersistent (S->x) ;
2098 mexMakeMemoryPersistent (S->z) ;
2099 k += 2 ;
2100 S->x = NULL ;
2101 S->z = NULL ;
2102 }
2103 else
2104 {
2105 /* copy S into a MATLAB sparse real double matrix */
2106 pargout [0] = mxCreateSparse (0, 0, 0, mxREAL) ;
2107 mxSetPr (pargout [0], S->x) ;
2108 mexMakeMemoryPersistent (S->x) ;
2109 k++ ;
2110 S->x = NULL ;
2111 }
2112
2113 mxSetM (pargout [0], S->nrow) ;
2114 mxSetN (pargout [0], S->ncol) ;
2115 mxSetNzmax (pargout [0], S->nzmax) ;
2116 MXFREE (mxGetJc (pargout [0])) ;
2117 MXFREE (mxGetIr (pargout [0])) ;
2118 mxSetJc (pargout [0], S->p) ;
2119 mxSetIr (pargout [0], S->i) ;
2120 mexMakeMemoryPersistent (S->p) ;
2121 mexMakeMemoryPersistent (S->i) ;
2122 k += 2 ;
2123
2124 /* free cholmod_sparse S, except for what has been given to MATLAB */
2125 S->p = NULL ;
2126 S->i = NULL ;
2127 cholmod_l_free_sparse (&S, cm) ;
2128 }
2129
2130 /* ---------------------------------------------------------------------- */
2131 /* return Z to MATLAB, if requested */
2132 /* ---------------------------------------------------------------------- */
2133
2134 if (nargout > 1)
2135 {
2136 if (Z == NULL)
2137 {
2138 /* Z not computed; return an empty matrix */
2139 Z = cholmod_l_spzeros (nrow, ncol, 0, CHOLMOD_REAL, cm) ;
2140 }
2141 pargout [1] = sputil_put_sparse (&Z, cm) ;
2142 }
2143
2144 cholmod_l_finish (cm) ;
2145 cholmod_l_print_common (" ", cm) ;
2146 /*
2147 if (cm->malloc_count != k) mexErrMsgTxt ("!") ;
2148 */
2149 }
2150