1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 PetscScalar diag; 6 } Mat_ConstantDiagonal; 7 8 static PetscErrorCode MatAXPY_ConstantDiagonal(Mat Y, PetscScalar a, Mat X, MatStructure str) 9 { 10 Mat_ConstantDiagonal *yctx = (Mat_ConstantDiagonal *)Y->data; 11 Mat_ConstantDiagonal *xctx = (Mat_ConstantDiagonal *)X->data; 12 13 PetscFunctionBegin; 14 yctx->diag += a * xctx->diag; 15 PetscFunctionReturn(PETSC_SUCCESS); 16 } 17 18 static PetscErrorCode MatEqual_ConstantDiagonal(Mat Y, Mat X, PetscBool *equal) 19 { 20 Mat_ConstantDiagonal *yctx = (Mat_ConstantDiagonal *)Y->data; 21 Mat_ConstantDiagonal *xctx = (Mat_ConstantDiagonal *)X->data; 22 23 PetscFunctionBegin; 24 *equal = (yctx->diag == xctx->diag) ? PETSC_TRUE : PETSC_FALSE; 25 PetscFunctionReturn(PETSC_SUCCESS); 26 } 27 28 static PetscErrorCode MatGetRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[]) 29 { 30 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data; 31 32 PetscFunctionBegin; 33 if (ncols) *ncols = 1; 34 if (cols) { 35 PetscCall(PetscMalloc1(1, cols)); 36 (*cols)[0] = row; 37 } 38 if (vals) { 39 PetscCall(PetscMalloc1(1, vals)); 40 (*vals)[0] = ctx->diag; 41 } 42 PetscFunctionReturn(PETSC_SUCCESS); 43 } 44 45 static PetscErrorCode MatRestoreRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[]) 46 { 47 PetscFunctionBegin; 48 if (cols) PetscCall(PetscFree(*cols)); 49 if (vals) PetscCall(PetscFree(*vals)); 50 PetscFunctionReturn(PETSC_SUCCESS); 51 } 52 53 static PetscErrorCode MatMultAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3) 54 { 55 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data; 56 57 PetscFunctionBegin; 58 if (v2 == v3) { 59 PetscCall(VecAXPBY(v3, ctx->diag, 1.0, v1)); 60 } else { 61 PetscCall(VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2)); 62 } 63 PetscFunctionReturn(PETSC_SUCCESS); 64 } 65 66 static PetscErrorCode MatMultHermitianTransposeAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3) 67 { 68 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data; 69 70 PetscFunctionBegin; 71 if (v2 == v3) { 72 PetscCall(VecAXPBY(v3, PetscConj(ctx->diag), 1.0, v1)); 73 } else { 74 PetscCall(VecAXPBYPCZ(v3, PetscConj(ctx->diag), 1.0, 0.0, v1, v2)); 75 } 76 PetscFunctionReturn(PETSC_SUCCESS); 77 } 78 79 static PetscErrorCode MatNorm_ConstantDiagonal(Mat A, NormType type, PetscReal *nrm) 80 { 81 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data; 82 83 PetscFunctionBegin; 84 if (type == NORM_FROBENIUS || type == NORM_2 || type == NORM_1 || type == NORM_INFINITY) *nrm = PetscAbsScalar(ctx->diag); 85 else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm"); 86 PetscFunctionReturn(PETSC_SUCCESS); 87 } 88 89 static PetscErrorCode MatCreateSubMatrices_ConstantDiagonal(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[]) 90 91 { 92 Mat B; 93 94 PetscFunctionBegin; 95 PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B)); 96 PetscCall(MatCreateSubMatrices(B, n, irow, icol, scall, submat)); 97 PetscCall(MatDestroy(&B)); 98 PetscFunctionReturn(PETSC_SUCCESS); 99 } 100 101 static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B) 102 { 103 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data; 104 105 PetscFunctionBegin; 106 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 107 PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 108 PetscCall(MatSetBlockSizesFromMats(*B, A, A)); 109 PetscCall(MatSetType(*B, MATCONSTANTDIAGONAL)); 110 PetscCall(PetscLayoutReference(A->rmap, &(*B)->rmap)); 111 PetscCall(PetscLayoutReference(A->cmap, &(*B)->cmap)); 112 if (op == MAT_COPY_VALUES) { 113 Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal *)(*B)->data; 114 bctx->diag = actx->diag; 115 } 116 PetscFunctionReturn(PETSC_SUCCESS); 117 } 118 119 static PetscErrorCode MatMissingDiagonal_ConstantDiagonal(Mat mat, PetscBool *missing, PetscInt *dd) 120 { 121 PetscFunctionBegin; 122 *missing = PETSC_FALSE; 123 PetscFunctionReturn(PETSC_SUCCESS); 124 } 125 126 static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat) 127 { 128 PetscFunctionBegin; 129 PetscCall(PetscFree(mat->data)); 130 mat->structural_symmetry_eternal = PETSC_FALSE; 131 mat->symmetry_eternal = PETSC_FALSE; 132 PetscFunctionReturn(PETSC_SUCCESS); 133 } 134 135 static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer) 136 { 137 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 138 PetscBool iascii; 139 140 PetscFunctionBegin; 141 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 142 if (iascii) { 143 PetscViewerFormat format; 144 145 PetscCall(PetscViewerGetFormat(viewer, &format)); 146 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS); 147 #if defined(PETSC_USE_COMPLEX) 148 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag))); 149 #else 150 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)(ctx->diag))); 151 #endif 152 } 153 PetscFunctionReturn(PETSC_SUCCESS); 154 } 155 156 static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J, MatAssemblyType mt) 157 { 158 PetscFunctionBegin; 159 PetscFunctionReturn(PETSC_SUCCESS); 160 } 161 162 static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y) 163 { 164 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 165 166 PetscFunctionBegin; 167 PetscCall(VecAXPBY(y, ctx->diag, 0.0, x)); 168 PetscFunctionReturn(PETSC_SUCCESS); 169 } 170 171 static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y) 172 { 173 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 174 175 PetscFunctionBegin; 176 PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x)); 177 PetscFunctionReturn(PETSC_SUCCESS); 178 } 179 180 static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x) 181 { 182 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 183 184 PetscFunctionBegin; 185 PetscCall(VecSet(x, ctx->diag)); 186 PetscFunctionReturn(PETSC_SUCCESS); 187 } 188 189 static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a) 190 { 191 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 192 193 PetscFunctionBegin; 194 ctx->diag += a; 195 PetscFunctionReturn(PETSC_SUCCESS); 196 } 197 198 static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a) 199 { 200 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 201 202 PetscFunctionBegin; 203 ctx->diag *= a; 204 PetscFunctionReturn(PETSC_SUCCESS); 205 } 206 207 static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y) 208 { 209 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 210 211 PetscFunctionBegin; 212 ctx->diag = 0.0; 213 PetscFunctionReturn(PETSC_SUCCESS); 214 } 215 216 static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x) 217 { 218 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data; 219 220 PetscFunctionBegin; 221 if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 222 else matin->factorerrortype = MAT_FACTOR_NOERROR; 223 PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b)); 224 PetscFunctionReturn(PETSC_SUCCESS); 225 } 226 227 static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y) 228 { 229 PetscFunctionBegin; 230 PetscCall(MatSolve_ConstantDiagonal(matin, x, y)); 231 PetscFunctionReturn(PETSC_SUCCESS); 232 } 233 234 static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info) 235 { 236 PetscFunctionBegin; 237 info->block_size = 1.0; 238 info->nz_allocated = 1.0; 239 info->nz_used = 1.0; 240 info->nz_unneeded = 0.0; 241 info->assemblies = A->num_ass; 242 info->mallocs = 0.0; 243 info->memory = 0; /* REVIEW ME */ 244 if (A->factortype) { 245 info->fill_ratio_given = 1.0; 246 info->fill_ratio_needed = 1.0; 247 info->factor_mallocs = 0.0; 248 } else { 249 info->fill_ratio_given = 0; 250 info->fill_ratio_needed = 0; 251 info->factor_mallocs = 0; 252 } 253 PetscFunctionReturn(PETSC_SUCCESS); 254 } 255 256 /*@ 257 MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal 258 259 Collective 260 261 Input Parameters: 262 + comm - MPI communicator 263 . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given) 264 This value should be the same as the local size used in creating the 265 y vector for the matrix-vector product y = Ax. 266 . n - This value should be the same as the local size used in creating the 267 x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have 268 calculated if `N` is given) For square matrices n is almost always `m`. 269 . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given) 270 . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given) 271 - diag - the diagonal value 272 273 Output Parameter: 274 . J - the diagonal matrix 275 276 Level: advanced 277 278 Notes: 279 Only supports square matrices with the same number of local rows and columns 280 281 .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()` 282 @*/ 283 PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J) 284 { 285 PetscFunctionBegin; 286 PetscCall(MatCreate(comm, J)); 287 PetscCall(MatSetSizes(*J, m, n, M, N)); 288 PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL)); 289 PetscCall(MatShift(*J, diag)); 290 PetscCall(MatSetUp(*J)); 291 PetscFunctionReturn(PETSC_SUCCESS); 292 } 293 294 PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A) 295 { 296 Mat_ConstantDiagonal *ctx; 297 298 PetscFunctionBegin; 299 PetscCall(PetscNew(&ctx)); 300 ctx->diag = 0.0; 301 A->data = (void *)ctx; 302 303 A->assembled = PETSC_TRUE; 304 A->preallocated = PETSC_TRUE; 305 A->structurally_symmetric = PETSC_BOOL3_TRUE; 306 A->structural_symmetry_eternal = PETSC_TRUE; 307 A->symmetric = PETSC_BOOL3_TRUE; 308 if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE; 309 A->symmetry_eternal = PETSC_TRUE; 310 311 A->ops->mult = MatMult_ConstantDiagonal; 312 A->ops->multadd = MatMultAdd_ConstantDiagonal; 313 A->ops->multtranspose = MatMult_ConstantDiagonal; 314 A->ops->multtransposeadd = MatMultAdd_ConstantDiagonal; 315 A->ops->multhermitiantranspose = MatMultHermitianTranspose_ConstantDiagonal; 316 A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal; 317 A->ops->solve = MatSolve_ConstantDiagonal; 318 A->ops->solvetranspose = MatSolve_ConstantDiagonal; 319 A->ops->norm = MatNorm_ConstantDiagonal; 320 A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal; 321 A->ops->duplicate = MatDuplicate_ConstantDiagonal; 322 A->ops->missingdiagonal = MatMissingDiagonal_ConstantDiagonal; 323 A->ops->getrow = MatGetRow_ConstantDiagonal; 324 A->ops->restorerow = MatRestoreRow_ConstantDiagonal; 325 A->ops->sor = MatSOR_ConstantDiagonal; 326 A->ops->shift = MatShift_ConstantDiagonal; 327 A->ops->scale = MatScale_ConstantDiagonal; 328 A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal; 329 A->ops->view = MatView_ConstantDiagonal; 330 A->ops->zeroentries = MatZeroEntries_ConstantDiagonal; 331 A->ops->assemblyend = MatAssemblyEnd_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