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