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