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