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 if (type == NORM_FROBENIUS || type == NORM_2 || type == NORM_1 || type == NORM_INFINITY) *nrm = PetscAbsScalar(ctx->diag); 84 else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm"); 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 PetscFunctionReturn(PETSC_SUCCESS); 132 } 133 134 static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer) 135 { 136 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 137 PetscBool iascii; 138 139 PetscFunctionBegin; 140 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 141 if (iascii) { 142 PetscViewerFormat format; 143 144 PetscCall(PetscViewerGetFormat(viewer, &format)); 145 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS); 146 #if defined(PETSC_USE_COMPLEX) 147 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag))); 148 #else 149 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)ctx->diag)); 150 #endif 151 } 152 PetscFunctionReturn(PETSC_SUCCESS); 153 } 154 155 static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y) 156 { 157 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 158 159 PetscFunctionBegin; 160 PetscCall(VecAXPBY(y, ctx->diag, 0.0, x)); 161 PetscFunctionReturn(PETSC_SUCCESS); 162 } 163 164 static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y) 165 { 166 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 167 168 PetscFunctionBegin; 169 PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x)); 170 PetscFunctionReturn(PETSC_SUCCESS); 171 } 172 173 static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x) 174 { 175 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 176 177 PetscFunctionBegin; 178 PetscCall(VecSet(x, ctx->diag)); 179 PetscFunctionReturn(PETSC_SUCCESS); 180 } 181 182 static PetscErrorCode MatShift_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 MatScale_ConstantDiagonal(Mat Y, PetscScalar a) 192 { 193 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 194 195 PetscFunctionBegin; 196 ctx->diag *= a; 197 PetscFunctionReturn(PETSC_SUCCESS); 198 } 199 200 static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y) 201 { 202 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 203 204 PetscFunctionBegin; 205 ctx->diag = 0.0; 206 PetscFunctionReturn(PETSC_SUCCESS); 207 } 208 209 static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x) 210 { 211 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data; 212 213 PetscFunctionBegin; 214 if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 215 else matin->factorerrortype = MAT_FACTOR_NOERROR; 216 PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b)); 217 PetscFunctionReturn(PETSC_SUCCESS); 218 } 219 220 static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y) 221 { 222 PetscFunctionBegin; 223 PetscCall(MatSolve_ConstantDiagonal(matin, x, y)); 224 PetscFunctionReturn(PETSC_SUCCESS); 225 } 226 227 static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info) 228 { 229 PetscFunctionBegin; 230 info->block_size = 1.0; 231 info->nz_allocated = 1.0; 232 info->nz_used = 1.0; 233 info->nz_unneeded = 0.0; 234 info->assemblies = A->num_ass; 235 info->mallocs = 0.0; 236 info->memory = 0; /* REVIEW ME */ 237 if (A->factortype) { 238 info->fill_ratio_given = 1.0; 239 info->fill_ratio_needed = 1.0; 240 info->factor_mallocs = 0.0; 241 } else { 242 info->fill_ratio_given = 0; 243 info->fill_ratio_needed = 0; 244 info->factor_mallocs = 0; 245 } 246 PetscFunctionReturn(PETSC_SUCCESS); 247 } 248 249 /*@ 250 MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal 251 252 Collective 253 254 Input Parameters: 255 + comm - MPI communicator 256 . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given) 257 This value should be the same as the local size used in creating the 258 y vector for the matrix-vector product y = Ax. 259 . n - This value should be the same as the local size used in creating the 260 x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have 261 calculated if `N` is given) For square matrices n is almost always `m`. 262 . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given) 263 . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given) 264 - diag - the diagonal value 265 266 Output Parameter: 267 . J - the diagonal matrix 268 269 Level: advanced 270 271 Notes: 272 Only supports square matrices with the same number of local rows and columns 273 274 .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()` 275 @*/ 276 PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J) 277 { 278 PetscFunctionBegin; 279 PetscCall(MatCreate(comm, J)); 280 PetscCall(MatSetSizes(*J, m, n, M, N)); 281 PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL)); 282 PetscCall(MatShift(*J, diag)); 283 PetscCall(MatSetUp(*J)); 284 PetscFunctionReturn(PETSC_SUCCESS); 285 } 286 287 /*MC 288 MATCONSTANTDIAGONAL - "constant-diagonal" - A diagonal matrix type with a uniform value 289 along the diagonal. 290 291 Level: advanced 292 293 .seealso: [](ch_matrices), `Mat`, `MatCreateConstantDiagonal()` 294 M*/ 295 PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A) 296 { 297 Mat_ConstantDiagonal *ctx; 298 299 PetscFunctionBegin; 300 PetscCall(PetscNew(&ctx)); 301 ctx->diag = 0.0; 302 A->data = (void *)ctx; 303 304 A->assembled = PETSC_TRUE; 305 A->preallocated = PETSC_TRUE; 306 A->structurally_symmetric = PETSC_BOOL3_TRUE; 307 A->structural_symmetry_eternal = PETSC_TRUE; 308 A->symmetric = PETSC_BOOL3_TRUE; 309 if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE; 310 A->symmetry_eternal = PETSC_TRUE; 311 312 A->ops->mult = MatMult_ConstantDiagonal; 313 A->ops->multadd = MatMultAdd_ConstantDiagonal; 314 A->ops->multtranspose = MatMult_ConstantDiagonal; 315 A->ops->multtransposeadd = MatMultAdd_ConstantDiagonal; 316 A->ops->multhermitiantranspose = MatMultHermitianTranspose_ConstantDiagonal; 317 A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal; 318 A->ops->solve = MatSolve_ConstantDiagonal; 319 A->ops->solvetranspose = MatSolve_ConstantDiagonal; 320 A->ops->norm = MatNorm_ConstantDiagonal; 321 A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal; 322 A->ops->duplicate = MatDuplicate_ConstantDiagonal; 323 A->ops->missingdiagonal = MatMissingDiagonal_ConstantDiagonal; 324 A->ops->getrow = MatGetRow_ConstantDiagonal; 325 A->ops->restorerow = MatRestoreRow_ConstantDiagonal; 326 A->ops->sor = MatSOR_ConstantDiagonal; 327 A->ops->shift = MatShift_ConstantDiagonal; 328 A->ops->scale = MatScale_ConstantDiagonal; 329 A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal; 330 A->ops->view = MatView_ConstantDiagonal; 331 A->ops->zeroentries = MatZeroEntries_ConstantDiagonal; 332 A->ops->destroy = MatDestroy_ConstantDiagonal; 333 A->ops->getinfo = MatGetInfo_ConstantDiagonal; 334 A->ops->equal = MatEqual_ConstantDiagonal; 335 A->ops->axpy = MatAXPY_ConstantDiagonal; 336 337 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL)); 338 PetscFunctionReturn(PETSC_SUCCESS); 339 } 340 341 static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info) 342 { 343 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data; 344 345 PetscFunctionBegin; 346 if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 347 else fact->factorerrortype = MAT_FACTOR_NOERROR; 348 fctx->diag = 1.0 / actx->diag; 349 fact->ops->solve = MatMult_ConstantDiagonal; 350 PetscFunctionReturn(PETSC_SUCCESS); 351 } 352 353 static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info) 354 { 355 PetscFunctionBegin; 356 fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal; 357 PetscFunctionReturn(PETSC_SUCCESS); 358 } 359 360 static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info) 361 { 362 PetscFunctionBegin; 363 fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal; 364 PetscFunctionReturn(PETSC_SUCCESS); 365 } 366 367 PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B) 368 { 369 PetscInt n = A->rmap->n, N = A->rmap->N; 370 371 PetscFunctionBegin; 372 PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B)); 373 374 (*B)->factortype = ftype; 375 (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 376 (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 377 (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 378 (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 379 380 (*B)->ops->shift = NULL; 381 (*B)->ops->scale = NULL; 382 (*B)->ops->mult = NULL; 383 (*B)->ops->sor = NULL; 384 (*B)->ops->zeroentries = NULL; 385 386 PetscCall(PetscFree((*B)->solvertype)); 387 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype)); 388 PetscFunctionReturn(PETSC_SUCCESS); 389 } 390