1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228
|
//------------------------------------------------------------------------------
// GB_Matrix_diag: construct a diagonal matrix from a vector
//------------------------------------------------------------------------------
// SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2022, All Rights Reserved.
// SPDX-License-Identifier: Apache-2.0
//------------------------------------------------------------------------------
#define GB_FREE_WORKSPACE \
{ \
GB_Matrix_free (&T) ; \
}
#define GB_FREE_ALL \
{ \
GB_FREE_WORKSPACE ; \
GB_phybix_free (C) ; \
}
#include "GB_diag.h"
#include "GB_unused.h"
GrB_Info GB_Matrix_diag // build a diagonal matrix from a vector
(
GrB_Matrix C, // output matrix
const GrB_Matrix V_in, // input vector (as an n-by-1 matrix)
int64_t k,
GB_Context Context
)
{
//--------------------------------------------------------------------------
// check inputs
//--------------------------------------------------------------------------
GrB_Info info ;
ASSERT_MATRIX_OK (C, "C input for GB_Matrix_diag", GB0) ;
ASSERT_MATRIX_OK (V_in, "V input for GB_Matrix_diag", GB0) ;
ASSERT (GB_VECTOR_OK (V_in)) ; // V_in is a vector on input
ASSERT (!GB_aliased (C, V_in)) ; // C and V_in cannot be aliased
ASSERT (!GB_IS_HYPERSPARSE (V_in)) ; // vectors cannot be hypersparse
struct GB_Matrix_opaque T_header ;
GrB_Matrix T = NULL ;
GrB_Type ctype = C->type ;
GrB_Type vtype = V_in->type ;
int64_t n = V_in->vlen + GB_IABS (k) ; // C must be n-by-n
ASSERT (GB_NROWS (C) == GB_NCOLS (C))
ASSERT (GB_NROWS (C) == n)
ASSERT (GB_Type_compatible (ctype, vtype)) ;
//--------------------------------------------------------------------------
// finish any pending work in V_in and clear the output matrix C
//--------------------------------------------------------------------------
GB_MATRIX_WAIT (V_in) ;
GB_phybix_free (C) ;
//--------------------------------------------------------------------------
// ensure V is not bitmap
//--------------------------------------------------------------------------
GrB_Matrix V ;
if (GB_IS_BITMAP (V_in))
{
// make a deep copy of V_in and convert to CSC
// set T->iso = V_in->iso OK
GB_CLEAR_STATIC_HEADER (T, &T_header) ;
GB_OK (GB_dup_worker (&T, V_in->iso, V_in, true, NULL, Context)) ;
GB_OK (GB_convert_bitmap_to_sparse (T, Context)) ;
V = T ;
}
else
{
// use V_in as-is
V = V_in ;
}
//--------------------------------------------------------------------------
// allocate C as sparse or hypersparse with vnz entries and vnz vectors
//--------------------------------------------------------------------------
// C is sparse if V is dense and k == 0, and hypersparse otherwise
const int64_t vnz = GB_nnz (V) ;
const bool V_is_full = GB_is_dense (V) ;
const int C_sparsity = (V_is_full && k == 0) ? GxB_SPARSE : GxB_HYPERSPARSE;
const bool C_iso = V->iso ;
if (C_iso)
{
GBURBLE ("(iso diag) ") ;
}
const bool csc = C->is_csc ;
const float bitmap_switch = C->bitmap_switch ;
const int sparsity_control = C->sparsity_control ;
// set C->iso = C_iso OK
GB_OK (GB_new_bix (&C, // existing header
ctype, n, n, GB_Ap_malloc, csc, C_sparsity, false,
C->hyper_switch, vnz, vnz, true, C_iso, Context)) ;
C->sparsity_control = sparsity_control ;
C->bitmap_switch = bitmap_switch ;
//--------------------------------------------------------------------------
// handle the CSR/CSC format of C and determine position of diagonal
//--------------------------------------------------------------------------
if (!csc)
{
// The kth diagonal of a CSC matrix is the same as the (-k)th diagonal
// of the CSR format, so if C is CSR, negate the value of k. Then
// treat C as if it were CSC in the rest of this method.
k = -k ;
}
int64_t kpositive, knegative ;
if (k >= 0)
{
kpositive = k ;
knegative = 0 ;
}
else
{
kpositive = 0 ;
knegative = -k ;
}
//--------------------------------------------------------------------------
// get the contents of C and determine # of threads to use
//--------------------------------------------------------------------------
GB_GET_NTHREADS_MAX (nthreads_max, chunk, Context) ;
int nthreads = GB_nthreads (vnz, chunk, nthreads_max) ;
int64_t *restrict Cp = C->p ;
int64_t *restrict Ch = C->h ;
int64_t *restrict Ci = C->i ;
//--------------------------------------------------------------------------
// copy the contents of V into the kth diagonal of C
//--------------------------------------------------------------------------
if (C_sparsity == GxB_SPARSE)
{
//----------------------------------------------------------------------
// V is full, or can be treated as full, and k == 0
//----------------------------------------------------------------------
// C->x = (ctype) V->x
GB_cast_matrix (C, V, Context) ;
// construct Cp and Ci
int64_t p ;
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (p = 0 ; p < vnz ; p++)
{
Cp [p] = p ;
Ci [p] = p ;
}
}
else if (V_is_full)
{
//----------------------------------------------------------------------
// V is full, or can be treated as full, and k != 0
//----------------------------------------------------------------------
// C->x = (ctype) V->x
GB_cast_matrix (C, V, Context) ;
// construct Cp, Ch, and Ci
int64_t p ;
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (p = 0 ; p < vnz ; p++)
{
Cp [p] = p ;
Ch [p] = p + kpositive ;
Ci [p] = p + knegative ;
}
}
else
{
//----------------------------------------------------------------------
// V is sparse
//----------------------------------------------------------------------
// C->x = (ctype) V->x
GB_cast_matrix (C, V, Context) ;
int64_t *restrict Vi = V->i ;
// construct Cp, Ch, and Ci
int64_t p ;
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (p = 0 ; p < vnz ; p++)
{
Cp [p] = p ;
Ch [p] = Vi [p] + kpositive ;
Ci [p] = Vi [p] + knegative ;
}
}
//--------------------------------------------------------------------------
// finalize the matrix C
//--------------------------------------------------------------------------
Cp [vnz] = vnz ;
C->nvals = vnz ;
C->nvec = vnz ;
C->nvec_nonempty = vnz ;
C->magic = GB_MAGIC ;
//--------------------------------------------------------------------------
// free workspace, conform C to its desired format, and return result
//--------------------------------------------------------------------------
GB_FREE_WORKSPACE ;
ASSERT_MATRIX_OK (C, "C before conform for GB_Matrix_diag", GB0) ;
GB_OK (GB_conform (C, Context)) ;
ASSERT_MATRIX_OK (C, "C output for GB_Matrix_diag", GB0) ;
return (GrB_SUCCESS) ;
}
|