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