1 /* ========================================================================== */
2 /* === Core/cholmod_complex ================================================= */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Core Module. Copyright (C) 2005-2006,
7 * Univ. of Florida. Author: Timothy A. Davis
8 * -------------------------------------------------------------------------- */
9
10 /* If you convert a matrix that contains uninitialized data, valgrind will
11 * complain. This can occur in a factor L which has gaps (a partial
12 * factorization, or after updates that change the nonzero pattern), an
13 * unpacked sparse matrix, a dense matrix with leading dimension d > # of rows,
14 * or any matrix (dense, sparse, triplet, or factor) with more space allocated
15 * than is used. You can safely ignore any of these complaints by valgrind. */
16
17 #include "cholmod_internal.h"
18 #include "cholmod_core.h"
19
20 /* ========================================================================== */
21 /* === cholmod_hypot ======================================================== */
22 /* ========================================================================== */
23
CHOLMOD(hypot)24 double CHOLMOD(hypot) (double x, double y)
25 {
26 return (SuiteSparse_config.hypot_func (x, y)) ;
27 }
28
29
30 /* ========================================================================== */
31 /* === cholmod_divcomplex =================================================== */
32 /* ========================================================================== */
33
34 /* c = a/b where c, a, and b are complex. The real and imaginary parts are
35 * passed as separate arguments to this routine. The NaN case is ignored
36 * for the double relop br >= bi. Returns 1 if the denominator is zero,
37 * 0 otherwise. Note that this return value is the single exception to the
38 * rule that all CHOLMOD routines that return int return TRUE if successful
39 * or FALSE otherise.
40 *
41 * This uses ACM Algo 116, by R. L. Smith, 1962, which tries to avoid
42 * underflow and overflow.
43 *
44 * c can be the same variable as a or b.
45 *
46 * Default value of the SuiteSparse_config.divcomplex_func pointer is
47 * SuiteSparse_divcomplex, located in SuiteSparse_config.c.
48 */
49
CHOLMOD(divcomplex)50 int CHOLMOD(divcomplex)
51 (
52 double ar, double ai, /* real and imaginary parts of a */
53 double br, double bi, /* real and imaginary parts of b */
54 double *cr, double *ci /* real and imaginary parts of c */
55 )
56 {
57 return (SuiteSparse_config.divcomplex_func (ar, ai, br, bi, cr, ci)) ;
58 }
59
60
61 /* ========================================================================== */
62 /* === change_complexity ==================================================== */
63 /* ========================================================================== */
64
65 /* X and Z represent an array of size nz, with numeric xtype given by xtype_in.
66 *
67 * If xtype_in is:
68 * CHOLMOD_PATTERN: X and Z must be NULL.
69 * CHOLMOD_REAL: X is of size nz, Z must be NULL.
70 * CHOLMOD_COMPLEX: X is of size 2*nz, Z must be NULL.
71 * CHOLMOD_ZOMPLEX: X is of size nz, Z is of size nz.
72 *
73 * The array is changed into the numeric xtype given by xtype_out, with the
74 * same definitions of X and Z above. Note that the input conditions, above,
75 * are not checked. These are checked in the caller routine.
76 *
77 * Returns TRUE if successful, FALSE otherwise. X and Z are not modified if
78 * not successful.
79 */
80
change_complexity(Int nz,int xtype_in,int xtype_out,int xtype1,int xtype2,void ** XX,void ** ZZ,cholmod_common * Common)81 static int change_complexity
82 (
83 /* ---- input ---- */
84 Int nz, /* size of X and/or Z */
85 int xtype_in, /* xtype of X and Z on input */
86 int xtype_out, /* requested xtype of X and Z on output */
87 int xtype1, /* xtype_out must be in the range [xtype1 .. xtype2] */
88 int xtype2,
89 /* ---- in/out --- */
90 void **XX, /* old X on input, new X on output */
91 void **ZZ, /* old Z on input, new Z on output */
92 /* --------------- */
93 cholmod_common *Common
94 )
95 {
96 double *Xold, *Zold, *Xnew, *Znew ;
97 Int k ;
98 size_t nz2 ;
99
100 if (xtype_out < xtype1 || xtype_out > xtype2)
101 {
102 ERROR (CHOLMOD_INVALID, "invalid xtype") ;
103 return (FALSE) ;
104 }
105
106 Common->status = CHOLMOD_OK ;
107 Xold = *XX ;
108 Zold = *ZZ ;
109
110 switch (xtype_in)
111 {
112
113 /* ------------------------------------------------------------------ */
114 /* converting from pattern */
115 /* ------------------------------------------------------------------ */
116
117 case CHOLMOD_PATTERN:
118
119 switch (xtype_out)
120 {
121
122 /* ---------------------------------------------------------- */
123 /* pattern -> real */
124 /* ---------------------------------------------------------- */
125
126 case CHOLMOD_REAL:
127 /* allocate X and set to all ones */
128 Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
129 if (Common->status < CHOLMOD_OK)
130 {
131 return (FALSE) ;
132 }
133 for (k = 0 ; k < nz ; k++)
134 {
135 Xnew [k] = 1 ;
136 }
137 *XX = Xnew ;
138 break ;
139
140 /* ---------------------------------------------------------- */
141 /* pattern -> complex */
142 /* ---------------------------------------------------------- */
143
144 case CHOLMOD_COMPLEX:
145 /* allocate X and set to all ones */
146 Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
147 if (Common->status < CHOLMOD_OK)
148 {
149 return (FALSE) ;
150 }
151 for (k = 0 ; k < nz ; k++)
152 {
153 Xnew [2*k ] = 1 ;
154 Xnew [2*k+1] = 0 ;
155 }
156 *XX = Xnew ;
157 break ;
158
159 /* ---------------------------------------------------------- */
160 /* pattern -> zomplex */
161 /* ---------------------------------------------------------- */
162
163 case CHOLMOD_ZOMPLEX:
164 /* allocate X and Z and set to all ones */
165 Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
166 Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
167 if (Common->status < CHOLMOD_OK)
168 {
169 CHOLMOD(free) (nz, sizeof (double), Xnew, Common) ;
170 CHOLMOD(free) (nz, sizeof (double), Znew, Common) ;
171 return (FALSE) ;
172 }
173 for (k = 0 ; k < nz ; k++)
174 {
175 Xnew [k] = 1 ;
176 Znew [k] = 0 ;
177 }
178 *XX = Xnew ;
179 *ZZ = Znew ;
180 break ;
181 }
182 break ;
183
184 /* ------------------------------------------------------------------ */
185 /* converting from real */
186 /* ------------------------------------------------------------------ */
187
188 case CHOLMOD_REAL:
189
190 switch (xtype_out)
191 {
192
193 /* ---------------------------------------------------------- */
194 /* real -> pattern */
195 /* ---------------------------------------------------------- */
196
197 case CHOLMOD_PATTERN:
198 /* free X */
199 *XX = CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
200 break ;
201
202 /* ---------------------------------------------------------- */
203 /* real -> complex */
204 /* ---------------------------------------------------------- */
205
206 case CHOLMOD_COMPLEX:
207 /* allocate a new X and copy the old X */
208 Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
209 if (Common->status < CHOLMOD_OK)
210 {
211 return (FALSE) ;
212 }
213 for (k = 0 ; k < nz ; k++)
214 {
215 Xnew [2*k ] = Xold [k] ;
216 Xnew [2*k+1] = 0 ;
217 }
218 CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
219 *XX = Xnew ;
220 break ;
221
222 /* ---------------------------------------------------------- */
223 /* real -> zomplex */
224 /* ---------------------------------------------------------- */
225
226 case CHOLMOD_ZOMPLEX:
227 /* allocate a new Z and set it to zero */
228 Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
229 if (Common->status < CHOLMOD_OK)
230 {
231 return (FALSE) ;
232 }
233 for (k = 0 ; k < nz ; k++)
234 {
235 Znew [k] = 0 ;
236 }
237 *ZZ = Znew ;
238 break ;
239 }
240 break ;
241
242 /* ------------------------------------------------------------------ */
243 /* converting from complex */
244 /* ------------------------------------------------------------------ */
245
246 case CHOLMOD_COMPLEX:
247
248 switch (xtype_out)
249 {
250
251 /* ---------------------------------------------------------- */
252 /* complex -> pattern */
253 /* ---------------------------------------------------------- */
254
255 case CHOLMOD_PATTERN:
256 /* free X */
257 *XX = CHOLMOD(free) (nz, 2*sizeof (double), *XX, Common) ;
258 break ;
259
260 /* ---------------------------------------------------------- */
261 /* complex -> real */
262 /* ---------------------------------------------------------- */
263
264 case CHOLMOD_REAL:
265 /* pack the real part of X, discarding the imaginary part */
266 for (k = 0 ; k < nz ; k++)
267 {
268 Xold [k] = Xold [2*k] ;
269 }
270 /* shrink X in half (this cannot fail) */
271 nz2 = 2*nz ;
272 *XX = CHOLMOD(realloc) (nz, sizeof (double), *XX, &nz2,
273 Common) ;
274 break ;
275
276 /* ---------------------------------------------------------- */
277 /* complex -> zomplex */
278 /* ---------------------------------------------------------- */
279
280 case CHOLMOD_ZOMPLEX:
281 /* allocate X and Z and copy the old X into them */
282 Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
283 Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
284 if (Common->status < CHOLMOD_OK)
285 {
286 CHOLMOD(free) (nz, sizeof (double), Xnew, Common) ;
287 CHOLMOD(free) (nz, sizeof (double), Znew, Common) ;
288 return (FALSE) ;
289 }
290 for (k = 0 ; k < nz ; k++)
291 {
292 Xnew [k] = Xold [2*k ] ;
293 Znew [k] = Xold [2*k+1] ;
294 }
295 CHOLMOD(free) (nz, 2*sizeof (double), *XX, Common) ;
296 *XX = Xnew ;
297 *ZZ = Znew ;
298 break ;
299 }
300 break ;
301
302 /* ------------------------------------------------------------------ */
303 /* converting from zomplex */
304 /* ------------------------------------------------------------------ */
305
306 case CHOLMOD_ZOMPLEX:
307
308 switch (xtype_out)
309 {
310
311 /* ---------------------------------------------------------- */
312 /* zomplex -> pattern */
313 /* ---------------------------------------------------------- */
314
315 case CHOLMOD_PATTERN:
316 /* free X and Z */
317 *XX = CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
318 *ZZ = CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
319 break ;
320
321 /* ---------------------------------------------------------- */
322 /* zomplex -> real */
323 /* ---------------------------------------------------------- */
324
325 case CHOLMOD_REAL:
326 /* free the imaginary part */
327 *ZZ = CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
328 break ;
329
330 /* ---------------------------------------------------------- */
331 /* zomplex -> complex */
332 /* ---------------------------------------------------------- */
333
334 case CHOLMOD_COMPLEX:
335 Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
336 if (Common->status < CHOLMOD_OK)
337 {
338 return (FALSE) ;
339 }
340 for (k = 0 ; k < nz ; k++)
341 {
342 Xnew [2*k ] = Xold [k] ;
343 Xnew [2*k+1] = Zold [k] ;
344 }
345 CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
346 CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
347 *XX = Xnew ;
348 *ZZ = NULL ;
349 break ;
350
351 }
352 break ;
353 }
354
355 return (TRUE) ;
356 }
357
358
359 /* ========================================================================== */
360 /* === cholmod_sparse_xtype ================================================= */
361 /* ========================================================================== */
362
363 /* Change the numeric xtype of a sparse matrix. Supports any type on input
364 * and output (pattern, real, complex, or zomplex). */
365
CHOLMOD(sparse_xtype)366 int CHOLMOD(sparse_xtype)
367 (
368 /* ---- input ---- */
369 int to_xtype, /* requested xtype */
370 /* ---- in/out --- */
371 cholmod_sparse *A, /* sparse matrix to change */
372 /* --------------- */
373 cholmod_common *Common
374 )
375 {
376 Int ok ;
377 RETURN_IF_NULL_COMMON (FALSE) ;
378 RETURN_IF_NULL (A, FALSE) ;
379 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
380
381 ok = change_complexity (A->nzmax, A->xtype, to_xtype,
382 CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, &(A->x), &(A->z), Common) ;
383 if (ok)
384 {
385 A->xtype = to_xtype ;
386 }
387 return (ok) ;
388 }
389
390
391 /* ========================================================================== */
392 /* === cholmod_triplet_xtype ================================================ */
393 /* ========================================================================== */
394
395 /* Change the numeric xtype of a triplet matrix. Supports any type on input
396 * and output (pattern, real, complex, or zomplex). */
397
CHOLMOD(triplet_xtype)398 int CHOLMOD(triplet_xtype)
399 (
400 /* ---- input ---- */
401 int to_xtype, /* requested xtype */
402 /* ---- in/out --- */
403 cholmod_triplet *T, /* triplet matrix to change */
404 /* --------------- */
405 cholmod_common *Common
406 )
407 {
408 Int ok ;
409 RETURN_IF_NULL_COMMON (FALSE) ;
410 RETURN_IF_NULL (T, FALSE) ;
411 RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
412 ok = change_complexity (T->nzmax, T->xtype, to_xtype,
413 CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, &(T->x), &(T->z), Common) ;
414 if (ok)
415 {
416 T->xtype = to_xtype ;
417 }
418 return (ok) ;
419 }
420
421
422 /* ========================================================================== */
423 /* === cholmod_dense_xtype ================================================= */
424 /* ========================================================================== */
425
426 /* Change the numeric xtype of a dense matrix. Supports real, complex or
427 * zomplex on input and output */
428
CHOLMOD(dense_xtype)429 int CHOLMOD(dense_xtype)
430 (
431 /* ---- input ---- */
432 int to_xtype, /* requested xtype */
433 /* ---- in/out --- */
434 cholmod_dense *X, /* dense matrix to change */
435 /* --------------- */
436 cholmod_common *Common
437 )
438 {
439 Int ok ;
440 RETURN_IF_NULL_COMMON (FALSE) ;
441 RETURN_IF_NULL (X, FALSE) ;
442 RETURN_IF_XTYPE_INVALID (X, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
443 ok = change_complexity (X->nzmax, X->xtype, to_xtype,
444 CHOLMOD_REAL, CHOLMOD_ZOMPLEX, &(X->x), &(X->z), Common) ;
445 if (ok)
446 {
447 X->xtype = to_xtype ;
448 }
449 return (ok) ;
450 }
451
452
453 /* ========================================================================== */
454 /* === cholmod_factor_xtype ================================================= */
455 /* ========================================================================== */
456
457 /* Change the numeric xtype of a factor. Supports real, complex or zomplex on
458 * input and output. Supernodal zomplex factors are not supported. */
459
CHOLMOD(factor_xtype)460 int CHOLMOD(factor_xtype)
461 (
462 /* ---- input ---- */
463 int to_xtype, /* requested xtype */
464 /* ---- in/out --- */
465 cholmod_factor *L, /* factor to change */
466 /* --------------- */
467 cholmod_common *Common
468 )
469 {
470 Int ok ;
471 RETURN_IF_NULL_COMMON (FALSE) ;
472 RETURN_IF_NULL (L, FALSE) ;
473 RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
474 if (L->is_super &&
475 (L->xtype == CHOLMOD_ZOMPLEX || to_xtype == CHOLMOD_ZOMPLEX))
476 {
477 ERROR (CHOLMOD_INVALID, "invalid xtype for supernodal L") ;
478 return (FALSE) ;
479 }
480 ok = change_complexity ((L->is_super ? L->xsize : L->nzmax), L->xtype,
481 to_xtype, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, &(L->x), &(L->z), Common) ;
482 if (ok)
483 {
484 L->xtype = to_xtype ;
485 }
486 return (ok) ;
487 }
488