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