Mercurial > octave-nkf

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