1 #include <../src/mat/impls/aij/mpi/mpiaij.h> /*I "petscmat.h" I*/ 2 #include <../src/mat/impls/sell/mpi/mpisell.h> /*I "petscmat.h" I*/ 3 #include <petsc/private/vecimpl.h> 4 #include <petsc/private/isimpl.h> 5 #include <petscblaslapack.h> 6 #include <petscsf.h> 7 8 /*MC 9 MATSELL - MATSELL = "sell" - A matrix type to be used for sparse matrices. 10 11 This matrix type is identical to `MATSEQSELL` when constructed with a single process communicator, 12 and `MATMPISELL` otherwise. As a result, for single process communicators, 13 `MatSeqSELLSetPreallocation()` is supported, and similarly `MatMPISELLSetPreallocation()` is supported 14 for communicators controlling multiple processes. It is recommended that you call both of 15 the above preallocation routines for simplicity. 16 17 Options Database Keys: 18 . -mat_type sell - sets the matrix type to `MATSELL` during a call to `MatSetFromOptions()` 19 20 Level: beginner 21 22 .seealso: `Mat`, `MATAIJ`, `MATBAIJ`, `MATSBAIJ`, `MatCreateSELL()`, `MatCreateSeqSELL()`, `MATSEQSELL`, `MATMPISELL` 23 M*/ 24 25 PetscErrorCode MatDiagonalSet_MPISELL(Mat Y, Vec D, InsertMode is) 26 { 27 Mat_MPISELL *sell = (Mat_MPISELL *)Y->data; 28 29 PetscFunctionBegin; 30 if (Y->assembled && Y->rmap->rstart == Y->cmap->rstart && Y->rmap->rend == Y->cmap->rend) { 31 PetscCall(MatDiagonalSet(sell->A, D, is)); 32 } else { 33 PetscCall(MatDiagonalSet_Default(Y, D, is)); 34 } 35 PetscFunctionReturn(PETSC_SUCCESS); 36 } 37 38 /* 39 Local utility routine that creates a mapping from the global column 40 number to the local number in the off-diagonal part of the local 41 storage of the matrix. When PETSC_USE_CTABLE is used this is scalable at 42 a slightly higher hash table cost; without it it is not scalable (each processor 43 has an order N integer array but is fast to access. 44 */ 45 PetscErrorCode MatCreateColmap_MPISELL_Private(Mat mat) 46 { 47 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 48 PetscInt n = sell->B->cmap->n, i; 49 50 PetscFunctionBegin; 51 PetscCheck(sell->garray, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MPISELL Matrix was assembled but is missing garray"); 52 #if defined(PETSC_USE_CTABLE) 53 PetscCall(PetscHMapICreateWithSize(n, &sell->colmap)); 54 for (i = 0; i < n; i++) PetscCall(PetscHMapISet(sell->colmap, sell->garray[i] + 1, i + 1)); 55 #else 56 PetscCall(PetscCalloc1(mat->cmap->N + 1, &sell->colmap)); 57 for (i = 0; i < n; i++) sell->colmap[sell->garray[i]] = i + 1; 58 #endif 59 PetscFunctionReturn(PETSC_SUCCESS); 60 } 61 62 #define MatSetValues_SeqSELL_A_Private(row, col, value, addv, orow, ocol) \ 63 { \ 64 if (col <= lastcol1) low1 = 0; \ 65 else high1 = nrow1; \ 66 lastcol1 = col; \ 67 while (high1 - low1 > 5) { \ 68 t = (low1 + high1) / 2; \ 69 if (cp1[sliceheight * t] > col) high1 = t; \ 70 else low1 = t; \ 71 } \ 72 for (_i = low1; _i < high1; _i++) { \ 73 if (cp1[sliceheight * _i] > col) break; \ 74 if (cp1[sliceheight * _i] == col) { \ 75 if (addv == ADD_VALUES) vp1[sliceheight * _i] += value; \ 76 else vp1[sliceheight * _i] = value; \ 77 inserted = PETSC_TRUE; \ 78 goto a_noinsert; \ 79 } \ 80 } \ 81 if (value == 0.0 && ignorezeroentries) { \ 82 low1 = 0; \ 83 high1 = nrow1; \ 84 goto a_noinsert; \ 85 } \ 86 if (nonew == 1) { \ 87 low1 = 0; \ 88 high1 = nrow1; \ 89 goto a_noinsert; \ 90 } \ 91 PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \ 92 MatSeqXSELLReallocateSELL(A, am, 1, nrow1, a->sliidx, a->sliceheight, row / sliceheight, row, col, a->colidx, a->val, cp1, vp1, nonew, MatScalar); \ 93 /* shift up all the later entries in this row */ \ 94 for (ii = nrow1 - 1; ii >= _i; ii--) { \ 95 cp1[sliceheight * (ii + 1)] = cp1[sliceheight * ii]; \ 96 vp1[sliceheight * (ii + 1)] = vp1[sliceheight * ii]; \ 97 } \ 98 cp1[sliceheight * _i] = col; \ 99 vp1[sliceheight * _i] = value; \ 100 a->nz++; \ 101 nrow1++; \ 102 A->nonzerostate++; \ 103 a_noinsert:; \ 104 a->rlen[row] = nrow1; \ 105 } 106 107 #define MatSetValues_SeqSELL_B_Private(row, col, value, addv, orow, ocol) \ 108 { \ 109 if (col <= lastcol2) low2 = 0; \ 110 else high2 = nrow2; \ 111 lastcol2 = col; \ 112 while (high2 - low2 > 5) { \ 113 t = (low2 + high2) / 2; \ 114 if (cp2[sliceheight * t] > col) high2 = t; \ 115 else low2 = t; \ 116 } \ 117 for (_i = low2; _i < high2; _i++) { \ 118 if (cp2[sliceheight * _i] > col) break; \ 119 if (cp2[sliceheight * _i] == col) { \ 120 if (addv == ADD_VALUES) vp2[sliceheight * _i] += value; \ 121 else vp2[sliceheight * _i] = value; \ 122 inserted = PETSC_TRUE; \ 123 goto b_noinsert; \ 124 } \ 125 } \ 126 if (value == 0.0 && ignorezeroentries) { \ 127 low2 = 0; \ 128 high2 = nrow2; \ 129 goto b_noinsert; \ 130 } \ 131 if (nonew == 1) { \ 132 low2 = 0; \ 133 high2 = nrow2; \ 134 goto b_noinsert; \ 135 } \ 136 PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \ 137 MatSeqXSELLReallocateSELL(B, bm, 1, nrow2, b->sliidx, b->sliceheight, row / sliceheight, row, col, b->colidx, b->val, cp2, vp2, nonew, MatScalar); \ 138 /* shift up all the later entries in this row */ \ 139 for (ii = nrow2 - 1; ii >= _i; ii--) { \ 140 cp2[sliceheight * (ii + 1)] = cp2[sliceheight * ii]; \ 141 vp2[sliceheight * (ii + 1)] = vp2[sliceheight * ii]; \ 142 } \ 143 cp2[sliceheight * _i] = col; \ 144 vp2[sliceheight * _i] = value; \ 145 b->nz++; \ 146 nrow2++; \ 147 B->nonzerostate++; \ 148 b_noinsert:; \ 149 b->rlen[row] = nrow2; \ 150 } 151 152 PetscErrorCode MatSetValues_MPISELL(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv) 153 { 154 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 155 PetscScalar value; 156 PetscInt i, j, rstart = mat->rmap->rstart, rend = mat->rmap->rend, shift1, shift2; 157 PetscInt cstart = mat->cmap->rstart, cend = mat->cmap->rend, row, col; 158 PetscBool roworiented = sell->roworiented; 159 160 /* Some Variables required in the macro */ 161 Mat A = sell->A; 162 Mat_SeqSELL *a = (Mat_SeqSELL *)A->data; 163 PetscBool ignorezeroentries = a->ignorezeroentries, found; 164 Mat B = sell->B; 165 Mat_SeqSELL *b = (Mat_SeqSELL *)B->data; 166 PetscInt *cp1, *cp2, ii, _i, nrow1, nrow2, low1, high1, low2, high2, t, lastcol1, lastcol2, sliceheight = a->sliceheight; 167 MatScalar *vp1, *vp2; 168 169 PetscFunctionBegin; 170 for (i = 0; i < m; i++) { 171 if (im[i] < 0) continue; 172 PetscCheck(im[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, im[i], mat->rmap->N - 1); 173 if (im[i] >= rstart && im[i] < rend) { 174 row = im[i] - rstart; 175 lastcol1 = -1; 176 shift1 = a->sliidx[row / sliceheight] + (row % sliceheight); /* starting index of the row */ 177 cp1 = a->colidx + shift1; 178 vp1 = a->val + shift1; 179 nrow1 = a->rlen[row]; 180 low1 = 0; 181 high1 = nrow1; 182 lastcol2 = -1; 183 shift2 = b->sliidx[row / sliceheight] + (row % sliceheight); /* starting index of the row */ 184 cp2 = b->colidx + shift2; 185 vp2 = b->val + shift2; 186 nrow2 = b->rlen[row]; 187 low2 = 0; 188 high2 = nrow2; 189 190 for (j = 0; j < n; j++) { 191 if (roworiented) value = v[i * n + j]; 192 else value = v[i + j * m]; 193 if (ignorezeroentries && value == 0.0 && (addv == ADD_VALUES)) continue; 194 if (in[j] >= cstart && in[j] < cend) { 195 col = in[j] - cstart; 196 MatSetValue_SeqSELL_Private(A, row, col, value, addv, im[i], in[j], cp1, vp1, lastcol1, low1, high1); /* set one value */ 197 #if defined(PETSC_HAVE_CUDA) 198 if (A->offloadmask != PETSC_OFFLOAD_UNALLOCATED && found) A->offloadmask = PETSC_OFFLOAD_CPU; 199 #endif 200 } else if (in[j] < 0) { 201 continue; 202 } else { 203 PetscCheck(in[j] < mat->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, in[j], mat->cmap->N - 1); 204 if (mat->was_assembled) { 205 if (!sell->colmap) PetscCall(MatCreateColmap_MPISELL_Private(mat)); 206 #if defined(PETSC_USE_CTABLE) 207 PetscCall(PetscHMapIGetWithDefault(sell->colmap, in[j] + 1, 0, &col)); 208 col--; 209 #else 210 col = sell->colmap[in[j]] - 1; 211 #endif 212 if (col < 0 && !((Mat_SeqSELL *)(sell->B->data))->nonew) { 213 PetscCall(MatDisAssemble_MPISELL(mat)); 214 col = in[j]; 215 /* Reinitialize the variables required by MatSetValues_SeqSELL_B_Private() */ 216 B = sell->B; 217 b = (Mat_SeqSELL *)B->data; 218 shift2 = b->sliidx[row / sliceheight] + (row % sliceheight); /* starting index of the row */ 219 cp2 = b->colidx + shift2; 220 vp2 = b->val + shift2; 221 nrow2 = b->rlen[row]; 222 low2 = 0; 223 high2 = nrow2; 224 found = PETSC_FALSE; 225 } else { 226 PetscCheck(col >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", im[i], in[j]); 227 } 228 } else col = in[j]; 229 MatSetValue_SeqSELL_Private(B, row, col, value, addv, im[i], in[j], cp2, vp2, lastcol2, low2, high2); /* set one value */ 230 #if defined(PETSC_HAVE_CUDA) 231 if (B->offloadmask != PETSC_OFFLOAD_UNALLOCATED && found) B->offloadmask = PETSC_OFFLOAD_CPU; 232 #endif 233 } 234 } 235 } else { 236 PetscCheck(!mat->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", im[i]); 237 if (!sell->donotstash) { 238 mat->assembled = PETSC_FALSE; 239 if (roworiented) { 240 PetscCall(MatStashValuesRow_Private(&mat->stash, im[i], n, in, v + i * n, (PetscBool)(ignorezeroentries && (addv == ADD_VALUES)))); 241 } else { 242 PetscCall(MatStashValuesCol_Private(&mat->stash, im[i], n, in, v + i, m, (PetscBool)(ignorezeroentries && (addv == ADD_VALUES)))); 243 } 244 } 245 } 246 } 247 PetscFunctionReturn(PETSC_SUCCESS); 248 } 249 250 PetscErrorCode MatGetValues_MPISELL(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[]) 251 { 252 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 253 PetscInt i, j, rstart = mat->rmap->rstart, rend = mat->rmap->rend; 254 PetscInt cstart = mat->cmap->rstart, cend = mat->cmap->rend, row, col; 255 256 PetscFunctionBegin; 257 for (i = 0; i < m; i++) { 258 if (idxm[i] < 0) continue; /* negative row */ 259 PetscCheck(idxm[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm[i], mat->rmap->N - 1); 260 if (idxm[i] >= rstart && idxm[i] < rend) { 261 row = idxm[i] - rstart; 262 for (j = 0; j < n; j++) { 263 if (idxn[j] < 0) continue; /* negative column */ 264 PetscCheck(idxn[j] < mat->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, idxn[j], mat->cmap->N - 1); 265 if (idxn[j] >= cstart && idxn[j] < cend) { 266 col = idxn[j] - cstart; 267 PetscCall(MatGetValues(sell->A, 1, &row, 1, &col, v + i * n + j)); 268 } else { 269 if (!sell->colmap) PetscCall(MatCreateColmap_MPISELL_Private(mat)); 270 #if defined(PETSC_USE_CTABLE) 271 PetscCall(PetscHMapIGetWithDefault(sell->colmap, idxn[j] + 1, 0, &col)); 272 col--; 273 #else 274 col = sell->colmap[idxn[j]] - 1; 275 #endif 276 if ((col < 0) || (sell->garray[col] != idxn[j])) *(v + i * n + j) = 0.0; 277 else PetscCall(MatGetValues(sell->B, 1, &row, 1, &col, v + i * n + j)); 278 } 279 } 280 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only local values currently supported"); 281 } 282 PetscFunctionReturn(PETSC_SUCCESS); 283 } 284 285 extern PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat, Vec, Vec); 286 287 PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat, MatAssemblyType mode) 288 { 289 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 290 PetscInt nstash, reallocs; 291 292 PetscFunctionBegin; 293 if (sell->donotstash || mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS); 294 295 PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range)); 296 PetscCall(MatStashGetInfo_Private(&mat->stash, &nstash, &reallocs)); 297 PetscCall(PetscInfo(sell->A, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs)); 298 PetscFunctionReturn(PETSC_SUCCESS); 299 } 300 301 PetscErrorCode MatAssemblyEnd_MPISELL(Mat mat, MatAssemblyType mode) 302 { 303 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 304 PetscMPIInt n; 305 PetscInt i, flg; 306 PetscInt *row, *col; 307 PetscScalar *val; 308 PetscBool other_disassembled; 309 /* do not use 'b = (Mat_SeqSELL*)sell->B->data' as B can be reset in disassembly */ 310 PetscFunctionBegin; 311 if (!sell->donotstash && !mat->nooffprocentries) { 312 while (1) { 313 PetscCall(MatStashScatterGetMesg_Private(&mat->stash, &n, &row, &col, &val, &flg)); 314 if (!flg) break; 315 316 for (i = 0; i < n; i++) { /* assemble one by one */ 317 PetscCall(MatSetValues_MPISELL(mat, 1, row + i, 1, col + i, val + i, mat->insertmode)); 318 } 319 } 320 PetscCall(MatStashScatterEnd_Private(&mat->stash)); 321 } 322 #if defined(PETSC_HAVE_CUDA) 323 if (mat->offloadmask == PETSC_OFFLOAD_CPU) sell->A->offloadmask = PETSC_OFFLOAD_CPU; 324 #endif 325 PetscCall(MatAssemblyBegin(sell->A, mode)); 326 PetscCall(MatAssemblyEnd(sell->A, mode)); 327 328 /* 329 determine if any processor has disassembled, if so we must 330 also disassemble ourselves, in order that we may reassemble. 331 */ 332 /* 333 if nonzero structure of submatrix B cannot change then we know that 334 no processor disassembled thus we can skip this stuff 335 */ 336 if (!((Mat_SeqSELL *)sell->B->data)->nonew) { 337 PetscCall(MPIU_Allreduce(&mat->was_assembled, &other_disassembled, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat))); 338 if (mat->was_assembled && !other_disassembled) PetscCall(MatDisAssemble_MPISELL(mat)); 339 } 340 if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) PetscCall(MatSetUpMultiply_MPISELL(mat)); 341 #if defined(PETSC_HAVE_CUDA) 342 if (mat->offloadmask == PETSC_OFFLOAD_CPU && sell->B->offloadmask != PETSC_OFFLOAD_UNALLOCATED) sell->B->offloadmask = PETSC_OFFLOAD_CPU; 343 #endif 344 PetscCall(MatAssemblyBegin(sell->B, mode)); 345 PetscCall(MatAssemblyEnd(sell->B, mode)); 346 PetscCall(PetscFree2(sell->rowvalues, sell->rowindices)); 347 sell->rowvalues = NULL; 348 PetscCall(VecDestroy(&sell->diag)); 349 350 /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */ 351 if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqSELL *)(sell->A->data))->nonew) { 352 PetscObjectState state = sell->A->nonzerostate + sell->B->nonzerostate; 353 PetscCall(MPIU_Allreduce(&state, &mat->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)mat))); 354 } 355 #if defined(PETSC_HAVE_CUDA) 356 mat->offloadmask = PETSC_OFFLOAD_BOTH; 357 #endif 358 PetscFunctionReturn(PETSC_SUCCESS); 359 } 360 361 PetscErrorCode MatZeroEntries_MPISELL(Mat A) 362 { 363 Mat_MPISELL *l = (Mat_MPISELL *)A->data; 364 365 PetscFunctionBegin; 366 PetscCall(MatZeroEntries(l->A)); 367 PetscCall(MatZeroEntries(l->B)); 368 PetscFunctionReturn(PETSC_SUCCESS); 369 } 370 371 PetscErrorCode MatMult_MPISELL(Mat A, Vec xx, Vec yy) 372 { 373 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 374 PetscInt nt; 375 376 PetscFunctionBegin; 377 PetscCall(VecGetLocalSize(xx, &nt)); 378 PetscCheck(nt == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible partition of A (%" PetscInt_FMT ") and xx (%" PetscInt_FMT ")", A->cmap->n, nt); 379 PetscCall(VecScatterBegin(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD)); 380 PetscCall((*a->A->ops->mult)(a->A, xx, yy)); 381 PetscCall(VecScatterEnd(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD)); 382 PetscCall((*a->B->ops->multadd)(a->B, a->lvec, yy, yy)); 383 PetscFunctionReturn(PETSC_SUCCESS); 384 } 385 386 PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat A, Vec bb, Vec xx) 387 { 388 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 389 390 PetscFunctionBegin; 391 PetscCall(MatMultDiagonalBlock(a->A, bb, xx)); 392 PetscFunctionReturn(PETSC_SUCCESS); 393 } 394 395 PetscErrorCode MatMultAdd_MPISELL(Mat A, Vec xx, Vec yy, Vec zz) 396 { 397 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 398 399 PetscFunctionBegin; 400 PetscCall(VecScatterBegin(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD)); 401 PetscCall((*a->A->ops->multadd)(a->A, xx, yy, zz)); 402 PetscCall(VecScatterEnd(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD)); 403 PetscCall((*a->B->ops->multadd)(a->B, a->lvec, zz, zz)); 404 PetscFunctionReturn(PETSC_SUCCESS); 405 } 406 407 PetscErrorCode MatMultTranspose_MPISELL(Mat A, Vec xx, Vec yy) 408 { 409 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 410 411 PetscFunctionBegin; 412 /* do nondiagonal part */ 413 PetscCall((*a->B->ops->multtranspose)(a->B, xx, a->lvec)); 414 /* do local part */ 415 PetscCall((*a->A->ops->multtranspose)(a->A, xx, yy)); 416 /* add partial results together */ 417 PetscCall(VecScatterBegin(a->Mvctx, a->lvec, yy, ADD_VALUES, SCATTER_REVERSE)); 418 PetscCall(VecScatterEnd(a->Mvctx, a->lvec, yy, ADD_VALUES, SCATTER_REVERSE)); 419 PetscFunctionReturn(PETSC_SUCCESS); 420 } 421 422 PetscErrorCode MatIsTranspose_MPISELL(Mat Amat, Mat Bmat, PetscReal tol, PetscBool *f) 423 { 424 MPI_Comm comm; 425 Mat_MPISELL *Asell = (Mat_MPISELL *)Amat->data, *Bsell; 426 Mat Adia = Asell->A, Bdia, Aoff, Boff, *Aoffs, *Boffs; 427 IS Me, Notme; 428 PetscInt M, N, first, last, *notme, i; 429 PetscMPIInt size; 430 431 PetscFunctionBegin; 432 /* Easy test: symmetric diagonal block */ 433 Bsell = (Mat_MPISELL *)Bmat->data; 434 Bdia = Bsell->A; 435 PetscCall(MatIsTranspose(Adia, Bdia, tol, f)); 436 if (!*f) PetscFunctionReturn(PETSC_SUCCESS); 437 PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm)); 438 PetscCallMPI(MPI_Comm_size(comm, &size)); 439 if (size == 1) PetscFunctionReturn(PETSC_SUCCESS); 440 441 /* Hard test: off-diagonal block. This takes a MatCreateSubMatrix. */ 442 PetscCall(MatGetSize(Amat, &M, &N)); 443 PetscCall(MatGetOwnershipRange(Amat, &first, &last)); 444 PetscCall(PetscMalloc1(N - last + first, ¬me)); 445 for (i = 0; i < first; i++) notme[i] = i; 446 for (i = last; i < M; i++) notme[i - last + first] = i; 447 PetscCall(ISCreateGeneral(MPI_COMM_SELF, N - last + first, notme, PETSC_COPY_VALUES, &Notme)); 448 PetscCall(ISCreateStride(MPI_COMM_SELF, last - first, first, 1, &Me)); 449 PetscCall(MatCreateSubMatrices(Amat, 1, &Me, &Notme, MAT_INITIAL_MATRIX, &Aoffs)); 450 Aoff = Aoffs[0]; 451 PetscCall(MatCreateSubMatrices(Bmat, 1, &Notme, &Me, MAT_INITIAL_MATRIX, &Boffs)); 452 Boff = Boffs[0]; 453 PetscCall(MatIsTranspose(Aoff, Boff, tol, f)); 454 PetscCall(MatDestroyMatrices(1, &Aoffs)); 455 PetscCall(MatDestroyMatrices(1, &Boffs)); 456 PetscCall(ISDestroy(&Me)); 457 PetscCall(ISDestroy(&Notme)); 458 PetscCall(PetscFree(notme)); 459 PetscFunctionReturn(PETSC_SUCCESS); 460 } 461 462 PetscErrorCode MatMultTransposeAdd_MPISELL(Mat A, Vec xx, Vec yy, Vec zz) 463 { 464 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 465 466 PetscFunctionBegin; 467 /* do nondiagonal part */ 468 PetscCall((*a->B->ops->multtranspose)(a->B, xx, a->lvec)); 469 /* do local part */ 470 PetscCall((*a->A->ops->multtransposeadd)(a->A, xx, yy, zz)); 471 /* add partial results together */ 472 PetscCall(VecScatterBegin(a->Mvctx, a->lvec, zz, ADD_VALUES, SCATTER_REVERSE)); 473 PetscCall(VecScatterEnd(a->Mvctx, a->lvec, zz, ADD_VALUES, SCATTER_REVERSE)); 474 PetscFunctionReturn(PETSC_SUCCESS); 475 } 476 477 /* 478 This only works correctly for square matrices where the subblock A->A is the 479 diagonal block 480 */ 481 PetscErrorCode MatGetDiagonal_MPISELL(Mat A, Vec v) 482 { 483 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 484 485 PetscFunctionBegin; 486 PetscCheck(A->rmap->N == A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Supports only square matrix where A->A is diag block"); 487 PetscCheck(A->rmap->rstart == A->cmap->rstart && A->rmap->rend == A->cmap->rend, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "row partition must equal col partition"); 488 PetscCall(MatGetDiagonal(a->A, v)); 489 PetscFunctionReturn(PETSC_SUCCESS); 490 } 491 492 PetscErrorCode MatScale_MPISELL(Mat A, PetscScalar aa) 493 { 494 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 495 496 PetscFunctionBegin; 497 PetscCall(MatScale(a->A, aa)); 498 PetscCall(MatScale(a->B, aa)); 499 PetscFunctionReturn(PETSC_SUCCESS); 500 } 501 502 PetscErrorCode MatDestroy_MPISELL(Mat mat) 503 { 504 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 505 506 PetscFunctionBegin; 507 #if defined(PETSC_USE_LOG) 508 PetscCall(PetscLogObjectState((PetscObject)mat, "Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT, mat->rmap->N, mat->cmap->N)); 509 #endif 510 PetscCall(MatStashDestroy_Private(&mat->stash)); 511 PetscCall(VecDestroy(&sell->diag)); 512 PetscCall(MatDestroy(&sell->A)); 513 PetscCall(MatDestroy(&sell->B)); 514 #if defined(PETSC_USE_CTABLE) 515 PetscCall(PetscHMapIDestroy(&sell->colmap)); 516 #else 517 PetscCall(PetscFree(sell->colmap)); 518 #endif 519 PetscCall(PetscFree(sell->garray)); 520 PetscCall(VecDestroy(&sell->lvec)); 521 PetscCall(VecScatterDestroy(&sell->Mvctx)); 522 PetscCall(PetscFree2(sell->rowvalues, sell->rowindices)); 523 PetscCall(PetscFree(sell->ld)); 524 PetscCall(PetscFree(mat->data)); 525 526 PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL)); 527 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatStoreValues_C", NULL)); 528 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatRetrieveValues_C", NULL)); 529 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatIsTranspose_C", NULL)); 530 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPISELLSetPreallocation_C", NULL)); 531 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisell_mpiaij_C", NULL)); 532 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalScaleLocal_C", NULL)); 533 PetscFunctionReturn(PETSC_SUCCESS); 534 } 535 536 #include <petscdraw.h> 537 PetscErrorCode MatView_MPISELL_ASCIIorDraworSocket(Mat mat, PetscViewer viewer) 538 { 539 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 540 PetscMPIInt rank = sell->rank, size = sell->size; 541 PetscBool isdraw, iascii, isbinary; 542 PetscViewer sviewer; 543 PetscViewerFormat format; 544 545 PetscFunctionBegin; 546 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 547 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 548 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 549 if (iascii) { 550 PetscCall(PetscViewerGetFormat(viewer, &format)); 551 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 552 MatInfo info; 553 PetscInt *inodes; 554 555 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat), &rank)); 556 PetscCall(MatGetInfo(mat, MAT_LOCAL, &info)); 557 PetscCall(MatInodeGetInodeSizes(sell->A, NULL, &inodes, NULL)); 558 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 559 if (!inodes) { 560 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", not using I-node routines\n", rank, mat->rmap->n, (PetscInt)info.nz_used, 561 (PetscInt)info.nz_allocated, (PetscInt)info.memory)); 562 } else { 563 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", using I-node routines\n", rank, mat->rmap->n, (PetscInt)info.nz_used, 564 (PetscInt)info.nz_allocated, (PetscInt)info.memory)); 565 } 566 PetscCall(MatGetInfo(sell->A, MAT_LOCAL, &info)); 567 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] on-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used)); 568 PetscCall(MatGetInfo(sell->B, MAT_LOCAL, &info)); 569 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] off-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used)); 570 PetscCall(PetscViewerFlush(viewer)); 571 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 572 PetscCall(PetscViewerASCIIPrintf(viewer, "Information on VecScatter used in matrix-vector product: \n")); 573 PetscCall(VecScatterView(sell->Mvctx, viewer)); 574 PetscFunctionReturn(PETSC_SUCCESS); 575 } else if (format == PETSC_VIEWER_ASCII_INFO) { 576 PetscInt inodecount, inodelimit, *inodes; 577 PetscCall(MatInodeGetInodeSizes(sell->A, &inodecount, &inodes, &inodelimit)); 578 if (inodes) { 579 PetscCall(PetscViewerASCIIPrintf(viewer, "using I-node (on process 0) routines: found %" PetscInt_FMT " nodes, limit used is %" PetscInt_FMT "\n", inodecount, inodelimit)); 580 } else { 581 PetscCall(PetscViewerASCIIPrintf(viewer, "not using I-node (on process 0) routines\n")); 582 } 583 PetscFunctionReturn(PETSC_SUCCESS); 584 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 585 PetscFunctionReturn(PETSC_SUCCESS); 586 } 587 } else if (isbinary) { 588 if (size == 1) { 589 PetscCall(PetscObjectSetName((PetscObject)sell->A, ((PetscObject)mat)->name)); 590 PetscCall(MatView(sell->A, viewer)); 591 } else { 592 /* PetscCall(MatView_MPISELL_Binary(mat,viewer)); */ 593 } 594 PetscFunctionReturn(PETSC_SUCCESS); 595 } else if (isdraw) { 596 PetscDraw draw; 597 PetscBool isnull; 598 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 599 PetscCall(PetscDrawIsNull(draw, &isnull)); 600 if (isnull) PetscFunctionReturn(PETSC_SUCCESS); 601 } 602 603 { 604 /* assemble the entire matrix onto first processor. */ 605 Mat A; 606 Mat_SeqSELL *Aloc; 607 PetscInt M = mat->rmap->N, N = mat->cmap->N, *acolidx, row, col, i, j; 608 MatScalar *aval; 609 PetscBool isnonzero; 610 611 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &A)); 612 if (rank == 0) { 613 PetscCall(MatSetSizes(A, M, N, M, N)); 614 } else { 615 PetscCall(MatSetSizes(A, 0, 0, M, N)); 616 } 617 /* This is just a temporary matrix, so explicitly using MATMPISELL is probably best */ 618 PetscCall(MatSetType(A, MATMPISELL)); 619 PetscCall(MatMPISELLSetPreallocation(A, 0, NULL, 0, NULL)); 620 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE)); 621 622 /* copy over the A part */ 623 Aloc = (Mat_SeqSELL *)sell->A->data; 624 acolidx = Aloc->colidx; 625 aval = Aloc->val; 626 for (i = 0; i < Aloc->totalslices; i++) { /* loop over slices */ 627 for (j = Aloc->sliidx[i]; j < Aloc->sliidx[i + 1]; j++) { 628 isnonzero = (PetscBool)((j - Aloc->sliidx[i]) / Aloc->sliceheight < Aloc->rlen[i * Aloc->sliceheight + j % Aloc->sliceheight]); 629 if (isnonzero) { /* check the mask bit */ 630 row = i * Aloc->sliceheight + j % Aloc->sliceheight + mat->rmap->rstart; 631 col = *acolidx + mat->rmap->rstart; 632 PetscCall(MatSetValues(A, 1, &row, 1, &col, aval, INSERT_VALUES)); 633 } 634 aval++; 635 acolidx++; 636 } 637 } 638 639 /* copy over the B part */ 640 Aloc = (Mat_SeqSELL *)sell->B->data; 641 acolidx = Aloc->colidx; 642 aval = Aloc->val; 643 for (i = 0; i < Aloc->totalslices; i++) { 644 for (j = Aloc->sliidx[i]; j < Aloc->sliidx[i + 1]; j++) { 645 isnonzero = (PetscBool)((j - Aloc->sliidx[i]) / Aloc->sliceheight < Aloc->rlen[i * Aloc->sliceheight + j % Aloc->sliceheight]); 646 if (isnonzero) { 647 row = i * Aloc->sliceheight + j % Aloc->sliceheight + mat->rmap->rstart; 648 col = sell->garray[*acolidx]; 649 PetscCall(MatSetValues(A, 1, &row, 1, &col, aval, INSERT_VALUES)); 650 } 651 aval++; 652 acolidx++; 653 } 654 } 655 656 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 657 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 658 /* 659 Everyone has to call to draw the matrix since the graphics waits are 660 synchronized across all processors that share the PetscDraw object 661 */ 662 PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer)); 663 if (rank == 0) { 664 PetscCall(PetscObjectSetName((PetscObject)((Mat_MPISELL *)(A->data))->A, ((PetscObject)mat)->name)); 665 PetscCall(MatView_SeqSELL(((Mat_MPISELL *)(A->data))->A, sviewer)); 666 } 667 PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer)); 668 PetscCall(PetscViewerFlush(viewer)); 669 PetscCall(MatDestroy(&A)); 670 } 671 PetscFunctionReturn(PETSC_SUCCESS); 672 } 673 674 PetscErrorCode MatView_MPISELL(Mat mat, PetscViewer viewer) 675 { 676 PetscBool iascii, isdraw, issocket, isbinary; 677 678 PetscFunctionBegin; 679 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 680 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 681 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 682 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSOCKET, &issocket)); 683 if (iascii || isdraw || isbinary || issocket) PetscCall(MatView_MPISELL_ASCIIorDraworSocket(mat, viewer)); 684 PetscFunctionReturn(PETSC_SUCCESS); 685 } 686 687 PetscErrorCode MatGetGhosts_MPISELL(Mat mat, PetscInt *nghosts, const PetscInt *ghosts[]) 688 { 689 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 690 691 PetscFunctionBegin; 692 PetscCall(MatGetSize(sell->B, NULL, nghosts)); 693 if (ghosts) *ghosts = sell->garray; 694 PetscFunctionReturn(PETSC_SUCCESS); 695 } 696 697 PetscErrorCode MatGetInfo_MPISELL(Mat matin, MatInfoType flag, MatInfo *info) 698 { 699 Mat_MPISELL *mat = (Mat_MPISELL *)matin->data; 700 Mat A = mat->A, B = mat->B; 701 PetscLogDouble isend[5], irecv[5]; 702 703 PetscFunctionBegin; 704 info->block_size = 1.0; 705 PetscCall(MatGetInfo(A, MAT_LOCAL, info)); 706 707 isend[0] = info->nz_used; 708 isend[1] = info->nz_allocated; 709 isend[2] = info->nz_unneeded; 710 isend[3] = info->memory; 711 isend[4] = info->mallocs; 712 713 PetscCall(MatGetInfo(B, MAT_LOCAL, info)); 714 715 isend[0] += info->nz_used; 716 isend[1] += info->nz_allocated; 717 isend[2] += info->nz_unneeded; 718 isend[3] += info->memory; 719 isend[4] += info->mallocs; 720 if (flag == MAT_LOCAL) { 721 info->nz_used = isend[0]; 722 info->nz_allocated = isend[1]; 723 info->nz_unneeded = isend[2]; 724 info->memory = isend[3]; 725 info->mallocs = isend[4]; 726 } else if (flag == MAT_GLOBAL_MAX) { 727 PetscCall(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)matin))); 728 729 info->nz_used = irecv[0]; 730 info->nz_allocated = irecv[1]; 731 info->nz_unneeded = irecv[2]; 732 info->memory = irecv[3]; 733 info->mallocs = irecv[4]; 734 } else if (flag == MAT_GLOBAL_SUM) { 735 PetscCall(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)matin))); 736 737 info->nz_used = irecv[0]; 738 info->nz_allocated = irecv[1]; 739 info->nz_unneeded = irecv[2]; 740 info->memory = irecv[3]; 741 info->mallocs = irecv[4]; 742 } 743 info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */ 744 info->fill_ratio_needed = 0; 745 info->factor_mallocs = 0; 746 PetscFunctionReturn(PETSC_SUCCESS); 747 } 748 749 PetscErrorCode MatSetOption_MPISELL(Mat A, MatOption op, PetscBool flg) 750 { 751 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 752 753 PetscFunctionBegin; 754 switch (op) { 755 case MAT_NEW_NONZERO_LOCATIONS: 756 case MAT_NEW_NONZERO_ALLOCATION_ERR: 757 case MAT_UNUSED_NONZERO_LOCATION_ERR: 758 case MAT_KEEP_NONZERO_PATTERN: 759 case MAT_NEW_NONZERO_LOCATION_ERR: 760 case MAT_USE_INODES: 761 case MAT_IGNORE_ZERO_ENTRIES: 762 MatCheckPreallocated(A, 1); 763 PetscCall(MatSetOption(a->A, op, flg)); 764 PetscCall(MatSetOption(a->B, op, flg)); 765 break; 766 case MAT_ROW_ORIENTED: 767 MatCheckPreallocated(A, 1); 768 a->roworiented = flg; 769 770 PetscCall(MatSetOption(a->A, op, flg)); 771 PetscCall(MatSetOption(a->B, op, flg)); 772 break; 773 case MAT_FORCE_DIAGONAL_ENTRIES: 774 case MAT_SORTED_FULL: 775 PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op])); 776 break; 777 case MAT_IGNORE_OFF_PROC_ENTRIES: 778 a->donotstash = flg; 779 break; 780 case MAT_SPD: 781 case MAT_SPD_ETERNAL: 782 break; 783 case MAT_SYMMETRIC: 784 MatCheckPreallocated(A, 1); 785 PetscCall(MatSetOption(a->A, op, flg)); 786 break; 787 case MAT_STRUCTURALLY_SYMMETRIC: 788 MatCheckPreallocated(A, 1); 789 PetscCall(MatSetOption(a->A, op, flg)); 790 break; 791 case MAT_HERMITIAN: 792 MatCheckPreallocated(A, 1); 793 PetscCall(MatSetOption(a->A, op, flg)); 794 break; 795 case MAT_SYMMETRY_ETERNAL: 796 MatCheckPreallocated(A, 1); 797 PetscCall(MatSetOption(a->A, op, flg)); 798 break; 799 case MAT_STRUCTURAL_SYMMETRY_ETERNAL: 800 MatCheckPreallocated(A, 1); 801 PetscCall(MatSetOption(a->A, op, flg)); 802 break; 803 default: 804 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %d", op); 805 } 806 PetscFunctionReturn(PETSC_SUCCESS); 807 } 808 809 PetscErrorCode MatDiagonalScale_MPISELL(Mat mat, Vec ll, Vec rr) 810 { 811 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 812 Mat a = sell->A, b = sell->B; 813 PetscInt s1, s2, s3; 814 815 PetscFunctionBegin; 816 PetscCall(MatGetLocalSize(mat, &s2, &s3)); 817 if (rr) { 818 PetscCall(VecGetLocalSize(rr, &s1)); 819 PetscCheck(s1 == s3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "right vector non-conforming local size"); 820 /* Overlap communication with computation. */ 821 PetscCall(VecScatterBegin(sell->Mvctx, rr, sell->lvec, INSERT_VALUES, SCATTER_FORWARD)); 822 } 823 if (ll) { 824 PetscCall(VecGetLocalSize(ll, &s1)); 825 PetscCheck(s1 == s2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "left vector non-conforming local size"); 826 PetscUseTypeMethod(b, diagonalscale, ll, NULL); 827 } 828 /* scale the diagonal block */ 829 PetscUseTypeMethod(a, diagonalscale, ll, rr); 830 831 if (rr) { 832 /* Do a scatter end and then right scale the off-diagonal block */ 833 PetscCall(VecScatterEnd(sell->Mvctx, rr, sell->lvec, INSERT_VALUES, SCATTER_FORWARD)); 834 PetscUseTypeMethod(b, diagonalscale, NULL, sell->lvec); 835 } 836 PetscFunctionReturn(PETSC_SUCCESS); 837 } 838 839 PetscErrorCode MatSetUnfactored_MPISELL(Mat A) 840 { 841 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 842 843 PetscFunctionBegin; 844 PetscCall(MatSetUnfactored(a->A)); 845 PetscFunctionReturn(PETSC_SUCCESS); 846 } 847 848 PetscErrorCode MatEqual_MPISELL(Mat A, Mat B, PetscBool *flag) 849 { 850 Mat_MPISELL *matB = (Mat_MPISELL *)B->data, *matA = (Mat_MPISELL *)A->data; 851 Mat a, b, c, d; 852 PetscBool flg; 853 854 PetscFunctionBegin; 855 a = matA->A; 856 b = matA->B; 857 c = matB->A; 858 d = matB->B; 859 860 PetscCall(MatEqual(a, c, &flg)); 861 if (flg) PetscCall(MatEqual(b, d, &flg)); 862 PetscCall(MPIU_Allreduce(&flg, flag, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 863 PetscFunctionReturn(PETSC_SUCCESS); 864 } 865 866 PetscErrorCode MatCopy_MPISELL(Mat A, Mat B, MatStructure str) 867 { 868 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 869 Mat_MPISELL *b = (Mat_MPISELL *)B->data; 870 871 PetscFunctionBegin; 872 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 873 if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) { 874 /* because of the column compression in the off-processor part of the matrix a->B, 875 the number of columns in a->B and b->B may be different, hence we cannot call 876 the MatCopy() directly on the two parts. If need be, we can provide a more 877 efficient copy than the MatCopy_Basic() by first uncompressing the a->B matrices 878 then copying the submatrices */ 879 PetscCall(MatCopy_Basic(A, B, str)); 880 } else { 881 PetscCall(MatCopy(a->A, b->A, str)); 882 PetscCall(MatCopy(a->B, b->B, str)); 883 } 884 PetscFunctionReturn(PETSC_SUCCESS); 885 } 886 887 PetscErrorCode MatSetUp_MPISELL(Mat A) 888 { 889 PetscFunctionBegin; 890 PetscCall(MatMPISELLSetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL)); 891 PetscFunctionReturn(PETSC_SUCCESS); 892 } 893 894 extern PetscErrorCode MatConjugate_SeqSELL(Mat); 895 896 PetscErrorCode MatConjugate_MPISELL(Mat mat) 897 { 898 PetscFunctionBegin; 899 if (PetscDefined(USE_COMPLEX)) { 900 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 901 902 PetscCall(MatConjugate_SeqSELL(sell->A)); 903 PetscCall(MatConjugate_SeqSELL(sell->B)); 904 } 905 PetscFunctionReturn(PETSC_SUCCESS); 906 } 907 908 PetscErrorCode MatRealPart_MPISELL(Mat A) 909 { 910 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 911 912 PetscFunctionBegin; 913 PetscCall(MatRealPart(a->A)); 914 PetscCall(MatRealPart(a->B)); 915 PetscFunctionReturn(PETSC_SUCCESS); 916 } 917 918 PetscErrorCode MatImaginaryPart_MPISELL(Mat A) 919 { 920 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 921 922 PetscFunctionBegin; 923 PetscCall(MatImaginaryPart(a->A)); 924 PetscCall(MatImaginaryPart(a->B)); 925 PetscFunctionReturn(PETSC_SUCCESS); 926 } 927 928 PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A, const PetscScalar **values) 929 { 930 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 931 932 PetscFunctionBegin; 933 PetscCall(MatInvertBlockDiagonal(a->A, values)); 934 A->factorerrortype = a->A->factorerrortype; 935 PetscFunctionReturn(PETSC_SUCCESS); 936 } 937 938 static PetscErrorCode MatSetRandom_MPISELL(Mat x, PetscRandom rctx) 939 { 940 Mat_MPISELL *sell = (Mat_MPISELL *)x->data; 941 942 PetscFunctionBegin; 943 PetscCall(MatSetRandom(sell->A, rctx)); 944 PetscCall(MatSetRandom(sell->B, rctx)); 945 PetscCall(MatAssemblyBegin(x, MAT_FINAL_ASSEMBLY)); 946 PetscCall(MatAssemblyEnd(x, MAT_FINAL_ASSEMBLY)); 947 PetscFunctionReturn(PETSC_SUCCESS); 948 } 949 950 PetscErrorCode MatSetFromOptions_MPISELL(Mat A, PetscOptionItems *PetscOptionsObject) 951 { 952 PetscFunctionBegin; 953 PetscOptionsHeadBegin(PetscOptionsObject, "MPISELL options"); 954 PetscOptionsHeadEnd(); 955 PetscFunctionReturn(PETSC_SUCCESS); 956 } 957 958 PetscErrorCode MatShift_MPISELL(Mat Y, PetscScalar a) 959 { 960 Mat_MPISELL *msell = (Mat_MPISELL *)Y->data; 961 Mat_SeqSELL *sell = (Mat_SeqSELL *)msell->A->data; 962 963 PetscFunctionBegin; 964 if (!Y->preallocated) { 965 PetscCall(MatMPISELLSetPreallocation(Y, 1, NULL, 0, NULL)); 966 } else if (!sell->nz) { 967 PetscInt nonew = sell->nonew; 968 PetscCall(MatSeqSELLSetPreallocation(msell->A, 1, NULL)); 969 sell->nonew = nonew; 970 } 971 PetscCall(MatShift_Basic(Y, a)); 972 PetscFunctionReturn(PETSC_SUCCESS); 973 } 974 975 PetscErrorCode MatMissingDiagonal_MPISELL(Mat A, PetscBool *missing, PetscInt *d) 976 { 977 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 978 979 PetscFunctionBegin; 980 PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only works for square matrices"); 981 PetscCall(MatMissingDiagonal(a->A, missing, d)); 982 if (d) { 983 PetscInt rstart; 984 PetscCall(MatGetOwnershipRange(A, &rstart, NULL)); 985 *d += rstart; 986 } 987 PetscFunctionReturn(PETSC_SUCCESS); 988 } 989 990 PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A, Mat *a) 991 { 992 PetscFunctionBegin; 993 *a = ((Mat_MPISELL *)A->data)->A; 994 PetscFunctionReturn(PETSC_SUCCESS); 995 } 996 997 static struct _MatOps MatOps_Values = {MatSetValues_MPISELL, 998 NULL, 999 NULL, 1000 MatMult_MPISELL, 1001 /* 4*/ MatMultAdd_MPISELL, 1002 MatMultTranspose_MPISELL, 1003 MatMultTransposeAdd_MPISELL, 1004 NULL, 1005 NULL, 1006 NULL, 1007 /*10*/ NULL, 1008 NULL, 1009 NULL, 1010 MatSOR_MPISELL, 1011 NULL, 1012 /*15*/ MatGetInfo_MPISELL, 1013 MatEqual_MPISELL, 1014 MatGetDiagonal_MPISELL, 1015 MatDiagonalScale_MPISELL, 1016 NULL, 1017 /*20*/ MatAssemblyBegin_MPISELL, 1018 MatAssemblyEnd_MPISELL, 1019 MatSetOption_MPISELL, 1020 MatZeroEntries_MPISELL, 1021 /*24*/ NULL, 1022 NULL, 1023 NULL, 1024 NULL, 1025 NULL, 1026 /*29*/ MatSetUp_MPISELL, 1027 NULL, 1028 NULL, 1029 MatGetDiagonalBlock_MPISELL, 1030 NULL, 1031 /*34*/ MatDuplicate_MPISELL, 1032 NULL, 1033 NULL, 1034 NULL, 1035 NULL, 1036 /*39*/ NULL, 1037 NULL, 1038 NULL, 1039 MatGetValues_MPISELL, 1040 MatCopy_MPISELL, 1041 /*44*/ NULL, 1042 MatScale_MPISELL, 1043 MatShift_MPISELL, 1044 MatDiagonalSet_MPISELL, 1045 NULL, 1046 /*49*/ MatSetRandom_MPISELL, 1047 NULL, 1048 NULL, 1049 NULL, 1050 NULL, 1051 /*54*/ MatFDColoringCreate_MPIXAIJ, 1052 NULL, 1053 MatSetUnfactored_MPISELL, 1054 NULL, 1055 NULL, 1056 /*59*/ NULL, 1057 MatDestroy_MPISELL, 1058 MatView_MPISELL, 1059 NULL, 1060 NULL, 1061 /*64*/ NULL, 1062 NULL, 1063 NULL, 1064 NULL, 1065 NULL, 1066 /*69*/ NULL, 1067 NULL, 1068 NULL, 1069 NULL, 1070 NULL, 1071 NULL, 1072 /*75*/ MatFDColoringApply_AIJ, /* reuse AIJ function */ 1073 MatSetFromOptions_MPISELL, 1074 NULL, 1075 NULL, 1076 NULL, 1077 /*80*/ NULL, 1078 NULL, 1079 NULL, 1080 /*83*/ NULL, 1081 NULL, 1082 NULL, 1083 NULL, 1084 NULL, 1085 NULL, 1086 /*89*/ NULL, 1087 NULL, 1088 NULL, 1089 NULL, 1090 NULL, 1091 /*94*/ NULL, 1092 NULL, 1093 NULL, 1094 NULL, 1095 NULL, 1096 /*99*/ NULL, 1097 NULL, 1098 NULL, 1099 MatConjugate_MPISELL, 1100 NULL, 1101 /*104*/ NULL, 1102 MatRealPart_MPISELL, 1103 MatImaginaryPart_MPISELL, 1104 NULL, 1105 NULL, 1106 /*109*/ NULL, 1107 NULL, 1108 NULL, 1109 NULL, 1110 MatMissingDiagonal_MPISELL, 1111 /*114*/ NULL, 1112 NULL, 1113 MatGetGhosts_MPISELL, 1114 NULL, 1115 NULL, 1116 /*119*/ NULL, 1117 NULL, 1118 NULL, 1119 NULL, 1120 NULL, 1121 /*124*/ NULL, 1122 NULL, 1123 MatInvertBlockDiagonal_MPISELL, 1124 NULL, 1125 NULL, 1126 /*129*/ NULL, 1127 NULL, 1128 NULL, 1129 NULL, 1130 NULL, 1131 /*134*/ NULL, 1132 NULL, 1133 NULL, 1134 NULL, 1135 NULL, 1136 /*139*/ NULL, 1137 NULL, 1138 NULL, 1139 MatFDColoringSetUp_MPIXAIJ, 1140 NULL, 1141 /*144*/ NULL, 1142 NULL, 1143 NULL, 1144 NULL, 1145 NULL, 1146 NULL, 1147 /*150*/ NULL, 1148 NULL}; 1149 1150 PetscErrorCode MatStoreValues_MPISELL(Mat mat) 1151 { 1152 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 1153 1154 PetscFunctionBegin; 1155 PetscCall(MatStoreValues(sell->A)); 1156 PetscCall(MatStoreValues(sell->B)); 1157 PetscFunctionReturn(PETSC_SUCCESS); 1158 } 1159 1160 PetscErrorCode MatRetrieveValues_MPISELL(Mat mat) 1161 { 1162 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 1163 1164 PetscFunctionBegin; 1165 PetscCall(MatRetrieveValues(sell->A)); 1166 PetscCall(MatRetrieveValues(sell->B)); 1167 PetscFunctionReturn(PETSC_SUCCESS); 1168 } 1169 1170 PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B, PetscInt d_rlenmax, const PetscInt d_rlen[], PetscInt o_rlenmax, const PetscInt o_rlen[]) 1171 { 1172 Mat_MPISELL *b; 1173 1174 PetscFunctionBegin; 1175 PetscCall(PetscLayoutSetUp(B->rmap)); 1176 PetscCall(PetscLayoutSetUp(B->cmap)); 1177 b = (Mat_MPISELL *)B->data; 1178 1179 if (!B->preallocated) { 1180 /* Explicitly create 2 MATSEQSELL matrices. */ 1181 PetscCall(MatCreate(PETSC_COMM_SELF, &b->A)); 1182 PetscCall(MatSetSizes(b->A, B->rmap->n, B->cmap->n, B->rmap->n, B->cmap->n)); 1183 PetscCall(MatSetBlockSizesFromMats(b->A, B, B)); 1184 PetscCall(MatSetType(b->A, MATSEQSELL)); 1185 PetscCall(MatCreate(PETSC_COMM_SELF, &b->B)); 1186 PetscCall(MatSetSizes(b->B, B->rmap->n, B->cmap->N, B->rmap->n, B->cmap->N)); 1187 PetscCall(MatSetBlockSizesFromMats(b->B, B, B)); 1188 PetscCall(MatSetType(b->B, MATSEQSELL)); 1189 } 1190 1191 PetscCall(MatSeqSELLSetPreallocation(b->A, d_rlenmax, d_rlen)); 1192 PetscCall(MatSeqSELLSetPreallocation(b->B, o_rlenmax, o_rlen)); 1193 B->preallocated = PETSC_TRUE; 1194 B->was_assembled = PETSC_FALSE; 1195 /* 1196 critical for MatAssemblyEnd to work. 1197 MatAssemblyBegin checks it to set up was_assembled 1198 and MatAssemblyEnd checks was_assembled to determine whether to build garray 1199 */ 1200 B->assembled = PETSC_FALSE; 1201 PetscFunctionReturn(PETSC_SUCCESS); 1202 } 1203 1204 PetscErrorCode MatDuplicate_MPISELL(Mat matin, MatDuplicateOption cpvalues, Mat *newmat) 1205 { 1206 Mat mat; 1207 Mat_MPISELL *a, *oldmat = (Mat_MPISELL *)matin->data; 1208 1209 PetscFunctionBegin; 1210 *newmat = NULL; 1211 PetscCall(MatCreate(PetscObjectComm((PetscObject)matin), &mat)); 1212 PetscCall(MatSetSizes(mat, matin->rmap->n, matin->cmap->n, matin->rmap->N, matin->cmap->N)); 1213 PetscCall(MatSetBlockSizesFromMats(mat, matin, matin)); 1214 PetscCall(MatSetType(mat, ((PetscObject)matin)->type_name)); 1215 a = (Mat_MPISELL *)mat->data; 1216 1217 mat->factortype = matin->factortype; 1218 mat->assembled = PETSC_TRUE; 1219 mat->insertmode = NOT_SET_VALUES; 1220 mat->preallocated = PETSC_TRUE; 1221 1222 a->size = oldmat->size; 1223 a->rank = oldmat->rank; 1224 a->donotstash = oldmat->donotstash; 1225 a->roworiented = oldmat->roworiented; 1226 a->rowindices = NULL; 1227 a->rowvalues = NULL; 1228 a->getrowactive = PETSC_FALSE; 1229 1230 PetscCall(PetscLayoutReference(matin->rmap, &mat->rmap)); 1231 PetscCall(PetscLayoutReference(matin->cmap, &mat->cmap)); 1232 1233 if (oldmat->colmap) { 1234 #if defined(PETSC_USE_CTABLE) 1235 PetscCall(PetscHMapIDuplicate(oldmat->colmap, &a->colmap)); 1236 #else 1237 PetscCall(PetscMalloc1(mat->cmap->N, &a->colmap)); 1238 PetscCall(PetscArraycpy(a->colmap, oldmat->colmap, mat->cmap->N)); 1239 #endif 1240 } else a->colmap = NULL; 1241 if (oldmat->garray) { 1242 PetscInt len; 1243 len = oldmat->B->cmap->n; 1244 PetscCall(PetscMalloc1(len + 1, &a->garray)); 1245 if (len) PetscCall(PetscArraycpy(a->garray, oldmat->garray, len)); 1246 } else a->garray = NULL; 1247 1248 PetscCall(VecDuplicate(oldmat->lvec, &a->lvec)); 1249 PetscCall(VecScatterCopy(oldmat->Mvctx, &a->Mvctx)); 1250 PetscCall(MatDuplicate(oldmat->A, cpvalues, &a->A)); 1251 PetscCall(MatDuplicate(oldmat->B, cpvalues, &a->B)); 1252 PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist, &((PetscObject)mat)->qlist)); 1253 *newmat = mat; 1254 PetscFunctionReturn(PETSC_SUCCESS); 1255 } 1256 1257 /*@C 1258 MatMPISELLSetPreallocation - Preallocates memory for a `MATMPISELL` sparse parallel matrix in sell format. 1259 For good matrix assembly performance the user should preallocate the matrix storage by 1260 setting the parameters `d_nz` (or `d_nnz`) and `o_nz` (or `o_nnz`). 1261 1262 Collective 1263 1264 Input Parameters: 1265 + B - the matrix 1266 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1267 (same value is used for all local rows) 1268 . d_nnz - array containing the number of nonzeros in the various rows of the 1269 DIAGONAL portion of the local submatrix (possibly different for each row) 1270 or NULL (`PETSC_NULL_INTEGER` in Fortran), if `d_nz` is used to specify the nonzero structure. 1271 The size of this array is equal to the number of local rows, i.e 'm'. 1272 For matrices that will be factored, you must leave room for (and set) 1273 the diagonal entry even if it is zero. 1274 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1275 submatrix (same value is used for all local rows). 1276 - o_nnz - array containing the number of nonzeros in the various rows of the 1277 OFF-DIAGONAL portion of the local submatrix (possibly different for 1278 each row) or NULL (`PETSC_NULL_INTEGER` in Fortran), if `o_nz` is used to specify the nonzero 1279 structure. The size of this array is equal to the number 1280 of local rows, i.e 'm'. 1281 1282 Example usage: 1283 Consider the following 8x8 matrix with 34 non-zero values, that is 1284 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1285 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1286 as follows 1287 1288 .vb 1289 1 2 0 | 0 3 0 | 0 4 1290 Proc0 0 5 6 | 7 0 0 | 8 0 1291 9 0 10 | 11 0 0 | 12 0 1292 ------------------------------------- 1293 13 0 14 | 15 16 17 | 0 0 1294 Proc1 0 18 0 | 19 20 21 | 0 0 1295 0 0 0 | 22 23 0 | 24 0 1296 ------------------------------------- 1297 Proc2 25 26 27 | 0 0 28 | 29 0 1298 30 0 0 | 31 32 33 | 0 34 1299 .ve 1300 1301 This can be represented as a collection of submatrices as 1302 1303 .vb 1304 A B C 1305 D E F 1306 G H I 1307 .ve 1308 1309 Where the submatrices A,B,C are owned by proc0, D,E,F are 1310 owned by proc1, G,H,I are owned by proc2. 1311 1312 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1313 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1314 The 'M','N' parameters are 8,8, and have the same values on all procs. 1315 1316 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1317 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1318 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1319 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1320 part as `MATSEQSELL` matrices. For example, proc1 will store [E] as a `MATSEQSELL` 1321 matrix, ans [DF] as another SeqSELL matrix. 1322 1323 When `d_nz`, `o_nz` parameters are specified, `d_nz` storage elements are 1324 allocated for every row of the local diagonal submatrix, and o_nz 1325 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1326 One way to choose `d_nz` and `o_nz` is to use the max nonzerors per local 1327 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1328 In this case, the values of d_nz,o_nz are 1329 .vb 1330 proc0 dnz = 2, o_nz = 2 1331 proc1 dnz = 3, o_nz = 2 1332 proc2 dnz = 1, o_nz = 4 1333 .ve 1334 We are allocating m*(d_nz+o_nz) storage locations for every proc. This 1335 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1336 for proc3. i.e we are using 12+15+10=37 storage locations to store 1337 34 values. 1338 1339 When `d_nnz`, `o_nnz` parameters are specified, the storage is specified 1340 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1341 In the above case the values for d_nnz,o_nnz are 1342 .vb 1343 proc0 d_nnz = [2,2,2] and o_nnz = [2,2,2] 1344 proc1 d_nnz = [3,3,2] and o_nnz = [2,1,1] 1345 proc2 d_nnz = [1,1] and o_nnz = [4,4] 1346 .ve 1347 Here the space allocated is according to nz (or maximum values in the nnz 1348 if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37 1349 1350 Level: intermediate 1351 1352 Notes: 1353 If the *_nnz parameter is given then the *_nz parameter is ignored 1354 1355 The stored row and column indices begin with zero. 1356 1357 The parallel matrix is partitioned such that the first m0 rows belong to 1358 process 0, the next m1 rows belong to process 1, the next m2 rows belong 1359 to process 2 etc.. where m0,m1,m2... are the input parameter 'm'. 1360 1361 The DIAGONAL portion of the local submatrix of a processor can be defined 1362 as the submatrix which is obtained by extraction the part corresponding to 1363 the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the 1364 first row that belongs to the processor, r2 is the last row belonging to 1365 the this processor, and c1-c2 is range of indices of the local part of a 1366 vector suitable for applying the matrix to. This is an mxn matrix. In the 1367 common case of a square matrix, the row and column ranges are the same and 1368 the DIAGONAL part is also square. The remaining portion of the local 1369 submatrix (mxN) constitute the OFF-DIAGONAL portion. 1370 1371 If `o_nnz`, `d_nnz` are specified, then `o_nz`, and `d_nz` are ignored. 1372 1373 You can call `MatGetInfo()` to get information on how effective the preallocation was; 1374 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1375 You can also run with the option -info and look for messages with the string 1376 malloc in them to see if additional memory allocation was needed. 1377 1378 .seealso: `Mat`, `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatCreatesell()`, 1379 `MATMPISELL`, `MatGetInfo()`, `PetscSplitOwnership()`, `MATSELL` 1380 @*/ 1381 PetscErrorCode MatMPISELLSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[]) 1382 { 1383 PetscFunctionBegin; 1384 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 1385 PetscValidType(B, 1); 1386 PetscTryMethod(B, "MatMPISELLSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz)); 1387 PetscFunctionReturn(PETSC_SUCCESS); 1388 } 1389 1390 /*MC 1391 MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices, 1392 based on the sliced Ellpack format 1393 1394 Options Database Key: 1395 . -mat_type sell - sets the matrix type to `MATSELL` during a call to `MatSetFromOptions()` 1396 1397 Level: beginner 1398 1399 .seealso: `Mat`, `MatCreateSell()`, `MATSEQSELL`, `MATSELL`, `MATSEQAIJ`, `MATAIJ`, `MATMPIAIJ` 1400 M*/ 1401 1402 /*@C 1403 MatCreateSELL - Creates a sparse parallel matrix in `MATSELL` format. 1404 1405 Collective 1406 1407 Input Parameters: 1408 + comm - MPI communicator 1409 . m - number of local rows (or `PETSC_DECIDE` to have calculated if M is given) 1410 This value should be the same as the local size used in creating the 1411 y vector for the matrix-vector product y = Ax. 1412 . n - This value should be the same as the local size used in creating the 1413 x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have 1414 calculated if `N` is given) For square matrices n is almost always `m`. 1415 . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given) 1416 . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given) 1417 . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix 1418 (same value is used for all local rows) 1419 . d_rlen - array containing the number of nonzeros in the various rows of the 1420 DIAGONAL portion of the local submatrix (possibly different for each row) 1421 or `NULL`, if d_rlenmax is used to specify the nonzero structure. 1422 The size of this array is equal to the number of local rows, i.e `m`. 1423 . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local 1424 submatrix (same value is used for all local rows). 1425 - o_rlen - array containing the number of nonzeros in the various rows of the 1426 OFF-DIAGONAL portion of the local submatrix (possibly different for 1427 each row) or `NULL`, if `o_rlenmax` is used to specify the nonzero 1428 structure. The size of this array is equal to the number 1429 of local rows, i.e `m`. 1430 1431 Output Parameter: 1432 . A - the matrix 1433 1434 Options Database Key: 1435 - -mat_sell_oneindex - Internally use indexing starting at 1 1436 rather than 0. When calling `MatSetValues()`, 1437 the user still MUST index entries starting at 0! 1438 1439 Example: 1440 Consider the following 8x8 matrix with 34 non-zero values, that is 1441 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1442 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1443 as follows 1444 1445 .vb 1446 1 2 0 | 0 3 0 | 0 4 1447 Proc0 0 5 6 | 7 0 0 | 8 0 1448 9 0 10 | 11 0 0 | 12 0 1449 ------------------------------------- 1450 13 0 14 | 15 16 17 | 0 0 1451 Proc1 0 18 0 | 19 20 21 | 0 0 1452 0 0 0 | 22 23 0 | 24 0 1453 ------------------------------------- 1454 Proc2 25 26 27 | 0 0 28 | 29 0 1455 30 0 0 | 31 32 33 | 0 34 1456 .ve 1457 1458 This can be represented as a collection of submatrices as 1459 .vb 1460 A B C 1461 D E F 1462 G H I 1463 .ve 1464 1465 Where the submatrices A,B,C are owned by proc0, D,E,F are 1466 owned by proc1, G,H,I are owned by proc2. 1467 1468 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1469 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1470 The 'M','N' parameters are 8,8, and have the same values on all procs. 1471 1472 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1473 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1474 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1475 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1476 part as `MATSEQSELL` matrices. For example, proc1 will store [E] as a `MATSEQSELL` 1477 matrix, ans [DF] as another `MATSEQSELL` matrix. 1478 1479 When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are 1480 allocated for every row of the local diagonal submatrix, and o_rlenmax 1481 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1482 One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local 1483 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1484 In this case, the values of d_rlenmax,o_rlenmax are 1485 .vb 1486 proc0 - d_rlenmax = 2, o_rlenmax = 2 1487 proc1 - d_rlenmax = 3, o_rlenmax = 2 1488 proc2 - d_rlenmax = 1, o_rlenmax = 4 1489 .ve 1490 We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This 1491 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1492 for proc3. i.e we are using 12+15+10=37 storage locations to store 1493 34 values. 1494 1495 When `d_rlen`, `o_rlen` parameters are specified, the storage is specified 1496 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1497 In the above case the values for `d_nnz`, `o_nnz` are 1498 .vb 1499 proc0 - d_nnz = [2,2,2] and o_nnz = [2,2,2] 1500 proc1 - d_nnz = [3,3,2] and o_nnz = [2,1,1] 1501 proc2 - d_nnz = [1,1] and o_nnz = [4,4] 1502 .ve 1503 Here the space allocated is still 37 though there are 34 nonzeros because 1504 the allocation is always done according to rlenmax. 1505 1506 Level: intermediate 1507 1508 Notes: 1509 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`, 1510 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1511 [MatXXXXSetPreallocation() is, for example, `MatSeqSELLSetPreallocation()`] 1512 1513 If the *_rlen parameter is given then the *_rlenmax parameter is ignored 1514 1515 `m`, `n`, `M`, `N` parameters specify the size of the matrix, and its partitioning across 1516 processors, while `d_rlenmax`, `d_rlen`, `o_rlenmax` , `o_rlen` parameters specify the approximate 1517 storage requirements for this matrix. 1518 1519 If `PETSC_DECIDE` or `PETSC_DETERMINE` is used for a particular argument on one 1520 processor than it must be used on all processors that share the object for 1521 that argument. 1522 1523 The user MUST specify either the local or global matrix dimensions 1524 (possibly both). 1525 1526 The parallel matrix is partitioned across processors such that the 1527 first m0 rows belong to process 0, the next m1 rows belong to 1528 process 1, the next m2 rows belong to process 2 etc.. where 1529 m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores 1530 values corresponding to [`m` x `N`] submatrix. 1531 1532 The columns are logically partitioned with the n0 columns belonging 1533 to 0th partition, the next n1 columns belonging to the next 1534 partition etc.. where n0,n1,n2... are the input parameter `n`. 1535 1536 The DIAGONAL portion of the local submatrix on any given processor 1537 is the submatrix corresponding to the rows and columns `m`, `n` 1538 corresponding to the given processor. i.e diagonal matrix on 1539 process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1] 1540 etc. The remaining portion of the local submatrix [m x (N-n)] 1541 constitute the OFF-DIAGONAL portion. The example below better 1542 illustrates this concept. 1543 1544 For a square global matrix we define each processor's diagonal portion 1545 to be its local rows and the corresponding columns (a square submatrix); 1546 each processor's off-diagonal portion encompasses the remainder of the 1547 local matrix (a rectangular submatrix). 1548 1549 If `o_rlen`, `d_rlen` are specified, then `o_rlenmax`, and `d_rlenmax` are ignored. 1550 1551 When calling this routine with a single process communicator, a matrix of 1552 type `MATSEQSELL` is returned. If a matrix of type `MATMPISELL` is desired for this 1553 type of communicator, use the construction mechanism 1554 .vb 1555 MatCreate(...,&A); 1556 MatSetType(A,MATMPISELL); 1557 MatSetSizes(A, m,n,M,N); 1558 MatMPISELLSetPreallocation(A,...); 1559 .ve 1560 1561 .seealso: `Mat`, `MATSELL`, `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatMPISELLSetPreallocation()`, `MatMPISELLSetPreallocationSELL()`, 1562 `MATMPISELL`, `MatCreateMPISELLWithArrays()` 1563 @*/ 1564 PetscErrorCode MatCreateSELL(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt d_rlenmax, const PetscInt d_rlen[], PetscInt o_rlenmax, const PetscInt o_rlen[], Mat *A) 1565 { 1566 PetscMPIInt size; 1567 1568 PetscFunctionBegin; 1569 PetscCall(MatCreate(comm, A)); 1570 PetscCall(MatSetSizes(*A, m, n, M, N)); 1571 PetscCallMPI(MPI_Comm_size(comm, &size)); 1572 if (size > 1) { 1573 PetscCall(MatSetType(*A, MATMPISELL)); 1574 PetscCall(MatMPISELLSetPreallocation(*A, d_rlenmax, d_rlen, o_rlenmax, o_rlen)); 1575 } else { 1576 PetscCall(MatSetType(*A, MATSEQSELL)); 1577 PetscCall(MatSeqSELLSetPreallocation(*A, d_rlenmax, d_rlen)); 1578 } 1579 PetscFunctionReturn(PETSC_SUCCESS); 1580 } 1581 1582 PetscErrorCode MatMPISELLGetSeqSELL(Mat A, Mat *Ad, Mat *Ao, const PetscInt *colmap[]) 1583 { 1584 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 1585 PetscBool flg; 1586 1587 PetscFunctionBegin; 1588 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISELL, &flg)); 1589 PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "This function requires a MATMPISELL matrix as input"); 1590 if (Ad) *Ad = a->A; 1591 if (Ao) *Ao = a->B; 1592 if (colmap) *colmap = a->garray; 1593 PetscFunctionReturn(PETSC_SUCCESS); 1594 } 1595 1596 /*@C 1597 MatMPISELLGetLocalMatCondensed - Creates a `MATSEQSELL` matrix from an `MATMPISELL` matrix by taking all its local rows and NON-ZERO columns 1598 1599 Not Collective 1600 1601 Input Parameters: 1602 + A - the matrix 1603 . scall - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX` 1604 . row - index sets of rows to extract (or `NULL`) 1605 - col - index sets of columns to extract (or `NULL`) 1606 1607 Output Parameter: 1608 . A_loc - the local sequential matrix generated 1609 1610 Level: developer 1611 1612 .seealso: `Mat`, `MATSEQSELL`, `MATMPISELL`, `MatGetOwnershipRange()`, `MatMPISELLGetLocalMat()` 1613 @*/ 1614 PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A, MatReuse scall, IS *row, IS *col, Mat *A_loc) 1615 { 1616 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 1617 PetscInt i, start, end, ncols, nzA, nzB, *cmap, imark, *idx; 1618 IS isrowa, iscola; 1619 Mat *aloc; 1620 PetscBool match; 1621 1622 PetscFunctionBegin; 1623 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISELL, &match)); 1624 PetscCheck(match, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Requires MATMPISELL matrix as input"); 1625 PetscCall(PetscLogEventBegin(MAT_Getlocalmatcondensed, A, 0, 0, 0)); 1626 if (!row) { 1627 start = A->rmap->rstart; 1628 end = A->rmap->rend; 1629 PetscCall(ISCreateStride(PETSC_COMM_SELF, end - start, start, 1, &isrowa)); 1630 } else { 1631 isrowa = *row; 1632 } 1633 if (!col) { 1634 start = A->cmap->rstart; 1635 cmap = a->garray; 1636 nzA = a->A->cmap->n; 1637 nzB = a->B->cmap->n; 1638 PetscCall(PetscMalloc1(nzA + nzB, &idx)); 1639 ncols = 0; 1640 for (i = 0; i < nzB; i++) { 1641 if (cmap[i] < start) idx[ncols++] = cmap[i]; 1642 else break; 1643 } 1644 imark = i; 1645 for (i = 0; i < nzA; i++) idx[ncols++] = start + i; 1646 for (i = imark; i < nzB; i++) idx[ncols++] = cmap[i]; 1647 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncols, idx, PETSC_OWN_POINTER, &iscola)); 1648 } else { 1649 iscola = *col; 1650 } 1651 if (scall != MAT_INITIAL_MATRIX) { 1652 PetscCall(PetscMalloc1(1, &aloc)); 1653 aloc[0] = *A_loc; 1654 } 1655 PetscCall(MatCreateSubMatrices(A, 1, &isrowa, &iscola, scall, &aloc)); 1656 *A_loc = aloc[0]; 1657 PetscCall(PetscFree(aloc)); 1658 if (!row) PetscCall(ISDestroy(&isrowa)); 1659 if (!col) PetscCall(ISDestroy(&iscola)); 1660 PetscCall(PetscLogEventEnd(MAT_Getlocalmatcondensed, A, 0, 0, 0)); 1661 PetscFunctionReturn(PETSC_SUCCESS); 1662 } 1663 1664 #include <../src/mat/impls/aij/mpi/mpiaij.h> 1665 1666 PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 1667 { 1668 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 1669 Mat B; 1670 Mat_MPIAIJ *b; 1671 1672 PetscFunctionBegin; 1673 PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled"); 1674 1675 if (reuse == MAT_REUSE_MATRIX) { 1676 B = *newmat; 1677 } else { 1678 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 1679 PetscCall(MatSetType(B, MATMPIAIJ)); 1680 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 1681 PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs)); 1682 PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL)); 1683 PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL)); 1684 } 1685 b = (Mat_MPIAIJ *)B->data; 1686 1687 if (reuse == MAT_REUSE_MATRIX) { 1688 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A)); 1689 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B)); 1690 } else { 1691 PetscCall(MatDestroy(&b->A)); 1692 PetscCall(MatDestroy(&b->B)); 1693 PetscCall(MatDisAssemble_MPISELL(A)); 1694 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A)); 1695 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B)); 1696 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 1697 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 1698 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 1699 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 1700 } 1701 1702 if (reuse == MAT_INPLACE_MATRIX) { 1703 PetscCall(MatHeaderReplace(A, &B)); 1704 } else { 1705 *newmat = B; 1706 } 1707 PetscFunctionReturn(PETSC_SUCCESS); 1708 } 1709 1710 PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 1711 { 1712 Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data; 1713 Mat B; 1714 Mat_MPISELL *b; 1715 1716 PetscFunctionBegin; 1717 PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled"); 1718 1719 if (reuse == MAT_REUSE_MATRIX) { 1720 B = *newmat; 1721 } else { 1722 Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)a->A->data, *Ba = (Mat_SeqAIJ *)a->B->data; 1723 PetscInt i, d_nz = 0, o_nz = 0, m = A->rmap->N, n = A->cmap->N, lm = A->rmap->n, ln = A->cmap->n; 1724 PetscInt *d_nnz, *o_nnz; 1725 PetscCall(PetscMalloc2(lm, &d_nnz, lm, &o_nnz)); 1726 for (i = 0; i < lm; i++) { 1727 d_nnz[i] = Aa->i[i + 1] - Aa->i[i]; 1728 o_nnz[i] = Ba->i[i + 1] - Ba->i[i]; 1729 if (d_nnz[i] > d_nz) d_nz = d_nnz[i]; 1730 if (o_nnz[i] > o_nz) o_nz = o_nnz[i]; 1731 } 1732 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 1733 PetscCall(MatSetType(B, MATMPISELL)); 1734 PetscCall(MatSetSizes(B, lm, ln, m, n)); 1735 PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs)); 1736 PetscCall(MatSeqSELLSetPreallocation(B, d_nz, d_nnz)); 1737 PetscCall(MatMPISELLSetPreallocation(B, d_nz, d_nnz, o_nz, o_nnz)); 1738 PetscCall(PetscFree2(d_nnz, o_nnz)); 1739 } 1740 b = (Mat_MPISELL *)B->data; 1741 1742 if (reuse == MAT_REUSE_MATRIX) { 1743 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A)); 1744 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B)); 1745 } else { 1746 PetscCall(MatDestroy(&b->A)); 1747 PetscCall(MatDestroy(&b->B)); 1748 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A)); 1749 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B)); 1750 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 1751 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 1752 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 1753 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 1754 } 1755 1756 if (reuse == MAT_INPLACE_MATRIX) { 1757 PetscCall(MatHeaderReplace(A, &B)); 1758 } else { 1759 *newmat = B; 1760 } 1761 PetscFunctionReturn(PETSC_SUCCESS); 1762 } 1763 1764 PetscErrorCode MatSOR_MPISELL(Mat matin, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx) 1765 { 1766 Mat_MPISELL *mat = (Mat_MPISELL *)matin->data; 1767 Vec bb1 = NULL; 1768 1769 PetscFunctionBegin; 1770 if (flag == SOR_APPLY_UPPER) { 1771 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1772 PetscFunctionReturn(PETSC_SUCCESS); 1773 } 1774 1775 if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) PetscCall(VecDuplicate(bb, &bb1)); 1776 1777 if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) { 1778 if (flag & SOR_ZERO_INITIAL_GUESS) { 1779 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1780 its--; 1781 } 1782 1783 while (its--) { 1784 PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1785 PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1786 1787 /* update rhs: bb1 = bb - B*x */ 1788 PetscCall(VecScale(mat->lvec, -1.0)); 1789 PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1)); 1790 1791 /* local sweep */ 1792 PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_SYMMETRIC_SWEEP, fshift, lits, 1, xx)); 1793 } 1794 } else if (flag & SOR_LOCAL_FORWARD_SWEEP) { 1795 if (flag & SOR_ZERO_INITIAL_GUESS) { 1796 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1797 its--; 1798 } 1799 while (its--) { 1800 PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1801 PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1802 1803 /* update rhs: bb1 = bb - B*x */ 1804 PetscCall(VecScale(mat->lvec, -1.0)); 1805 PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1)); 1806 1807 /* local sweep */ 1808 PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_FORWARD_SWEEP, fshift, lits, 1, xx)); 1809 } 1810 } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) { 1811 if (flag & SOR_ZERO_INITIAL_GUESS) { 1812 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1813 its--; 1814 } 1815 while (its--) { 1816 PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1817 PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1818 1819 /* update rhs: bb1 = bb - B*x */ 1820 PetscCall(VecScale(mat->lvec, -1.0)); 1821 PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1)); 1822 1823 /* local sweep */ 1824 PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_BACKWARD_SWEEP, fshift, lits, 1, xx)); 1825 } 1826 } else SETERRQ(PetscObjectComm((PetscObject)matin), PETSC_ERR_SUP, "Parallel SOR not supported"); 1827 1828 PetscCall(VecDestroy(&bb1)); 1829 1830 matin->factorerrortype = mat->A->factorerrortype; 1831 PetscFunctionReturn(PETSC_SUCCESS); 1832 } 1833 1834 /*MC 1835 MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices. 1836 1837 Options Database Keys: 1838 . -mat_type mpisell - sets the matrix type to `MATMPISELL` during a call to `MatSetFromOptions()` 1839 1840 Level: beginner 1841 1842 .seealso: `Mat`, `MATSELL`, `MATSEQSELL` `MatCreateSELL()` 1843 M*/ 1844 PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B) 1845 { 1846 Mat_MPISELL *b; 1847 PetscMPIInt size; 1848 1849 PetscFunctionBegin; 1850 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 1851 PetscCall(PetscNew(&b)); 1852 B->data = (void *)b; 1853 PetscCall(PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps))); 1854 B->assembled = PETSC_FALSE; 1855 B->insertmode = NOT_SET_VALUES; 1856 b->size = size; 1857 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B), &b->rank)); 1858 /* build cache for off array entries formed */ 1859 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash)); 1860 1861 b->donotstash = PETSC_FALSE; 1862 b->colmap = NULL; 1863 b->garray = NULL; 1864 b->roworiented = PETSC_TRUE; 1865 1866 /* stuff used for matrix vector multiply */ 1867 b->lvec = NULL; 1868 b->Mvctx = NULL; 1869 1870 /* stuff for MatGetRow() */ 1871 b->rowindices = NULL; 1872 b->rowvalues = NULL; 1873 b->getrowactive = PETSC_FALSE; 1874 1875 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_MPISELL)); 1876 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_MPISELL)); 1877 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsTranspose_C", MatIsTranspose_MPISELL)); 1878 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPISELLSetPreallocation_C", MatMPISELLSetPreallocation_MPISELL)); 1879 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisell_mpiaij_C", MatConvert_MPISELL_MPIAIJ)); 1880 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDiagonalScaleLocal_C", MatDiagonalScaleLocal_MPISELL)); 1881 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATMPISELL)); 1882 PetscFunctionReturn(PETSC_SUCCESS); 1883 } 1884