1 //------------------------------------------------------------------------------ 2 // GB_matrix.h: definitions for GrB_Matrix and GrB_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 // The GrB_Matrix and GrB_Vector objects are different names for the same 11 // content. A GrB_Vector is held as an m-by-1 non-hypersparse CSC matrix. 12 // This file is #include'd in GB.h to define the GB_Matrix_opaque and 13 // GB_Vector_opaque structs. It would be cleaner to define just one opaque 14 // struct, and then GrB_Matrix and GrB_Vector would be typedef'd as pointers to 15 // the same struct, but then the compiler gets confused with Generic(x). 16 17 // For a GrB_Vector object, as an m-by-1 non-hypersparse CSC matrix: 18 // bool is_csc ; // always true 19 // int64_t plen ; // always 1, so A->p always has length 2, and 20 // // contains [0 k] if the vector has k entries; 21 // // A->p is NULL if the GrB_Vector is bitmap. 22 // int64_t vdim ; // always 1 23 // int64_t nvec ; // always 1 24 // int64_t *h ; // always NULL 25 26 //------------------------------------------------------------------------------ 27 // basic information: magic, error logger, and type 28 //------------------------------------------------------------------------------ 29 30 // The first four items exactly match the first four items in the 31 // GrB_Descriptor struct. 32 33 int64_t magic ; // for detecting uninitialized objects 34 size_t header_size ; // size of the malloc'd block for this struct, or 0 35 char *logger ; // error logger string 36 size_t logger_size ; // size of the malloc'd block for logger, or 0 37 38 // The remaining items are specific the GrB_Matrix, GrB_Vector and GxB_Scalar 39 // structs, and do not appear in the GrB_Descriptor struct: 40 GrB_Type type ; // the type of each numerical entry 41 42 //------------------------------------------------------------------------------ 43 // compressed sparse vector data structure 44 //------------------------------------------------------------------------------ 45 46 // The matrix can be held in one of 8 formats, each one consisting of a set of 47 // vectors. The vector "names" are in the range 0 to A->vdim-1. Each 48 // vector has length A->vlen. These two values define the dimension of the 49 // matrix, where A is m-by-n. The m and n dimenions are vlen and vdim for the 50 // CSC formats, and reversed for the CSR formats. 51 52 // Ap, Ai, Ax, Ah, and Ab are abbreviations for A->p, A->i, A->x, A->h, and 53 // A->b, respectively. 54 55 // For the sparse and hypersparse formats, Ap is an integer array of size 56 // A->plen+1, with Ap [0] always zero. The matrix contains A->nvec sparse 57 // vectors, where A->nvec <= A->plen <= A->vdim. The arrays Ai and Ax are both 58 // of size A->nzmax, and define the indices and values in each sparse vector. 59 // The total number of entries in the matrix is Ap [nvec] <= A->nzmax. 60 // For the bitmap and full sparsity structures, Ap and Ai are NULL. 61 62 // For both hypersparse and non-hypersparse matrices, if A->nvec_nonempty is 63 // computed, it is the number of vectors that contain at least one entry, where 64 // 0 <= A->nvec_nonempty <= A->nvec always holds. If not computed, 65 // A->nvec_nonempty is equal to -1. 66 67 //------------------------------------------------------------------------------ 68 // The 8 formats: (hypersparse, sparse, bitmap, full) x (CSR or CSC) 69 //------------------------------------------------------------------------------ 70 71 // -------------------------------------- 72 // Full structure: 73 // -------------------------------------- 74 75 // Ah, Ap, Ai, and Ab are all NULL. 76 // A->nvec == A->vdim. A->plen is not needed (set to -1) 77 78 // -------------------------------------- 79 // A->is_csc is true: full CSC format 80 // -------------------------------------- 81 82 // A is m-by-n: where A->vdim = n, and A->vlen = m 83 84 // Column A(:,j) is held in Ax [p1:p2-1] where p1 = k*m, p2 = (k+1)*m. 85 // A(i,j) at position p has row index i = p%m and value Ax [p] 86 87 // -------------------------------------- 88 // A->is_csc is false: full CSR format 89 // -------------------------------------- 90 91 // A is m-by-n: where A->vdim = m, and A->vlen = n 92 93 // Row A(i,:) is held in Ax [p1:p2-1] where p1 = k*n, p2 = (k+1)*n. 94 // A(i,j) at position p has column index j = p%n and value Ax [p] 95 96 // -------------------------------------- 97 // Bitmap structure: 98 // -------------------------------------- 99 100 // Ah, Ap, and Ai are NULL. Ab is an int8_t array of size m*n. 101 // A->nvec == A->vdim. A->plen is not needed (set to -1) 102 103 // The bitmap structure is identical to the full structure, except for the 104 // addition of the bitmap array A->b. 105 106 // -------------------------------------- 107 // A->is_csc is true: bitmap CSC format 108 // -------------------------------------- 109 110 // A is m-by-n: where A->vdim = n, and A->vlen = m 111 112 // Column A(:,j) is held in Ax [p1:p2-1] where p1 = k*m, p2 = (k+1)*m. 113 // A(i,j) at position p has row index i = p%m and value Ax [p]. 114 // The entry A(i,j) is present if Ab [p] == 1, and not present if 115 // Ab [p] == 0. 116 117 // -------------------------------------- 118 // A->is_csc is false: bitmap CSR format 119 // -------------------------------------- 120 121 // A is m-by-n: where A->vdim = m, and A->vlen = n 122 123 // Row A(i,:) is held in Ax [p1:p2-1] where p1 = k*n, p2 = (k+1)*n. 124 // A(i,j) at position p has column index j = p%n and value Ax [p] 125 // The entry A(i,j) is present if Ab [p] == 1, and not present if 126 // Ab [p] == 0. 127 128 // -------------------------------------- 129 // Sparse structure: 130 // -------------------------------------- 131 132 // Ah and Ab are NULL 133 // A->nvec == A->plen == A->vdim 134 135 // -------------------------------------- 136 // A->is_csc is true: sparse CSC format 137 // -------------------------------------- 138 139 // Ap, Ai, and Ax store a sparse matrix in the a very similar style 140 // as MATLAB and CSparse, as a collection of sparse column vectors. 141 142 // Column A(:,j) is held in two parts: the row indices are in 143 // Ai [Ap [j]...Ap [j+1]-1], and the numerical values are in the 144 // same positions in Ax. 145 146 // A is m-by-n: where A->vdim = n, and A->vlen = m 147 148 // -------------------------------------- 149 // A->is_csc is false: sparse CSR format 150 // -------------------------------------- 151 152 // Ap, Ai, and Ax store a sparse matrix in CSR format, as a collection 153 // of sparse row vectors. 154 155 // Row A(i,:) is held in two parts: the column indices are in 156 // Ai [Ap [i]...Ap [i+1]-1], and the numerical values are in the 157 // same positions in Ax. 158 159 // A is m-by-n: where A->vdim = m, and A->vlen = n 160 161 // -------------------------------------- 162 // Hypersparse structure: 163 // -------------------------------------- 164 165 // Ab is NULL 166 // Ah is non-NULL and has size A->plen; it is always kept sorted, 167 // A->nvec <= A->plen <= A->vdim 168 169 // -------------------------------------- 170 // A->is_csc is true: hypersparse CSC format 171 // -------------------------------------- 172 173 // A is held as a set of A->nvec sparse column vectors, but not all 174 // columns 0 to n-1 are present. 175 176 // If column A(:,j) has any entries, then j = Ah [k] for some 177 // k in the range 0 to A->nvec-1. 178 179 // Column A(:,j) is held in two parts: the row indices are in Ai [Ap 180 // [k]...Ap [k+1]-1], and the numerical values are in the same 181 // positions in Ax. 182 183 // A is m-by-n: where A->vdim = n, and A->vlen = m 184 185 // -------------------------------------- 186 // A->is_csc is false: hypersparse CSR format 187 // -------------------------------------- 188 189 // A is held as a set of A->nvec sparse row vectors, but not all 190 // row 0 to m-1 are present. 191 192 // If row A(i,:) has any entries, then i = Ah [k] for some 193 // k in the range 0 to A->nvec-1. 194 195 // Row A(i,:) is held in two parts: the column indices are in Ai 196 // [Ap [k]...Ap [k+1]-1], and the numerical values are in the same 197 // positions in Ax. 198 199 // A is m-by-n: where A->vdim = n, and A->vlen = m 200 201 //------------------------------------------------------------------------------ 202 // GraphBLAS vs MATLAB vs CSparse 203 //------------------------------------------------------------------------------ 204 205 // Like MATLAB, the indices in a completed GraphBLAS matrix (as implemented 206 // here) are always kept sorted. If all vectors in a matrix have row indices 207 // in strictly ascending order, the matrix is called "unjumbled" in this code. 208 // A matrix with one or more unsorted vectors is "jumbled". 209 210 // GraphBLAS allows for pending operations, in a matrix with pending work. 211 // Pending work includes one or more of the following (1) the presence of 212 // zombies, (2) pending tuples, and (3) the matrix is jumbled. 213 214 // Unlike MATLAB, explicit zeros are never dropped in a GraphBLAS matrix. They 215 // cannot be since the semiring "zero" might be something else, like -Infinity 216 // for a max-plus semiring. However, dropping zeros is a minor nuance in the 217 // data structure. 218 219 // Like GraphBLAS, CSparse also keeps explicit zeros. CSparse allows its 220 // matrices to be jumbled at any time, and this is not considered an unfinished 221 // GraphBLAS matrix. 222 223 // Finally, MATLAB only allows for boolean ("logical" class) and double 224 // precision sparse matrices (complex and real). CSparse only supports double. 225 // By contrast, GraphBLAS supports any type, including types defined at run 226 // time by the user application. In the GraphBLAS code, the term "nonzero" is 227 // sometimes used in the comments, but this is short-hand for the phrase "an 228 // entry A(i,j) whose value is explicity held in the matrix and which appears 229 // in the pattern; its value can be anything". Entries not in the pattern are 230 // simply "not there"; see for example GrB_*_extractElement. The actual 231 // numerical value of these implicit entries is dependent upon the identity 232 // value of the semiring's monoid operation used on the matrix. The actual 233 // semiring is not held in the matrix itself, and there are no restrictions on 234 // using a matrix in multiple semirings. 235 236 //------------------------------------------------------------------------------ 237 // primary matrix content 238 //------------------------------------------------------------------------------ 239 240 int64_t plen ; // size of A->p and A->h 241 int64_t vlen ; // length of each sparse vector 242 int64_t vdim ; // number of vectors in the matrix 243 int64_t nvec ; // number of non-empty vectors for hypersparse form, 244 // or same as vdim otherwise. nvec <= plen. 245 // some of these vectors in Ah may actually be empty. 246 int64_t nvec_nonempty ; // the actual number of non-empty vectors, or -1 if 247 // not known 248 249 int64_t *h ; // list of non-empty vectors: h_size >= 8*max(plen,1) 250 int64_t *p ; // pointers: p_size >= 8*(plen+1) 251 int64_t *i ; // indices: i_size >= 8*max(nzmax,1) 252 void *x ; // values: x_size >= max(nzmax*A->type->size,1) 253 int8_t *b ; // bitmap: b_size >= max(nzmax,1) 254 int64_t nzmax ; // max possible # of entries 255 int64_t nvals ; // nvals(A) if A is bitmap 256 257 size_t p_size ; // exact size of A->p in bytes, zero if A->p is NULL 258 size_t h_size ; // exact size of A->h in bytes, zero if A->h is NULL 259 size_t b_size ; // exact size of A->b in bytes, zero if A->b is NULL 260 size_t i_size ; // exact size of A->i in bytes, zero if A->i is NULL 261 size_t x_size ; // exact size of A->x in bytes, zero if A->x is NULL 262 263 //------------------------------------------------------------------------------ 264 // pending tuples 265 //------------------------------------------------------------------------------ 266 267 // The list of pending tuples is a feature that does not appear in MATLAB or 268 // CSparse, although something like it appears in CHOLMOD as the "unpacked" 269 // matrix format, which allows the pattern of a matrix to be modified when 270 // updating/downdating a Cholesky factorization. 271 272 // If an entry A(i,j) does not appear in the data structure, assigning A(i,j)=x 273 // requires all entries in vectors j to the end of matrix to be shifted down by 274 // one, taking up to O(nnz(A)) time to do so. This very slow, and it is why in 275 // both MATLAB and CSparse, the recommendation is to create a list of tuples, 276 // and to build a sparse matrix all at once. This is done by the MATLAB 277 // "sparse" function, the CSparse "cs_compress", and the GraphBLAS 278 // GrB_Matrix_build and GrB_Vector_build functions. 279 280 // MATLAB does not have a "non-blocking" mode of operation, so A(i,j)=x can be 281 // very slow for a single scalar x. With GraphBLAS' non-blocking mode, tuples 282 // from GrB_setElement and GrB_*assign can be held in another format that is 283 // easy to update: a conventional list of tuples, held inside the matrix 284 // itself. A list of tuples is easy to update but hard to work with in most 285 // operations, so whenever another GraphBLAS method or operation needs to 286 // access the matrix, the matrix is "flattened" by applying all the pending 287 // tuples. 288 289 // When a new entry is added that does not exist in the matrix, it is added to 290 // this list of pending tuples. Only when the matrix is needed in another 291 // operation are the pending tuples assembled into the compressed sparse vector 292 // form, A->h, A->p, A->i, and A->x. 293 294 // The type of the list of pending tuples (Pending->type) need not be the same 295 // as the type of the matrix. The GrB_assign and GxB_subassign operations can 296 // leave pending tuples in the matrix. The accum operator, if not NULL, 297 // becomes the pending operator for assembling the pending tuples and adding 298 // them to the matrix. For typecasting z=accum(x,y), the pending tuples are 299 // typecasted to the type of y. 300 // 301 // Let aij by the value of the pending tuple of a matrix C. There are up to 5 302 // different types to consider: Pending->type (the type of aij), ztype, xtype, 303 // ytype, and ctype = C->type, (the type of the matrix C with pending tuples). 304 // 305 // If this is the first update to C(i,j), or if there is no accum operator, 306 // for for GrB_setElement: 307 // 308 // aij of Pending->type 309 // cij = (ctype) aij 310 // 311 // For subsequent tuples with GrB_assign or GxB_subassign, when accum is 312 // present: 313 // 314 // y = (ytype) aij 315 // x = (xtype) cij 316 // z = accum (x,y), result of ztype 317 // cij = (ctype) z ; 318 // 319 // Since the pending tuple must be typecasted to either ctype or ytype, 320 // depending on where it appears, it must be stored in its original type. 321 322 GB_Pending Pending ; // list of pending tuples 323 324 //----------------------------------------------------------------------------- 325 // zombies 326 //----------------------------------------------------------------------------- 327 328 // A "zombie" is the opposite of a pending tuple. It is an entry A(i,j) that 329 // has been marked for deletion, but has not been deleted yet because it is more 330 // efficient to delete all zombies all at once, rather than one (or a few) at a 331 // time. An entry A(i,j) is marked as a zombie by 'flipping' its index via 332 // GB_FLIP(i). A flipped index is negative, and the actual index can be 333 // obtained by GB_UNFLIP(i). GB_FLIP(i) is a function that is its own inverse: 334 // GB_FLIP(GB_FLIP(x))=x for all x. 335 336 // Using zombies allows entries to be marked for deletion. Their index is 337 // still important, for two reasons: (1) the indices in each vector of the 338 // matrix are kept sorted to enable the use of binary search, (2) a zombie may 339 // be restored as a regular entry by a subsequent update, via setElement, 340 // subassign, or assign. In this case its index is unflipped and its value 341 // modified. Had the zombie not been there, the update would have to be placed 342 // in the pending tuple list. It is more efficient to keep the pending tuple 343 // lists as short as possible, so zombies are kept as long as possible to 344 // facilitate faster subsequent updates. 345 346 // Unlike pending tuples, no list of zombies is needed since they are already 347 // in the right place in the matrix. However, methods and operations in 348 // GraphBLAS that cannot tolerate zombies in their input matries can check the 349 // condition (A->nzombies > 0), and then delete all of them if they appear, via 350 // GB_Matrix_wait. 351 352 uint64_t nzombies ; // number of zombies marked for deletion 353 354 //------------------------------------------------------------------------------ 355 // sparsity control 356 //------------------------------------------------------------------------------ 357 358 // The hyper_switch determines how the matrix is converted between the 359 // hypersparse and non-hypersparse formats. Let n = A->vdim and let k be the 360 // actual number of non-empty vectors. If A is hypersparse, k can be less than 361 // A->nvec since the latter can include vectors that appear in A->h but are 362 // actually empty. 363 364 // If a matrix is currently hypersparse, it can be converted to non-hypersparse 365 // if the condition (n <= 1 || k > n*hyper_switch*2) holds. Otherwise, it 366 // stays hypersparse. Note that if n <= 1 the matrix is always stored as 367 // non-hypersparse. 368 369 // If currently non-hypersparse, it can be converted to hypersparse if the 370 // condition (n > 1 && k <= n*hyper_switch) holds. Otherwise, it stays 371 // non-hypersparse. Note that if n <= 1 the matrix remains non-hypersparse. 372 373 // The default value of hyper_switch is assigned to be GxB_HYPER_DEFAULT at 374 // startup by GrB_init, and can then be modified globally with 375 // GxB_Global_Option_set. All new matrices are created with the same 376 // hyper_switch. Once a particular matrix has been constructed, its 377 // hyper_switch can be modified from the default with GxB_Matrix_Option_set. 378 // GrB_Vectors are never stored as hypersparse. 379 380 // A new matrix created via GrB_Matrix_new starts with k=0 and is created in 381 // hypersparse form unless (n <= 1 || 0 > hyper_switch) holds, where 382 // hyper_switch is the global default value. GrB_Vectors are always 383 // non-hypersparse. 384 385 // To force a matrix to always stay non-hypersparse, use hyper_switch = -1 (or 386 // any negative number). To force a matrix to always stay hypersparse, use 387 // hyper_switch = 1 or more. For code readability, these values are also 388 // predefined for the user application as GxB_ALWAYS_HYPER and GxB_NEVER_HYPER. 389 390 // Summary for switching between formats: 391 392 // (1) by-row and by-column: there is no automatic switch between CSR/CSC. 393 // By default, all GrB_Matrices are held in CSR form, unless they are 394 // n-by-1 (then they are CSC). The GrB_vector is always CSC. 395 396 // (2) If A->sparsity is GxB_AUTO_SPARSITY (15), then the following rules are 397 // used to control the sparsity structure: 398 // 399 // (a) When a matrix is created, it is empty and starts as hypersparse, 400 // except that a GrB_Vector is never hypersparse. 401 // 402 // (b) A hypersparse matrix with k non-empty vectors and 403 // k > 2*n*A->hyper_switch is changed to sparse, and a sparse matrix 404 // with k <= 1*n*A->hyper_switch is changed to hypersparse. 405 // A->hyper_switch = (1/16) by default. See GB_convert*test. 406 // 407 // (c) A matrix with all entries present is converted to full (anz = 408 // GB_NNZ(A) = anz_dense = (A->vlen)*(A->vdim)). 409 // 410 // (d) A matrix with anz = GB_NNZ(A) entries and dimension A->vlen by 411 // A->vdim can have at most anz_dense = (A->vlen)*(A->vdim) entries. 412 // If A is sparse/hypersparse with anz > A->bitmap_switch * anz_dense, 413 // then it switches to bitmap. If A is bitmap and anz = 414 // (A->bitmap_switch / 2) * anz_dense, it switches to sparse. In 415 // between those two regions, the sparsity structure is unchanged. 416 417 float hyper_switch ; // controls conversion hyper to/from sparse 418 float bitmap_switch ; // controls conversion sparse to/from bitmap 419 int sparsity ; // controls sparsity structure: hypersparse, 420 // sparse, bitmap, or full, or any combination. 421 422 //------------------------------------------------------------------------------ 423 // shallow matrices 424 //------------------------------------------------------------------------------ 425 426 // Internal matrices in this implementation of GraphBLAS may have "shallow" 427 // components. These are pointers A->p, A->h, A->i, A->b, and A->x that point 428 // to the content of another matrix. Using shallow components speeds up 429 // computations and saves memory, but shallow matrices are never passed back to 430 // the user application. 431 432 // If the following are true, then the corresponding component of the 433 // object is a pointer into components of another object. They must not 434 // be freed when freeing this object. 435 436 bool p_shallow ; // true if p is a shallow copy 437 bool h_shallow ; // true if h is a shallow copy 438 bool b_shallow ; // true if b is a shallow copy 439 bool i_shallow ; // true if i is a shallow copy 440 bool x_shallow ; // true if x is a shallow copy 441 bool static_header ; // true if this struct is statically allocated 442 443 //------------------------------------------------------------------------------ 444 // other bool content 445 //------------------------------------------------------------------------------ 446 447 bool is_csc ; // true if stored by column, false if by row 448 bool jumbled ; // true if the matrix may be jumbled. bitmap and full 449 // matrices are never jumbled. 450 451 //------------------------------------------------------------------------------ 452 // iterating through a matrix 453 //------------------------------------------------------------------------------ 454 455 // The matrix can be held in 8 formats: (hypersparse, sparse, bitmap, full) x 456 // (CSR, CSC). The comments below assume A is in CSC format but the code works 457 // for both CSR and CSC. 458 459 #ifdef for_comments_only // only so vim will add color to the code below: 460 461 // for reference: 462 #define GBI(Ai,p,avlen) ((Ai == NULL) ? ((p) % (avlen)) : Ai [p]) 463 #define GBB(Ab,p) ((Ab == NULL) ? 1 : Ab [p]) 464 #define GBP(Ap,k,avlen) ((Ap == NULL) ? ((k) * (avlen)) : Ap [k]) 465 #define GBH(Ah,k) ((Ah == NULL) ? (k) : Ah [k]) 466 467 // A->vdim: the vector dimension of A (ncols(A)) 468 // A->nvec: # of vectors that appear in A. For the hypersparse case, 469 // these are the number of column indices in Ah [0..nvec-1], since 470 // A is CSC. For all cases, Ap [0...nvec] are the pointers. 471 472 //-------------------- 473 // (1) full // A->h, A->p, A->i, A->b are NULL, A->nvec == A->vdim 474 475 int64_t vlen = A->vlen ; 476 for (k = 0 ; k < A->nvec ; k++) 477 { 478 j = k ; 479 // operate on column A(:,j) 480 int64_t pA_start = k * vlen ; 481 int64_t pA_end = (k+1) * vlen ; 482 for (p = pA_start ; p < pA_end ; p++) 483 { 484 // A(i,j) has row i = (p % vlen), value aij = Ax [p] 485 } 486 } 487 488 //-------------------- 489 // (2) bitmap // A->h, A->p, A->i are NULL, A->nvec == A->vdim 490 491 int64_t vlen = A->vlen ; 492 for (k = 0 ; k < A->nvec ; k++) 493 { 494 j = k ; 495 // operate on column A(:,j) 496 int64_t pA_start = k * vlen ; 497 int64_t pA_end = (k+1) * vlen ; 498 for (p = pA_start ; p < pA_end ; p++) 499 { 500 if (Ab [p] != 0) 501 { 502 // A(i,j) has row i = (p % vlen), value aij = Ax [p] 503 } 504 else 505 { 506 // A(i,j) is not present 507 } 508 } 509 } 510 511 //-------------------- 512 // (3) sparse // A->h is NULL, A->nvec == A->vdim 513 514 for (k = 0 ; k < A->nvec ; k++) 515 { 516 j = k ; 517 // operate on column A(:,j) 518 for (p = Ap [k] ; p < Ap [k+1] ; p++) 519 { 520 // A(i,j) has row i = Ai [p], value aij = Ax [p] 521 } 522 } 523 524 //-------------------- 525 // (4) hypersparse // A->h is non-NULL, A->nvec <= A->dim 526 527 for (k = 0 ; k < A->nvec ; k++) 528 { 529 j = A->h [k] 530 // operate on column A(:,j) 531 for (p = Ap [k] ; p < Ap [k+1] ; p++) 532 { 533 // A(i,j) has row i = Ai [p], value aij = Ax [p] 534 } 535 } 536 537 //-------------------- 538 // generic: for any matrix 539 540 int64_t vlen = A->vlen ; 541 for (k = 0 ; k < A->nvec ; k++) 542 { 543 j = GBH (Ah, k) ; 544 // operate on column A(:,j) 545 int64_t pA_start = GBP (Ap, k, vlen) ; 546 int64_t pA_end = GBP (Ap, k+1, vlen) ; 547 for (p = pA_start ; p < pA_end ; p++) 548 { 549 // A(i,j) has row index i, value aij = Ax [p] 550 if (!GBB (Ab, p)) continue ; 551 int64_t i = GBI (Ai, p, vlen) ; 552 double aij = Ax [p] ; 553 } 554 } 555 556 #endif 557 558