1 /* ========================================================================== */
2 /* === CHOLMOD/MATLAB/ldlsolve mexFunction ================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/MATLAB Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * http://www.suitesparse.com
8 * MATLAB(tm) is a Trademark of The MathWorks, Inc.
9 * -------------------------------------------------------------------------- */
10
11 /* Solve LDL'x=b given an LDL' factorization computed by ldlchol.
12 *
13 * Usage:
14 *
15 * x = ldlsolve (LD,b)
16 *
17 * b can be dense or sparse.
18 */
19
20 #include "cholmod_matlab.h"
21
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])22 void mexFunction
23 (
24 int nargout,
25 mxArray *pargout [ ],
26 int nargin,
27 const mxArray *pargin [ ]
28 )
29 {
30 double dummy = 0, rcond ;
31 Long *Lp, *Lnz, *Lprev, *Lnext ;
32 cholmod_sparse *Bs, Bspmatrix, *Xs ;
33 cholmod_dense *B, Bmatrix, *X ;
34 cholmod_factor *L ;
35 cholmod_common Common, *cm ;
36 Long j, k, n, B_is_sparse, head, tail ;
37
38 /* ---------------------------------------------------------------------- */
39 /* start CHOLMOD and set parameters */
40 /* ---------------------------------------------------------------------- */
41
42 cm = &Common ;
43 cholmod_l_start (cm) ;
44 sputil_config (SPUMONI, cm) ;
45
46 /* ---------------------------------------------------------------------- */
47 /* check inputs */
48 /* ---------------------------------------------------------------------- */
49
50 if (nargout > 1 || nargin != 2)
51 {
52 mexErrMsgTxt ("Usage: x = ldlsolve (LD, b)") ;
53 }
54
55 n = mxGetN (pargin [0]) ;
56 k = mxGetN (pargin [1]) ;
57
58 if (!mxIsSparse (pargin [0]) || n != mxGetM (pargin [0]))
59 {
60 mexErrMsgTxt ("ldlsolve: LD must be sparse and square") ;
61 }
62 if (n != mxGetM (pargin [1]))
63 {
64 mexErrMsgTxt ("ldlsolve: b wrong dimension") ;
65 }
66
67 /* ---------------------------------------------------------------------- */
68 /* get b */
69 /* ---------------------------------------------------------------------- */
70
71 /* get sparse or dense matrix B */
72 B = NULL ;
73 Bs = NULL ;
74 B_is_sparse = mxIsSparse (pargin [1]) ;
75 if (B_is_sparse)
76 {
77 /* get sparse matrix B (unsymmetric) */
78 Bs = sputil_get_sparse (pargin [1], &Bspmatrix, &dummy, 0) ;
79 }
80 else
81 {
82 /* get dense matrix B */
83 B = sputil_get_dense (pargin [1], &Bmatrix, &dummy) ;
84 }
85
86 /* ---------------------------------------------------------------------- */
87 /* construct a shallow copy of the input sparse matrix L */
88 /* ---------------------------------------------------------------------- */
89
90 /* the construction of the CHOLMOD takes O(n) time and memory */
91
92 /* allocate the CHOLMOD symbolic L */
93 L = cholmod_l_allocate_factor (n, cm) ;
94 L->ordering = CHOLMOD_NATURAL ;
95
96 /* get the MATLAB L */
97 L->p = mxGetJc (pargin [0]) ;
98 L->i = mxGetIr (pargin [0]) ;
99 L->x = mxGetPr (pargin [0]) ;
100 L->z = mxGetPi (pargin [0]) ;
101
102 /* allocate and initialize the rest of L */
103 L->nz = cholmod_l_malloc (n, sizeof (Long), cm) ;
104 Lp = L->p ;
105 Lnz = L->nz ;
106 for (j = 0 ; j < n ; j++)
107 {
108 Lnz [j] = Lp [j+1] - Lp [j] ;
109 }
110 L->prev = cholmod_l_malloc (n+2, sizeof (Long), cm) ;
111 L->next = cholmod_l_malloc (n+2, sizeof (Long), cm) ;
112 Lprev = L->prev ;
113 Lnext = L->next ;
114
115 head = n+1 ;
116 tail = n ;
117 Lnext [head] = 0 ;
118 Lprev [head] = -1 ;
119 Lnext [tail] = -1 ;
120 Lprev [tail] = n-1 ;
121 for (j = 0 ; j < n ; j++)
122 {
123 Lnext [j] = j+1 ;
124 Lprev [j] = j-1 ;
125 }
126 Lprev [0] = head ;
127
128 L->xtype = (mxIsComplex (pargin [0])) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
129 L->nzmax = Lp [n] ;
130
131 /* ---------------------------------------------------------------------- */
132 /* solve and return solution to MATLAB */
133 /* ---------------------------------------------------------------------- */
134
135 if (B_is_sparse)
136 {
137 /* solve LDL'X=B with sparse X and B; return sparse X to MATLAB */
138 Xs = cholmod_l_spsolve (CHOLMOD_LDLt, L, Bs, cm) ;
139 pargout [0] = sputil_put_sparse (&Xs, cm) ;
140 }
141 else
142 {
143 /* solve AX=B with dense X and B; return dense X to MATLAB */
144 X = cholmod_l_solve (CHOLMOD_LDLt, L, B, cm) ;
145 pargout [0] = sputil_put_dense (&X, cm) ;
146 }
147
148 rcond = cholmod_l_rcond (L, cm) ;
149 if (rcond == 0)
150 {
151 mexWarnMsgTxt ("Matrix is indefinite or singular to working precision");
152 }
153 else if (rcond < DBL_EPSILON)
154 {
155 mexWarnMsgTxt ("Matrix is close to singular or badly scaled.") ;
156 mexPrintf (" Results may be inaccurate. RCOND = %g.\n", rcond) ;
157 }
158
159 /* ---------------------------------------------------------------------- */
160 /* free workspace and the CHOLMOD L, except for what is copied to MATLAB */
161 /* ---------------------------------------------------------------------- */
162
163 L->p = NULL ;
164 L->i = NULL ;
165 L->x = NULL ;
166 L->z = NULL ;
167 cholmod_l_free_factor (&L, cm) ;
168 cholmod_l_finish (cm) ;
169 cholmod_l_print_common (" ", cm) ;
170 /*
171 if (cm->malloc_count !=
172 (mxIsComplex (pargout [0]) + (mxIsSparse (pargout[0]) ? 3:1)))
173 mexErrMsgTxt ("!") ;
174 */
175 }
176