1 /* ========================================================================== */
2 /* === SuiteSparse_config =================================================== */
3 /* ========================================================================== */
4
5 /* SuiteSparse configuration : memory manager and printf functions. */
6
7 /* Copyright (c) 2013-2018, Timothy A. Davis. No licensing restrictions
8 * apply to this file or to the SuiteSparse_config directory.
9 * Author: Timothy A. Davis.
10 */
11
12 #include <math.h>
13 #include <stdlib.h>
14
15 #ifndef NPRINT
16 #include <stdio.h>
17 #endif
18
19 #ifdef MATLAB_MEX_FILE
20 #include "mex.h"
21 #include "matrix.h"
22 #endif
23
24 #ifndef NULL
25 #define NULL ((void *) 0)
26 #endif
27
28 #include "SuiteSparse_config.h"
29
30 /* -------------------------------------------------------------------------- */
31 /* SuiteSparse_config : a global extern struct */
32 /* -------------------------------------------------------------------------- */
33
34 /* The SuiteSparse_config struct is available to all SuiteSparse functions and
35 to all applications that use those functions. It must be modified with
36 care, particularly in a multithreaded context. Normally, the application
37 will initialize this object once, via SuiteSparse_start, possibily followed
38 by application-specific modifications if the applications wants to use
39 alternative memory manager functions.
40
41 The user can redefine these global pointers at run-time to change the
42 memory manager and printf function used by SuiteSparse.
43
44 If -DNMALLOC is defined at compile-time, then no memory-manager is
45 specified. You must define them at run-time, after calling
46 SuiteSparse_start.
47
48 If -DPRINT is defined a compile time, then printf is disabled, and
49 SuiteSparse will not use printf.
50 */
51
52 struct SuiteSparse_config_struct SuiteSparse_config =
53 {
54
55 /* memory management functions */
56 #ifndef NMALLOC
57 #ifdef MATLAB_MEX_FILE
58 /* MATLAB mexFunction: */
59 mxMalloc, mxCalloc, mxRealloc, mxFree,
60 #else
61 /* standard ANSI C: */
62 malloc, calloc, realloc, free,
63 #endif
64 #else
65 /* no memory manager defined; you must define one at run-time: */
66 NULL, NULL, NULL, NULL,
67 #endif
68
69 /* printf function */
70 #ifndef NPRINT
71 #ifdef MATLAB_MEX_FILE
72 /* MATLAB mexFunction: */
73 mexPrintf,
74 #else
75 /* standard ANSI C: */
76 printf,
77 #endif
78 #else
79 /* printf is disabled */
80 NULL,
81 #endif
82
83 SuiteSparse_hypot,
84 SuiteSparse_divcomplex
85
86 } ;
87
88 /* -------------------------------------------------------------------------- */
89 /* SuiteSparse_start */
90 /* -------------------------------------------------------------------------- */
91
92 /* All applications that use SuiteSparse should call SuiteSparse_start prior
93 to using any SuiteSparse function. Only a single thread should call this
94 function, in a multithreaded application. Currently, this function is
95 optional, since all this function currently does is to set the four memory
96 function pointers to NULL (which tells SuiteSparse to use the default
97 functions). In a multi- threaded application, only a single thread should
98 call this function.
99
100 Future releases of SuiteSparse might enforce a requirement that
101 SuiteSparse_start be called prior to calling any SuiteSparse function.
102 */
103
SuiteSparse_start(void)104 void SuiteSparse_start ( void )
105 {
106
107 /* memory management functions */
108 #ifndef NMALLOC
109 #ifdef MATLAB_MEX_FILE
110 /* MATLAB mexFunction: */
111 SuiteSparse_config.malloc_func = mxMalloc ;
112 SuiteSparse_config.calloc_func = mxCalloc ;
113 SuiteSparse_config.realloc_func = mxRealloc ;
114 SuiteSparse_config.free_func = mxFree ;
115 #else
116 /* standard ANSI C: */
117 SuiteSparse_config.malloc_func = malloc ;
118 SuiteSparse_config.calloc_func = calloc ;
119 SuiteSparse_config.realloc_func = realloc ;
120 SuiteSparse_config.free_func = free ;
121 #endif
122 #else
123 /* no memory manager defined; you must define one after calling
124 SuiteSparse_start */
125 SuiteSparse_config.malloc_func = NULL ;
126 SuiteSparse_config.calloc_func = NULL ;
127 SuiteSparse_config.realloc_func = NULL ;
128 SuiteSparse_config.free_func = NULL ;
129 #endif
130
131 /* printf function */
132 #ifndef NPRINT
133 #ifdef MATLAB_MEX_FILE
134 /* MATLAB mexFunction: */
135 SuiteSparse_config.printf_func = mexPrintf ;
136 #else
137 /* standard ANSI C: */
138 SuiteSparse_config.printf_func = printf ;
139 #endif
140 #else
141 /* printf is disabled */
142 SuiteSparse_config.printf_func = NULL ;
143 #endif
144
145 /* math functions */
146 SuiteSparse_config.hypot_func = SuiteSparse_hypot ;
147 SuiteSparse_config.divcomplex_func = SuiteSparse_divcomplex ;
148 }
149
150 /* -------------------------------------------------------------------------- */
151 /* SuiteSparse_finish */
152 /* -------------------------------------------------------------------------- */
153
154 /* This currently does nothing, but in the future, applications should call
155 SuiteSparse_start before calling any SuiteSparse function, and then
156 SuiteSparse_finish after calling the last SuiteSparse function, just before
157 exiting. In a multithreaded application, only a single thread should call
158 this function.
159
160 Future releases of SuiteSparse might use this function for any
161 SuiteSparse-wide cleanup operations or finalization of statistics.
162 */
163
SuiteSparse_finish(void)164 void SuiteSparse_finish ( void )
165 {
166 /* do nothing */ ;
167 }
168
169 /* -------------------------------------------------------------------------- */
170 /* SuiteSparse_malloc: malloc wrapper */
171 /* -------------------------------------------------------------------------- */
172
SuiteSparse_malloc(size_t nitems,size_t size_of_item)173 void *SuiteSparse_malloc /* pointer to allocated block of memory */
174 (
175 size_t nitems, /* number of items to malloc */
176 size_t size_of_item /* sizeof each item */
177 )
178 {
179 void *p ;
180 size_t size ;
181 if (nitems < 1) nitems = 1 ;
182 if (size_of_item < 1) size_of_item = 1 ;
183 size = nitems * size_of_item ;
184
185 if (size != ((double) nitems) * size_of_item)
186 {
187 /* size_t overflow */
188 p = NULL ;
189 }
190 else
191 {
192 p = (void *) (SuiteSparse_config.malloc_func) (size) ;
193 }
194 return (p) ;
195 }
196
197
198 /* -------------------------------------------------------------------------- */
199 /* SuiteSparse_calloc: calloc wrapper */
200 /* -------------------------------------------------------------------------- */
201
SuiteSparse_calloc(size_t nitems,size_t size_of_item)202 void *SuiteSparse_calloc /* pointer to allocated block of memory */
203 (
204 size_t nitems, /* number of items to calloc */
205 size_t size_of_item /* sizeof each item */
206 )
207 {
208 void *p ;
209 size_t size ;
210 if (nitems < 1) nitems = 1 ;
211 if (size_of_item < 1) size_of_item = 1 ;
212 size = nitems * size_of_item ;
213
214 if (size != ((double) nitems) * size_of_item)
215 {
216 /* size_t overflow */
217 p = NULL ;
218 }
219 else
220 {
221 p = (void *) (SuiteSparse_config.calloc_func) (nitems, size_of_item) ;
222 }
223 return (p) ;
224 }
225
226 /* -------------------------------------------------------------------------- */
227 /* SuiteSparse_realloc: realloc wrapper */
228 /* -------------------------------------------------------------------------- */
229
230 /* If p is non-NULL on input, it points to a previously allocated object of
231 size nitems_old * size_of_item. The object is reallocated to be of size
232 nitems_new * size_of_item. If p is NULL on input, then a new object of that
233 size is allocated. On success, a pointer to the new object is returned,
234 and ok is returned as 1. If the allocation fails, ok is set to 0 and a
235 pointer to the old (unmodified) object is returned.
236 */
237
SuiteSparse_realloc(size_t nitems_new,size_t nitems_old,size_t size_of_item,void * p,int * ok)238 void *SuiteSparse_realloc /* pointer to reallocated block of memory, or
239 to original block if the realloc failed. */
240 (
241 size_t nitems_new, /* new number of items in the object */
242 size_t nitems_old, /* old number of items in the object */
243 size_t size_of_item, /* sizeof each item */
244 void *p, /* old object to reallocate */
245 int *ok /* 1 if successful, 0 otherwise */
246 )
247 {
248 size_t size ;
249 if (nitems_old < 1) nitems_old = 1 ;
250 if (nitems_new < 1) nitems_new = 1 ;
251 if (size_of_item < 1) size_of_item = 1 ;
252 size = nitems_new * size_of_item ;
253
254 if (size != ((double) nitems_new) * size_of_item)
255 {
256 /* size_t overflow */
257 (*ok) = 0 ;
258 }
259 else if (p == NULL)
260 {
261 /* a fresh object is being allocated */
262 p = SuiteSparse_malloc (nitems_new, size_of_item) ;
263 (*ok) = (p != NULL) ;
264 }
265 else if (nitems_old == nitems_new)
266 {
267 /* the object does not change; do nothing */
268 (*ok) = 1 ;
269 }
270 else
271 {
272 /* change the size of the object from nitems_old to nitems_new */
273 void *pnew ;
274 pnew = (void *) (SuiteSparse_config.realloc_func) (p, size) ;
275 if (pnew == NULL)
276 {
277 if (nitems_new < nitems_old)
278 {
279 /* the attempt to reduce the size of the block failed, but
280 the old block is unchanged. So pretend to succeed. */
281 (*ok) = 1 ;
282 }
283 else
284 {
285 /* out of memory */
286 (*ok) = 0 ;
287 }
288 }
289 else
290 {
291 /* success */
292 p = pnew ;
293 (*ok) = 1 ;
294 }
295 }
296 return (p) ;
297 }
298
299 /* -------------------------------------------------------------------------- */
300 /* SuiteSparse_free: free wrapper */
301 /* -------------------------------------------------------------------------- */
302
SuiteSparse_free(void * p)303 void *SuiteSparse_free /* always returns NULL */
304 (
305 void *p /* block to free */
306 )
307 {
308 if (p)
309 {
310 (SuiteSparse_config.free_func) (p) ;
311 }
312 return (NULL) ;
313 }
314
315
316 /* -------------------------------------------------------------------------- */
317 /* SuiteSparse_tic: return current wall clock time */
318 /* -------------------------------------------------------------------------- */
319
320 /* Returns the number of seconds (tic [0]) and nanoseconds (tic [1]) since some
321 * unspecified but fixed time in the past. If no timer is installed, zero is
322 * returned. A scalar double precision value for 'tic' could be used, but this
323 * might cause loss of precision because clock_getttime returns the time from
324 * some distant time in the past. Thus, an array of size 2 is used.
325 *
326 * The timer is enabled by default. To disable the timer, compile with
327 * -DNTIMER. If enabled on a POSIX C 1993 system, the timer requires linking
328 * with the -lrt library.
329 *
330 * example:
331 *
332 * double tic [2], r, s, t ;
333 * SuiteSparse_tic (tic) ; // start the timer
334 * // do some work A
335 * t = SuiteSparse_toc (tic) ; // t is time for work A, in seconds
336 * // do some work B
337 * s = SuiteSparse_toc (tic) ; // s is time for work A and B, in seconds
338 * SuiteSparse_tic (tic) ; // restart the timer
339 * // do some work C
340 * r = SuiteSparse_toc (tic) ; // s is time for work C, in seconds
341 *
342 * A double array of size 2 is used so that this routine can be more easily
343 * ported to non-POSIX systems. The caller does not rely on the POSIX
344 * <time.h> include file.
345 */
346
347 #ifdef SUITESPARSE_TIMER_ENABLED
348
349 #include <time.h>
350
SuiteSparse_tic(double tic[2])351 void SuiteSparse_tic
352 (
353 double tic [2] /* output, contents undefined on input */
354 )
355 {
356 /* POSIX C 1993 timer, requires -librt */
357 struct timespec t ;
358 clock_gettime (CLOCK_MONOTONIC, &t) ;
359 tic [0] = (double) (t.tv_sec) ;
360 tic [1] = (double) (t.tv_nsec) ;
361 }
362
363 #else
364
SuiteSparse_tic(double tic[2])365 void SuiteSparse_tic
366 (
367 double tic [2] /* output, contents undefined on input */
368 )
369 {
370 /* no timer installed */
371 tic [0] = 0 ;
372 tic [1] = 0 ;
373 }
374
375 #endif
376
377
378 /* -------------------------------------------------------------------------- */
379 /* SuiteSparse_toc: return time since last tic */
380 /* -------------------------------------------------------------------------- */
381
382 /* Assuming SuiteSparse_tic is accurate to the nanosecond, this function is
383 * accurate down to the nanosecond for 2^53 nanoseconds since the last call to
384 * SuiteSparse_tic, which is sufficient for SuiteSparse (about 104 days). If
385 * additional accuracy is required, the caller can use two calls to
386 * SuiteSparse_tic and do the calculations differently.
387 */
388
SuiteSparse_toc(double tic[2])389 double SuiteSparse_toc /* returns time in seconds since last tic */
390 (
391 double tic [2] /* input, not modified from last call to SuiteSparse_tic */
392 )
393 {
394 double toc [2] ;
395 SuiteSparse_tic (toc) ;
396 return ((toc [0] - tic [0]) + 1e-9 * (toc [1] - tic [1])) ;
397 }
398
399
400 /* -------------------------------------------------------------------------- */
401 /* SuiteSparse_time: return current wallclock time in seconds */
402 /* -------------------------------------------------------------------------- */
403
404 /* This function might not be accurate down to the nanosecond. */
405
SuiteSparse_time(void)406 double SuiteSparse_time /* returns current wall clock time in seconds */
407 (
408 void
409 )
410 {
411 double toc [2] ;
412 SuiteSparse_tic (toc) ;
413 return (toc [0] + 1e-9 * toc [1]) ;
414 }
415
416
417 /* -------------------------------------------------------------------------- */
418 /* SuiteSparse_version: return the current version of SuiteSparse */
419 /* -------------------------------------------------------------------------- */
420
SuiteSparse_version(int version[3])421 int SuiteSparse_version
422 (
423 int version [3]
424 )
425 {
426 if (version != NULL)
427 {
428 version [0] = SUITESPARSE_MAIN_VERSION ;
429 version [1] = SUITESPARSE_SUB_VERSION ;
430 version [2] = SUITESPARSE_SUBSUB_VERSION ;
431 }
432 return (SUITESPARSE_VERSION) ;
433 }
434
435 /* -------------------------------------------------------------------------- */
436 /* SuiteSparse_hypot */
437 /* -------------------------------------------------------------------------- */
438
439 /* There is an equivalent routine called hypot in <math.h>, which conforms
440 * to ANSI C99. However, SuiteSparse does not assume that ANSI C99 is
441 * available. You can use the ANSI C99 hypot routine with:
442 *
443 * #include <math.h>
444 *i SuiteSparse_config.hypot_func = hypot ;
445 *
446 * Default value of the SuiteSparse_config.hypot_func pointer is
447 * SuiteSparse_hypot, defined below.
448 *
449 * s = hypot (x,y) computes s = sqrt (x*x + y*y) but does so more accurately.
450 * The NaN cases for the double relops x >= y and x+y == x are safely ignored.
451 *
452 * Source: Algorithm 312, "Absolute value and square root of a complex number,"
453 * P. Friedland, Comm. ACM, vol 10, no 10, October 1967, page 665.
454 */
455
SuiteSparse_hypot(double x,double y)456 double SuiteSparse_hypot (double x, double y)
457 {
458 double s, r ;
459 x = fabs (x) ;
460 y = fabs (y) ;
461 if (x >= y)
462 {
463 if (x + y == x)
464 {
465 s = x ;
466 }
467 else
468 {
469 r = y / x ;
470 s = x * sqrt (1.0 + r*r) ;
471 }
472 }
473 else
474 {
475 if (y + x == y)
476 {
477 s = y ;
478 }
479 else
480 {
481 r = x / y ;
482 s = y * sqrt (1.0 + r*r) ;
483 }
484 }
485 return (s) ;
486 }
487
488 /* -------------------------------------------------------------------------- */
489 /* SuiteSparse_divcomplex */
490 /* -------------------------------------------------------------------------- */
491
492 /* c = a/b where c, a, and b are complex. The real and imaginary parts are
493 * passed as separate arguments to this routine. The NaN case is ignored
494 * for the double relop br >= bi. Returns 1 if the denominator is zero,
495 * 0 otherwise.
496 *
497 * This uses ACM Algo 116, by R. L. Smith, 1962, which tries to avoid
498 * underflow and overflow.
499 *
500 * c can be the same variable as a or b.
501 *
502 * Default value of the SuiteSparse_config.divcomplex_func pointer is
503 * SuiteSparse_divcomplex.
504 */
505
SuiteSparse_divcomplex(double ar,double ai,double br,double bi,double * cr,double * ci)506 int SuiteSparse_divcomplex
507 (
508 double ar, double ai, /* real and imaginary parts of a */
509 double br, double bi, /* real and imaginary parts of b */
510 double *cr, double *ci /* real and imaginary parts of c */
511 )
512 {
513 double tr, ti, r, den ;
514 if (fabs (br) >= fabs (bi))
515 {
516 r = bi / br ;
517 den = br + r * bi ;
518 tr = (ar + ai * r) / den ;
519 ti = (ai - ar * r) / den ;
520 }
521 else
522 {
523 r = br / bi ;
524 den = r * br + bi ;
525 tr = (ar * r + ai) / den ;
526 ti = (ai * r - ar) / den ;
527 }
528 *cr = tr ;
529 *ci = ti ;
530 return (den == 0.) ;
531 }
532