1 //------------------------------------------------------------------------------
2 // GB_subassign_symbolic: S = C(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 #include "GB_subassign_methods.h"
11 #include "GB_subref.h"
12
13 #undef GB_FREE_ALL
14 #define GB_FREE_ALL GB_phbix_free (S) ;
15
GB_subassign_symbolic(GrB_Matrix S,const GrB_Matrix C,const GrB_Index * I,const int64_t ni,const GrB_Index * J,const int64_t nj,const bool S_must_not_be_jumbled,GB_Context Context)16 GrB_Info GB_subassign_symbolic
17 (
18 // output
19 GrB_Matrix S, // S = symbolic(C(I,J)), static header
20 // inputs, not modified:
21 const GrB_Matrix C, // matrix to extract the pattern of
22 const GrB_Index *I, // index list for S = C(I,J), or GrB_ALL, etc.
23 const int64_t ni, // length of I, or special
24 const GrB_Index *J, // index list for S = C(I,J), or GrB_ALL, etc.
25 const int64_t nj, // length of J, or special
26 const bool S_must_not_be_jumbled,
27 GB_Context Context
28 )
29 {
30
31 //--------------------------------------------------------------------------
32 // check inputs
33 //--------------------------------------------------------------------------
34
35 GrB_Info info ;
36 ASSERT (!GB_IS_BITMAP (C)) ; // the caller cannot tolerate C bitmap
37 ASSERT (S != NULL && S->static_header) ;
38
39 //--------------------------------------------------------------------------
40 // extract the pattern: S = C(I,J) for S_Extraction method, and quick mask
41 //--------------------------------------------------------------------------
42
43 // S is a sparse int64_t matrix. Its "values" are not numerical, but
44 // indices into C. For example, suppose 100 = I [5] and 200 = J [7]. Then
45 // S(5,7) is the entry C(I(5),J(7)), and the value of S(5,7) is the
46 // position in C that holds that particular entry C(100,200):
47 // pC = S->x [...] gives the location of the value C->x [pC] and row index
48 // 100 = C->i [pC], and pC will be between C->p [200] ... C->p [200+1]-1
49 // if C is non-hypersparse. If C is hyperparse then pC will be still
50 // reside inside the vector jC, in the range C->p [k] ... C->p [k+1]-1,
51 // if jC is the kth non-empty vector in the hyperlist of C.
52
53 //--------------------------------------------------------------------------
54 // extract symbolic structure S=C(I,J)
55 //--------------------------------------------------------------------------
56
57 // FUTURE::: if whole_C_matrix is true, then C(:,:) = ... and S == C,
58 // except that S is zombie-free, read-only; and C collects zombies.
59
60 // FUTURE:: the properties of I and J are already known, and thus do
61 // not need to be recomputed by GB_subref.
62
63 // S and C have the same CSR/CSC format. S can be jumbled. It is in
64 // in the same hypersparse form as C (unless S is empty, in which case
65 // it is always returned as hypersparse). This also checks I and J.
66 GB_OK (GB_subref (S, C->is_csc, C, I, ni, J, nj, true, Context)) ;
67 ASSERT (GB_JUMBLED_OK (S)) ; // GB_subref can return S as unsorted
68
69 //--------------------------------------------------------------------------
70 // sort S if requested
71 //--------------------------------------------------------------------------
72
73 if (S_must_not_be_jumbled)
74 {
75 GB_MATRIX_WAIT_IF_JUMBLED (S) ; // but the caller requires S sorted
76 ASSERT (!GB_JUMBLED (S)) ;
77 }
78
79 //--------------------------------------------------------------------------
80 // check the result of S=C(I,J)
81 //--------------------------------------------------------------------------
82
83 #ifdef GB_DEBUG
84 ASSERT_MATRIX_OK (C, "C for subref extraction", GB0) ;
85 ASSERT_MATRIX_OK (S, "S for subref extraction", GB0) ;
86
87 // since C is not bitmap, neither is S
88 ASSERT (!GB_IS_BITMAP (S)) ;
89
90 // GB_subref sorts its input matrix, so C is no longer jumbled
91 ASSERT (!GB_JUMBLED (C)) ;
92
93 // this body of code explains what S contains.
94 // S is nI-by-nJ where nI = length (I) and nJ = length (J)
95
96 int64_t nI, Icolon [3], nJ, Jcolon [3] ;
97 int Ikind, Jkind ;
98 GB_ijlength (I, ni, C->vlen, &nI, &Ikind, Icolon) ;
99 GB_ijlength (J, nj, C->vdim, &nJ, &Jkind, Jcolon) ;
100
101 // get S
102 const int64_t *restrict Sp = S->p ;
103 const int64_t *restrict Sh = S->h ;
104 const int64_t *restrict Si = S->i ;
105 const int64_t *restrict Sx = (int64_t *) S->x ;
106 // for each vector of S
107 for (int64_t k = 0 ; k < S->nvec ; k++)
108 {
109 // prepare to iterate over the entries of vector S(:,jnew)
110 int64_t jnew = GBH (Sh, k) ;
111 int64_t pS_start = GBP (Sp, k, S->vlen) ;
112 int64_t pS_end = GBP (Sp, k+1, S->vlen) ;
113 // S (inew,jnew) corresponds to C (iC, jC) ;
114 // jC = J [j] ; or J is a colon expression
115 int64_t jC = GB_ijlist (J, jnew, Jkind, Jcolon) ;
116 for (int64_t pS = pS_start ; pS < pS_end ; pS++)
117 {
118 // S (inew,jnew) is a pointer back into C (I(inew), J(jnew))
119 int64_t inew = GBI (Si, pS, S->vlen) ;
120 ASSERT (inew >= 0 && inew < nI) ;
121 // iC = I [iA] ; or I is a colon expression
122 int64_t iC = GB_ijlist (I, inew, Ikind, Icolon) ;
123 int64_t p = Sx [pS] ;
124 ASSERT (p >= 0 && p < GB_NNZ (C)) ;
125 int64_t pC_start, pC_end, pleft = 0, pright = C->nvec-1 ;
126 bool found = GB_lookup (C->h != NULL, C->h, C->p, C->vlen,
127 &pleft, pright, jC, &pC_start, &pC_end) ;
128 ASSERT (found) ;
129 // If iC == I [inew] and jC == J [jnew], (or the equivaleent
130 // for GB_ALL, GB_RANGE, GB_STRIDE) then A(inew,jnew) will be
131 // assigned to C(iC,jC), and p = S(inew,jnew) gives the pointer
132 // into C to where the entry (C(iC,jC) appears in C:
133 ASSERT (pC_start <= p && p < pC_end) ;
134 ASSERT (iC == GB_UNFLIP (GBI (C->i, p, C->vlen))) ;
135 }
136 }
137 #endif
138
139 //--------------------------------------------------------------------------
140 // return result
141 //--------------------------------------------------------------------------
142
143 return (GrB_SUCCESS) ;
144 }
145
146