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