1 //------------------------------------------------------------------------------
2 // SLIP_LU/slip_cast_array: ` and typecast an array
3 //------------------------------------------------------------------------------
4
5 // SLIP_LU: (c) 2019-2020, Chris Lourenco, Jinhao Chen, Erick Moreno-Centeno,
6 // Timothy A. Davis, Texas A&M University All Rights Reserved. See
7 // SLIP_LU/License for the license.
8
9 //------------------------------------------------------------------------------
10
11 // Scales and typecasts an input array X, into the output array Y.
12
13 // X: an array of type xtype, of size n.
14 // Y: an array of type ytype, of size n.
15
16 // Note about the scaling factors:
17 //
18 // This function copies the scaled values of X into the array Y.
19 // If Y is mpz_t, the values in X must be scaled to be integral
20 // if they are not already integral. As a result, this function
21 // will expand the values and set y_scale = this factor.
22 //`
23 // Conversely, if Y is not mpz_t and X is, we apply X's scaling
24 // factor here to get the final values of Y. For instance,
25 // if Y is FP64 and X is mpz_t, the values of Y are obtained
26 // as Y = X / x_scale.
27 //
28 // The final value of x_scale is not modified.
29 // The final value of y_scale is set as follows.
30 // If Y is mpz_t, y_scale is set to the appropriate value
31 // in order to make all of its entries integral.
32 //
33 // If Y is any other data type, this function always sets
34 // y_scale = 1.
35 //
36
37 #define SLIP_FREE_ALL \
38 SLIP_MPQ_CLEAR(temp); \
39
40 #include "slip_internal.h"
41 #pragma GCC diagnostic ignored "-Wunused-variable"
42
slip_cast_array(void * Y,SLIP_type ytype,void * X,SLIP_type xtype,int64_t n,mpq_t y_scale,mpq_t x_scale,const SLIP_options * option)43 SLIP_info slip_cast_array
44 (
45 void *Y, // output array, of size n
46 SLIP_type ytype, // type of Y
47 void *X, // input array, of size n
48 SLIP_type xtype, // type of X
49 int64_t n, // size of Y and X
50 mpq_t y_scale, // scale factor applied if Y is mpz_t
51 mpq_t x_scale, // scale factor applied if x is mpz_t
52 const SLIP_options *option// Command options. If NULL, set to default values
53 )
54 {
55
56 //--------------------------------------------------------------------------
57 // check inputs
58 // xtype and ytype are checked in SLIP_matrix_copy
59 //--------------------------------------------------------------------------
60
61 if (Y == NULL || X == NULL)
62 {
63 return (SLIP_INCORRECT_INPUT) ;
64 }
65 SLIP_info info ;
66 int r;
67 mpq_t temp; SLIP_MPQ_SET_NULL(temp);
68
69 mpfr_rnd_t round = SLIP_OPTION_ROUND (option) ;
70
71 //--------------------------------------------------------------------------
72 // Y [0:n-1] = (ytype) X [0:n-1]
73 //--------------------------------------------------------------------------
74
75 switch (ytype)
76 {
77
78 //----------------------------------------------------------------------
79 // output array Y is mpz_t
80 // If X is not mpz_t or int64, the values of X are scaled and y_scale is
81 // set to be this scaling factor.
82 //----------------------------------------------------------------------
83
84 case SLIP_MPZ:
85 {
86 mpz_t *y = (mpz_t *) Y ;
87 switch (xtype)
88 {
89
90 case SLIP_MPZ: // mpz_t to mpz_t
91 {
92 mpz_t *x = (mpz_t *) X ;
93 for (int64_t k = 0 ; k < n ; k++)
94 {
95 SLIP_CHECK (SLIP_mpz_set (y [k], x[k])) ;
96 }
97 // y is a direct copy of x. Set y_scale = x_scale
98 SLIP_CHECK(SLIP_mpq_set(y_scale, x_scale));
99 }
100 break ;
101
102 case SLIP_MPQ: // mpq_t to mpz_t
103 {
104 mpq_t *x = (mpq_t *) X ;
105 SLIP_CHECK (slip_expand_mpq_array(Y, X, y_scale, n,option));
106 }
107 break ;
108
109 case SLIP_MPFR: // mpfr_t to mpz_t
110 {
111 mpfr_t *x = (mpfr_t *) X ;
112 SLIP_CHECK (slip_expand_mpfr_array (Y, X, y_scale, n,
113 option)) ;
114 }
115 break ;
116
117 case SLIP_INT64: // int64_t to mpz_t
118 {
119 int64_t *x = (int64_t *) X ;
120 for (int64_t k = 0 ; k < n ; k++)
121 {
122 SLIP_CHECK (SLIP_mpz_set_si (y [k], x [k])) ;
123 }
124 SLIP_CHECK (SLIP_mpq_set_ui (y_scale, 1, 1)) ;
125 }
126 break ;
127
128 case SLIP_FP64: // double to mpz_t
129 {
130 double *x = (double *) X ;
131 SLIP_CHECK (slip_expand_double_array (y, x, y_scale, n,
132 option)) ;
133 }
134 break ;
135
136 }
137 }
138 break ;
139
140 //----------------------------------------------------------------------
141 // output array Y is mpq_t
142 //----------------------------------------------------------------------
143
144 case SLIP_MPQ:
145 {
146 mpq_t *y = (mpq_t *) Y ;
147 switch (xtype)
148 {
149
150 case SLIP_MPZ: // mpz_t to mpq_t
151 {
152 // In this case, x is mpz_t and y is mpq_t. the scaling
153 // factor x_scale must be used. If x_scale is not equal to
154 // 1, each value in y is divided by x_scale
155
156 // Check if x_scale == 1
157 SLIP_CHECK(SLIP_mpq_cmp_ui(&r, x_scale, 1, 1));
158 mpz_t *x = (mpz_t *) X ;
159
160 if (r == 0)
161 {
162 // x_scale = 1. Simply do a direct copy.
163 for (int64_t k = 0 ; k < n ; k++)
164 {
165 SLIP_CHECK (SLIP_mpq_set_z (y [k], x [k])) ;
166 }
167 }
168 else
169 {
170 // x_scale != 1. In this case, we divide each entry
171 // of Y by x_scale
172 for (int64_t k = 0 ; k < n ; k++)
173 {
174 SLIP_CHECK (SLIP_mpq_set_z (y [k], x [k])) ;
175 SLIP_CHECK (SLIP_mpq_div(y[k], y[k], x_scale));
176 }
177 }
178 }
179 break ;
180
181 case SLIP_MPQ: // mpq_t to mpq_t
182 {
183 mpq_t *x = (mpq_t *) X ;
184 for (int64_t k = 0 ; k < n ; k++)
185 {
186 SLIP_CHECK (SLIP_mpq_set (y [k], x [k])) ;
187 }
188 }
189 break ;
190
191 case SLIP_MPFR: // mpfr_t to mpq_t
192 {
193 mpfr_t *x = (mpfr_t *) X ;
194 for (int64_t k = 0 ; k < n ; k++)
195 {
196 SLIP_CHECK (SLIP_mpfr_get_q( y[k], x[k], round));
197 }
198 }
199 break ;
200
201 case SLIP_INT64: // int64 to mpq_t
202 {
203 int64_t *x = (int64_t *) X ;
204 for (int64_t k = 0 ; k < n ; k++)
205 {
206 SLIP_CHECK (SLIP_mpq_set_si (y [k], x [k], 1)) ;
207 }
208 }
209 break ;
210
211 case SLIP_FP64: // double to mpq_t
212 {
213 double *x = (double *) X ;
214 for (int64_t k = 0 ; k < n ; k++)
215 {
216 SLIP_CHECK (SLIP_mpq_set_d (y [k], x [k])) ;
217 }
218 }
219 break ;
220
221 }
222 }
223 break ;
224
225 //----------------------------------------------------------------------
226 // output array Y is mpfr_t
227 //----------------------------------------------------------------------
228
229 case SLIP_MPFR:
230 {
231 mpfr_t *y = (mpfr_t *) Y ;
232 switch (xtype)
233 {
234
235 case SLIP_MPZ: // mpz_t to mpfr_t
236 {
237 // x is mpz_t and y is mpfr_t. Like in the above mpq_t
238 // case, if the scaling factor of x is not equal to 1, the
239 // values of y must be scaled.
240 mpz_t *x = (mpz_t *) X ;
241 SLIP_CHECK(SLIP_mpq_cmp_ui(&r, x_scale, 1, 1));
242
243 if (r == 0)
244 {
245 // x_scale = 1. Simply do a direct copy.
246 for (int64_t k = 0 ; k < n ; k++)
247 {
248 SLIP_CHECK (SLIP_mpfr_set_z (y [k], x [k], round)) ;
249 }
250 }
251 else
252 {
253 // x_scale != 1. In this case, we divide each entry
254 // of Y by x_scale. To do this, we will cast each
255 // x_k to mpq_t, then divide by the scale, then
256 // cast the result to mpfr_t
257 SLIP_CHECK(SLIP_mpq_init(temp));
258 for (int64_t k = 0 ; k < n ; k++)
259 {
260 SLIP_CHECK( SLIP_mpq_set_z( temp, x[k]));
261 SLIP_CHECK( SLIP_mpq_div(temp, temp, x_scale));
262 SLIP_CHECK(SLIP_mpfr_set_q(y[k], temp, round));
263 }
264 }
265 }
266 break ;
267
268 case SLIP_MPQ: // mpq_t to mpfr_t
269 {
270 mpq_t *x = (mpq_t *) X ;
271 for (int64_t k = 0 ; k < n ; k++)
272 {
273 SLIP_CHECK (SLIP_mpfr_set_q (y [k], x [k], round)) ;
274 }
275 }
276 break ;
277
278 case SLIP_MPFR: // mpfr_t to mpfr_t
279 {
280 mpfr_t *x = (mpfr_t *) X ;
281 for (int64_t k = 0 ; k < n ; k++)
282 {
283 SLIP_CHECK (SLIP_mpfr_set (y [k], x [k], round)) ;
284 }
285 }
286 break ;
287
288 case SLIP_INT64: // int64 to mpfr_t
289 {
290 int64_t *x = (int64_t *) X ;
291 for (int64_t k = 0 ; k < n ; k++)
292 {
293 SLIP_CHECK(SLIP_mpfr_set_si(y[k], x[k], round));
294 }
295 }
296 break ;
297
298 case SLIP_FP64: // double to mpfr_t
299 {
300 double *x = (double *) X ;
301 for (int64_t k = 0 ; k < n ; k++)
302 {
303 SLIP_CHECK (SLIP_mpfr_set_d (y [k], x [k], round)) ;
304 }
305 }
306 break ;
307
308 }
309 }
310 break ;
311
312 //----------------------------------------------------------------------
313 // output array Y is int64_t
314 //----------------------------------------------------------------------
315
316 case SLIP_INT64:
317 {
318 int64_t *y = (int64_t *) Y ;
319 switch (xtype)
320 {
321
322 case SLIP_MPZ: // mpz_t to int64_t
323 {
324 // x is mpz_t and y is int64_t. Same as above,
325 // if x_scale > 1 it is applied
326 mpz_t *x = (mpz_t *) X ;
327 SLIP_CHECK(SLIP_mpq_cmp_ui(&r, x_scale, 1, 1));
328
329 if (r == 0)
330 {
331 // x_scale = 1. Simply do a direct copy.
332 for (int64_t k = 0 ; k < n ; k++)
333 {
334 SLIP_CHECK(SLIP_mpz_get_si( &(y[k]), x[k]));
335 }
336 }
337 else
338 {
339 // x_scale != 1. In this case, we divide each entry
340 // of Y by x_scale. To do this, we will cast each
341 // x_k to mpq_t, then divide by the scale, then
342 // cast the result to double and cast the double to int
343 SLIP_CHECK(SLIP_mpq_init(temp));
344 for (int64_t k = 0 ; k < n ; k++)
345 {
346 SLIP_CHECK( SLIP_mpq_set_z( temp, x[k]));
347 SLIP_CHECK( SLIP_mpq_div(temp, temp, x_scale));
348 double temp2;
349 SLIP_CHECK(SLIP_mpq_get_d(&temp2, temp));
350 y[k] = slip_cast_double_to_int64(temp2);
351 }
352 }
353 }
354 break ;
355
356 case SLIP_MPQ: // mpq_t to int64_t
357 {
358 mpq_t *x = (mpq_t *) X ;
359 for (int64_t k = 0 ; k < n ; k++)
360 {
361 double t ;
362 SLIP_CHECK (SLIP_mpq_get_d (&t, x [k])) ;
363 y [k] = slip_cast_double_to_int64 (t) ;
364 }
365 }
366 break ;
367
368 case SLIP_MPFR: // mpfr_t to int64_t
369 {
370 mpfr_t *x = (mpfr_t *) X ;
371 for (int64_t k = 0 ; k < n ; k++)
372 {
373 SLIP_CHECK( SLIP_mpfr_get_si( &(y[k]),x[k], round));
374 }
375 }
376 break ;
377
378 case SLIP_INT64: // int64_t to int64_t
379 {
380 memcpy (Y, X, n * sizeof (int64_t)) ;
381 }
382 break ;
383
384 case SLIP_FP64: // double to int64_t
385 {
386 double *x = (double *) X ;
387 for (int64_t k = 0 ; k < n ; k++)
388 {
389 y [k] = slip_cast_double_to_int64 (x [k]) ;
390 }
391 }
392 break ;
393
394 }
395 }
396 break ;
397
398 //----------------------------------------------------------------------
399 // output array Y is double
400 //----------------------------------------------------------------------
401
402 case SLIP_FP64:
403 {
404 double *y = (double *) Y ;
405 switch (xtype)
406 {
407
408 case SLIP_MPZ: // mpz_t to double
409 {
410 // Same as above, x is mpz_t, y is double. Must
411 // divide by x_scale if x_scale != 1.
412 mpz_t *x = (mpz_t *) X ;
413 SLIP_CHECK(SLIP_mpq_cmp_ui(&r, x_scale, 1, 1));
414
415 if (r == 0)
416 {
417 // x_scale = 1. Simply do a direct copy.
418 for (int64_t k = 0 ; k < n ; k++)
419 {
420 SLIP_CHECK(SLIP_mpz_get_d( &(y[k]), x[k]));
421 }
422 }
423 else
424 {
425 // x_scale != 1. In this case, we divide each entry
426 // of Y by x_scale. To do this, we will cast each
427 // x_k to mpq_t, then divide by the scale, then
428 // cast the result to double
429 SLIP_CHECK(SLIP_mpq_init(temp));
430 for (int64_t k = 0 ; k < n ; k++)
431 {
432 SLIP_CHECK( SLIP_mpq_set_z( temp, x[k]));
433 SLIP_CHECK( SLIP_mpq_div(temp, temp, x_scale));
434 SLIP_CHECK(SLIP_mpq_get_d(&(y[k]), temp));
435 }
436 }
437 }
438 break ;
439
440 case SLIP_MPQ: // mpq_t to double
441 {
442 mpq_t *x = (mpq_t *) X ;
443 for (int64_t k = 0 ; k < n ; k++)
444 {
445 SLIP_CHECK (SLIP_mpq_get_d (&(y [k]), x [k])) ;
446 }
447 }
448 break ;
449
450 case SLIP_MPFR: // mpfr_t to double
451 {
452 mpfr_t *x = (mpfr_t *) X ;
453 for (int64_t k = 0 ; k < n ; k++)
454 {
455 SLIP_CHECK (SLIP_mpfr_get_d (&(y [k]), x [k],
456 round));
457 }
458 }
459 break ;
460
461 case SLIP_INT64: // int64_t to double
462 {
463 int64_t *x = (int64_t *) X ;
464 for (int64_t k = 0 ; k < n ; k++)
465 {
466 y [k] = (double) (x [k]) ;
467 }
468 }
469 break ;
470
471 case SLIP_FP64: // double to double
472 {
473 memcpy (Y, X, n * sizeof (double)) ;
474 }
475 break ;
476
477 }
478 }
479 break ;
480
481 }
482 SLIP_FREE_ALL
483 return (SLIP_OK) ;
484 }
485