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