1 //------------------------------------------------------------------------------
2 // GB_Vector_diag: extract a diagonal from a matrix, as a vector
3 //------------------------------------------------------------------------------
4
5 // SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
6 // SPDX-License-Identifier: Apache-2.0
7
8 //------------------------------------------------------------------------------
9
10 #define GB_FREE_WORK \
11 GB_Matrix_free (&T) ;
12
13 #define GB_FREE_ALL \
14 GB_FREE_WORK ; \
15 GB_phbix_free (V) ;
16
17 #include "GB_diag.h"
18 #include "GB_select.h"
19
GB_Vector_diag(GrB_Matrix V,const GrB_Matrix A,int64_t k,GB_Context Context)20 GrB_Info GB_Vector_diag // extract a diagonal from a matrix, as a vector
21 (
22 GrB_Matrix V, // output vector (as an n-by-1 matrix)
23 const GrB_Matrix A, // input matrix
24 int64_t k,
25 GB_Context Context
26 )
27 {
28
29 //--------------------------------------------------------------------------
30 // check inputs
31 //--------------------------------------------------------------------------
32
33 GrB_Info info ;
34 ASSERT_MATRIX_OK (A, "A input for GB_Vector_diag", GB0) ;
35 ASSERT_MATRIX_OK (V, "V input for GB_Vector_diag", GB0) ;
36 ASSERT (GB_VECTOR_OK (V)) ; // V is a vector on input
37 ASSERT (!GB_aliased (A, V)) ; // A and V cannot be aliased
38 ASSERT (!GB_IS_HYPERSPARSE (V)) ; // vectors cannot be hypersparse
39
40 struct GB_Matrix_opaque T_header ;
41 GrB_Matrix T = GB_clear_static_header (&T_header) ;
42
43 GrB_Type atype = A->type ;
44 GrB_Type vtype = V->type ;
45 int64_t nrows = GB_NROWS (A) ;
46 int64_t ncols = GB_NCOLS (A) ;
47 int64_t n ;
48 if (k >= ncols || k <= -nrows)
49 {
50 // output vector V must have zero length
51 n = 0 ;
52 }
53 else if (k >= 0)
54 {
55 // if k is in range 0 to n-1, V must have length min (m,n-k)
56 n = GB_IMIN (nrows, ncols - k) ;
57 }
58 else
59 {
60 // if k is in range -1 to -m+1, V must have length min (m+k,n)
61 n = GB_IMIN (nrows + k, ncols) ;
62 }
63
64 if (n != V->vlen)
65 {
66 GB_ERROR (GrB_DIMENSION_MISMATCH,
67 "Input vector must have size " GBd "\n", n) ;
68 }
69
70 if (!GB_Type_compatible (atype, vtype))
71 {
72 GB_ERROR (GrB_DOMAIN_MISMATCH, "Input matrix of type [%s] "
73 "cannot be typecast to output of type [%s]\n",
74 atype->name, vtype->name) ;
75 }
76
77 //--------------------------------------------------------------------------
78 // finish any pending work in A and clear the output vector V
79 //--------------------------------------------------------------------------
80
81 GB_MATRIX_WAIT (A) ;
82 GB_phbix_free (V) ;
83
84 //--------------------------------------------------------------------------
85 // handle the CSR/CSC format of A
86 //--------------------------------------------------------------------------
87
88 bool csc = A->is_csc ;
89 if (!csc)
90 {
91 // The kth diagonal of a CSC matrix is the same as the (-k)th diagonal
92 // of the CSR format, so if A is CSR, negate the value of k. Then
93 // treat A as if it were CSC in the rest of this method.
94 k = -k ;
95 }
96
97 //--------------------------------------------------------------------------
98 // extract the kth diagonal of A into the temporary hypersparse matrix T
99 //--------------------------------------------------------------------------
100
101 GB_OK (GB_selector (T, GB_DIAG_opcode, NULL, false, A, k, NULL, Context)) ;
102 GB_OK (GB_convert_any_to_hyper (T, Context)) ;
103 GB_MATRIX_WAIT (T) ;
104 ASSERT_MATRIX_OK (T, "T = diag (A,k)", GB0) ;
105
106 //--------------------------------------------------------------------------
107 // transplant the pattern of T into the sparse vector V
108 //--------------------------------------------------------------------------
109
110 int64_t vnz = GB_NNZ (T) ;
111 float bitmap_switch = V->bitmap_switch ;
112 int sparsity_control = V->sparsity ;
113 bool static_header = V->static_header ;
114
115 GB_OK (GB_new (&V, static_header, // prior static or dynamic header
116 vtype, n, 1, GB_Ap_malloc, true, GxB_SPARSE,
117 GxB_NEVER_HYPER, 1, Context)) ;
118 V->sparsity = sparsity_control ;
119 V->bitmap_switch = bitmap_switch ;
120
121 V->p [0] = 0 ;
122 V->p [1] = vnz ;
123 if (k >= 0)
124 {
125 // transplant T->i into V->i
126 V->i = T->i ;
127 V->i_size = T->i_size ;
128 T->i = NULL ;
129 }
130 else
131 {
132 // transplant T->h into V->i
133 V->i = T->h ;
134 V->i_size = T->h_size ;
135 T->h = NULL ;
136 }
137
138 //--------------------------------------------------------------------------
139 // transplant or typecast the values from T to V
140 //--------------------------------------------------------------------------
141
142 GB_Type_code vcode = vtype->code ;
143 GB_Type_code acode = atype->code ;
144 if (vcode == acode)
145 {
146 // transplant T->x into V->x
147 V->x = T->x ;
148 V->x_size = T->x_size ;
149 T->x = NULL ;
150 }
151 else
152 {
153 // V->x = (vtype) T->x
154 GB_GET_NTHREADS_MAX (nthreads_max, chunk, Context) ;
155 int nthreads = GB_nthreads (vnz, chunk, nthreads_max) ;
156 size_t vsize = vtype->size ;
157 size_t asize = atype->size ;
158 V->x = GB_MALLOC (vnz * vsize, GB_void, &(V->x_size)) ;
159 if (V->x == NULL)
160 {
161 // out of memory
162 GB_FREE_ALL ;
163 return (GrB_OUT_OF_MEMORY) ;
164 }
165 GB_cast_array ((GB_void *) V->x, vcode, (GB_void *) T->x, acode,
166 NULL, asize, vnz, nthreads) ;
167 }
168
169 //--------------------------------------------------------------------------
170 // finalize the vector V
171 //--------------------------------------------------------------------------
172
173 V->nzmax = T->nzmax ;
174 V->jumbled = T->jumbled ;
175 V->nvec_nonempty = (vnz == 0) ? 0 : 1 ;
176 V->magic = GB_MAGIC ;
177
178 //--------------------------------------------------------------------------
179 // free workspace, conform V to its desired format, and return result
180 //--------------------------------------------------------------------------
181
182 GB_FREE_WORK ;
183 ASSERT_MATRIX_OK (V, "V before conform for GB_Vector_diag", GB0) ;
184 GB_OK (GB_conform (V, Context)) ;
185 ASSERT_MATRIX_OK (V, "V output for GB_Vector_diag", GB0) ;
186 return (GrB_SUCCESS) ;
187 }
188
189