1 /* ========================================================================== */
2 /* === UMFPACK_get_determinant ============================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* UMFPACK Copyright (c) Timothy A. Davis, CISE, */
7 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */
8 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */
9 /* UMFPACK_get_determinant contributed by David Bateman, Motorola, Paris. */
10 /* -------------------------------------------------------------------------- */
11
12 /*
13 User-callable. From the LU factors, scale factor, and permutation vectors
14 held in the Numeric object, calculates the determinant of the matrix A.
15 See umfpack_get_determinant.h for a more detailed description.
16
17 Dynamic memory usage: calls UMF_malloc once, for a total space of
18 n integers, and then frees all of it via UMF_free when done.
19
20 Contributed by David Bateman, Motorola, Nov. 2004.
21 Modified for V4.4, Jan. 2005.
22 */
23
24 #include "umf_internal.h"
25 #include "umf_valid_numeric.h"
26 #include "umf_malloc.h"
27 #include "umf_free.h"
28
29 /* ========================================================================== */
30 /* === rescale_determinant ================================================== */
31 /* ========================================================================== */
32
33 /* If the mantissa is too big or too small, rescale it and change exponent */
34
rescale_determinant(Entry * d_mantissa,double * d_exponent)35 PRIVATE Int rescale_determinant
36 (
37 Entry *d_mantissa,
38 double *d_exponent
39 )
40 {
41 double d_abs ;
42
43 ABS (d_abs, *d_mantissa) ;
44
45 if (SCALAR_IS_ZERO (d_abs))
46 {
47 /* the determinant is zero */
48 *d_exponent = 0 ;
49 return (FALSE) ;
50 }
51
52 if (SCALAR_IS_NAN (d_abs))
53 {
54 /* the determinant is NaN */
55 return (FALSE) ;
56 }
57
58 while (d_abs < 1.)
59 {
60 SCALE (*d_mantissa, 10.0) ;
61 *d_exponent = *d_exponent - 1.0 ;
62 ABS (d_abs, *d_mantissa) ;
63 }
64
65 while (d_abs >= 10.)
66 {
67 SCALE (*d_mantissa, 0.1) ;
68 *d_exponent = *d_exponent + 1.0 ;
69 ABS (d_abs, *d_mantissa) ;
70 }
71
72 return (TRUE) ;
73 }
74
75 /* ========================================================================== */
76 /* === UMFPACK_get_determinant ============================================== */
77 /* ========================================================================== */
78
UMFPACK_get_determinant(double * Mx,double * Mz,double * Ex,void * NumericHandle,double User_Info[UMFPACK_INFO])79 GLOBAL Int UMFPACK_get_determinant
80 (
81 double *Mx,
82 #ifdef COMPLEX
83 double *Mz,
84 #endif
85 double *Ex,
86 void *NumericHandle,
87 double User_Info [UMFPACK_INFO]
88 )
89 {
90 /* ---------------------------------------------------------------------- */
91 /* local variables */
92 /* ---------------------------------------------------------------------- */
93
94 Entry d_mantissa, d_tmp ;
95 double d_exponent, Info2 [UMFPACK_INFO], one [2] = {1.0, 0.0}, d_sign ;
96 Entry *D ;
97 double *Info, *Rs ;
98 NumericType *Numeric ;
99 Int i, n, itmp, npiv, *Wi, *Rperm, *Cperm, do_scale ;
100
101 #ifndef NRECIPROCAL
102 Int do_recip ;
103 #endif
104
105 /* ---------------------------------------------------------------------- */
106 /* check input parameters */
107 /* ---------------------------------------------------------------------- */
108
109 if (User_Info != (double *) NULL)
110 {
111 /* return Info in user's array */
112 Info = User_Info ;
113 }
114 else
115 {
116 /* no Info array passed - use local one instead */
117 Info = Info2 ;
118 for (i = 0 ; i < UMFPACK_INFO ; i++)
119 {
120 Info [i] = EMPTY ;
121 }
122 }
123
124 Info [UMFPACK_STATUS] = UMFPACK_OK ;
125
126 Numeric = (NumericType *) NumericHandle ;
127 if (!UMF_valid_numeric (Numeric))
128 {
129 Info [UMFPACK_STATUS] = UMFPACK_ERROR_invalid_Numeric_object ;
130 return (UMFPACK_ERROR_invalid_Numeric_object) ;
131 }
132
133 if (Numeric->n_row != Numeric->n_col)
134 {
135 /* only square systems can be handled */
136 Info [UMFPACK_STATUS] = UMFPACK_ERROR_invalid_system ;
137 return (UMFPACK_ERROR_invalid_system) ;
138 }
139
140 if (Mx == (double *) NULL)
141 {
142 Info [UMFPACK_STATUS] = UMFPACK_ERROR_argument_missing ;
143 return (UMFPACK_ERROR_argument_missing) ;
144 }
145
146 n = Numeric->n_row ;
147
148 /* ---------------------------------------------------------------------- */
149 /* allocate workspace */
150 /* ---------------------------------------------------------------------- */
151
152 Wi = (Int *) UMF_malloc (n, sizeof (Int)) ;
153
154 if (!Wi)
155 {
156 DEBUGm4 (("out of memory: get determinant\n")) ;
157 Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
158 return (UMFPACK_ERROR_out_of_memory) ;
159 }
160
161 /* ---------------------------------------------------------------------- */
162 /* compute the determinant */
163 /* ---------------------------------------------------------------------- */
164
165 Rs = Numeric->Rs ; /* row scale factors */
166 do_scale = (Rs != (double *) NULL) ;
167
168 #ifndef NRECIPROCAL
169 do_recip = Numeric->do_recip ;
170 #endif
171
172 d_mantissa = ((Entry *) one) [0] ;
173 d_exponent = 0.0 ;
174 D = Numeric->D ;
175
176 /* compute product of diagonal entries of U */
177 for (i = 0 ; i < n ; i++)
178 {
179 MULT (d_tmp, d_mantissa, D [i]) ;
180 d_mantissa = d_tmp ;
181
182 if (!rescale_determinant (&d_mantissa, &d_exponent))
183 {
184 /* the determinant is zero or NaN */
185 Info [UMFPACK_STATUS] = UMFPACK_WARNING_singular_matrix ;
186 /* no need to compute the determinant of R */
187 do_scale = FALSE ;
188 break ;
189 }
190 }
191
192 /* compute product of diagonal entries of R (or its inverse) */
193 if (do_scale)
194 {
195 for (i = 0 ; i < n ; i++)
196 {
197 #ifndef NRECIPROCAL
198 if (do_recip)
199 {
200 /* compute determinant of R inverse */
201 SCALE_DIV (d_mantissa, Rs [i]) ;
202 }
203 else
204 #endif
205 {
206 /* compute determinant of R */
207 SCALE (d_mantissa, Rs [i]) ;
208 }
209 if (!rescale_determinant (&d_mantissa, &d_exponent))
210 {
211 /* the determinant is zero or NaN. This is very unlikey to
212 * occur here, since the scale factors for a tiny or zero row
213 * are set to 1. */
214 Info [UMFPACK_STATUS] = UMFPACK_WARNING_singular_matrix ;
215 break ;
216 }
217 }
218 }
219
220 /* ---------------------------------------------------------------------- */
221 /* determine if P and Q are odd or even permutations */
222 /* ---------------------------------------------------------------------- */
223
224 npiv = 0 ;
225 Rperm = Numeric->Rperm ;
226
227 for (i = 0 ; i < n ; i++)
228 {
229 Wi [i] = Rperm [i] ;
230 }
231
232 for (i = 0 ; i < n ; i++)
233 {
234 while (Wi [i] != i)
235 {
236 itmp = Wi [Wi [i]] ;
237 Wi [Wi [i]] = Wi [i] ;
238 Wi [i] = itmp ;
239 npiv++ ;
240 }
241 }
242
243 Cperm = Numeric->Cperm ;
244
245 for (i = 0 ; i < n ; i++)
246 {
247 Wi [i] = Cperm [i] ;
248 }
249
250 for (i = 0 ; i < n ; i++)
251 {
252 while (Wi [i] != i)
253 {
254 itmp = Wi [Wi [i]] ;
255 Wi [Wi [i]] = Wi [i] ;
256 Wi [i] = itmp ;
257 npiv++ ;
258 }
259 }
260
261 /* if npiv is odd, the sign is -1. if it is even, the sign is +1 */
262 d_sign = (npiv % 2) ? -1. : 1. ;
263
264 /* ---------------------------------------------------------------------- */
265 /* free workspace */
266 /* ---------------------------------------------------------------------- */
267
268 (void) UMF_free ((void *) Wi) ;
269
270 /* ---------------------------------------------------------------------- */
271 /* compute the magnitude and exponent of the determinant */
272 /* ---------------------------------------------------------------------- */
273
274 if (Ex == (double *) NULL)
275 {
276 /* Ex is not provided, so return the entire determinant in d_mantissa */
277 SCALE (d_mantissa, pow (10.0, d_exponent)) ;
278 }
279 else
280 {
281 Ex [0] = d_exponent ;
282 }
283
284 Mx [0] = d_sign * REAL_COMPONENT (d_mantissa) ;
285
286 #ifdef COMPLEX
287 if (SPLIT (Mz))
288 {
289 Mz [0] = d_sign * IMAG_COMPONENT (d_mantissa) ;
290 }
291 else
292 {
293 Mx [1] = d_sign * IMAG_COMPONENT (d_mantissa) ;
294 }
295 #endif
296
297 /* determine if the determinant has (or will) overflow or underflow */
298 if (d_exponent + 1.0 > log10 (DBL_MAX))
299 {
300 Info [UMFPACK_STATUS] = UMFPACK_WARNING_determinant_overflow ;
301 }
302 else if (d_exponent - 1.0 < log10 (DBL_MIN))
303 {
304 Info [UMFPACK_STATUS] = UMFPACK_WARNING_determinant_underflow ;
305 }
306
307 return (UMFPACK_OK) ;
308 }
309