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