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 MatDestroy_ConstantDiagonal(Mat mat) 119 { 120 PetscFunctionBegin; 121 PetscCall(PetscFree(mat->data)); 122 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConstantDiagonalGetConstant_C", NULL)); 123 PetscFunctionReturn(PETSC_SUCCESS); 124 } 125 126 static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer) 127 { 128 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 129 PetscBool isascii; 130 131 PetscFunctionBegin; 132 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 133 if (isascii) { 134 PetscViewerFormat format; 135 136 PetscCall(PetscViewerGetFormat(viewer, &format)); 137 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS); 138 if (PetscImaginaryPart(ctx->diag) == 0) { 139 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)PetscRealPart(ctx->diag))); 140 } else { 141 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag))); 142 } 143 } 144 PetscFunctionReturn(PETSC_SUCCESS); 145 } 146 147 static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y) 148 { 149 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 150 151 PetscFunctionBegin; 152 PetscCall(VecAXPBY(y, ctx->diag, 0.0, x)); 153 PetscFunctionReturn(PETSC_SUCCESS); 154 } 155 156 static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y) 157 { 158 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 159 160 PetscFunctionBegin; 161 PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x)); 162 PetscFunctionReturn(PETSC_SUCCESS); 163 } 164 165 static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x) 166 { 167 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 168 169 PetscFunctionBegin; 170 PetscCall(VecSet(x, ctx->diag)); 171 PetscFunctionReturn(PETSC_SUCCESS); 172 } 173 174 static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a) 175 { 176 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 177 178 PetscFunctionBegin; 179 ctx->diag += a; 180 PetscFunctionReturn(PETSC_SUCCESS); 181 } 182 183 static PetscErrorCode MatScale_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 MatZeroEntries_ConstantDiagonal(Mat Y) 193 { 194 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 195 196 PetscFunctionBegin; 197 ctx->diag = 0.0; 198 PetscFunctionReturn(PETSC_SUCCESS); 199 } 200 201 static PetscErrorCode MatConjugate_ConstantDiagonal(Mat Y) 202 { 203 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 204 205 PetscFunctionBegin; 206 ctx->diag = PetscConj(ctx->diag); 207 PetscFunctionReturn(PETSC_SUCCESS); 208 } 209 210 static PetscErrorCode MatTranspose_ConstantDiagonal(Mat A, MatReuse reuse, Mat *matout) 211 { 212 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data; 213 214 PetscFunctionBegin; 215 if (reuse == MAT_INPLACE_MATRIX) { 216 PetscLayout tmplayout = A->rmap; 217 218 A->rmap = A->cmap; 219 A->cmap = tmplayout; 220 } else { 221 if (reuse == MAT_INITIAL_MATRIX) { 222 PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N, ctx->diag, matout)); 223 } else { 224 PetscCall(MatZeroEntries(*matout)); 225 PetscCall(MatShift(*matout, ctx->diag)); 226 } 227 } 228 PetscFunctionReturn(PETSC_SUCCESS); 229 } 230 231 static PetscErrorCode MatSetRandom_ConstantDiagonal(Mat A, PetscRandom rand) 232 { 233 PetscMPIInt rank; 234 MPI_Comm comm; 235 PetscScalar v = 0.0; 236 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data; 237 238 PetscFunctionBegin; 239 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 240 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 241 if (!rank) PetscCall(PetscRandomGetValue(rand, &v)); 242 PetscCallMPI(MPI_Bcast(&v, 1, MPIU_SCALAR, 0, comm)); 243 ctx->diag = v; 244 PetscFunctionReturn(PETSC_SUCCESS); 245 } 246 247 static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x) 248 { 249 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data; 250 251 PetscFunctionBegin; 252 if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 253 else matin->factorerrortype = MAT_FACTOR_NOERROR; 254 PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b)); 255 PetscFunctionReturn(PETSC_SUCCESS); 256 } 257 258 static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y) 259 { 260 PetscFunctionBegin; 261 PetscCall(MatSolve_ConstantDiagonal(matin, x, y)); 262 PetscFunctionReturn(PETSC_SUCCESS); 263 } 264 265 static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info) 266 { 267 PetscFunctionBegin; 268 info->block_size = 1.0; 269 info->nz_allocated = 1.0; 270 info->nz_used = 1.0; 271 info->nz_unneeded = 0.0; 272 info->assemblies = A->num_ass; 273 info->mallocs = 0.0; 274 info->memory = 0; /* REVIEW ME */ 275 if (A->factortype) { 276 info->fill_ratio_given = 1.0; 277 info->fill_ratio_needed = 1.0; 278 info->factor_mallocs = 0.0; 279 } else { 280 info->fill_ratio_given = 0; 281 info->fill_ratio_needed = 0; 282 info->factor_mallocs = 0; 283 } 284 PetscFunctionReturn(PETSC_SUCCESS); 285 } 286 287 /*@ 288 MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal 289 290 Collective 291 292 Input Parameters: 293 + comm - MPI communicator 294 . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given) 295 This value should be the same as the local size used in creating the 296 y vector for the matrix-vector product y = Ax. 297 . n - This value should be the same as the local size used in creating the 298 x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have 299 calculated if `N` is given) For square matrices n is almost always `m`. 300 . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given) 301 . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given) 302 - diag - the diagonal value 303 304 Output Parameter: 305 . J - the diagonal matrix 306 307 Level: advanced 308 309 Notes: 310 Only supports square matrices with the same number of local rows and columns 311 312 .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()` 313 @*/ 314 PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J) 315 { 316 PetscFunctionBegin; 317 PetscCall(MatCreate(comm, J)); 318 PetscCall(MatSetSizes(*J, m, n, M, N)); 319 PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL)); 320 PetscCall(MatShift(*J, diag)); 321 PetscCall(MatSetUp(*J)); 322 PetscFunctionReturn(PETSC_SUCCESS); 323 } 324 325 /*@ 326 MatConstantDiagonalGetConstant - Get the scalar constant of a constant diagonal matrix 327 328 Not collective 329 330 Input Parameter: 331 . mat - a `MATCONSTANTDIAGONAL` 332 333 Output Parameter: 334 . value - the scalar value 335 336 Level: developer 337 338 .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL` 339 @*/ 340 PetscErrorCode MatConstantDiagonalGetConstant(Mat mat, PetscScalar *value) 341 { 342 PetscFunctionBegin; 343 PetscUseMethod(mat, "MatConstantDiagonalGetConstant_C", (Mat, PetscScalar *), (mat, value)); 344 PetscFunctionReturn(PETSC_SUCCESS); 345 } 346 347 static PetscErrorCode MatConstantDiagonalGetConstant_ConstantDiagonal(Mat mat, PetscScalar *value) 348 { 349 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data; 350 351 PetscFunctionBegin; 352 *value = ctx->diag; 353 PetscFunctionReturn(PETSC_SUCCESS); 354 } 355 356 /*MC 357 MATCONSTANTDIAGONAL - "constant-diagonal" - A diagonal matrix type with a uniform value 358 along the diagonal. 359 360 Level: advanced 361 362 .seealso: [](ch_matrices), `Mat`, `MatCreateConstantDiagonal()` 363 M*/ 364 PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A) 365 { 366 Mat_ConstantDiagonal *ctx; 367 368 PetscFunctionBegin; 369 PetscCall(PetscNew(&ctx)); 370 ctx->diag = 0.0; 371 A->data = (void *)ctx; 372 373 A->assembled = PETSC_TRUE; 374 A->preallocated = PETSC_TRUE; 375 A->structurally_symmetric = PETSC_BOOL3_TRUE; 376 A->structural_symmetry_eternal = PETSC_TRUE; 377 A->symmetric = PETSC_BOOL3_TRUE; 378 if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE; 379 A->symmetry_eternal = PETSC_TRUE; 380 381 A->ops->mult = MatMult_ConstantDiagonal; 382 A->ops->multadd = MatMultAdd_ConstantDiagonal; 383 A->ops->multtranspose = MatMult_ConstantDiagonal; 384 A->ops->multtransposeadd = MatMultAdd_ConstantDiagonal; 385 A->ops->multhermitiantranspose = MatMultHermitianTranspose_ConstantDiagonal; 386 A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal; 387 A->ops->solve = MatSolve_ConstantDiagonal; 388 A->ops->solvetranspose = MatSolve_ConstantDiagonal; 389 A->ops->norm = MatNorm_ConstantDiagonal; 390 A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal; 391 A->ops->duplicate = MatDuplicate_ConstantDiagonal; 392 A->ops->getrow = MatGetRow_ConstantDiagonal; 393 A->ops->restorerow = MatRestoreRow_ConstantDiagonal; 394 A->ops->sor = MatSOR_ConstantDiagonal; 395 A->ops->shift = MatShift_ConstantDiagonal; 396 A->ops->scale = MatScale_ConstantDiagonal; 397 A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal; 398 A->ops->view = MatView_ConstantDiagonal; 399 A->ops->zeroentries = MatZeroEntries_ConstantDiagonal; 400 A->ops->destroy = MatDestroy_ConstantDiagonal; 401 A->ops->getinfo = MatGetInfo_ConstantDiagonal; 402 A->ops->equal = MatEqual_ConstantDiagonal; 403 A->ops->axpy = MatAXPY_ConstantDiagonal; 404 A->ops->setrandom = MatSetRandom_ConstantDiagonal; 405 A->ops->conjugate = MatConjugate_ConstantDiagonal; 406 A->ops->transpose = MatTranspose_ConstantDiagonal; 407 408 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL)); 409 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConstantDiagonalGetConstant_C", MatConstantDiagonalGetConstant_ConstantDiagonal)); 410 PetscFunctionReturn(PETSC_SUCCESS); 411 } 412 413 static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info) 414 { 415 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data; 416 417 PetscFunctionBegin; 418 if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 419 else fact->factorerrortype = MAT_FACTOR_NOERROR; 420 fctx->diag = 1.0 / actx->diag; 421 fact->ops->solve = MatMult_ConstantDiagonal; 422 PetscFunctionReturn(PETSC_SUCCESS); 423 } 424 425 static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info) 426 { 427 PetscFunctionBegin; 428 fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal; 429 PetscFunctionReturn(PETSC_SUCCESS); 430 } 431 432 static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info) 433 { 434 PetscFunctionBegin; 435 fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal; 436 PetscFunctionReturn(PETSC_SUCCESS); 437 } 438 439 PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B) 440 { 441 PetscInt n = A->rmap->n, N = A->rmap->N; 442 443 PetscFunctionBegin; 444 PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B)); 445 446 (*B)->factortype = ftype; 447 (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 448 (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 449 (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 450 (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 451 452 (*B)->ops->shift = NULL; 453 (*B)->ops->scale = NULL; 454 (*B)->ops->mult = NULL; 455 (*B)->ops->sor = NULL; 456 (*B)->ops->zeroentries = NULL; 457 458 PetscCall(PetscFree((*B)->solvertype)); 459 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype)); 460 PetscFunctionReturn(PETSC_SUCCESS); 461 } 462