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 acess. 44 */ 45 PetscErrorCode MatCreateColmap_MPISELL_Private(Mat mat) 46 { 47 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 48 PetscInt n = sell->B->cmap->n, i; 49 50 PetscFunctionBegin; 51 PetscCheck(sell->garray, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MPISELL Matrix was assembled but is missing garray"); 52 #if defined(PETSC_USE_CTABLE) 53 PetscCall(PetscHMapICreateWithSize(n, &sell->colmap)); 54 for (i = 0; i < n; i++) PetscCall(PetscHMapISet(sell->colmap, sell->garray[i] + 1, i + 1)); 55 #else 56 PetscCall(PetscCalloc1(mat->cmap->N + 1, &sell->colmap)); 57 for (i = 0; i < n; i++) sell->colmap[sell->garray[i]] = i + 1; 58 #endif 59 PetscFunctionReturn(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(PetscHMapIGetWithDefault(sell->colmap, in[j] + 1, 0, &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(PetscHMapIGetWithDefault(sell->colmap, idxn[j] + 1, 0, &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(PetscHMapIDestroy(&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 NULL}; 1135 1136 /* ----------------------------------------------------------------------------------------*/ 1137 1138 PetscErrorCode MatStoreValues_MPISELL(Mat mat) 1139 { 1140 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 1141 1142 PetscFunctionBegin; 1143 PetscCall(MatStoreValues(sell->A)); 1144 PetscCall(MatStoreValues(sell->B)); 1145 PetscFunctionReturn(0); 1146 } 1147 1148 PetscErrorCode MatRetrieveValues_MPISELL(Mat mat) 1149 { 1150 Mat_MPISELL *sell = (Mat_MPISELL *)mat->data; 1151 1152 PetscFunctionBegin; 1153 PetscCall(MatRetrieveValues(sell->A)); 1154 PetscCall(MatRetrieveValues(sell->B)); 1155 PetscFunctionReturn(0); 1156 } 1157 1158 PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B, PetscInt d_rlenmax, const PetscInt d_rlen[], PetscInt o_rlenmax, const PetscInt o_rlen[]) 1159 { 1160 Mat_MPISELL *b; 1161 1162 PetscFunctionBegin; 1163 PetscCall(PetscLayoutSetUp(B->rmap)); 1164 PetscCall(PetscLayoutSetUp(B->cmap)); 1165 b = (Mat_MPISELL *)B->data; 1166 1167 if (!B->preallocated) { 1168 /* Explicitly create 2 MATSEQSELL matrices. */ 1169 PetscCall(MatCreate(PETSC_COMM_SELF, &b->A)); 1170 PetscCall(MatSetSizes(b->A, B->rmap->n, B->cmap->n, B->rmap->n, B->cmap->n)); 1171 PetscCall(MatSetBlockSizesFromMats(b->A, B, B)); 1172 PetscCall(MatSetType(b->A, MATSEQSELL)); 1173 PetscCall(MatCreate(PETSC_COMM_SELF, &b->B)); 1174 PetscCall(MatSetSizes(b->B, B->rmap->n, B->cmap->N, B->rmap->n, B->cmap->N)); 1175 PetscCall(MatSetBlockSizesFromMats(b->B, B, B)); 1176 PetscCall(MatSetType(b->B, MATSEQSELL)); 1177 } 1178 1179 PetscCall(MatSeqSELLSetPreallocation(b->A, d_rlenmax, d_rlen)); 1180 PetscCall(MatSeqSELLSetPreallocation(b->B, o_rlenmax, o_rlen)); 1181 B->preallocated = PETSC_TRUE; 1182 B->was_assembled = PETSC_FALSE; 1183 /* 1184 critical for MatAssemblyEnd to work. 1185 MatAssemblyBegin checks it to set up was_assembled 1186 and MatAssemblyEnd checks was_assembled to determine whether to build garray 1187 */ 1188 B->assembled = PETSC_FALSE; 1189 PetscFunctionReturn(0); 1190 } 1191 1192 PetscErrorCode MatDuplicate_MPISELL(Mat matin, MatDuplicateOption cpvalues, Mat *newmat) 1193 { 1194 Mat mat; 1195 Mat_MPISELL *a, *oldmat = (Mat_MPISELL *)matin->data; 1196 1197 PetscFunctionBegin; 1198 *newmat = NULL; 1199 PetscCall(MatCreate(PetscObjectComm((PetscObject)matin), &mat)); 1200 PetscCall(MatSetSizes(mat, matin->rmap->n, matin->cmap->n, matin->rmap->N, matin->cmap->N)); 1201 PetscCall(MatSetBlockSizesFromMats(mat, matin, matin)); 1202 PetscCall(MatSetType(mat, ((PetscObject)matin)->type_name)); 1203 a = (Mat_MPISELL *)mat->data; 1204 1205 mat->factortype = matin->factortype; 1206 mat->assembled = PETSC_TRUE; 1207 mat->insertmode = NOT_SET_VALUES; 1208 mat->preallocated = PETSC_TRUE; 1209 1210 a->size = oldmat->size; 1211 a->rank = oldmat->rank; 1212 a->donotstash = oldmat->donotstash; 1213 a->roworiented = oldmat->roworiented; 1214 a->rowindices = NULL; 1215 a->rowvalues = NULL; 1216 a->getrowactive = PETSC_FALSE; 1217 1218 PetscCall(PetscLayoutReference(matin->rmap, &mat->rmap)); 1219 PetscCall(PetscLayoutReference(matin->cmap, &mat->cmap)); 1220 1221 if (oldmat->colmap) { 1222 #if defined(PETSC_USE_CTABLE) 1223 PetscCall(PetscHMapIDuplicate(oldmat->colmap, &a->colmap)); 1224 #else 1225 PetscCall(PetscMalloc1(mat->cmap->N, &a->colmap)); 1226 PetscCall(PetscArraycpy(a->colmap, oldmat->colmap, mat->cmap->N)); 1227 #endif 1228 } else a->colmap = NULL; 1229 if (oldmat->garray) { 1230 PetscInt len; 1231 len = oldmat->B->cmap->n; 1232 PetscCall(PetscMalloc1(len + 1, &a->garray)); 1233 if (len) PetscCall(PetscArraycpy(a->garray, oldmat->garray, len)); 1234 } else a->garray = NULL; 1235 1236 PetscCall(VecDuplicate(oldmat->lvec, &a->lvec)); 1237 PetscCall(VecScatterCopy(oldmat->Mvctx, &a->Mvctx)); 1238 PetscCall(MatDuplicate(oldmat->A, cpvalues, &a->A)); 1239 PetscCall(MatDuplicate(oldmat->B, cpvalues, &a->B)); 1240 PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist, &((PetscObject)mat)->qlist)); 1241 *newmat = mat; 1242 PetscFunctionReturn(0); 1243 } 1244 1245 /*@C 1246 MatMPISELLSetPreallocation - Preallocates memory for a `MATMPISELL` sparse parallel matrix in sell format. 1247 For good matrix assembly performance the user should preallocate the matrix storage by 1248 setting the parameters d_nz (or d_nnz) and o_nz (or o_nnz). 1249 1250 Collective 1251 1252 Input Parameters: 1253 + B - the matrix 1254 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1255 (same value is used for all local rows) 1256 . d_nnz - array containing the number of nonzeros in the various rows of the 1257 DIAGONAL portion of the local submatrix (possibly different for each row) 1258 or NULL (`PETSC_NULL_INTEGER` in Fortran), if d_nz is used to specify the nonzero structure. 1259 The size of this array is equal to the number of local rows, i.e 'm'. 1260 For matrices that will be factored, you must leave room for (and set) 1261 the diagonal entry even if it is zero. 1262 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1263 submatrix (same value is used for all local rows). 1264 - o_nnz - array containing the number of nonzeros in the various rows of the 1265 OFF-DIAGONAL portion of the local submatrix (possibly different for 1266 each row) or NULL (`PETSC_NULL_INTEGER` in Fortran), if o_nz is used to specify the nonzero 1267 structure. The size of this array is equal to the number 1268 of local rows, i.e 'm'. 1269 1270 If the *_nnz parameter is given then the *_nz parameter is ignored 1271 1272 The stored row and column indices begin with zero. 1273 1274 The parallel matrix is partitioned such that the first m0 rows belong to 1275 process 0, the next m1 rows belong to process 1, the next m2 rows belong 1276 to process 2 etc.. where m0,m1,m2... are the input parameter 'm'. 1277 1278 The DIAGONAL portion of the local submatrix of a processor can be defined 1279 as the submatrix which is obtained by extraction the part corresponding to 1280 the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the 1281 first row that belongs to the processor, r2 is the last row belonging to 1282 the this processor, and c1-c2 is range of indices of the local part of a 1283 vector suitable for applying the matrix to. This is an mxn matrix. In the 1284 common case of a square matrix, the row and column ranges are the same and 1285 the DIAGONAL part is also square. The remaining portion of the local 1286 submatrix (mxN) constitute the OFF-DIAGONAL portion. 1287 1288 If o_nnz, d_nnz are specified, then o_nz, and d_nz are ignored. 1289 1290 You can call `MatGetInfo()` to get information on how effective the preallocation was; 1291 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1292 You can also run with the option -info and look for messages with the string 1293 malloc in them to see if additional memory allocation was needed. 1294 1295 Example usage: 1296 1297 Consider the following 8x8 matrix with 34 non-zero values, that is 1298 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1299 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1300 as follows: 1301 1302 .vb 1303 1 2 0 | 0 3 0 | 0 4 1304 Proc0 0 5 6 | 7 0 0 | 8 0 1305 9 0 10 | 11 0 0 | 12 0 1306 ------------------------------------- 1307 13 0 14 | 15 16 17 | 0 0 1308 Proc1 0 18 0 | 19 20 21 | 0 0 1309 0 0 0 | 22 23 0 | 24 0 1310 ------------------------------------- 1311 Proc2 25 26 27 | 0 0 28 | 29 0 1312 30 0 0 | 31 32 33 | 0 34 1313 .ve 1314 1315 This can be represented as a collection of submatrices as: 1316 1317 .vb 1318 A B C 1319 D E F 1320 G H I 1321 .ve 1322 1323 Where the submatrices A,B,C are owned by proc0, D,E,F are 1324 owned by proc1, G,H,I are owned by proc2. 1325 1326 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1327 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1328 The 'M','N' parameters are 8,8, and have the same values on all procs. 1329 1330 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1331 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1332 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1333 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1334 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1335 matrix, ans [DF] as another SeqSELL matrix. 1336 1337 When d_nz, o_nz parameters are specified, d_nz storage elements are 1338 allocated for every row of the local diagonal submatrix, and o_nz 1339 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1340 One way to choose d_nz and o_nz is to use the max nonzerors per local 1341 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1342 In this case, the values of d_nz,o_nz are: 1343 .vb 1344 proc0 : dnz = 2, o_nz = 2 1345 proc1 : dnz = 3, o_nz = 2 1346 proc2 : dnz = 1, o_nz = 4 1347 .ve 1348 We are allocating m*(d_nz+o_nz) storage locations for every proc. This 1349 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1350 for proc3. i.e we are using 12+15+10=37 storage locations to store 1351 34 values. 1352 1353 When d_nnz, o_nnz parameters are specified, the storage is specified 1354 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1355 In the above case the values for d_nnz,o_nnz are: 1356 .vb 1357 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1358 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1359 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1360 .ve 1361 Here the space allocated is according to nz (or maximum values in the nnz 1362 if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37 1363 1364 Level: intermediate 1365 1366 .seealso: `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatCreatesell()`, 1367 `MATMPISELL`, `MatGetInfo()`, `PetscSplitOwnership()`, `MATSELL` 1368 @*/ 1369 PetscErrorCode MatMPISELLSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[]) 1370 { 1371 PetscFunctionBegin; 1372 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 1373 PetscValidType(B, 1); 1374 PetscTryMethod(B, "MatMPISELLSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz)); 1375 PetscFunctionReturn(0); 1376 } 1377 1378 /*MC 1379 MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices, 1380 based on the sliced Ellpack format 1381 1382 Options Database Keys: 1383 . -mat_type sell - sets the matrix type to "seqsell" during a call to MatSetFromOptions() 1384 1385 Level: beginner 1386 1387 .seealso: `MatCreateSell()`, `MATSEQSELL`, `MATSELL`, `MATSEQAIJ`, `MATAIJ`, `MATMPIAIJ` 1388 M*/ 1389 1390 /*@C 1391 MatCreateSELL - Creates a sparse parallel matrix in `MATSELL` format. 1392 1393 Collective 1394 1395 Input Parameters: 1396 + comm - MPI communicator 1397 . m - number of local rows (or `PETSC_DECIDE` to have calculated if M is given) 1398 This value should be the same as the local size used in creating the 1399 y vector for the matrix-vector product y = Ax. 1400 . n - This value should be the same as the local size used in creating the 1401 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 1402 calculated if N is given) For square matrices n is almost always m. 1403 . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given) 1404 . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given) 1405 . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix 1406 (same value is used for all local rows) 1407 . d_rlen - array containing the number of nonzeros in the various rows of the 1408 DIAGONAL portion of the local submatrix (possibly different for each row) 1409 or NULL, if d_rlenmax is used to specify the nonzero structure. 1410 The size of this array is equal to the number of local rows, i.e 'm'. 1411 . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local 1412 submatrix (same value is used for all local rows). 1413 - o_rlen - array containing the number of nonzeros in the various rows of the 1414 OFF-DIAGONAL portion of the local submatrix (possibly different for 1415 each row) or NULL, if o_rlenmax is used to specify the nonzero 1416 structure. The size of this array is equal to the number 1417 of local rows, i.e 'm'. 1418 1419 Output Parameter: 1420 . A - the matrix 1421 1422 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`, 1423 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1424 [MatXXXXSetPreallocation() is, for example, `MatSeqSELLSetPreallocation()`] 1425 1426 Notes: 1427 If the *_rlen parameter is given then the *_rlenmax parameter is ignored 1428 1429 m,n,M,N parameters specify the size of the matrix, and its partitioning across 1430 processors, while d_rlenmax,d_rlen,o_rlenmax,o_rlen parameters specify the approximate 1431 storage requirements for this matrix. 1432 1433 If `PETSC_DECIDE` or `PETSC_DETERMINE` is used for a particular argument on one 1434 processor than it must be used on all processors that share the object for 1435 that argument. 1436 1437 The user MUST specify either the local or global matrix dimensions 1438 (possibly both). 1439 1440 The parallel matrix is partitioned across processors such that the 1441 first m0 rows belong to process 0, the next m1 rows belong to 1442 process 1, the next m2 rows belong to process 2 etc.. where 1443 m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores 1444 values corresponding to [m x N] submatrix. 1445 1446 The columns are logically partitioned with the n0 columns belonging 1447 to 0th partition, the next n1 columns belonging to the next 1448 partition etc.. where n0,n1,n2... are the input parameter 'n'. 1449 1450 The DIAGONAL portion of the local submatrix on any given processor 1451 is the submatrix corresponding to the rows and columns m,n 1452 corresponding to the given processor. i.e diagonal matrix on 1453 process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1] 1454 etc. The remaining portion of the local submatrix [m x (N-n)] 1455 constitute the OFF-DIAGONAL portion. The example below better 1456 illustrates this concept. 1457 1458 For a square global matrix we define each processor's diagonal portion 1459 to be its local rows and the corresponding columns (a square submatrix); 1460 each processor's off-diagonal portion encompasses the remainder of the 1461 local matrix (a rectangular submatrix). 1462 1463 If o_rlen, d_rlen are specified, then o_rlenmax, and d_rlenmax are ignored. 1464 1465 When calling this routine with a single process communicator, a matrix of 1466 type `MATSEQSELL` is returned. If a matrix of type `MATMPISELL` is desired for this 1467 type of communicator, use the construction mechanism: 1468 `MatCreate`(...,&A); `MatSetType`(A,`MATMPISELL`); `MatSetSizes`(A, m,n,M,N); `MatMPISELLSetPreallocation`(A,...); 1469 1470 Options Database Keys: 1471 - -mat_sell_oneindex - Internally use indexing starting at 1 1472 rather than 0. Note that when calling MatSetValues(), 1473 the user still MUST index entries starting at 0! 1474 1475 Example usage: 1476 1477 Consider the following 8x8 matrix with 34 non-zero values, that is 1478 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1479 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1480 as follows: 1481 1482 .vb 1483 1 2 0 | 0 3 0 | 0 4 1484 Proc0 0 5 6 | 7 0 0 | 8 0 1485 9 0 10 | 11 0 0 | 12 0 1486 ------------------------------------- 1487 13 0 14 | 15 16 17 | 0 0 1488 Proc1 0 18 0 | 19 20 21 | 0 0 1489 0 0 0 | 22 23 0 | 24 0 1490 ------------------------------------- 1491 Proc2 25 26 27 | 0 0 28 | 29 0 1492 30 0 0 | 31 32 33 | 0 34 1493 .ve 1494 1495 This can be represented as a collection of submatrices as: 1496 1497 .vb 1498 A B C 1499 D E F 1500 G H I 1501 .ve 1502 1503 Where the submatrices A,B,C are owned by proc0, D,E,F are 1504 owned by proc1, G,H,I are owned by proc2. 1505 1506 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1507 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1508 The 'M','N' parameters are 8,8, and have the same values on all procs. 1509 1510 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1511 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1512 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1513 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1514 part as `MATSSESELL` matrices. for eg: proc1 will store [E] as a `MATSEQSELL` 1515 matrix, ans [DF] as another `MATSEQSELL` matrix. 1516 1517 When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are 1518 allocated for every row of the local diagonal submatrix, and o_rlenmax 1519 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1520 One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local 1521 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1522 In this case, the values of d_rlenmax,o_rlenmax are: 1523 .vb 1524 proc0 : d_rlenmax = 2, o_rlenmax = 2 1525 proc1 : d_rlenmax = 3, o_rlenmax = 2 1526 proc2 : d_rlenmax = 1, o_rlenmax = 4 1527 .ve 1528 We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This 1529 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1530 for proc3. i.e we are using 12+15+10=37 storage locations to store 1531 34 values. 1532 1533 When d_rlen, o_rlen parameters are specified, the storage is specified 1534 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1535 In the above case the values for d_nnz,o_nnz are: 1536 .vb 1537 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1538 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1539 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1540 .ve 1541 Here the space allocated is still 37 though there are 34 nonzeros because 1542 the allocation is always done according to rlenmax. 1543 1544 Level: intermediate 1545 1546 .seealso: `MATSELL`, `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatMPISELLSetPreallocation()`, `MatMPISELLSetPreallocationSELL()`, 1547 `MATMPISELL`, `MatCreateMPISELLWithArrays()` 1548 @*/ 1549 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) 1550 { 1551 PetscMPIInt size; 1552 1553 PetscFunctionBegin; 1554 PetscCall(MatCreate(comm, A)); 1555 PetscCall(MatSetSizes(*A, m, n, M, N)); 1556 PetscCallMPI(MPI_Comm_size(comm, &size)); 1557 if (size > 1) { 1558 PetscCall(MatSetType(*A, MATMPISELL)); 1559 PetscCall(MatMPISELLSetPreallocation(*A, d_rlenmax, d_rlen, o_rlenmax, o_rlen)); 1560 } else { 1561 PetscCall(MatSetType(*A, MATSEQSELL)); 1562 PetscCall(MatSeqSELLSetPreallocation(*A, d_rlenmax, d_rlen)); 1563 } 1564 PetscFunctionReturn(0); 1565 } 1566 1567 PetscErrorCode MatMPISELLGetSeqSELL(Mat A, Mat *Ad, Mat *Ao, const PetscInt *colmap[]) 1568 { 1569 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 1570 PetscBool flg; 1571 1572 PetscFunctionBegin; 1573 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISELL, &flg)); 1574 PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "This function requires a MATMPISELL matrix as input"); 1575 if (Ad) *Ad = a->A; 1576 if (Ao) *Ao = a->B; 1577 if (colmap) *colmap = a->garray; 1578 PetscFunctionReturn(0); 1579 } 1580 1581 /*@C 1582 MatMPISELLGetLocalMatCondensed - Creates a `MATSEQSELL` matrix from an `MATMPISELL` matrix by taking all its local rows and NON-ZERO columns 1583 1584 Not Collective 1585 1586 Input Parameters: 1587 + A - the matrix 1588 . scall - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX` 1589 - row, col - index sets of rows and columns to extract (or NULL) 1590 1591 Output Parameter: 1592 . A_loc - the local sequential matrix generated 1593 1594 Level: developer 1595 1596 .seealso: `MATSEQSELL`, `MATMPISELL`, `MatGetOwnershipRange()`, `MatMPISELLGetLocalMat()` 1597 @*/ 1598 PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A, MatReuse scall, IS *row, IS *col, Mat *A_loc) 1599 { 1600 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 1601 PetscInt i, start, end, ncols, nzA, nzB, *cmap, imark, *idx; 1602 IS isrowa, iscola; 1603 Mat *aloc; 1604 PetscBool match; 1605 1606 PetscFunctionBegin; 1607 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISELL, &match)); 1608 PetscCheck(match, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Requires MATMPISELL matrix as input"); 1609 PetscCall(PetscLogEventBegin(MAT_Getlocalmatcondensed, A, 0, 0, 0)); 1610 if (!row) { 1611 start = A->rmap->rstart; 1612 end = A->rmap->rend; 1613 PetscCall(ISCreateStride(PETSC_COMM_SELF, end - start, start, 1, &isrowa)); 1614 } else { 1615 isrowa = *row; 1616 } 1617 if (!col) { 1618 start = A->cmap->rstart; 1619 cmap = a->garray; 1620 nzA = a->A->cmap->n; 1621 nzB = a->B->cmap->n; 1622 PetscCall(PetscMalloc1(nzA + nzB, &idx)); 1623 ncols = 0; 1624 for (i = 0; i < nzB; i++) { 1625 if (cmap[i] < start) idx[ncols++] = cmap[i]; 1626 else break; 1627 } 1628 imark = i; 1629 for (i = 0; i < nzA; i++) idx[ncols++] = start + i; 1630 for (i = imark; i < nzB; i++) idx[ncols++] = cmap[i]; 1631 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncols, idx, PETSC_OWN_POINTER, &iscola)); 1632 } else { 1633 iscola = *col; 1634 } 1635 if (scall != MAT_INITIAL_MATRIX) { 1636 PetscCall(PetscMalloc1(1, &aloc)); 1637 aloc[0] = *A_loc; 1638 } 1639 PetscCall(MatCreateSubMatrices(A, 1, &isrowa, &iscola, scall, &aloc)); 1640 *A_loc = aloc[0]; 1641 PetscCall(PetscFree(aloc)); 1642 if (!row) PetscCall(ISDestroy(&isrowa)); 1643 if (!col) PetscCall(ISDestroy(&iscola)); 1644 PetscCall(PetscLogEventEnd(MAT_Getlocalmatcondensed, A, 0, 0, 0)); 1645 PetscFunctionReturn(0); 1646 } 1647 1648 #include <../src/mat/impls/aij/mpi/mpiaij.h> 1649 1650 PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 1651 { 1652 Mat_MPISELL *a = (Mat_MPISELL *)A->data; 1653 Mat B; 1654 Mat_MPIAIJ *b; 1655 1656 PetscFunctionBegin; 1657 PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled"); 1658 1659 if (reuse == MAT_REUSE_MATRIX) { 1660 B = *newmat; 1661 } else { 1662 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 1663 PetscCall(MatSetType(B, MATMPIAIJ)); 1664 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 1665 PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs)); 1666 PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL)); 1667 PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL)); 1668 } 1669 b = (Mat_MPIAIJ *)B->data; 1670 1671 if (reuse == MAT_REUSE_MATRIX) { 1672 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A)); 1673 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B)); 1674 } else { 1675 PetscCall(MatDestroy(&b->A)); 1676 PetscCall(MatDestroy(&b->B)); 1677 PetscCall(MatDisAssemble_MPISELL(A)); 1678 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A)); 1679 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B)); 1680 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 1681 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 1682 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 1683 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 1684 } 1685 1686 if (reuse == MAT_INPLACE_MATRIX) { 1687 PetscCall(MatHeaderReplace(A, &B)); 1688 } else { 1689 *newmat = B; 1690 } 1691 PetscFunctionReturn(0); 1692 } 1693 1694 PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 1695 { 1696 Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data; 1697 Mat B; 1698 Mat_MPISELL *b; 1699 1700 PetscFunctionBegin; 1701 PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled"); 1702 1703 if (reuse == MAT_REUSE_MATRIX) { 1704 B = *newmat; 1705 } else { 1706 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 1707 PetscCall(MatSetType(B, MATMPISELL)); 1708 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 1709 PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs)); 1710 PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL)); 1711 PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL)); 1712 } 1713 b = (Mat_MPISELL *)B->data; 1714 1715 if (reuse == MAT_REUSE_MATRIX) { 1716 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A)); 1717 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B)); 1718 } else { 1719 PetscCall(MatDestroy(&b->A)); 1720 PetscCall(MatDestroy(&b->B)); 1721 PetscCall(MatDisAssemble_MPIAIJ(A)); 1722 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A)); 1723 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B)); 1724 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 1725 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 1726 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 1727 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 1728 } 1729 1730 if (reuse == MAT_INPLACE_MATRIX) { 1731 PetscCall(MatHeaderReplace(A, &B)); 1732 } else { 1733 *newmat = B; 1734 } 1735 PetscFunctionReturn(0); 1736 } 1737 1738 PetscErrorCode MatSOR_MPISELL(Mat matin, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx) 1739 { 1740 Mat_MPISELL *mat = (Mat_MPISELL *)matin->data; 1741 Vec bb1 = NULL; 1742 1743 PetscFunctionBegin; 1744 if (flag == SOR_APPLY_UPPER) { 1745 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1746 PetscFunctionReturn(0); 1747 } 1748 1749 if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) PetscCall(VecDuplicate(bb, &bb1)); 1750 1751 if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) { 1752 if (flag & SOR_ZERO_INITIAL_GUESS) { 1753 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1754 its--; 1755 } 1756 1757 while (its--) { 1758 PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1759 PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1760 1761 /* update rhs: bb1 = bb - B*x */ 1762 PetscCall(VecScale(mat->lvec, -1.0)); 1763 PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1)); 1764 1765 /* local sweep */ 1766 PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_SYMMETRIC_SWEEP, fshift, lits, 1, xx)); 1767 } 1768 } else if (flag & SOR_LOCAL_FORWARD_SWEEP) { 1769 if (flag & SOR_ZERO_INITIAL_GUESS) { 1770 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1771 its--; 1772 } 1773 while (its--) { 1774 PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1775 PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1776 1777 /* update rhs: bb1 = bb - B*x */ 1778 PetscCall(VecScale(mat->lvec, -1.0)); 1779 PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1)); 1780 1781 /* local sweep */ 1782 PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_FORWARD_SWEEP, fshift, lits, 1, xx)); 1783 } 1784 } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) { 1785 if (flag & SOR_ZERO_INITIAL_GUESS) { 1786 PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx)); 1787 its--; 1788 } 1789 while (its--) { 1790 PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1791 PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD)); 1792 1793 /* update rhs: bb1 = bb - B*x */ 1794 PetscCall(VecScale(mat->lvec, -1.0)); 1795 PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1)); 1796 1797 /* local sweep */ 1798 PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_BACKWARD_SWEEP, fshift, lits, 1, xx)); 1799 } 1800 } else SETERRQ(PetscObjectComm((PetscObject)matin), PETSC_ERR_SUP, "Parallel SOR not supported"); 1801 1802 PetscCall(VecDestroy(&bb1)); 1803 1804 matin->factorerrortype = mat->A->factorerrortype; 1805 PetscFunctionReturn(0); 1806 } 1807 1808 /*MC 1809 MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices. 1810 1811 Options Database Keys: 1812 . -mat_type mpisell - sets the matrix type to `MATMPISELL` during a call to `MatSetFromOptions()` 1813 1814 Level: beginner 1815 1816 .seealso: `MATSELL`, `MATSEQSELL` `MatCreateSELL()` 1817 M*/ 1818 PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B) 1819 { 1820 Mat_MPISELL *b; 1821 PetscMPIInt size; 1822 1823 PetscFunctionBegin; 1824 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 1825 PetscCall(PetscNew(&b)); 1826 B->data = (void *)b; 1827 PetscCall(PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps))); 1828 B->assembled = PETSC_FALSE; 1829 B->insertmode = NOT_SET_VALUES; 1830 b->size = size; 1831 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B), &b->rank)); 1832 /* build cache for off array entries formed */ 1833 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash)); 1834 1835 b->donotstash = PETSC_FALSE; 1836 b->colmap = NULL; 1837 b->garray = NULL; 1838 b->roworiented = PETSC_TRUE; 1839 1840 /* stuff used for matrix vector multiply */ 1841 b->lvec = NULL; 1842 b->Mvctx = NULL; 1843 1844 /* stuff for MatGetRow() */ 1845 b->rowindices = NULL; 1846 b->rowvalues = NULL; 1847 b->getrowactive = PETSC_FALSE; 1848 1849 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_MPISELL)); 1850 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_MPISELL)); 1851 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsTranspose_C", MatIsTranspose_MPISELL)); 1852 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPISELLSetPreallocation_C", MatMPISELLSetPreallocation_MPISELL)); 1853 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisell_mpiaij_C", MatConvert_MPISELL_MPIAIJ)); 1854 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDiagonalScaleLocal_C", MatDiagonalScaleLocal_MPISELL)); 1855 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATMPISELL)); 1856 PetscFunctionReturn(0); 1857 } 1858