1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 2 3 typedef struct { 4 PetscScalar diag; 5 } Mat_ConstantDiagonal; 6 7 static PetscErrorCode MatAXPY_ConstantDiagonal(Mat Y, PetscScalar a, Mat X, MatStructure str) 8 { 9 Mat_ConstantDiagonal *yctx = (Mat_ConstantDiagonal *)Y->data; 10 Mat_ConstantDiagonal *xctx = (Mat_ConstantDiagonal *)X->data; 11 12 PetscFunctionBegin; 13 yctx->diag += a * xctx->diag; 14 PetscFunctionReturn(PETSC_SUCCESS); 15 } 16 17 static PetscErrorCode MatEqual_ConstantDiagonal(Mat Y, Mat X, PetscBool *equal) 18 { 19 Mat_ConstantDiagonal *yctx = (Mat_ConstantDiagonal *)Y->data; 20 Mat_ConstantDiagonal *xctx = (Mat_ConstantDiagonal *)X->data; 21 22 PetscFunctionBegin; 23 *equal = (yctx->diag == xctx->diag) ? PETSC_TRUE : PETSC_FALSE; 24 PetscFunctionReturn(PETSC_SUCCESS); 25 } 26 27 static PetscErrorCode MatGetRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[]) 28 { 29 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data; 30 31 PetscFunctionBegin; 32 if (ncols) *ncols = 1; 33 if (cols) { 34 PetscCall(PetscMalloc1(1, cols)); 35 (*cols)[0] = row; 36 } 37 if (vals) { 38 PetscCall(PetscMalloc1(1, vals)); 39 (*vals)[0] = ctx->diag; 40 } 41 PetscFunctionReturn(PETSC_SUCCESS); 42 } 43 44 static PetscErrorCode MatRestoreRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[]) 45 { 46 PetscFunctionBegin; 47 if (cols) PetscCall(PetscFree(*cols)); 48 if (vals) PetscCall(PetscFree(*vals)); 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 static PetscErrorCode MatMultAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3) 53 { 54 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data; 55 56 PetscFunctionBegin; 57 if (v2 == v3) { 58 PetscCall(VecAXPBY(v3, ctx->diag, 1.0, v1)); 59 } else { 60 PetscCall(VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2)); 61 } 62 PetscFunctionReturn(PETSC_SUCCESS); 63 } 64 65 static PetscErrorCode MatMultHermitianTransposeAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3) 66 { 67 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data; 68 69 PetscFunctionBegin; 70 if (v2 == v3) { 71 PetscCall(VecAXPBY(v3, PetscConj(ctx->diag), 1.0, v1)); 72 } else { 73 PetscCall(VecAXPBYPCZ(v3, PetscConj(ctx->diag), 1.0, 0.0, v1, v2)); 74 } 75 PetscFunctionReturn(PETSC_SUCCESS); 76 } 77 78 static PetscErrorCode MatNorm_ConstantDiagonal(Mat A, NormType type, PetscReal *nrm) 79 { 80 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data; 81 82 PetscFunctionBegin; 83 PetscCheck(type == NORM_FROBENIUS || type == NORM_2 || type == NORM_1 || type == NORM_INFINITY, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm"); 84 *nrm = PetscAbsScalar(ctx->diag); 85 PetscFunctionReturn(PETSC_SUCCESS); 86 } 87 88 static PetscErrorCode MatCreateSubMatrices_ConstantDiagonal(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[]) 89 90 { 91 Mat B; 92 93 PetscFunctionBegin; 94 PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B)); 95 PetscCall(MatCreateSubMatrices(B, n, irow, icol, scall, submat)); 96 PetscCall(MatDestroy(&B)); 97 PetscFunctionReturn(PETSC_SUCCESS); 98 } 99 100 static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B) 101 { 102 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data; 103 104 PetscFunctionBegin; 105 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 106 PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 107 PetscCall(MatSetBlockSizesFromMats(*B, A, A)); 108 PetscCall(MatSetType(*B, MATCONSTANTDIAGONAL)); 109 PetscCall(PetscLayoutReference(A->rmap, &(*B)->rmap)); 110 PetscCall(PetscLayoutReference(A->cmap, &(*B)->cmap)); 111 if (op == MAT_COPY_VALUES) { 112 Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal *)(*B)->data; 113 bctx->diag = actx->diag; 114 } 115 PetscFunctionReturn(PETSC_SUCCESS); 116 } 117 118 static PetscErrorCode MatMissingDiagonal_ConstantDiagonal(Mat mat, PetscBool *missing, PetscInt *dd) 119 { 120 PetscFunctionBegin; 121 *missing = PETSC_FALSE; 122 PetscFunctionReturn(PETSC_SUCCESS); 123 } 124 125 static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat) 126 { 127 PetscFunctionBegin; 128 PetscCall(PetscFree(mat->data)); 129 mat->structural_symmetry_eternal = PETSC_FALSE; 130 mat->symmetry_eternal = PETSC_FALSE; 131 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConstantDiagonalGetConstant_C", NULL)); 132 PetscFunctionReturn(PETSC_SUCCESS); 133 } 134 135 static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer) 136 { 137 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 138 PetscBool isascii; 139 140 PetscFunctionBegin; 141 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 142 if (isascii) { 143 PetscViewerFormat format; 144 145 PetscCall(PetscViewerGetFormat(viewer, &format)); 146 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS); 147 if (PetscImaginaryPart(ctx->diag) == 0) { 148 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)PetscRealPart(ctx->diag))); 149 } else { 150 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag))); 151 } 152 } 153 PetscFunctionReturn(PETSC_SUCCESS); 154 } 155 156 static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y) 157 { 158 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 159 160 PetscFunctionBegin; 161 PetscCall(VecAXPBY(y, ctx->diag, 0.0, x)); 162 PetscFunctionReturn(PETSC_SUCCESS); 163 } 164 165 static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y) 166 { 167 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 168 169 PetscFunctionBegin; 170 PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x)); 171 PetscFunctionReturn(PETSC_SUCCESS); 172 } 173 174 static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x) 175 { 176 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 177 178 PetscFunctionBegin; 179 PetscCall(VecSet(x, ctx->diag)); 180 PetscFunctionReturn(PETSC_SUCCESS); 181 } 182 183 static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a) 184 { 185 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 186 187 PetscFunctionBegin; 188 ctx->diag += a; 189 PetscFunctionReturn(PETSC_SUCCESS); 190 } 191 192 static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a) 193 { 194 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 195 196 PetscFunctionBegin; 197 ctx->diag *= a; 198 PetscFunctionReturn(PETSC_SUCCESS); 199 } 200 201 static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y) 202 { 203 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 204 205 PetscFunctionBegin; 206 ctx->diag = 0.0; 207 PetscFunctionReturn(PETSC_SUCCESS); 208 } 209 210 static PetscErrorCode MatConjugate_ConstantDiagonal(Mat Y) 211 { 212 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 213 214 PetscFunctionBegin; 215 ctx->diag = PetscConj(ctx->diag); 216 PetscFunctionReturn(PETSC_SUCCESS); 217 } 218 219 static PetscErrorCode MatTranspose_ConstantDiagonal(Mat A, MatReuse reuse, Mat *matout) 220 { 221 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data; 222 223 PetscFunctionBegin; 224 if (reuse == MAT_INPLACE_MATRIX) { 225 PetscLayout tmplayout = A->rmap; 226 227 A->rmap = A->cmap; 228 A->cmap = tmplayout; 229 } else { 230 if (reuse == MAT_INITIAL_MATRIX) { 231 PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N, ctx->diag, matout)); 232 } else { 233 PetscCall(MatZeroEntries(*matout)); 234 PetscCall(MatShift(*matout, ctx->diag)); 235 } 236 } 237 PetscFunctionReturn(PETSC_SUCCESS); 238 } 239 240 static PetscErrorCode MatSetRandom_ConstantDiagonal(Mat A, PetscRandom rand) 241 { 242 PetscMPIInt rank; 243 MPI_Comm comm; 244 PetscScalar v = 0.0; 245 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data; 246 247 PetscFunctionBegin; 248 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 249 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 250 if (!rank) PetscCall(PetscRandomGetValue(rand, &v)); 251 PetscCallMPI(MPI_Bcast(&v, 1, MPIU_SCALAR, 0, comm)); 252 ctx->diag = v; 253 PetscFunctionReturn(PETSC_SUCCESS); 254 } 255 256 static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x) 257 { 258 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data; 259 260 PetscFunctionBegin; 261 if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 262 else matin->factorerrortype = MAT_FACTOR_NOERROR; 263 PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b)); 264 PetscFunctionReturn(PETSC_SUCCESS); 265 } 266 267 static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y) 268 { 269 PetscFunctionBegin; 270 PetscCall(MatSolve_ConstantDiagonal(matin, x, y)); 271 PetscFunctionReturn(PETSC_SUCCESS); 272 } 273 274 static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info) 275 { 276 PetscFunctionBegin; 277 info->block_size = 1.0; 278 info->nz_allocated = 1.0; 279 info->nz_used = 1.0; 280 info->nz_unneeded = 0.0; 281 info->assemblies = A->num_ass; 282 info->mallocs = 0.0; 283 info->memory = 0; /* REVIEW ME */ 284 if (A->factortype) { 285 info->fill_ratio_given = 1.0; 286 info->fill_ratio_needed = 1.0; 287 info->factor_mallocs = 0.0; 288 } else { 289 info->fill_ratio_given = 0; 290 info->fill_ratio_needed = 0; 291 info->factor_mallocs = 0; 292 } 293 PetscFunctionReturn(PETSC_SUCCESS); 294 } 295 296 /*@ 297 MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal 298 299 Collective 300 301 Input Parameters: 302 + comm - MPI communicator 303 . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given) 304 This value should be the same as the local size used in creating the 305 y vector for the matrix-vector product y = Ax. 306 . n - This value should be the same as the local size used in creating the 307 x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have 308 calculated if `N` is given) For square matrices n is almost always `m`. 309 . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given) 310 . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given) 311 - diag - the diagonal value 312 313 Output Parameter: 314 . J - the diagonal matrix 315 316 Level: advanced 317 318 Notes: 319 Only supports square matrices with the same number of local rows and columns 320 321 .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()` 322 @*/ 323 PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J) 324 { 325 PetscFunctionBegin; 326 PetscCall(MatCreate(comm, J)); 327 PetscCall(MatSetSizes(*J, m, n, M, N)); 328 PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL)); 329 PetscCall(MatShift(*J, diag)); 330 PetscCall(MatSetUp(*J)); 331 PetscFunctionReturn(PETSC_SUCCESS); 332 } 333 334 /*@ 335 MatConstantDiagonalGetConstant - Get the scalar constant of a constant diagonal matrix 336 337 Not collective 338 339 Input Parameter: 340 . mat - a `MATCONSTANTDIAGONAL` 341 342 Output Parameter: 343 . value - the scalar value 344 345 Level: developer 346 347 .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL` 348 @*/ 349 PetscErrorCode MatConstantDiagonalGetConstant(Mat mat, PetscScalar *value) 350 { 351 PetscFunctionBegin; 352 PetscUseMethod(mat, "MatConstantDiagonalGetConstant_C", (Mat, PetscScalar *), (mat, value)); 353 PetscFunctionReturn(PETSC_SUCCESS); 354 } 355 356 static PetscErrorCode MatConstantDiagonalGetConstant_ConstantDiagonal(Mat mat, PetscScalar *value) 357 { 358 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data; 359 360 PetscFunctionBegin; 361 *value = ctx->diag; 362 PetscFunctionReturn(PETSC_SUCCESS); 363 } 364 365 /*MC 366 MATCONSTANTDIAGONAL - "constant-diagonal" - A diagonal matrix type with a uniform value 367 along the diagonal. 368 369 Level: advanced 370 371 .seealso: [](ch_matrices), `Mat`, `MatCreateConstantDiagonal()` 372 M*/ 373 PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A) 374 { 375 Mat_ConstantDiagonal *ctx; 376 377 PetscFunctionBegin; 378 PetscCall(PetscNew(&ctx)); 379 ctx->diag = 0.0; 380 A->data = (void *)ctx; 381 382 A->assembled = PETSC_TRUE; 383 A->preallocated = PETSC_TRUE; 384 A->structurally_symmetric = PETSC_BOOL3_TRUE; 385 A->structural_symmetry_eternal = PETSC_TRUE; 386 A->symmetric = PETSC_BOOL3_TRUE; 387 if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE; 388 A->symmetry_eternal = PETSC_TRUE; 389 390 A->ops->mult = MatMult_ConstantDiagonal; 391 A->ops->multadd = MatMultAdd_ConstantDiagonal; 392 A->ops->multtranspose = MatMult_ConstantDiagonal; 393 A->ops->multtransposeadd = MatMultAdd_ConstantDiagonal; 394 A->ops->multhermitiantranspose = MatMultHermitianTranspose_ConstantDiagonal; 395 A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal; 396 A->ops->solve = MatSolve_ConstantDiagonal; 397 A->ops->solvetranspose = MatSolve_ConstantDiagonal; 398 A->ops->norm = MatNorm_ConstantDiagonal; 399 A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal; 400 A->ops->duplicate = MatDuplicate_ConstantDiagonal; 401 A->ops->missingdiagonal = MatMissingDiagonal_ConstantDiagonal; 402 A->ops->getrow = MatGetRow_ConstantDiagonal; 403 A->ops->restorerow = MatRestoreRow_ConstantDiagonal; 404 A->ops->sor = MatSOR_ConstantDiagonal; 405 A->ops->shift = MatShift_ConstantDiagonal; 406 A->ops->scale = MatScale_ConstantDiagonal; 407 A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal; 408 A->ops->view = MatView_ConstantDiagonal; 409 A->ops->zeroentries = MatZeroEntries_ConstantDiagonal; 410 A->ops->destroy = MatDestroy_ConstantDiagonal; 411 A->ops->getinfo = MatGetInfo_ConstantDiagonal; 412 A->ops->equal = MatEqual_ConstantDiagonal; 413 A->ops->axpy = MatAXPY_ConstantDiagonal; 414 A->ops->setrandom = MatSetRandom_ConstantDiagonal; 415 A->ops->conjugate = MatConjugate_ConstantDiagonal; 416 A->ops->transpose = MatTranspose_ConstantDiagonal; 417 418 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL)); 419 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConstantDiagonalGetConstant_C", MatConstantDiagonalGetConstant_ConstantDiagonal)); 420 PetscFunctionReturn(PETSC_SUCCESS); 421 } 422 423 static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info) 424 { 425 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data; 426 427 PetscFunctionBegin; 428 if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 429 else fact->factorerrortype = MAT_FACTOR_NOERROR; 430 fctx->diag = 1.0 / actx->diag; 431 fact->ops->solve = MatMult_ConstantDiagonal; 432 PetscFunctionReturn(PETSC_SUCCESS); 433 } 434 435 static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info) 436 { 437 PetscFunctionBegin; 438 fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal; 439 PetscFunctionReturn(PETSC_SUCCESS); 440 } 441 442 static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info) 443 { 444 PetscFunctionBegin; 445 fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal; 446 PetscFunctionReturn(PETSC_SUCCESS); 447 } 448 449 PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B) 450 { 451 PetscInt n = A->rmap->n, N = A->rmap->N; 452 453 PetscFunctionBegin; 454 PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B)); 455 456 (*B)->factortype = ftype; 457 (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 458 (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 459 (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 460 (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 461 462 (*B)->ops->shift = NULL; 463 (*B)->ops->scale = NULL; 464 (*B)->ops->mult = NULL; 465 (*B)->ops->sor = NULL; 466 (*B)->ops->zeroentries = NULL; 467 468 PetscCall(PetscFree((*B)->solvertype)); 469 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype)); 470 PetscFunctionReturn(PETSC_SUCCESS); 471 } 472