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 PetscCall(PetscLogObjectState((PetscObject)mat, "Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT, mat->rmap->N, mat->cmap->N)); 508 PetscCall(MatStashDestroy_Private(&mat->stash)); 509 PetscCall(VecDestroy(&sell->diag)); 510 PetscCall(MatDestroy(&sell->A)); 511 PetscCall(MatDestroy(&sell->B)); 512 #if defined(PETSC_USE_CTABLE) 513 PetscCall(PetscHMapIDestroy(&sell->colmap)); 514 #else 515 PetscCall(PetscFree(sell->colmap)); 516 #endif 517 PetscCall(PetscFree(sell->garray)); 518 PetscCall(VecDestroy(&sell->lvec)); 519 PetscCall(VecScatterDestroy(&sell->Mvctx)); 520 PetscCall(PetscFree2(sell->rowvalues, sell->rowindices)); 521 PetscCall(PetscFree(sell->ld)); 522 PetscCall(PetscFree(mat->data)); 523 524 PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL)); 525 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatStoreValues_C", NULL)); 526 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatRetrieveValues_C", NULL)); 527 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatIsTranspose_C", NULL)); 528 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPISELLSetPreallocation_C", NULL)); 529 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisell_mpiaij_C", NULL)); 530 #if defined(PETSC_HAVE_CUDA) 531 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisell_mpisellcuda_C", NULL)); 532 #endif 533 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalScaleLocal_C", NULL)); 534 PetscFunctionReturn(PETSC_SUCCESS); 535 } 536 537 #include <petscdraw.h> 538 PetscErrorCode MatView_MPISELL_ASCIIorDraworSocket(Mat mat, PetscViewer viewer) 539 { 540 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 541 PetscMPIInt rank = sell->rank, size = sell->size; 542 PetscBool isdraw, iascii, isbinary; 543 PetscViewer sviewer; 544 PetscViewerFormat format; 545 546 PetscFunctionBegin; 547 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 548 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 549 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 550 if (iascii) { 551 PetscCall(PetscViewerGetFormat(viewer, &format)); 552 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 553 MatInfo info; 554 PetscInt *inodes; 555 556 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat), &rank)); 557 PetscCall(MatGetInfo(mat, MAT_LOCAL, &info)); 558 PetscCall(MatInodeGetInodeSizes(sell->A, NULL, &inodes, NULL)); 559 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 560 if (!inodes) { 561 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, 562 (PetscInt)info.nz_allocated, (PetscInt)info.memory)); 563 } else { 564 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, 565 (PetscInt)info.nz_allocated, (PetscInt)info.memory)); 566 } 567 PetscCall(MatGetInfo(sell->A, MAT_LOCAL, &info)); 568 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] on-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used)); 569 PetscCall(MatGetInfo(sell->B, MAT_LOCAL, &info)); 570 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] off-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used)); 571 PetscCall(PetscViewerFlush(viewer)); 572 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 573 PetscCall(PetscViewerASCIIPrintf(viewer, "Information on VecScatter used in matrix-vector product: \n")); 574 PetscCall(VecScatterView(sell->Mvctx, viewer)); 575 PetscFunctionReturn(PETSC_SUCCESS); 576 } else if (format == PETSC_VIEWER_ASCII_INFO) { 577 PetscInt inodecount, inodelimit, *inodes; 578 PetscCall(MatInodeGetInodeSizes(sell->A, &inodecount, &inodes, &inodelimit)); 579 if (inodes) { 580 PetscCall(PetscViewerASCIIPrintf(viewer, "using I-node (on process 0) routines: found %" PetscInt_FMT " nodes, limit used is %" PetscInt_FMT "\n", inodecount, inodelimit)); 581 } else { 582 PetscCall(PetscViewerASCIIPrintf(viewer, "not using I-node (on process 0) routines\n")); 583 } 584 PetscFunctionReturn(PETSC_SUCCESS); 585 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 586 PetscFunctionReturn(PETSC_SUCCESS); 587 } 588 } else if (isbinary) { 589 if (size == 1) { 590 PetscCall(PetscObjectSetName((PetscObject)sell->A, ((PetscObject)mat)->name)); 591 PetscCall(MatView(sell->A, viewer)); 592 } else { 593 /* PetscCall(MatView_MPISELL_Binary(mat,viewer)); */ 594 } 595 PetscFunctionReturn(PETSC_SUCCESS); 596 } else if (isdraw) { 597 PetscDraw draw; 598 PetscBool isnull; 599 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 600 PetscCall(PetscDrawIsNull(draw, &isnull)); 601 if (isnull) PetscFunctionReturn(PETSC_SUCCESS); 602 } 603 604 { 605 /* assemble the entire matrix onto first processor. */ 606 Mat A; 607 Mat_SeqSELL *Aloc; 608 PetscInt M = mat->rmap->N, N = mat->cmap->N, *acolidx, row, col, i, j; 609 MatScalar *aval; 610 PetscBool isnonzero; 611 612 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &A)); 613 if (rank == 0) { 614 PetscCall(MatSetSizes(A, M, N, M, N)); 615 } else { 616 PetscCall(MatSetSizes(A, 0, 0, M, N)); 617 } 618 /* This is just a temporary matrix, so explicitly using MATMPISELL is probably best */ 619 PetscCall(MatSetType(A, MATMPISELL)); 620 PetscCall(MatMPISELLSetPreallocation(A, 0, NULL, 0, NULL)); 621 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE)); 622 623 /* copy over the A part */ 624 Aloc = (Mat_SeqSELL *)sell->A->data; 625 acolidx = Aloc->colidx; 626 aval = Aloc->val; 627 for (i = 0; i < Aloc->totalslices; i++) { /* loop over slices */ 628 for (j = Aloc->sliidx[i]; j < Aloc->sliidx[i + 1]; j++) { 629 isnonzero = (PetscBool)((j - Aloc->sliidx[i]) / Aloc->sliceheight < Aloc->rlen[i * Aloc->sliceheight + j % Aloc->sliceheight]); 630 if (isnonzero) { /* check the mask bit */ 631 row = i * Aloc->sliceheight + j % Aloc->sliceheight + mat->rmap->rstart; 632 col = *acolidx + mat->rmap->rstart; 633 PetscCall(MatSetValues(A, 1, &row, 1, &col, aval, INSERT_VALUES)); 634 } 635 aval++; 636 acolidx++; 637 } 638 } 639 640 /* copy over the B part */ 641 Aloc = (Mat_SeqSELL *)sell->B->data; 642 acolidx = Aloc->colidx; 643 aval = Aloc->val; 644 for (i = 0; i < Aloc->totalslices; i++) { 645 for (j = Aloc->sliidx[i]; j < Aloc->sliidx[i + 1]; j++) { 646 isnonzero = (PetscBool)((j - Aloc->sliidx[i]) / Aloc->sliceheight < Aloc->rlen[i * Aloc->sliceheight + j % Aloc->sliceheight]); 647 if (isnonzero) { 648 row = i * Aloc->sliceheight + j % Aloc->sliceheight + mat->rmap->rstart; 649 col = sell->garray[*acolidx]; 650 PetscCall(MatSetValues(A, 1, &row, 1, &col, aval, INSERT_VALUES)); 651 } 652 aval++; 653 acolidx++; 654 } 655 } 656 657 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 658 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 659 /* 660 Everyone has to call to draw the matrix since the graphics waits are 661 synchronized across all processors that share the PetscDraw object 662 */ 663 PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer)); 664 if (rank == 0) { 665 PetscCall(PetscObjectSetName((PetscObject)((Mat_MPISELL *)(A->data))->A, ((PetscObject)mat)->name)); 666 PetscCall(MatView_SeqSELL(((Mat_MPISELL *)(A->data))->A, sviewer)); 667 } 668 PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer)); 669 PetscCall(PetscViewerFlush(viewer)); 670 PetscCall(MatDestroy(&A)); 671 } 672 PetscFunctionReturn(PETSC_SUCCESS); 673 } 674 675 PetscErrorCode MatView_MPISELL(Mat mat, PetscViewer viewer) 676 { 677 PetscBool iascii, isdraw, issocket, isbinary; 678 679 PetscFunctionBegin; 680 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 681 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 682 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 683 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSOCKET, &issocket)); 684 if (iascii || isdraw || isbinary || issocket) PetscCall(MatView_MPISELL_ASCIIorDraworSocket(mat, viewer)); 685 PetscFunctionReturn(PETSC_SUCCESS); 686 } 687 688 PetscErrorCode MatGetGhosts_MPISELL(Mat mat, PetscInt *nghosts, const PetscInt *ghosts[]) 689 { 690 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 691 692 PetscFunctionBegin; 693 PetscCall(MatGetSize(sell->B, NULL, nghosts)); 694 if (ghosts) *ghosts = sell->garray; 695 PetscFunctionReturn(PETSC_SUCCESS); 696 } 697 698 PetscErrorCode MatGetInfo_MPISELL(Mat matin, MatInfoType flag, MatInfo *info) 699 { 700 Mat_MPISELL *mat = (Mat_MPISELL *)matin->data; 701 Mat A = mat->A, B = mat->B; 702 PetscLogDouble isend[5], irecv[5]; 703 704 PetscFunctionBegin; 705 info->block_size = 1.0; 706 PetscCall(MatGetInfo(A, MAT_LOCAL, info)); 707 708 isend[0] = info->nz_used; 709 isend[1] = info->nz_allocated; 710 isend[2] = info->nz_unneeded; 711 isend[3] = info->memory; 712 isend[4] = info->mallocs; 713 714 PetscCall(MatGetInfo(B, MAT_LOCAL, info)); 715 716 isend[0] += info->nz_used; 717 isend[1] += info->nz_allocated; 718 isend[2] += info->nz_unneeded; 719 isend[3] += info->memory; 720 isend[4] += info->mallocs; 721 if (flag == MAT_LOCAL) { 722 info->nz_used = isend[0]; 723 info->nz_allocated = isend[1]; 724 info->nz_unneeded = isend[2]; 725 info->memory = isend[3]; 726 info->mallocs = isend[4]; 727 } else if (flag == MAT_GLOBAL_MAX) { 728 PetscCall(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)matin))); 729 730 info->nz_used = irecv[0]; 731 info->nz_allocated = irecv[1]; 732 info->nz_unneeded = irecv[2]; 733 info->memory = irecv[3]; 734 info->mallocs = irecv[4]; 735 } else if (flag == MAT_GLOBAL_SUM) { 736 PetscCall(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)matin))); 737 738 info->nz_used = irecv[0]; 739 info->nz_allocated = irecv[1]; 740 info->nz_unneeded = irecv[2]; 741 info->memory = irecv[3]; 742 info->mallocs = irecv[4]; 743 } 744 info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */ 745 info->fill_ratio_needed = 0; 746 info->factor_mallocs = 0; 747 PetscFunctionReturn(PETSC_SUCCESS); 748 } 749 750 PetscErrorCode MatSetOption_MPISELL(Mat A, MatOption op, PetscBool flg) 751 { 752 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 753 754 PetscFunctionBegin; 755 switch (op) { 756 case MAT_NEW_NONZERO_LOCATIONS: 757 case MAT_NEW_NONZERO_ALLOCATION_ERR: 758 case MAT_UNUSED_NONZERO_LOCATION_ERR: 759 case MAT_KEEP_NONZERO_PATTERN: 760 case MAT_NEW_NONZERO_LOCATION_ERR: 761 case MAT_USE_INODES: 762 case MAT_IGNORE_ZERO_ENTRIES: 763 MatCheckPreallocated(A, 1); 764 PetscCall(MatSetOption(a->A, op, flg)); 765 PetscCall(MatSetOption(a->B, op, flg)); 766 break; 767 case MAT_ROW_ORIENTED: 768 MatCheckPreallocated(A, 1); 769 a->roworiented = flg; 770 771 PetscCall(MatSetOption(a->A, op, flg)); 772 PetscCall(MatSetOption(a->B, op, flg)); 773 break; 774 case MAT_FORCE_DIAGONAL_ENTRIES: 775 case MAT_SORTED_FULL: 776 PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op])); 777 break; 778 case MAT_IGNORE_OFF_PROC_ENTRIES: 779 a->donotstash = flg; 780 break; 781 case MAT_SPD: 782 case MAT_SPD_ETERNAL: 783 break; 784 case MAT_SYMMETRIC: 785 MatCheckPreallocated(A, 1); 786 PetscCall(MatSetOption(a->A, op, flg)); 787 break; 788 case MAT_STRUCTURALLY_SYMMETRIC: 789 MatCheckPreallocated(A, 1); 790 PetscCall(MatSetOption(a->A, op, flg)); 791 break; 792 case MAT_HERMITIAN: 793 MatCheckPreallocated(A, 1); 794 PetscCall(MatSetOption(a->A, op, flg)); 795 break; 796 case MAT_SYMMETRY_ETERNAL: 797 MatCheckPreallocated(A, 1); 798 PetscCall(MatSetOption(a->A, op, flg)); 799 break; 800 case MAT_STRUCTURAL_SYMMETRY_ETERNAL: 801 MatCheckPreallocated(A, 1); 802 PetscCall(MatSetOption(a->A, op, flg)); 803 break; 804 default: 805 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %d", op); 806 } 807 PetscFunctionReturn(PETSC_SUCCESS); 808 } 809 810 PetscErrorCode MatDiagonalScale_MPISELL(Mat mat, Vec ll, Vec rr) 811 { 812 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 813 Mat a = sell->A, b = sell->B; 814 PetscInt s1, s2, s3; 815 816 PetscFunctionBegin; 817 PetscCall(MatGetLocalSize(mat, &s2, &s3)); 818 if (rr) { 819 PetscCall(VecGetLocalSize(rr, &s1)); 820 PetscCheck(s1 == s3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "right vector non-conforming local size"); 821 /* Overlap communication with computation. */ 822 PetscCall(VecScatterBegin(sell->Mvctx, rr, sell->lvec, INSERT_VALUES, SCATTER_FORWARD)); 823 } 824 if (ll) { 825 PetscCall(VecGetLocalSize(ll, &s1)); 826 PetscCheck(s1 == s2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "left vector non-conforming local size"); 827 PetscUseTypeMethod(b, diagonalscale, ll, NULL); 828 } 829 /* scale the diagonal block */ 830 PetscUseTypeMethod(a, diagonalscale, ll, rr); 831 832 if (rr) { 833 /* Do a scatter end and then right scale the off-diagonal block */ 834 PetscCall(VecScatterEnd(sell->Mvctx, rr, sell->lvec, INSERT_VALUES, SCATTER_FORWARD)); 835 PetscUseTypeMethod(b, diagonalscale, NULL, sell->lvec); 836 } 837 PetscFunctionReturn(PETSC_SUCCESS); 838 } 839 840 PetscErrorCode MatSetUnfactored_MPISELL(Mat A) 841 { 842 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 843 844 PetscFunctionBegin; 845 PetscCall(MatSetUnfactored(a->A)); 846 PetscFunctionReturn(PETSC_SUCCESS); 847 } 848 849 PetscErrorCode MatEqual_MPISELL(Mat A, Mat B, PetscBool *flag) 850 { 851 Mat_MPISELL *matB = (Mat_MPISELL *)B->data, *matA = (Mat_MPISELL *)A->data; 852 Mat a, b, c, d; 853 PetscBool flg; 854 855 PetscFunctionBegin; 856 a = matA->A; 857 b = matA->B; 858 c = matB->A; 859 d = matB->B; 860 861 PetscCall(MatEqual(a, c, &flg)); 862 if (flg) PetscCall(MatEqual(b, d, &flg)); 863 PetscCall(MPIU_Allreduce(&flg, flag, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 864 PetscFunctionReturn(PETSC_SUCCESS); 865 } 866 867 PetscErrorCode MatCopy_MPISELL(Mat A, Mat B, MatStructure str) 868 { 869 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 870 Mat_MPISELL *b = (Mat_MPISELL *)B->data; 871 872 PetscFunctionBegin; 873 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 874 if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) { 875 /* because of the column compression in the off-processor part of the matrix a->B, 876 the number of columns in a->B and b->B may be different, hence we cannot call 877 the MatCopy() directly on the two parts. If need be, we can provide a more 878 efficient copy than the MatCopy_Basic() by first uncompressing the a->B matrices 879 then copying the submatrices */ 880 PetscCall(MatCopy_Basic(A, B, str)); 881 } else { 882 PetscCall(MatCopy(a->A, b->A, str)); 883 PetscCall(MatCopy(a->B, b->B, str)); 884 } 885 PetscFunctionReturn(PETSC_SUCCESS); 886 } 887 888 PetscErrorCode MatSetUp_MPISELL(Mat A) 889 { 890 PetscFunctionBegin; 891 PetscCall(MatMPISELLSetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL)); 892 PetscFunctionReturn(PETSC_SUCCESS); 893 } 894 895 extern PetscErrorCode MatConjugate_SeqSELL(Mat); 896 897 PetscErrorCode MatConjugate_MPISELL(Mat mat) 898 { 899 PetscFunctionBegin; 900 if (PetscDefined(USE_COMPLEX)) { 901 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 902 903 PetscCall(MatConjugate_SeqSELL(sell->A)); 904 PetscCall(MatConjugate_SeqSELL(sell->B)); 905 } 906 PetscFunctionReturn(PETSC_SUCCESS); 907 } 908 909 PetscErrorCode MatRealPart_MPISELL(Mat A) 910 { 911 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 912 913 PetscFunctionBegin; 914 PetscCall(MatRealPart(a->A)); 915 PetscCall(MatRealPart(a->B)); 916 PetscFunctionReturn(PETSC_SUCCESS); 917 } 918 919 PetscErrorCode MatImaginaryPart_MPISELL(Mat A) 920 { 921 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 922 923 PetscFunctionBegin; 924 PetscCall(MatImaginaryPart(a->A)); 925 PetscCall(MatImaginaryPart(a->B)); 926 PetscFunctionReturn(PETSC_SUCCESS); 927 } 928 929 PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A, const PetscScalar **values) 930 { 931 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 932 933 PetscFunctionBegin; 934 PetscCall(MatInvertBlockDiagonal(a->A, values)); 935 A->factorerrortype = a->A->factorerrortype; 936 PetscFunctionReturn(PETSC_SUCCESS); 937 } 938 939 static PetscErrorCode MatSetRandom_MPISELL(Mat x, PetscRandom rctx) 940 { 941 Mat_MPISELL *sell = (Mat_MPISELL *)x->data; 942 943 PetscFunctionBegin; 944 PetscCall(MatSetRandom(sell->A, rctx)); 945 PetscCall(MatSetRandom(sell->B, rctx)); 946 PetscCall(MatAssemblyBegin(x, MAT_FINAL_ASSEMBLY)); 947 PetscCall(MatAssemblyEnd(x, MAT_FINAL_ASSEMBLY)); 948 PetscFunctionReturn(PETSC_SUCCESS); 949 } 950 951 PetscErrorCode MatSetFromOptions_MPISELL(Mat A, PetscOptionItems *PetscOptionsObject) 952 { 953 PetscFunctionBegin; 954 PetscOptionsHeadBegin(PetscOptionsObject, "MPISELL options"); 955 PetscOptionsHeadEnd(); 956 PetscFunctionReturn(PETSC_SUCCESS); 957 } 958 959 PetscErrorCode MatShift_MPISELL(Mat Y, PetscScalar a) 960 { 961 Mat_MPISELL *msell = (Mat_MPISELL *)Y->data; 962 Mat_SeqSELL *sell = (Mat_SeqSELL *)msell->A->data; 963 964 PetscFunctionBegin; 965 if (!Y->preallocated) { 966 PetscCall(MatMPISELLSetPreallocation(Y, 1, NULL, 0, NULL)); 967 } else if (!sell->nz) { 968 PetscInt nonew = sell->nonew; 969 PetscCall(MatSeqSELLSetPreallocation(msell->A, 1, NULL)); 970 sell->nonew = nonew; 971 } 972 PetscCall(MatShift_Basic(Y, a)); 973 PetscFunctionReturn(PETSC_SUCCESS); 974 } 975 976 PetscErrorCode MatMissingDiagonal_MPISELL(Mat A, PetscBool *missing, PetscInt *d) 977 { 978 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 979 980 PetscFunctionBegin; 981 PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only works for square matrices"); 982 PetscCall(MatMissingDiagonal(a->A, missing, d)); 983 if (d) { 984 PetscInt rstart; 985 PetscCall(MatGetOwnershipRange(A, &rstart, NULL)); 986 *d += rstart; 987 } 988 PetscFunctionReturn(PETSC_SUCCESS); 989 } 990 991 PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A, Mat *a) 992 { 993 PetscFunctionBegin; 994 *a = ((Mat_MPISELL *)A->data)->A; 995 PetscFunctionReturn(PETSC_SUCCESS); 996 } 997 998 static struct _MatOps MatOps_Values = {MatSetValues_MPISELL, 999 NULL, 1000 NULL, 1001 MatMult_MPISELL, 1002 /* 4*/ MatMultAdd_MPISELL, 1003 MatMultTranspose_MPISELL, 1004 MatMultTransposeAdd_MPISELL, 1005 NULL, 1006 NULL, 1007 NULL, 1008 /*10*/ NULL, 1009 NULL, 1010 NULL, 1011 MatSOR_MPISELL, 1012 NULL, 1013 /*15*/ MatGetInfo_MPISELL, 1014 MatEqual_MPISELL, 1015 MatGetDiagonal_MPISELL, 1016 MatDiagonalScale_MPISELL, 1017 NULL, 1018 /*20*/ MatAssemblyBegin_MPISELL, 1019 MatAssemblyEnd_MPISELL, 1020 MatSetOption_MPISELL, 1021 MatZeroEntries_MPISELL, 1022 /*24*/ NULL, 1023 NULL, 1024 NULL, 1025 NULL, 1026 NULL, 1027 /*29*/ MatSetUp_MPISELL, 1028 NULL, 1029 NULL, 1030 MatGetDiagonalBlock_MPISELL, 1031 NULL, 1032 /*34*/ MatDuplicate_MPISELL, 1033 NULL, 1034 NULL, 1035 NULL, 1036 NULL, 1037 /*39*/ NULL, 1038 NULL, 1039 NULL, 1040 MatGetValues_MPISELL, 1041 MatCopy_MPISELL, 1042 /*44*/ NULL, 1043 MatScale_MPISELL, 1044 MatShift_MPISELL, 1045 MatDiagonalSet_MPISELL, 1046 NULL, 1047 /*49*/ MatSetRandom_MPISELL, 1048 NULL, 1049 NULL, 1050 NULL, 1051 NULL, 1052 /*54*/ MatFDColoringCreate_MPIXAIJ, 1053 NULL, 1054 MatSetUnfactored_MPISELL, 1055 NULL, 1056 NULL, 1057 /*59*/ NULL, 1058 MatDestroy_MPISELL, 1059 MatView_MPISELL, 1060 NULL, 1061 NULL, 1062 /*64*/ NULL, 1063 NULL, 1064 NULL, 1065 NULL, 1066 NULL, 1067 /*69*/ NULL, 1068 NULL, 1069 NULL, 1070 NULL, 1071 NULL, 1072 NULL, 1073 /*75*/ MatFDColoringApply_AIJ, /* reuse AIJ function */ 1074 MatSetFromOptions_MPISELL, 1075 NULL, 1076 NULL, 1077 NULL, 1078 /*80*/ NULL, 1079 NULL, 1080 NULL, 1081 /*83*/ NULL, 1082 NULL, 1083 NULL, 1084 NULL, 1085 NULL, 1086 NULL, 1087 /*89*/ NULL, 1088 NULL, 1089 NULL, 1090 NULL, 1091 NULL, 1092 /*94*/ NULL, 1093 NULL, 1094 NULL, 1095 NULL, 1096 NULL, 1097 /*99*/ NULL, 1098 NULL, 1099 NULL, 1100 MatConjugate_MPISELL, 1101 NULL, 1102 /*104*/ NULL, 1103 MatRealPart_MPISELL, 1104 MatImaginaryPart_MPISELL, 1105 NULL, 1106 NULL, 1107 /*109*/ NULL, 1108 NULL, 1109 NULL, 1110 NULL, 1111 MatMissingDiagonal_MPISELL, 1112 /*114*/ NULL, 1113 NULL, 1114 MatGetGhosts_MPISELL, 1115 NULL, 1116 NULL, 1117 /*119*/ NULL, 1118 NULL, 1119 NULL, 1120 NULL, 1121 NULL, 1122 /*124*/ NULL, 1123 NULL, 1124 MatInvertBlockDiagonal_MPISELL, 1125 NULL, 1126 NULL, 1127 /*129*/ NULL, 1128 NULL, 1129 NULL, 1130 NULL, 1131 NULL, 1132 /*134*/ NULL, 1133 NULL, 1134 NULL, 1135 NULL, 1136 NULL, 1137 /*139*/ NULL, 1138 NULL, 1139 NULL, 1140 MatFDColoringSetUp_MPIXAIJ, 1141 NULL, 1142 /*144*/ NULL, 1143 NULL, 1144 NULL, 1145 NULL, 1146 NULL, 1147 NULL, 1148 /*150*/ NULL, 1149 NULL}; 1150 1151 PetscErrorCode MatStoreValues_MPISELL(Mat mat) 1152 { 1153 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 1154 1155 PetscFunctionBegin; 1156 PetscCall(MatStoreValues(sell->A)); 1157 PetscCall(MatStoreValues(sell->B)); 1158 PetscFunctionReturn(PETSC_SUCCESS); 1159 } 1160 1161 PetscErrorCode MatRetrieveValues_MPISELL(Mat mat) 1162 { 1163 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 1164 1165 PetscFunctionBegin; 1166 PetscCall(MatRetrieveValues(sell->A)); 1167 PetscCall(MatRetrieveValues(sell->B)); 1168 PetscFunctionReturn(PETSC_SUCCESS); 1169 } 1170 1171 PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B, PetscInt d_rlenmax, const PetscInt d_rlen[], PetscInt o_rlenmax, const PetscInt o_rlen[]) 1172 { 1173 Mat_MPISELL *b; 1174 1175 PetscFunctionBegin; 1176 PetscCall(PetscLayoutSetUp(B->rmap)); 1177 PetscCall(PetscLayoutSetUp(B->cmap)); 1178 b = (Mat_MPISELL *)B->data; 1179 1180 if (!B->preallocated) { 1181 /* Explicitly create 2 MATSEQSELL matrices. */ 1182 PetscCall(MatCreate(PETSC_COMM_SELF, &b->A)); 1183 PetscCall(MatSetSizes(b->A, B->rmap->n, B->cmap->n, B->rmap->n, B->cmap->n)); 1184 PetscCall(MatSetBlockSizesFromMats(b->A, B, B)); 1185 PetscCall(MatSetType(b->A, MATSEQSELL)); 1186 PetscCall(MatCreate(PETSC_COMM_SELF, &b->B)); 1187 PetscCall(MatSetSizes(b->B, B->rmap->n, B->cmap->N, B->rmap->n, B->cmap->N)); 1188 PetscCall(MatSetBlockSizesFromMats(b->B, B, B)); 1189 PetscCall(MatSetType(b->B, MATSEQSELL)); 1190 } 1191 1192 PetscCall(MatSeqSELLSetPreallocation(b->A, d_rlenmax, d_rlen)); 1193 PetscCall(MatSeqSELLSetPreallocation(b->B, o_rlenmax, o_rlen)); 1194 B->preallocated = PETSC_TRUE; 1195 B->was_assembled = PETSC_FALSE; 1196 /* 1197 critical for MatAssemblyEnd to work. 1198 MatAssemblyBegin checks it to set up was_assembled 1199 and MatAssemblyEnd checks was_assembled to determine whether to build garray 1200 */ 1201 B->assembled = PETSC_FALSE; 1202 PetscFunctionReturn(PETSC_SUCCESS); 1203 } 1204 1205 PetscErrorCode MatDuplicate_MPISELL(Mat matin, MatDuplicateOption cpvalues, Mat *newmat) 1206 { 1207 Mat mat; 1208 Mat_MPISELL *a, *oldmat = (Mat_MPISELL *)matin->data; 1209 1210 PetscFunctionBegin; 1211 *newmat = NULL; 1212 PetscCall(MatCreate(PetscObjectComm((PetscObject)matin), &mat)); 1213 PetscCall(MatSetSizes(mat, matin->rmap->n, matin->cmap->n, matin->rmap->N, matin->cmap->N)); 1214 PetscCall(MatSetBlockSizesFromMats(mat, matin, matin)); 1215 PetscCall(MatSetType(mat, ((PetscObject)matin)->type_name)); 1216 a = (Mat_MPISELL *)mat->data; 1217 1218 mat->factortype = matin->factortype; 1219 mat->assembled = PETSC_TRUE; 1220 mat->insertmode = NOT_SET_VALUES; 1221 mat->preallocated = PETSC_TRUE; 1222 1223 a->size = oldmat->size; 1224 a->rank = oldmat->rank; 1225 a->donotstash = oldmat->donotstash; 1226 a->roworiented = oldmat->roworiented; 1227 a->rowindices = NULL; 1228 a->rowvalues = NULL; 1229 a->getrowactive = PETSC_FALSE; 1230 1231 PetscCall(PetscLayoutReference(matin->rmap, &mat->rmap)); 1232 PetscCall(PetscLayoutReference(matin->cmap, &mat->cmap)); 1233 1234 if (oldmat->colmap) { 1235 #if defined(PETSC_USE_CTABLE) 1236 PetscCall(PetscHMapIDuplicate(oldmat->colmap, &a->colmap)); 1237 #else 1238 PetscCall(PetscMalloc1(mat->cmap->N, &a->colmap)); 1239 PetscCall(PetscArraycpy(a->colmap, oldmat->colmap, mat->cmap->N)); 1240 #endif 1241 } else a->colmap = NULL; 1242 if (oldmat->garray) { 1243 PetscInt len; 1244 len = oldmat->B->cmap->n; 1245 PetscCall(PetscMalloc1(len + 1, &a->garray)); 1246 if (len) PetscCall(PetscArraycpy(a->garray, oldmat->garray, len)); 1247 } else a->garray = NULL; 1248 1249 PetscCall(VecDuplicate(oldmat->lvec, &a->lvec)); 1250 PetscCall(VecScatterCopy(oldmat->Mvctx, &a->Mvctx)); 1251 PetscCall(MatDuplicate(oldmat->A, cpvalues, &a->A)); 1252 PetscCall(MatDuplicate(oldmat->B, cpvalues, &a->B)); 1253 PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist, &((PetscObject)mat)->qlist)); 1254 *newmat = mat; 1255 PetscFunctionReturn(PETSC_SUCCESS); 1256 } 1257 1258 /*@C 1259 MatMPISELLSetPreallocation - Preallocates memory for a `MATMPISELL` sparse parallel matrix in sell format. 1260 For good matrix assembly performance the user should preallocate the matrix storage by 1261 setting the parameters `d_nz` (or `d_nnz`) and `o_nz` (or `o_nnz`). 1262 1263 Collective 1264 1265 Input Parameters: 1266 + B - the matrix 1267 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1268 (same value is used for all local rows) 1269 . d_nnz - array containing the number of nonzeros in the various rows of the 1270 DIAGONAL portion of the local submatrix (possibly different for each row) 1271 or NULL (`PETSC_NULL_INTEGER` in Fortran), if `d_nz` is used to specify the nonzero structure. 1272 The size of this array is equal to the number of local rows, i.e 'm'. 1273 For matrices that will be factored, you must leave room for (and set) 1274 the diagonal entry even if it is zero. 1275 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1276 submatrix (same value is used for all local rows). 1277 - o_nnz - array containing the number of nonzeros in the various rows of the 1278 OFF-DIAGONAL portion of the local submatrix (possibly different for 1279 each row) or NULL (`PETSC_NULL_INTEGER` in Fortran), if `o_nz` is used to specify the nonzero 1280 structure. The size of this array is equal to the number 1281 of local rows, i.e 'm'. 1282 1283 Example usage: 1284 Consider the following 8x8 matrix with 34 non-zero values, that is 1285 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1286 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1287 as follows 1288 1289 .vb 1290 1 2 0 | 0 3 0 | 0 4 1291 Proc0 0 5 6 | 7 0 0 | 8 0 1292 9 0 10 | 11 0 0 | 12 0 1293 ------------------------------------- 1294 13 0 14 | 15 16 17 | 0 0 1295 Proc1 0 18 0 | 19 20 21 | 0 0 1296 0 0 0 | 22 23 0 | 24 0 1297 ------------------------------------- 1298 Proc2 25 26 27 | 0 0 28 | 29 0 1299 30 0 0 | 31 32 33 | 0 34 1300 .ve 1301 1302 This can be represented as a collection of submatrices as 1303 1304 .vb 1305 A B C 1306 D E F 1307 G H I 1308 .ve 1309 1310 Where the submatrices A,B,C are owned by proc0, D,E,F are 1311 owned by proc1, G,H,I are owned by proc2. 1312 1313 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1314 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1315 The 'M','N' parameters are 8,8, and have the same values on all procs. 1316 1317 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1318 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1319 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1320 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1321 part as `MATSEQSELL` matrices. For example, proc1 will store [E] as a `MATSEQSELL` 1322 matrix, ans [DF] as another SeqSELL matrix. 1323 1324 When `d_nz`, `o_nz` parameters are specified, `d_nz` storage elements are 1325 allocated for every row of the local diagonal submatrix, and o_nz 1326 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1327 One way to choose `d_nz` and `o_nz` is to use the max nonzerors per local 1328 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1329 In this case, the values of d_nz,o_nz are 1330 .vb 1331 proc0 dnz = 2, o_nz = 2 1332 proc1 dnz = 3, o_nz = 2 1333 proc2 dnz = 1, o_nz = 4 1334 .ve 1335 We are allocating m*(d_nz+o_nz) storage locations for every proc. This 1336 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1337 for proc3. i.e we are using 12+15+10=37 storage locations to store 1338 34 values. 1339 1340 When `d_nnz`, `o_nnz` parameters are specified, the storage is specified 1341 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1342 In the above case the values for d_nnz,o_nnz are 1343 .vb 1344 proc0 d_nnz = [2,2,2] and o_nnz = [2,2,2] 1345 proc1 d_nnz = [3,3,2] and o_nnz = [2,1,1] 1346 proc2 d_nnz = [1,1] and o_nnz = [4,4] 1347 .ve 1348 Here the space allocated is according to nz (or maximum values in the nnz 1349 if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37 1350 1351 Level: intermediate 1352 1353 Notes: 1354 If the *_nnz parameter is given then the *_nz parameter is ignored 1355 1356 The stored row and column indices begin with zero. 1357 1358 The parallel matrix is partitioned such that the first m0 rows belong to 1359 process 0, the next m1 rows belong to process 1, the next m2 rows belong 1360 to process 2 etc.. where m0,m1,m2... are the input parameter 'm'. 1361 1362 The DIAGONAL portion of the local submatrix of a processor can be defined 1363 as the submatrix which is obtained by extraction the part corresponding to 1364 the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the 1365 first row that belongs to the processor, r2 is the last row belonging to 1366 the this processor, and c1-c2 is range of indices of the local part of a 1367 vector suitable for applying the matrix to. This is an mxn matrix. In the 1368 common case of a square matrix, the row and column ranges are the same and 1369 the DIAGONAL part is also square. The remaining portion of the local 1370 submatrix (mxN) constitute the OFF-DIAGONAL portion. 1371 1372 If `o_nnz`, `d_nnz` are specified, then `o_nz`, and `d_nz` are ignored. 1373 1374 You can call `MatGetInfo()` to get information on how effective the preallocation was; 1375 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1376 You can also run with the option -info and look for messages with the string 1377 malloc in them to see if additional memory allocation was needed. 1378 1379 .seealso: `Mat`, `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatCreatesell()`, 1380 `MATMPISELL`, `MatGetInfo()`, `PetscSplitOwnership()`, `MATSELL` 1381 @*/ 1382 PetscErrorCode MatMPISELLSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[]) 1383 { 1384 PetscFunctionBegin; 1385 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 1386 PetscValidType(B, 1); 1387 PetscTryMethod(B, "MatMPISELLSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz)); 1388 PetscFunctionReturn(PETSC_SUCCESS); 1389 } 1390 1391 /*MC 1392 MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices, 1393 based on the sliced Ellpack format 1394 1395 Options Database Key: 1396 . -mat_type sell - sets the matrix type to `MATSELL` during a call to `MatSetFromOptions()` 1397 1398 Level: beginner 1399 1400 .seealso: `Mat`, `MatCreateSell()`, `MATSEQSELL`, `MATSELL`, `MATSEQAIJ`, `MATAIJ`, `MATMPIAIJ` 1401 M*/ 1402 1403 /*@C 1404 MatCreateSELL - Creates a sparse parallel matrix in `MATSELL` format. 1405 1406 Collective 1407 1408 Input Parameters: 1409 + comm - MPI communicator 1410 . m - number of local rows (or `PETSC_DECIDE` to have calculated if M is given) 1411 This value should be the same as the local size used in creating the 1412 y vector for the matrix-vector product y = Ax. 1413 . n - This value should be the same as the local size used in creating the 1414 x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have 1415 calculated if `N` is given) For square matrices n is almost always `m`. 1416 . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given) 1417 . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given) 1418 . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix 1419 (same value is used for all local rows) 1420 . d_rlen - array containing the number of nonzeros in the various rows of the 1421 DIAGONAL portion of the local submatrix (possibly different for each row) 1422 or `NULL`, if d_rlenmax is used to specify the nonzero structure. 1423 The size of this array is equal to the number of local rows, i.e `m`. 1424 . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local 1425 submatrix (same value is used for all local rows). 1426 - o_rlen - array containing the number of nonzeros in the various rows of the 1427 OFF-DIAGONAL portion of the local submatrix (possibly different for 1428 each row) or `NULL`, if `o_rlenmax` is used to specify the nonzero 1429 structure. The size of this array is equal to the number 1430 of local rows, i.e `m`. 1431 1432 Output Parameter: 1433 . A - the matrix 1434 1435 Options Database Key: 1436 . -mat_sell_oneindex - Internally use indexing starting at 1 1437 rather than 0. When calling `MatSetValues()`, 1438 the user still MUST index entries starting at 0! 1439 1440 Example: 1441 Consider the following 8x8 matrix with 34 non-zero values, that is 1442 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1443 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1444 as follows 1445 1446 .vb 1447 1 2 0 | 0 3 0 | 0 4 1448 Proc0 0 5 6 | 7 0 0 | 8 0 1449 9 0 10 | 11 0 0 | 12 0 1450 ------------------------------------- 1451 13 0 14 | 15 16 17 | 0 0 1452 Proc1 0 18 0 | 19 20 21 | 0 0 1453 0 0 0 | 22 23 0 | 24 0 1454 ------------------------------------- 1455 Proc2 25 26 27 | 0 0 28 | 29 0 1456 30 0 0 | 31 32 33 | 0 34 1457 .ve 1458 1459 This can be represented as a collection of submatrices as 1460 .vb 1461 A B C 1462 D E F 1463 G H I 1464 .ve 1465 1466 Where the submatrices A,B,C are owned by proc0, D,E,F are 1467 owned by proc1, G,H,I are owned by proc2. 1468 1469 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1470 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1471 The 'M','N' parameters are 8,8, and have the same values on all procs. 1472 1473 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1474 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1475 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1476 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1477 part as `MATSEQSELL` matrices. For example, proc1 will store [E] as a `MATSEQSELL` 1478 matrix, ans [DF] as another `MATSEQSELL` matrix. 1479 1480 When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are 1481 allocated for every row of the local diagonal submatrix, and o_rlenmax 1482 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1483 One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local 1484 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1485 In this case, the values of d_rlenmax,o_rlenmax are 1486 .vb 1487 proc0 - d_rlenmax = 2, o_rlenmax = 2 1488 proc1 - d_rlenmax = 3, o_rlenmax = 2 1489 proc2 - d_rlenmax = 1, o_rlenmax = 4 1490 .ve 1491 We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This 1492 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1493 for proc3. i.e we are using 12+15+10=37 storage locations to store 1494 34 values. 1495 1496 When `d_rlen`, `o_rlen` parameters are specified, the storage is specified 1497 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1498 In the above case the values for `d_nnz`, `o_nnz` are 1499 .vb 1500 proc0 - d_nnz = [2,2,2] and o_nnz = [2,2,2] 1501 proc1 - d_nnz = [3,3,2] and o_nnz = [2,1,1] 1502 proc2 - d_nnz = [1,1] and o_nnz = [4,4] 1503 .ve 1504 Here the space allocated is still 37 though there are 34 nonzeros because 1505 the allocation is always done according to rlenmax. 1506 1507 Level: intermediate 1508 1509 Notes: 1510 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`, 1511 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1512 [MatXXXXSetPreallocation() is, for example, `MatSeqSELLSetPreallocation()`] 1513 1514 If the *_rlen parameter is given then the *_rlenmax parameter is ignored 1515 1516 `m`, `n`, `M`, `N` parameters specify the size of the matrix, and its partitioning across 1517 processors, while `d_rlenmax`, `d_rlen`, `o_rlenmax` , `o_rlen` parameters specify the approximate 1518 storage requirements for this matrix. 1519 1520 If `PETSC_DECIDE` or `PETSC_DETERMINE` is used for a particular argument on one 1521 processor than it must be used on all processors that share the object for 1522 that argument. 1523 1524 The user MUST specify either the local or global matrix dimensions 1525 (possibly both). 1526 1527 The parallel matrix is partitioned across processors such that the 1528 first m0 rows belong to process 0, the next m1 rows belong to 1529 process 1, the next m2 rows belong to process 2 etc.. where 1530 m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores 1531 values corresponding to [`m` x `N`] submatrix. 1532 1533 The columns are logically partitioned with the n0 columns belonging 1534 to 0th partition, the next n1 columns belonging to the next 1535 partition etc.. where n0,n1,n2... are the input parameter `n`. 1536 1537 The DIAGONAL portion of the local submatrix on any given processor 1538 is the submatrix corresponding to the rows and columns `m`, `n` 1539 corresponding to the given processor. i.e diagonal matrix on 1540 process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1] 1541 etc. The remaining portion of the local submatrix [m x (N-n)] 1542 constitute the OFF-DIAGONAL portion. The example below better 1543 illustrates this concept. 1544 1545 For a square global matrix we define each processor's diagonal portion 1546 to be its local rows and the corresponding columns (a square submatrix); 1547 each processor's off-diagonal portion encompasses the remainder of the 1548 local matrix (a rectangular submatrix). 1549 1550 If `o_rlen`, `d_rlen` are specified, then `o_rlenmax`, and `d_rlenmax` are ignored. 1551 1552 When calling this routine with a single process communicator, a matrix of 1553 type `MATSEQSELL` is returned. If a matrix of type `MATMPISELL` is desired for this 1554 type of communicator, use the construction mechanism 1555 .vb 1556 MatCreate(...,&A); 1557 MatSetType(A,MATMPISELL); 1558 MatSetSizes(A, m,n,M,N); 1559 MatMPISELLSetPreallocation(A,...); 1560 .ve 1561 1562 .seealso: `Mat`, `MATSELL`, `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatMPISELLSetPreallocation()`, `MatMPISELLSetPreallocationSELL()`, 1563 `MATMPISELL`, `MatCreateMPISELLWithArrays()` 1564 @*/ 1565 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) 1566 { 1567 PetscMPIInt size; 1568 1569 PetscFunctionBegin; 1570 PetscCall(MatCreate(comm, A)); 1571 PetscCall(MatSetSizes(*A, m, n, M, N)); 1572 PetscCallMPI(MPI_Comm_size(comm, &size)); 1573 if (size > 1) { 1574 PetscCall(MatSetType(*A, MATMPISELL)); 1575 PetscCall(MatMPISELLSetPreallocation(*A, d_rlenmax, d_rlen, o_rlenmax, o_rlen)); 1576 } else { 1577 PetscCall(MatSetType(*A, MATSEQSELL)); 1578 PetscCall(MatSeqSELLSetPreallocation(*A, d_rlenmax, d_rlen)); 1579 } 1580 PetscFunctionReturn(PETSC_SUCCESS); 1581 } 1582 1583 PetscErrorCode MatMPISELLGetSeqSELL(Mat A, Mat *Ad, Mat *Ao, const PetscInt *colmap[]) 1584 { 1585 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 1586 PetscBool flg; 1587 1588 PetscFunctionBegin; 1589 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISELL, &flg)); 1590 PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "This function requires a MATMPISELL matrix as input"); 1591 if (Ad) *Ad = a->A; 1592 if (Ao) *Ao = a->B; 1593 if (colmap) *colmap = a->garray; 1594 PetscFunctionReturn(PETSC_SUCCESS); 1595 } 1596 1597 /*@C 1598 MatMPISELLGetLocalMatCondensed - Creates a `MATSEQSELL` matrix from an `MATMPISELL` matrix by taking all its local rows and NON-ZERO columns 1599 1600 Not Collective 1601 1602 Input Parameters: 1603 + A - the matrix 1604 . scall - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX` 1605 . row - index sets of rows to extract (or `NULL`) 1606 - col - index sets of columns to extract (or `NULL`) 1607 1608 Output Parameter: 1609 . A_loc - the local sequential matrix generated 1610 1611 Level: developer 1612 1613 .seealso: `Mat`, `MATSEQSELL`, `MATMPISELL`, `MatGetOwnershipRange()`, `MatMPISELLGetLocalMat()` 1614 @*/ 1615 PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A, MatReuse scall, IS *row, IS *col, Mat *A_loc) 1616 { 1617 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 1618 PetscInt i, start, end, ncols, nzA, nzB, *cmap, imark, *idx; 1619 IS isrowa, iscola; 1620 Mat *aloc; 1621 PetscBool match; 1622 1623 PetscFunctionBegin; 1624 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISELL, &match)); 1625 PetscCheck(match, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Requires MATMPISELL matrix as input"); 1626 PetscCall(PetscLogEventBegin(MAT_Getlocalmatcondensed, A, 0, 0, 0)); 1627 if (!row) { 1628 start = A->rmap->rstart; 1629 end = A->rmap->rend; 1630 PetscCall(ISCreateStride(PETSC_COMM_SELF, end - start, start, 1, &isrowa)); 1631 } else { 1632 isrowa = *row; 1633 } 1634 if (!col) { 1635 start = A->cmap->rstart; 1636 cmap = a->garray; 1637 nzA = a->A->cmap->n; 1638 nzB = a->B->cmap->n; 1639 PetscCall(PetscMalloc1(nzA + nzB, &idx)); 1640 ncols = 0; 1641 for (i = 0; i < nzB; i++) { 1642 if (cmap[i] < start) idx[ncols++] = cmap[i]; 1643 else break; 1644 } 1645 imark = i; 1646 for (i = 0; i < nzA; i++) idx[ncols++] = start + i; 1647 for (i = imark; i < nzB; i++) idx[ncols++] = cmap[i]; 1648 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncols, idx, PETSC_OWN_POINTER, &iscola)); 1649 } else { 1650 iscola = *col; 1651 } 1652 if (scall != MAT_INITIAL_MATRIX) { 1653 PetscCall(PetscMalloc1(1, &aloc)); 1654 aloc[0] = *A_loc; 1655 } 1656 PetscCall(MatCreateSubMatrices(A, 1, &isrowa, &iscola, scall, &aloc)); 1657 *A_loc = aloc[0]; 1658 PetscCall(PetscFree(aloc)); 1659 if (!row) PetscCall(ISDestroy(&isrowa)); 1660 if (!col) PetscCall(ISDestroy(&iscola)); 1661 PetscCall(PetscLogEventEnd(MAT_Getlocalmatcondensed, A, 0, 0, 0)); 1662 PetscFunctionReturn(PETSC_SUCCESS); 1663 } 1664 1665 #include <../src/mat/impls/aij/mpi/mpiaij.h> 1666 1667 PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 1668 { 1669 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 1670 Mat B; 1671 Mat_MPIAIJ *b; 1672 1673 PetscFunctionBegin; 1674 PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled"); 1675 1676 if (reuse == MAT_REUSE_MATRIX) { 1677 B = *newmat; 1678 } else { 1679 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 1680 PetscCall(MatSetType(B, MATMPIAIJ)); 1681 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 1682 PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs)); 1683 PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL)); 1684 PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL)); 1685 } 1686 b = (Mat_MPIAIJ *)B->data; 1687 1688 if (reuse == MAT_REUSE_MATRIX) { 1689 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A)); 1690 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B)); 1691 } else { 1692 PetscCall(MatDestroy(&b->A)); 1693 PetscCall(MatDestroy(&b->B)); 1694 PetscCall(MatDisAssemble_MPISELL(A)); 1695 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A)); 1696 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B)); 1697 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 1698 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 1699 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 1700 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 1701 } 1702 1703 if (reuse == MAT_INPLACE_MATRIX) { 1704 PetscCall(MatHeaderReplace(A, &B)); 1705 } else { 1706 *newmat = B; 1707 } 1708 PetscFunctionReturn(PETSC_SUCCESS); 1709 } 1710 1711 PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 1712 { 1713 Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data; 1714 Mat B; 1715 Mat_MPISELL *b; 1716 1717 PetscFunctionBegin; 1718 PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled"); 1719 1720 if (reuse == MAT_REUSE_MATRIX) { 1721 B = *newmat; 1722 } else { 1723 Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)a->A->data, *Ba = (Mat_SeqAIJ *)a->B->data; 1724 PetscInt i, d_nz = 0, o_nz = 0, m = A->rmap->N, n = A->cmap->N, lm = A->rmap->n, ln = A->cmap->n; 1725 PetscInt *d_nnz, *o_nnz; 1726 PetscCall(PetscMalloc2(lm, &d_nnz, lm, &o_nnz)); 1727 for (i = 0; i < lm; i++) { 1728 d_nnz[i] = Aa->i[i + 1] - Aa->i[i]; 1729 o_nnz[i] = Ba->i[i + 1] - Ba->i[i]; 1730 if (d_nnz[i] > d_nz) d_nz = d_nnz[i]; 1731 if (o_nnz[i] > o_nz) o_nz = o_nnz[i]; 1732 } 1733 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 1734 PetscCall(MatSetType(B, MATMPISELL)); 1735 PetscCall(MatSetSizes(B, lm, ln, m, n)); 1736 PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs)); 1737 PetscCall(MatSeqSELLSetPreallocation(B, d_nz, d_nnz)); 1738 PetscCall(MatMPISELLSetPreallocation(B, d_nz, d_nnz, o_nz, o_nnz)); 1739 PetscCall(PetscFree2(d_nnz, o_nnz)); 1740 } 1741 b = (Mat_MPISELL *)B->data; 1742 1743 if (reuse == MAT_REUSE_MATRIX) { 1744 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A)); 1745 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B)); 1746 } else { 1747 PetscCall(MatDestroy(&b->A)); 1748 PetscCall(MatDestroy(&b->B)); 1749 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A)); 1750 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B)); 1751 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 1752 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 1753 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 1754 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 1755 } 1756 1757 if (reuse == MAT_INPLACE_MATRIX) { 1758 PetscCall(MatHeaderReplace(A, &B)); 1759 } else { 1760 *newmat = B; 1761 } 1762 PetscFunctionReturn(PETSC_SUCCESS); 1763 } 1764 1765 PetscErrorCode MatSOR_MPISELL(Mat matin, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx) 1766 { 1767 Mat_MPISELL *mat = (Mat_MPISELL *)matin->data; 1768 Vec bb1 = NULL; 1769 1770 PetscFunctionBegin; 1771 if (flag == SOR_APPLY_UPPER) { 1772 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1773 PetscFunctionReturn(PETSC_SUCCESS); 1774 } 1775 1776 if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) PetscCall(VecDuplicate(bb, &bb1)); 1777 1778 if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) { 1779 if (flag & SOR_ZERO_INITIAL_GUESS) { 1780 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1781 its--; 1782 } 1783 1784 while (its--) { 1785 PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1786 PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1787 1788 /* update rhs: bb1 = bb - B*x */ 1789 PetscCall(VecScale(mat->lvec, -1.0)); 1790 PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1)); 1791 1792 /* local sweep */ 1793 PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_SYMMETRIC_SWEEP, fshift, lits, 1, xx)); 1794 } 1795 } else if (flag & SOR_LOCAL_FORWARD_SWEEP) { 1796 if (flag & SOR_ZERO_INITIAL_GUESS) { 1797 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1798 its--; 1799 } 1800 while (its--) { 1801 PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1802 PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1803 1804 /* update rhs: bb1 = bb - B*x */ 1805 PetscCall(VecScale(mat->lvec, -1.0)); 1806 PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1)); 1807 1808 /* local sweep */ 1809 PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_FORWARD_SWEEP, fshift, lits, 1, xx)); 1810 } 1811 } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) { 1812 if (flag & SOR_ZERO_INITIAL_GUESS) { 1813 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1814 its--; 1815 } 1816 while (its--) { 1817 PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1818 PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1819 1820 /* update rhs: bb1 = bb - B*x */ 1821 PetscCall(VecScale(mat->lvec, -1.0)); 1822 PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1)); 1823 1824 /* local sweep */ 1825 PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_BACKWARD_SWEEP, fshift, lits, 1, xx)); 1826 } 1827 } else SETERRQ(PetscObjectComm((PetscObject)matin), PETSC_ERR_SUP, "Parallel SOR not supported"); 1828 1829 PetscCall(VecDestroy(&bb1)); 1830 1831 matin->factorerrortype = mat->A->factorerrortype; 1832 PetscFunctionReturn(PETSC_SUCCESS); 1833 } 1834 1835 #if defined(PETSC_HAVE_CUDA) 1836 PETSC_INTERN PetscErrorCode MatConvert_MPISELL_MPISELLCUDA(Mat, MatType, MatReuse, Mat *); 1837 #endif 1838 1839 /*MC 1840 MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices. 1841 1842 Options Database Keys: 1843 . -mat_type mpisell - sets the matrix type to `MATMPISELL` during a call to `MatSetFromOptions()` 1844 1845 Level: beginner 1846 1847 .seealso: `Mat`, `MATSELL`, `MATSEQSELL` `MatCreateSELL()` 1848 M*/ 1849 PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B) 1850 { 1851 Mat_MPISELL *b; 1852 PetscMPIInt size; 1853 1854 PetscFunctionBegin; 1855 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 1856 PetscCall(PetscNew(&b)); 1857 B->data = (void *)b; 1858 B->ops[0] = MatOps_Values; 1859 B->assembled = PETSC_FALSE; 1860 B->insertmode = NOT_SET_VALUES; 1861 b->size = size; 1862 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B), &b->rank)); 1863 /* build cache for off array entries formed */ 1864 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash)); 1865 1866 b->donotstash = PETSC_FALSE; 1867 b->colmap = NULL; 1868 b->garray = NULL; 1869 b->roworiented = PETSC_TRUE; 1870 1871 /* stuff used for matrix vector multiply */ 1872 b->lvec = NULL; 1873 b->Mvctx = NULL; 1874 1875 /* stuff for MatGetRow() */ 1876 b->rowindices = NULL; 1877 b->rowvalues = NULL; 1878 b->getrowactive = PETSC_FALSE; 1879 1880 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_MPISELL)); 1881 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_MPISELL)); 1882 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsTranspose_C", MatIsTranspose_MPISELL)); 1883 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPISELLSetPreallocation_C", MatMPISELLSetPreallocation_MPISELL)); 1884 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisell_mpiaij_C", MatConvert_MPISELL_MPIAIJ)); 1885 #if defined(PETSC_HAVE_CUDA) 1886 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisell_mpisellcuda_C", MatConvert_MPISELL_MPISELLCUDA)); 1887 #endif 1888 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDiagonalScaleLocal_C", MatDiagonalScaleLocal_MPISELL)); 1889 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATMPISELL)); 1890 PetscFunctionReturn(PETSC_SUCCESS); 1891 } 1892