1 #include "sparseinv.h"
2
3 /*
4 Z = sparseinv_mex (L, d, UT, Zpattern)
5
6 Given (L+I)*D*(UT+I)' = A, and the symbolic Cholesky factorization of A+A',
7 compute the sparse inverse subset, Z. UT is stored by column, so U = UT'
8 is implicitly stored by row, and is implicitly unit-diagonal. The diagonal
9 is not present. L is stored by column, and is also unit-diagonal. The
10 diagonal is not present in L, either. d is a full vector of size n.
11
12 This mexFunction is only meant to be called from the sparsinv m-file.
13 An optional 2nd output argument returns the flop count.
14
15 Copyright 2011, Timothy A. Davis, http://www.suitesparse.com
16 */
17
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])18 void mexFunction
19 (
20 int nargout,
21 mxArray *pargout [ ],
22 int nargin,
23 const mxArray *pargin [ ]
24 )
25 {
26 Int *Zp, *Zi, *Lp, *Li, *Up, *Uj, *Zpatp, *Zpati, n, *Zdiagp, *Lmunch,
27 znz, j, p, flops ;
28 double *Zx, *Lx, *Ux, *z, *d ;
29
30 /* check inputs */
31 if (nargin != 4 || nargout > 2)
32 {
33 mexErrMsgTxt ("Usage: [Z flops] = sparseinv_mex (L, d, UT, Zpattern)") ;
34 }
35 n = mxGetN (pargin [0]) ;
36 for (j = 0 ; j < 4 ; j++)
37 {
38 if (j == 1) continue ;
39 if (!mxIsSparse (pargin [j]) || mxIsComplex (pargin [j]) ||
40 (mxGetM (pargin [j]) != n) || (mxGetN (pargin [j]) != n))
41 {
42 mexErrMsgTxt ("Matrices must be sparse, real, square, & same size");
43 }
44 }
45 if (mxIsSparse (pargin [1]) || mxIsComplex (pargin [1]) ||
46 (mxGetM (pargin [1]) != n) || (mxGetN (pargin [1]) != 1))
47 {
48 mexErrMsgTxt ("Input d must be a real dense vector of right size") ;
49 }
50
51 /* get inputs */
52 Lp = (Int *) mxGetJc (pargin [0]) ;
53 Li = (Int *) mxGetIr (pargin [0]) ;
54 Lx = mxGetPr (pargin [0]) ;
55
56 d = mxGetPr (pargin [1]) ;
57
58 Up = (Int *) mxGetJc (pargin [2]) ;
59 Uj = (Int *) mxGetIr (pargin [2]) ;
60 Ux = mxGetPr (pargin [2]) ;
61
62 Zpatp = (Int *) mxGetJc (pargin [3]) ;
63 Zpati = (Int *) mxGetIr (pargin [3]) ;
64 znz = Zpatp [n] ;
65
66 /* create output */
67 pargout [0] = mxCreateSparse (n, n, znz, mxREAL) ;
68 Zx = mxGetPr (pargout [0]) ;
69
70 /* get workspace */
71 z = mxCalloc (n, sizeof (double)) ;
72 Zdiagp = mxMalloc (n * sizeof (Int)) ;
73 Lmunch = mxMalloc (n * sizeof (Int)) ;
74
75 /* do the work */
76 flops = sparseinv (n, Lp, Li, Lx, d, Up, Uj, Ux, Zpatp, Zpati, Zx,
77 z, Zdiagp, Lmunch) ;
78
79 /* free workspace */
80 mxFree (z) ;
81 mxFree (Zdiagp) ;
82 mxFree (Lmunch) ;
83
84 /* return results to MATLAB */
85 Zp = (Int *) mxGetJc (pargout [0]) ;
86 Zi = (Int *) mxGetIr (pargout [0]) ;
87 for (j = 0 ; j <= n ; j++)
88 {
89 Zp [j] = Zpatp [j] ;
90 }
91 for (p = 0 ; p < znz ; p++)
92 {
93 Zi [p] = Zpati [p] ;
94 }
95
96 /* drop explicit zeros from the output Z matrix */
97 znz = 0 ;
98 for (j = 0 ; j < n ; j++)
99 {
100 p = Zp [j] ; /* get current location of col j */
101 Zp [j] = znz ; /* record new location of col j */
102 for ( ; p < Zp [j+1] ; p++)
103 {
104 if (Zx [p] != 0)
105 {
106 Zx [znz] = Zx [p] ; /* keep Z(i,j) */
107 Zi [znz++] = Zi [p] ;
108 }
109 }
110 }
111 Zp [n] = znz ; /* finalize Z */
112
113 if (nargout > 1)
114 {
115 pargout [1] = mxCreateDoubleScalar ((double) flops) ;
116 }
117 }
118