//------------------------------------------------------------------------------ // SPEX_Utilities/SPEX_matrix_copy: create a copy of a matrix //------------------------------------------------------------------------------ // SPEX_Utilities: (c) 2019-2024, Christopher Lourenco, Jinhao Chen, // Lorena Mejia Domenzain, Erick Moreno-Centeno, and Timothy A. Davis. // All Rights Reserved. // SPDX-License-Identifier: GPL-2.0-or-later or LGPL-3.0-or-later //------------------------------------------------------------------------------ // SPEX_matrix_copy creates a SPEX_matrix C that is a modified copy of a // SPEX_matrix A. The new matrix C can have a different kind and type // than A. // The input matrix A is assumed to be valid. It can be checked first with // SPEX_matrix_check, if desired. If the input matrix A is not valid, results // are undefined. // SPEX supports 16 matrix formats: 15 of them are all combinations of // (CSC, triplet, dense) x (mpz, mpq, mpfr, int64, double). // This function can convert an input // matrix A in any of these 16 formats, into an output matrix C in any of the // 16 supported formats. #define SPEX_FREE_WORK \ SPEX_matrix_free(&T, option); \ SPEX_matrix_free(&Y, option); \ SPEX_FREE(W); #define SPEX_FREE_ALL \ SPEX_FREE_WORK; \ SPEX_matrix_free(&C, option); #include "spex_util_internal.h" SPEX_info SPEX_matrix_copy ( SPEX_matrix *C_handle, // matrix to create (never shallow) // inputs, not modified: SPEX_kind C_kind, // C->kind: CSC, triplet, dense SPEX_type C_type, // C->type: mpz_t, mpq_t, mpfr_t, int64_t, or double const SPEX_matrix A, // matrix to make a copy of (may be shallow) const SPEX_options option ) { //-------------------------------------------------------------------------- // check inputs //-------------------------------------------------------------------------- SPEX_info info ; if (!spex_initialized ( )) return (SPEX_PANIC); int64_t nz; SPEX_matrix C = NULL ; SPEX_matrix Y = NULL ; SPEX_matrix T = NULL ; int64_t *W = NULL ; SPEX_CHECK (SPEX_matrix_nnz (&nz, A, option)); ASSERT( nz >= 0); if (C_handle == NULL || nz < 0 || // checked in SPEX_matrix_nnz: //A == NULL || A->kind < SPEX_CSC A->type < SPEX_MPZ || A->type > SPEX_FP64 || C_kind < SPEX_CSC || C_kind > SPEX_DENSE || C_type < SPEX_MPZ || C_type > SPEX_FP64 ) { return (SPEX_INCORRECT_INPUT); } (*C_handle) = NULL ; int64_t m = A->m ; int64_t n = A->n ; mpfr_rnd_t round = SPEX_OPTION_ROUND (option); //-------------------------------------------------------------------------- // copy and convert A into C //-------------------------------------------------------------------------- switch (C_kind) { //---------------------------------------------------------------------- // C is CSC //---------------------------------------------------------------------- case SPEX_CSC: { switch (A->kind) { //-------------------------------------------------------------- // A is CSC, C is CSC //-------------------------------------------------------------- case SPEX_CSC: { // allocate C SPEX_CHECK (SPEX_matrix_allocate (&C, SPEX_CSC, C_type, m, n, nz, false, true, option)); // copy the pattern of A into C memcpy (C->p, A->p, (n+1) * sizeof (int64_t)); memcpy (C->i, A->i, nz * sizeof (int64_t)); // copy and typecast A->x into C->x SPEX_CHECK (spex_cast_array (SPEX_X (C), C->type, SPEX_X (A), A->type, nz, C->scale, A->scale, option)); } break ; //-------------------------------------------------------------- // A is triplet, C is CSC //-------------------------------------------------------------- case SPEX_TRIPLET: { // Y = typecast the values of A into the type of C // (not the pattern; Y is SPEX_DENSE) SPEX_CHECK (spex_cast_matrix (&Y, C_type, A, option)); // allocate workspace W = (int64_t *) SPEX_calloc (n, sizeof (int64_t)); if (W == NULL) { SPEX_FREE_ALL; return (SPEX_OUT_OF_MEMORY); } // allocate C SPEX_CHECK (SPEX_matrix_allocate (&C, SPEX_CSC, C_type, m, n, nz, false, true, option)); // Scaling factor of C is currently in Y, set it // here SPEX_mpq_set(C->scale, Y->scale); // count the # of entries in each column for (int64_t k = 0 ; k < nz ; k++) { W [A->j [k]]++ ; } // C->p = cumulative sum of W spex_cumsum (C->p, W, n); // build the matrix switch (C->type) { case SPEX_MPZ: for (int64_t k = 0 ; k < nz ; k++) { int64_t p = W [A->j [k]]++ ; C->i [p] = A->i [k] ; SPEX_MPZ_SET ( SPEX_1D (C, p, mpz), SPEX_1D (Y, k, mpz)); } break ; case SPEX_MPQ: for (int64_t k = 0 ; k < nz ; k++) { int64_t p = W [A->j [k]]++ ; C->i [p] = A->i [k] ; SPEX_MPQ_SET ( SPEX_1D (C, p, mpq), SPEX_1D (Y, k, mpq)); } break ; case SPEX_MPFR: for (int64_t k = 0 ; k < nz ; k++) { int64_t p = W [A->j [k]]++ ; C->i [p] = A->i [k] ; SPEX_MPFR_SET ( SPEX_1D (C, p, mpfr), SPEX_1D (Y, k, mpfr), round); } break ; case SPEX_INT64: for (int64_t k = 0 ; k < nz ; k++) { int64_t p = W [A->j [k]]++ ; C->i [p] = A->i [k] ; SPEX_1D (C, p, int64) = SPEX_1D (Y, k, int64); } break ; case SPEX_FP64: for (int64_t k = 0 ; k < nz ; k++) { int64_t p = W [A->j [k]]++ ; C->i [p] = A->i [k] ; SPEX_1D (C, p, fp64) = SPEX_1D (Y, k, fp64); } break ; } } break ; //-------------------------------------------------------------- // A is dense, C is CSC //-------------------------------------------------------------- case SPEX_DENSE: { // Y = typecast the values of A into the type of C SPEX_CHECK( spex_cast_matrix(&Y, C_type, A, option) ); int s ; // count the actual nonzeros in Y int64_t actual = 0 ; switch (Y->type) { case SPEX_MPZ: for (int64_t k = 0 ; k < nz ; k++) { SPEX_MPZ_SGN (&s, SPEX_1D (Y, k, mpz)); if (s != 0) actual++ ; } break ; case SPEX_MPQ: for (int64_t k = 0 ; k < nz ; k++) { SPEX_MPQ_SGN (&s, SPEX_1D (Y, k, mpq)); if (s != 0) actual++ ; } break ; case SPEX_MPFR: for (int64_t k = 0 ; k < nz ; k++) { SPEX_MPFR_SGN (&s, SPEX_1D (Y, k, mpfr)); if (s != 0) actual++ ; } break ; case SPEX_INT64: for (int64_t k = 0 ; k < nz ; k++) { if (SPEX_1D (Y, k, int64) != 0) actual++ ; } break ; case SPEX_FP64: for (int64_t k = 0 ; k < nz ; k++) { if (SPEX_1D (Y, k, fp64) != 0) actual++ ; } break ; } // allocate C SPEX_CHECK (SPEX_matrix_allocate (&C, SPEX_CSC, C_type, m, n, actual, false, true, option)); // C's scaling factor is currently in Y. Set it here SPEX_mpq_set(C->scale, Y->scale); // Construct C nz = 0 ; switch (C->type) { case SPEX_MPZ: for (int64_t j = 0 ; j < n ; j++) { C->p [j] = nz ; for (int64_t i = 0 ; i < m ; i++) { SPEX_MPZ_SGN(&s, Y->x.mpz[ i + j*A->m]); if (s != 0) { C->i [nz] = i ; SPEX_MPZ_SET ( SPEX_1D (C, nz, mpz), Y->x.mpz[ i + j*A->m] ); nz++ ; } } } break ; case SPEX_MPQ: for (int64_t j = 0 ; j < n ; j++) { C->p [j] = nz ; for (int64_t i = 0 ; i < m ; i++) { SPEX_MPQ_SGN (&s, Y->x.mpq[ i + j*A->m]); if (s != 0) { C->i [nz] = i ; SPEX_MPQ_SET ( SPEX_1D(C, nz, mpq), Y->x.mpq[ i + j*A->m]); nz++ ; } } } break ; case SPEX_MPFR: for (int64_t j = 0 ; j < n ; j++) { C->p [j] = nz ; for (int64_t i = 0 ; i < m ; i++) { SPEX_MPFR_SGN (&s, Y->x.mpfr[i + j*A->m]); if (s != 0) { C->i [nz] = i ; SPEX_MPFR_SET ( SPEX_1D (C, nz, mpfr), Y->x.mpfr[i + j*A->m], round); nz++ ; } } } break ; case SPEX_INT64: for (int64_t j = 0 ; j < n ; j++) { C->p [j] = nz ; for (int64_t i = 0 ; i < m ; i++) { if ( Y->x.int64[i +j*A->m] != 0) { C->i [nz] = i ; SPEX_1D (C, nz, int64) = Y->x.int64[i +j*A->m] ; nz++ ; } } } break ; case SPEX_FP64: for (int64_t j = 0 ; j < n ; j++) { C->p [j] = nz ; for (int64_t i = 0 ; i < m ; i++) { if ( Y->x.fp64[i +j*A->m] != 0) { C->i [nz] = i ; SPEX_1D (C, nz, fp64) = Y->x.fp64[i +j*A->m]; nz++ ; } } } break ; } C->p [n] = nz ; } break ; } } break ; //---------------------------------------------------------------------- // C is triplet //---------------------------------------------------------------------- case SPEX_TRIPLET: { switch (A->kind) { //-------------------------------------------------------------- // A is CSC, C is triplet //-------------------------------------------------------------- case SPEX_CSC: { // allocate C SPEX_CHECK (SPEX_matrix_allocate (&C, SPEX_TRIPLET, C_type, m, n, nz, false, true, option)); // copy and typecast A->x into C->x SPEX_CHECK (spex_cast_array (SPEX_X (C), C->type, SPEX_X (A), A->type, nz, C->scale, A->scale, option)); // copy the row indices A->i into C->i memcpy (C->i, A->i, nz * sizeof (int64_t)); // construct C->j for (int64_t j = 0 ; j < n ; j++) { for (int64_t p = A->p [j] ; p < A->p [j+1] ; p++) { C->j [p] = j ; } } // set C->nz C->nz = nz; } break ; //-------------------------------------------------------------- // A is triplet, C is triplet //-------------------------------------------------------------- case SPEX_TRIPLET: { // allocate C SPEX_CHECK (SPEX_matrix_allocate (&C, SPEX_TRIPLET, C_type, m, n, nz, false, true, option)); // copy the pattern of A into C memcpy (C->j, A->j, nz * sizeof (int64_t)); memcpy (C->i, A->i, nz * sizeof (int64_t)); // copy and typecast A->x into C->x SPEX_CHECK (spex_cast_array (SPEX_X (C), C->type, SPEX_X (A), A->type, nz, C->scale, A->scale, option)); // set C->nz C->nz = nz; } break ; //-------------------------------------------------------------- // A is dense, C is triplet //-------------------------------------------------------------- case SPEX_DENSE: { // convert A to a temporary CSC matrix SPEX_CHECK (SPEX_matrix_copy (&T, SPEX_CSC, C_type, A, option)); // convert T from CSC to triplet SPEX_CHECK (SPEX_matrix_copy (&C, SPEX_TRIPLET, C_type, T, option)); SPEX_matrix_free (&T, option); // set C->nz C->nz = nz; } break ; } } break ; //---------------------------------------------------------------------- // C is dense //---------------------------------------------------------------------- case SPEX_DENSE: { switch (A->kind) { //-------------------------------------------------------------- // A is CSC, C is dense //-------------------------------------------------------------- case SPEX_CSC: { // allocate C SPEX_CHECK (SPEX_matrix_allocate (&C, SPEX_DENSE, C_type, m, n, nz, false, true, option)); // Y = typecast the values of A into the type of C SPEX_CHECK (spex_cast_matrix (&Y, C->type, A, option)); // Set C's scaling factor SPEX_mpq_set(C->scale, Y->scale); switch (C->type) { case SPEX_MPZ: for (int64_t j = 0 ; j < n ; j++) { for (int64_t p = A->p [j] ; p < A->p [j+1] ;p++) { int64_t i = A->i [p] ; SPEX_MPZ_SET ( SPEX_2D (C, i, j, mpz), SPEX_1D (Y, p, mpz)); } } break ; case SPEX_MPQ: for (int64_t j = 0 ; j < n ; j++) { for (int64_t p = A->p [j] ; p < A->p [j+1] ;p++) { int64_t i = A->i [p] ; SPEX_MPQ_SET ( SPEX_2D (C, i, j, mpq), SPEX_1D (Y, p, mpq)); } } break ; case SPEX_MPFR: for (int64_t j = 0 ; j < n ; j++) { for (int64_t p = A->p [j] ; p < A->p [j+1] ;p++) { int64_t i = A->i [p] ; SPEX_MPFR_SET ( SPEX_2D (C, i, j, mpfr), SPEX_1D (Y, p, mpfr), round); } } break ; case SPEX_INT64: for (int64_t j = 0 ; j < n ; j++) { for (int64_t p = A->p [j] ; p < A->p [j+1] ;p++) { int64_t i = A->i [p] ; SPEX_2D (C, i, j, int64) = SPEX_1D (Y, p, int64); } } break ; case SPEX_FP64: for (int64_t j = 0 ; j < n ; j++) { for (int64_t p = A->p [j] ; p < A->p [j+1] ;p++) { int64_t i = A->i [p] ; SPEX_2D (C, i, j, fp64) = SPEX_1D (Y, p, fp64); } } break ; } } break ; //-------------------------------------------------------------- // A is triplet, C is dense //-------------------------------------------------------------- case SPEX_TRIPLET: { // allocate C SPEX_CHECK (SPEX_matrix_allocate (&C, SPEX_DENSE, C_type, m, n, nz, false, true, option)); // Y = typecast the values of A into the type of C SPEX_CHECK (spex_cast_matrix (&Y, C->type, A, option)); // Set C's scaling factor SPEX_mpq_set(C->scale, Y->scale); switch (C->type) { case SPEX_MPZ: for (int64_t k = 0 ; k < nz ; k++) { int64_t i = A->i [k] ; int64_t j = A->j [k] ; SPEX_MPZ_SET ( SPEX_2D (C, i, j, mpz), SPEX_1D (Y, k, mpz)); } break ; case SPEX_MPQ: for (int64_t k = 0 ; k < nz ; k++) { int64_t i = A->i [k] ; int64_t j = A->j [k] ; SPEX_MPQ_SET ( SPEX_2D (C, i, j, mpq), SPEX_1D (Y, k, mpq)); } break ; case SPEX_MPFR: for (int64_t k = 0 ; k < nz ; k++) { int64_t i = A->i [k] ; int64_t j = A->j [k] ; SPEX_MPFR_SET ( SPEX_2D (C, i, j, mpfr), SPEX_1D (Y, k, mpfr), round); } break ; case SPEX_INT64: for (int64_t k = 0 ; k < nz ; k++) { int64_t i = A->i [k] ; int64_t j = A->j [k] ; SPEX_2D (C, i, j, int64) = SPEX_1D (Y, k, int64); } break ; case SPEX_FP64: for (int64_t k = 0 ; k < nz ; k++) { int64_t i = A->i [k] ; int64_t j = A->j [k] ; SPEX_2D (C, i, j, fp64) = SPEX_1D (Y, k, fp64); } break ; } } break ; //-------------------------------------------------------------- // A is dense, C is dense //-------------------------------------------------------------- case SPEX_DENSE: { // allocate C SPEX_CHECK (SPEX_matrix_allocate (&C, SPEX_DENSE, C_type, m, n, nz, false, true, option)); // copy and typecast A->x into C->x SPEX_CHECK (spex_cast_array (SPEX_X (C), C->type, SPEX_X (A), A->type, nz, C->scale, A->scale, option)); } break ; } } break ; } //-------------------------------------------------------------------------- // free workspace and return result //-------------------------------------------------------------------------- SPEX_FREE_WORK; (*C_handle) = C; return (SPEX_OK); }