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