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