1 //------------------------------------------------------------------------------
2 // SLIP_LU/SLIP_matrix_copy: create a copy of a matrix
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 // SLIP_matrix_copy creates a SLIP_matrix C that is a modified copy of a
12 // SLIP_matrix A. The new matrix C can have a different kind and type
13 // than A.
14
15 // The input matrix A is assumed to be valid. It can be checked first with
16 // SLIP_matrix_check, if desired. If the input matrix A is not valid, results
17 // are undefined.
18
19 #define SLIP_FREE_WORK \
20 SLIP_matrix_free (&T, option) ; \
21 SLIP_matrix_free (&Y, option) ; \
22 SLIP_FREE (W) ;
23
24 #define SLIP_FREE_ALL \
25 SLIP_FREE_WORK ; \
26 SLIP_matrix_free (&C, option) ;
27
28 #include "slip_internal.h"
29
SLIP_matrix_copy(SLIP_matrix ** C_handle,SLIP_kind C_kind,SLIP_type C_type,SLIP_matrix * A,const SLIP_options * option)30 SLIP_info SLIP_matrix_copy
31 (
32 SLIP_matrix **C_handle, // matrix to create (never shallow)
33 // inputs, not modified:
34 SLIP_kind C_kind, // C->kind: CSC, triplet, or dense
35 SLIP_type C_type, // C->type: mpz_t, mpq_t, mpfr_t, int64_t, or double
36 SLIP_matrix *A, // matrix to make a copy of (may be shallow)
37 const SLIP_options *option
38 )
39 {
40
41 //--------------------------------------------------------------------------
42 // check inputs
43 //--------------------------------------------------------------------------
44
45 SLIP_info info ;
46 if (!slip_initialized ( )) return (SLIP_PANIC) ;
47
48 int64_t nz = SLIP_matrix_nnz (A, option) ;
49 if (C_handle == NULL || nz < 0 ||
50 //checked in SLIP_matrix_nnz
51 //A == NULL || A->kind < SLIP_CSC || A->kind > SLIP_DENSE ||
52 A->type < SLIP_MPZ || A->type > SLIP_FP64 ||
53 C_kind < SLIP_CSC || C_kind > SLIP_DENSE ||
54 C_type < SLIP_MPZ || C_type > SLIP_FP64)
55 {
56 return (SLIP_INCORRECT_INPUT) ;
57 }
58 (*C_handle) = NULL ;
59 SLIP_matrix *C = NULL ;
60 SLIP_matrix *Y = NULL ;
61 SLIP_matrix *T = NULL ;
62 int64_t *W = NULL ;
63 int64_t m = A->m ;
64 int64_t n = A->n ;
65 mpfr_rnd_t round = SLIP_OPTION_ROUND (option) ;
66
67 //--------------------------------------------------------------------------
68 // copy and convert A into C
69 //--------------------------------------------------------------------------
70
71 switch (C_kind)
72 {
73
74 //----------------------------------------------------------------------
75 // C is CSC
76 //----------------------------------------------------------------------
77
78 case SLIP_CSC:
79 {
80
81 switch (A->kind)
82 {
83
84 //--------------------------------------------------------------
85 // A is CSC, C is CSC
86 //--------------------------------------------------------------
87
88 case SLIP_CSC:
89 {
90 // allocate C
91 SLIP_CHECK (SLIP_matrix_allocate (&C, SLIP_CSC, C_type,
92 m, n, nz, false, true, option)) ;
93 // copy the pattern of A into C
94 memcpy (C->p, A->p, (n+1) * sizeof (int64_t)) ;
95 memcpy (C->i, A->i, nz * sizeof (int64_t)) ;
96 // copy and typecast A->x into C->x
97 SLIP_CHECK (slip_cast_array (SLIP_X (C), C->type,
98 SLIP_X (A), A->type, nz, C->scale, A->scale, option)) ;
99 }
100 break ;
101
102 //--------------------------------------------------------------
103 // A is triplet, C is CSC
104 //--------------------------------------------------------------
105
106 case SLIP_TRIPLET:
107 {
108
109 // Y = typecast the values of A into the type of C
110 // (not the pattern; Y is SLIP_DENSE)
111 SLIP_CHECK (slip_cast_matrix (&Y, C_type, A, option)) ;
112
113 // allocate workspace
114 W = (int64_t *) SLIP_calloc (n, sizeof (int64_t)) ;
115 if (W == NULL)
116 {
117 SLIP_FREE_ALL ;
118 return (SLIP_OUT_OF_MEMORY) ;
119 }
120
121 // allocate C
122 SLIP_CHECK (SLIP_matrix_allocate (&C, SLIP_CSC,
123 C_type, m, n, nz, false, true, option)) ;
124
125 // Scaling factor of C is currently in Y, set it
126 // here
127 SLIP_mpq_set(C->scale, Y->scale);
128
129 // count the # of entries in each column
130 for (int64_t k = 0 ; k < nz ; k++)
131 {
132 W [A->j [k]]++ ;
133 }
134
135 // C->p = cumulative sum of W
136 slip_cumsum (C->p, W, n) ;
137
138 // build the matrix
139 switch (C->type)
140 {
141 case SLIP_MPZ:
142 for (int64_t k = 0 ; k < nz ; k++)
143 {
144 int64_t p = W [A->j [k]]++ ;
145 C->i [p] = A->i [k] ;
146 SLIP_CHECK (SLIP_mpz_set (
147 SLIP_1D (C, p, mpz),
148 SLIP_1D (Y, k, mpz))) ;
149 }
150 break ;
151
152 case SLIP_MPQ:
153 for (int64_t k = 0 ; k < nz ; k++)
154 {
155 int64_t p = W [A->j [k]]++ ;
156 C->i [p] = A->i [k] ;
157 SLIP_CHECK (SLIP_mpq_set (
158 SLIP_1D (C, p, mpq),
159 SLIP_1D (Y, k, mpq))) ;
160 }
161 break ;
162
163 case SLIP_MPFR:
164 for (int64_t k = 0 ; k < nz ; k++)
165 {
166 int64_t p = W [A->j [k]]++ ;
167 C->i [p] = A->i [k] ;
168 SLIP_CHECK (SLIP_mpfr_set (
169 SLIP_1D (C, p, mpfr),
170 SLIP_1D (Y, k, mpfr),
171 round)) ;
172 }
173 break ;
174
175 case SLIP_INT64:
176 for (int64_t k = 0 ; k < nz ; k++)
177 {
178 int64_t p = W [A->j [k]]++ ;
179 C->i [p] = A->i [k] ;
180 SLIP_1D (C, p, int64) =
181 SLIP_1D (Y, k, int64) ;
182 }
183 break ;
184
185 case SLIP_FP64:
186 for (int64_t k = 0 ; k < nz ; k++)
187 {
188 int64_t p = W [A->j [k]]++ ;
189 C->i [p] = A->i [k] ;
190 SLIP_1D (C, p, fp64) =
191 SLIP_1D (Y, k, fp64) ;
192 }
193 break ;
194
195 }
196
197 }
198 break ;
199
200 //--------------------------------------------------------------
201 // A is dense, C is CSC
202 //--------------------------------------------------------------
203
204 case SLIP_DENSE:
205 {
206 // Y = typecast the values of A into the type of C
207 SLIP_CHECK (slip_cast_matrix (&Y, C_type, A, option)) ;
208 int s ;
209
210 // count the actual nonzeros in Y
211 int64_t actual = 0 ;
212 switch (Y->type)
213 {
214
215 case SLIP_MPZ:
216 for (int64_t k = 0 ; k < nz ; k++)
217 {
218 SLIP_CHECK (SLIP_mpz_sgn (&s,
219 SLIP_1D (Y, k, mpz))) ;
220 if (s != 0) actual++ ;
221 }
222 break ;
223
224 case SLIP_MPQ:
225 for (int64_t k = 0 ; k < nz ; k++)
226 {
227 SLIP_CHECK (SLIP_mpq_sgn (&s,
228 SLIP_1D (Y, k, mpq))) ;
229 if (s != 0) actual++ ;
230 }
231 break ;
232
233 case SLIP_MPFR:
234 for (int64_t k = 0 ; k < nz ; k++)
235 {
236 SLIP_CHECK (SLIP_mpfr_sgn (&s,
237 SLIP_1D (Y, k, mpfr))) ;
238 if (s != 0) actual++ ;
239 }
240 break ;
241
242 case SLIP_INT64:
243 for (int64_t k = 0 ; k < nz ; k++)
244 {
245 if (SLIP_1D (Y, k, int64) != 0) actual++ ;
246 }
247 break ;
248
249 case SLIP_FP64:
250 for (int64_t k = 0 ; k < nz ; k++)
251 {
252 if (SLIP_1D (Y, k, fp64) != 0) actual++ ;
253 }
254 break ;
255
256 }
257 // allocate C
258 SLIP_CHECK (SLIP_matrix_allocate (&C, SLIP_CSC, C_type,
259 m, n, actual, false, true, option)) ;
260
261 // C's scaling factor is currently in Y. Set it here
262 SLIP_mpq_set(C->scale, Y->scale);
263
264 // Construct C
265 nz = 0 ;
266 switch (C->type)
267 {
268
269 case SLIP_MPZ:
270 for (int64_t j = 0 ; j < n ; j++)
271 {
272 C->p [j] = nz ;
273 for (int64_t i = 0 ; i < m ; i++)
274 {
275 SLIP_CHECK( SLIP_mpz_sgn( &s, Y->x.mpz[ i + j*A->m]));
276 if (s != 0)
277 {
278 C->i [nz] = i ;
279 SLIP_CHECK( SLIP_mpz_set ( SLIP_1D (C, nz, mpz),
280 Y->x.mpz[ i + j*A->m] ));
281 nz++ ;
282 }
283 }
284 }
285 break ;
286
287 case SLIP_MPQ:
288 for (int64_t j = 0 ; j < n ; j++)
289 {
290 C->p [j] = nz ;
291 for (int64_t i = 0 ; i < m ; i++)
292 {
293 SLIP_CHECK (SLIP_mpq_sgn (&s,
294 Y->x.mpq[ i + j*A->m])) ;
295 if (s != 0)
296 {
297 C->i [nz] = i ;
298 SLIP_CHECK(SLIP_mpq_set (
299 SLIP_1D(C, nz, mpq),
300 Y->x.mpq[ i + j*A->m]));
301 nz++ ;
302 }
303 }
304 }
305 break ;
306
307 case SLIP_MPFR:
308 for (int64_t j = 0 ; j < n ; j++)
309 {
310 C->p [j] = nz ;
311 for (int64_t i = 0 ; i < m ; i++)
312 {
313 SLIP_CHECK (SLIP_mpfr_sgn (&s,
314 Y->x.mpfr[i + j*A->m])) ;
315 if (s != 0)
316 {
317 C->i [nz] = i ;
318 SLIP_CHECK (SLIP_mpfr_set (
319 SLIP_1D (C, nz, mpfr),
320 Y->x.mpfr[i + j*A->m],
321 round)) ;
322
323 nz++ ;
324 }
325 }
326 }
327 break ;
328
329 case SLIP_INT64:
330 for (int64_t j = 0 ; j < n ; j++)
331 {
332 C->p [j] = nz ;
333 for (int64_t i = 0 ; i < m ; i++)
334 {
335 if ( Y->x.int64[i +j*A->m] != 0)
336 {
337 C->i [nz] = i ;
338 SLIP_1D (C, nz, int64) =
339 Y->x.int64[i +j*A->m] ;
340 nz++ ;
341 }
342 }
343 }
344 break ;
345
346 case SLIP_FP64:
347 for (int64_t j = 0 ; j < n ; j++)
348 {
349 C->p [j] = nz ;
350 for (int64_t i = 0 ; i < m ; i++)
351 {
352 if ( Y->x.fp64[i +j*A->m] != 0)
353 {
354 C->i [nz] = i ;
355 SLIP_1D (C, nz, fp64) =
356 Y->x.fp64[i +j*A->m];
357 nz++ ;
358 }
359 }
360 }
361 break ;
362 }
363 C->p [n] = nz ;
364 }
365 break ;
366
367 }
368
369 }
370 break ;
371
372 //----------------------------------------------------------------------
373 // C is triplet
374 //----------------------------------------------------------------------
375
376 case SLIP_TRIPLET:
377 {
378
379 switch (A->kind)
380 {
381
382 //--------------------------------------------------------------
383 // A is CSC, C is triplet
384 //--------------------------------------------------------------
385
386 case SLIP_CSC:
387 {
388 // allocate C
389 SLIP_CHECK (SLIP_matrix_allocate (&C, SLIP_TRIPLET, C_type,
390 m, n, nz, false, true, option)) ;
391 // copy and typecast A->x into C->x
392 SLIP_CHECK (slip_cast_array (SLIP_X (C), C->type,
393 SLIP_X (A), A->type, nz, C->scale, A->scale, option)) ;
394 // copy the row indices A->i into C->i
395 memcpy (C->i, A->i, nz * sizeof (int64_t)) ;
396 // construct C->j
397 for (int64_t j = 0 ; j < n ; j++)
398 {
399 for (int64_t p = A->p [j] ; p < A->p [j+1] ; p++)
400 {
401 C->j [p] = j ;
402 }
403 }
404 // set C->nz
405 C->nz = nz;
406 }
407 break ;
408
409 //--------------------------------------------------------------
410 // A is triplet, C is triplet
411 //--------------------------------------------------------------
412
413 case SLIP_TRIPLET:
414 {
415 // allocate C
416 SLIP_CHECK (SLIP_matrix_allocate (&C, SLIP_TRIPLET, C_type,
417 m, n, nz, false, true, option)) ;
418 // copy the pattern of A into C
419 memcpy (C->j, A->j, nz * sizeof (int64_t)) ;
420 memcpy (C->i, A->i, nz * sizeof (int64_t)) ;
421 // copy and typecast A->x into C->x
422 SLIP_CHECK (slip_cast_array (SLIP_X (C), C->type,
423 SLIP_X (A), A->type, nz, C->scale, A->scale, option)) ;
424 // set C->nz
425 C->nz = nz;
426 }
427 break ;
428
429 //--------------------------------------------------------------
430 // A is dense, C is triplet
431 //--------------------------------------------------------------
432
433 case SLIP_DENSE:
434 {
435 // convert A to a temporary CSC matrix
436 SLIP_CHECK (SLIP_matrix_copy (&T, SLIP_CSC, C_type,
437 A, option)) ;
438 // convert T from CSC to triplet
439 SLIP_CHECK (SLIP_matrix_copy (&C, SLIP_TRIPLET, C_type,
440 T, option)) ;
441 SLIP_matrix_free (&T, option) ;
442 // set C->nz
443 C->nz = nz;
444 }
445 break ;
446
447 }
448
449 }
450 break ;
451
452 //----------------------------------------------------------------------
453 // C is dense
454 //----------------------------------------------------------------------
455
456 case SLIP_DENSE:
457 {
458
459 // allocate C
460 SLIP_CHECK (SLIP_matrix_allocate (&C, SLIP_DENSE, C_type,
461 m, n, nz, false, true, option)) ;
462
463 switch (A->kind)
464 {
465
466 //--------------------------------------------------------------
467 // A is CSC, C is dense
468 //--------------------------------------------------------------
469
470 case SLIP_CSC:
471 {
472 // Y = typecast the values of A into the type of C
473 SLIP_CHECK (slip_cast_matrix (&Y, C->type, A, option)) ;
474
475 // Set C's scaling factor
476 SLIP_mpq_set(C->scale, Y->scale);
477
478 switch (C->type)
479 {
480
481 case SLIP_MPZ:
482 for (int64_t j = 0 ; j < n ; j++)
483 {
484 for (int64_t p = A->p [j] ; p < A->p [j+1] ;p++)
485 {
486 int64_t i = A->i [p] ;
487 SLIP_CHECK (SLIP_mpz_set (
488 SLIP_2D (C, i, j, mpz),
489 SLIP_1D (Y, p, mpz))) ;
490 }
491 }
492 break ;
493
494 case SLIP_MPQ:
495 for (int64_t j = 0 ; j < n ; j++)
496 {
497 for (int64_t p = A->p [j] ; p < A->p [j+1] ;p++)
498 {
499 int64_t i = A->i [p] ;
500 SLIP_CHECK (SLIP_mpq_set (
501 SLIP_2D (C, i, j, mpq),
502 SLIP_1D (Y, p, mpq))) ;
503 }
504 }
505 break ;
506
507 case SLIP_MPFR:
508 for (int64_t j = 0 ; j < n ; j++)
509 {
510 for (int64_t p = A->p [j] ; p < A->p [j+1] ;p++)
511 {
512 int64_t i = A->i [p] ;
513 SLIP_CHECK (SLIP_mpfr_set (
514 SLIP_2D (C, i, j, mpfr),
515 SLIP_1D (Y, p, mpfr),
516 round)) ;
517 }
518 }
519 break ;
520
521 case SLIP_INT64:
522 for (int64_t j = 0 ; j < n ; j++)
523 {
524 for (int64_t p = A->p [j] ; p < A->p [j+1] ;p++)
525 {
526 int64_t i = A->i [p] ;
527 SLIP_2D (C, i, j, int64) =
528 SLIP_1D (Y, p, int64) ;
529 }
530 }
531 break ;
532
533 case SLIP_FP64:
534 for (int64_t j = 0 ; j < n ; j++)
535 {
536 for (int64_t p = A->p [j] ; p < A->p [j+1] ;p++)
537 {
538 int64_t i = A->i [p] ;
539 SLIP_2D (C, i, j, fp64) =
540 SLIP_1D (Y, p, fp64) ;
541 }
542 }
543 break ;
544
545 }
546
547 }
548 break ;
549
550 //--------------------------------------------------------------
551 // A is triplet, C is dense
552 //--------------------------------------------------------------
553
554 case SLIP_TRIPLET:
555 {
556 // Y = typecast the values of A into the type of C
557 SLIP_CHECK (slip_cast_matrix (&Y, C->type, A, option)) ;
558
559 // Set C's scaling factor
560 SLIP_mpq_set(C->scale, Y->scale);
561 switch (C->type)
562 {
563
564 case SLIP_MPZ:
565 for (int64_t k = 0 ; k < nz ; k++)
566 {
567 int64_t i = A->i [k] ;
568 int64_t j = A->j [k] ;
569 SLIP_CHECK (SLIP_mpz_set (
570 SLIP_2D (C, i, j, mpz),
571 SLIP_1D (Y, k, mpz))) ;
572 }
573 break ;
574
575 case SLIP_MPQ:
576 for (int64_t k = 0 ; k < nz ; k++)
577 {
578 int64_t i = A->i [k] ;
579 int64_t j = A->j [k] ;
580 SLIP_CHECK (SLIP_mpq_set (
581 SLIP_2D (C, i, j, mpq),
582 SLIP_1D (Y, k, mpq))) ;
583 }
584 break ;
585
586 case SLIP_MPFR:
587 for (int64_t k = 0 ; k < nz ; k++)
588 {
589 int64_t i = A->i [k] ;
590 int64_t j = A->j [k] ;
591 SLIP_CHECK (SLIP_mpfr_set (
592 SLIP_2D (C, i, j, mpfr),
593 SLIP_1D (Y, k, mpfr),
594 round)) ;
595 }
596 break ;
597
598 case SLIP_INT64:
599 for (int64_t k = 0 ; k < nz ; k++)
600 {
601 int64_t i = A->i [k] ;
602 int64_t j = A->j [k] ;
603 SLIP_2D (C, i, j, int64) =
604 SLIP_1D (Y, k, int64) ;
605 }
606 break ;
607
608 case SLIP_FP64:
609 for (int64_t k = 0 ; k < nz ; k++)
610 {
611 int64_t i = A->i [k] ;
612 int64_t j = A->j [k] ;
613 SLIP_2D (C, i, j, fp64) =
614 SLIP_1D (Y, k, fp64) ;
615 }
616 break ;
617
618 }
619 }
620 break ;
621
622 //--------------------------------------------------------------
623 // A is dense, C is dense
624 //--------------------------------------------------------------
625
626 case SLIP_DENSE:
627 {
628 // copy and typecast A->x into C->x
629 SLIP_CHECK (slip_cast_array (SLIP_X (C), C->type,
630 SLIP_X (A), A->type, nz, C->scale, A->scale, option)) ;
631 }
632 break ;
633
634 }
635
636 }
637 break ;
638
639 }
640
641 //--------------------------------------------------------------------------
642 // free workspace and return result
643 //--------------------------------------------------------------------------
644
645 SLIP_FREE_WORK ;
646 (*C_handle) = C ;
647
648 return (SLIP_OK) ;
649 }
650
651