1 //------------------------------------------------------------------------------
2 // GB_subref: C = A(I,J)
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 // C=A(I,J), either symbolic or numeric. In a symbolic extraction, Cx [p] is
11 // not the value of A(i,j), but its position in Ai,Ax. That is, pA = Cx [p]
12 // means that the entry at position p in C is the same as the entry in A at
13 // position pA. In this case, Cx has a type of int64_t.
14
15 // Numeric extraction:
16
17 // Sparse submatrix reference, C = A(I,J), extracting the values. This is
18 // an internal function called by GB_extract with symbolic==false, which
19 // does the work of the user-callable GrB_*_extract methods. It is also
20 // called by GB_assign to extract the submask. No pending tuples or
21 // zombies appear in A.
22
23 // Symbolic extraction:
24
25 // Sparse submatrix reference, C = A(I,J), extracting the pattern, not the
26 // values. For the symbolic case, this function is called only by
27 // GB_subassign_symbolic. Symbolic extraction creates a matrix C with the
28 // same pattern (C->p and C->i) as numeric extraction, but with different
29 // values, C->x. For numeric extracion if C(inew,jnew) = A(i,j), the
30 // value of A(i,j) is copied into C(i,j). For symbolic extraction, its
31 // *pointer* is copied into C(i,j). Suppose an entry A(i,j) is held in Ai
32 // [pa] and Ax [pa], and it appears in the output matrix C in Ci [pc] and
33 // Cx [pc]. Then the two methods differ as follows:
34
35 // this is the same:
36
37 // i = Ai [pa] ; // index i of entry A(i,j)
38
39 // aij = Ax [pa] ; // value of the entry A(i,j)
40
41 // Ci [pc] = inew ; // index inew of C(inew,jnew)
42
43 // this is different:
44
45 // Cx [pc] = aij ; // for numeric extraction
46
47 // Cx [pc] = pa ; // for symbolic extraction
48
49 // This function is called with symbolic==true by only by
50 // GB_subassign_symbolic, which uses it to extract the pattern of C(I,J),
51 // for the submatrix assignment C(I,J)=A. In this case, this function
52 // needs to deal with zombie entries. GB_subassign_symbolic uses this
53 // function on its C matrix, which is called A here because it is not
54 // modified here.
55
56 // Reading a zombie entry: A zombie entry A(i,j) has been marked by
57 // flipping its index. The value of a zombie is not important, just its
58 // presence in the pattern. All zombies have been flipped (i < 0), and
59 // all regular entries are not flipped (i >= 0). Zombies are entries that
60 // have been marked for deletion but have not been removed from the matrix
61 // yet, since it's more efficient to delete zombies all at once rather
62 // than one at a time.
63
64 // The symbolic case is zombie-agnostic, in the sense that it does not
65 // delete them. It treats them like regular entries. However, their
66 // normal index must be used, not their flipped indices. The output
67 // matrix C contains all unflipped indices, and its references to zombies
68 // and regular entries are identical. Zombies in A are dealt with later.
69 // They cannot be detected in the output C matrix, but they can be
70 // detected in A. Since pa = Cx [pc] holds the position of the entry in
71 // A, the entry is a zombie if Ai [pa] has been flipped.
72
73 #define GB_FREE_WORK \
74 { \
75 GB_FREE_WERK (&TaskList, TaskList_size) ; \
76 GB_FREE_WERK (&Ap_start, Ap_start_size) ; \
77 GB_FREE_WERK (&Ap_end, Ap_end_size) ; \
78 GB_FREE_WERK (&Mark, Mark_size) ; \
79 GB_FREE_WERK (&Inext, Inext_size) ; \
80 }
81
82 #define GB_FREE_ALL \
83 { \
84 GB_FREE (&Cp, Cp_size) ; \
85 GB_FREE (&Ch, Ch_size) ; \
86 GB_FREE_WORK ; \
87 }
88
89 #include "GB_subref.h"
90
91 GB_PUBLIC // accessed by the MATLAB tests in GraphBLAS/Test only
GB_subref(GrB_Matrix C,const bool C_is_csc,const GrB_Matrix A,const GrB_Index * I,const int64_t ni,const GrB_Index * J,const int64_t nj,const bool symbolic,GB_Context Context)92 GrB_Info GB_subref // C = A(I,J): either symbolic or numeric
93 (
94 // output
95 GrB_Matrix C, // output matrix, static header
96 // input, not modified
97 const bool C_is_csc, // requested format of C
98 const GrB_Matrix A,
99 const GrB_Index *I, // index list for C = A(I,J), or GrB_ALL, etc.
100 const int64_t ni, // length of I, or special
101 const GrB_Index *J, // index list for C = A(I,J), or GrB_ALL, etc.
102 const int64_t nj, // length of J, or special
103 const bool symbolic, // if true, construct C as symbolic
104 GB_Context Context
105 )
106 {
107
108 //--------------------------------------------------------------------------
109 // check inputs
110 //--------------------------------------------------------------------------
111
112 GrB_Info info ;
113 ASSERT (C != NULL && C->static_header) ;
114 ASSERT_MATRIX_OK (A, "A for C=A(I,J) subref", GB0) ;
115 ASSERT (GB_ZOMBIES_OK (A)) ;
116 ASSERT (GB_JUMBLED_OK (A)) ; // A is sorted, below, if jumbled on input
117 ASSERT (GB_PENDING_OK (A)) ;
118
119 //--------------------------------------------------------------------------
120 // handle bitmap and full cases
121 //--------------------------------------------------------------------------
122
123 if (GB_IS_BITMAP (A) || GB_IS_FULL (A))
124 {
125 // C is constructed with same sparsity as A (bitmap or full)
126 return (GB_bitmap_subref (C, C_is_csc, A, I, ni, J, nj, symbolic,
127 Context)) ;
128 }
129
130 //--------------------------------------------------------------------------
131 // initializations
132 //--------------------------------------------------------------------------
133
134 int64_t *Cp = NULL ; size_t Cp_size = 0 ;
135 int64_t *Ch = NULL ; size_t Ch_size = 0 ;
136 int64_t *Ap_start = NULL ; size_t Ap_start_size = 0 ;
137 int64_t *Ap_end = NULL ; size_t Ap_end_size = 0 ;
138 int64_t *Mark = NULL ; size_t Mark_size = 0 ;
139 int64_t *Inext = NULL ; size_t Inext_size = 0 ;
140 GB_task_struct *TaskList = NULL ; size_t TaskList_size = 0 ;
141
142 int64_t Cnvec = 0, nI = 0, nJ, Icolon [3], Cnvec_nonempty, ndupl ;
143 bool post_sort, need_qsort ;
144 int Ikind, ntasks, nthreads ;
145
146 //--------------------------------------------------------------------------
147 // ensure A is unjumbled
148 //--------------------------------------------------------------------------
149
150 // ensure input matrix is not jumbled. Zombies are OK.
151 GB_MATRIX_WAIT_IF_JUMBLED (A) ;
152
153 //--------------------------------------------------------------------------
154 // phase0: find vectors for C=A(I,J), and I,J properties
155 //--------------------------------------------------------------------------
156
157 GB_OK (GB_subref_phase0 (
158 // computed by phase0:
159 &Ch, &Ch_size, &Ap_start, &Ap_start_size, &Ap_end, &Ap_end_size,
160 &Cnvec, &need_qsort, &Ikind, &nI, Icolon, &nJ,
161 // original input:
162 A, I, ni, J, nj, Context)) ;
163
164 //--------------------------------------------------------------------------
165 // phase0b: split C=A(I,J) into tasks for phase1 and phase2
166 //--------------------------------------------------------------------------
167
168 // This phase also inverts I if needed.
169
170 GB_OK (GB_subref_slice (
171 // computed by phase0b:
172 &TaskList, &TaskList_size, &ntasks, &nthreads, &post_sort,
173 &Mark, &Mark_size, &Inext, &Inext_size, &ndupl,
174 // computed by phase0:
175 Ap_start, Ap_end, Cnvec, need_qsort, Ikind, nI, Icolon,
176 // original input:
177 A->vlen, GB_NNZ (A), I, Context)) ;
178
179 //--------------------------------------------------------------------------
180 // phase1: count the number of entries in each vector of C
181 //--------------------------------------------------------------------------
182
183 GB_OK (GB_subref_phase1 (
184 // computed by phase1:
185 &Cp, &Cp_size, &Cnvec_nonempty,
186 // computed by phase0b:
187 TaskList, ntasks, nthreads, Mark, Inext, ndupl,
188 // computed by phase0:
189 Ap_start, Ap_end, Cnvec, need_qsort, Ikind, nI, Icolon,
190 // original input:
191 A, I, symbolic, Context)) ;
192
193 //--------------------------------------------------------------------------
194 // phase2: compute the entries (indices and values) in each vector of C
195 //--------------------------------------------------------------------------
196
197 GB_OK (GB_subref_phase2 (
198 // computed by phase2:
199 C,
200 // from phase1:
201 &Cp, Cp_size, Cnvec_nonempty,
202 // from phase0b:
203 TaskList, ntasks, nthreads, post_sort, Mark, Inext, ndupl,
204 // from phase0:
205 &Ch, Ch_size, Ap_start, Ap_end, Cnvec, need_qsort,
206 Ikind, nI, Icolon, nJ,
207 // original input:
208 C_is_csc, A, I, symbolic, Context)) ;
209
210 // Cp and Ch have been imported into C->p and C->h, or freed if phase2
211 // fails. Either way, Cp and Ch are set to NULL so that they cannot be
212 // freed here (except by freeing C itself).
213
214 // free workspace
215 GB_FREE_WORK ;
216
217 //--------------------------------------------------------------------------
218 // return result
219 //--------------------------------------------------------------------------
220
221 // C can be returned jumbled, even if A is not jumbled
222 ASSERT_MATRIX_OK (C, "C output for C=A(I,J)", GB0) ;
223 ASSERT (GB_ZOMBIES_OK (C)) ;
224 ASSERT (GB_JUMBLED_OK (C)) ;
225 ASSERT (GB_IS_SPARSE (A) || GB_IS_HYPERSPARSE (A)) ;
226 return (GrB_SUCCESS) ;
227 }
228
229