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 (ncols) *ncols = 0; 49 if (cols) PetscCall(PetscFree(*cols)); 50 if (vals) PetscCall(PetscFree(*vals)); 51 PetscFunctionReturn(PETSC_SUCCESS); 52 } 53 54 static PetscErrorCode MatMultAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3) 55 { 56 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data; 57 58 PetscFunctionBegin; 59 if (v2 == v3) { 60 PetscCall(VecAXPBY(v3, ctx->diag, 1.0, v1)); 61 } else { 62 PetscCall(VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2)); 63 } 64 PetscFunctionReturn(PETSC_SUCCESS); 65 } 66 67 static PetscErrorCode MatMultHermitianTransposeAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3) 68 { 69 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data; 70 71 PetscFunctionBegin; 72 if (v2 == v3) { 73 PetscCall(VecAXPBY(v3, PetscConj(ctx->diag), 1.0, v1)); 74 } else { 75 PetscCall(VecAXPBYPCZ(v3, PetscConj(ctx->diag), 1.0, 0.0, v1, v2)); 76 } 77 PetscFunctionReturn(PETSC_SUCCESS); 78 } 79 80 static PetscErrorCode MatNorm_ConstantDiagonal(Mat A, NormType type, PetscReal *nrm) 81 { 82 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data; 83 84 PetscFunctionBegin; 85 if (type == NORM_FROBENIUS || type == NORM_2 || type == NORM_1 || type == NORM_INFINITY) *nrm = PetscAbsScalar(ctx->diag); 86 else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm"); 87 PetscFunctionReturn(PETSC_SUCCESS); 88 } 89 90 static PetscErrorCode MatCreateSubMatrices_ConstantDiagonal(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[]) 91 92 { 93 Mat B; 94 95 PetscFunctionBegin; 96 PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B)); 97 PetscCall(MatCreateSubMatrices(B, n, irow, icol, scall, submat)); 98 PetscCall(MatDestroy(&B)); 99 PetscFunctionReturn(PETSC_SUCCESS); 100 } 101 102 static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B) 103 { 104 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data; 105 106 PetscFunctionBegin; 107 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 108 PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 109 PetscCall(MatSetBlockSizesFromMats(*B, A, A)); 110 PetscCall(MatSetType(*B, MATCONSTANTDIAGONAL)); 111 PetscCall(PetscLayoutReference(A->rmap, &(*B)->rmap)); 112 PetscCall(PetscLayoutReference(A->cmap, &(*B)->cmap)); 113 if (op == MAT_COPY_VALUES) { 114 Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal *)(*B)->data; 115 bctx->diag = actx->diag; 116 } 117 PetscFunctionReturn(PETSC_SUCCESS); 118 } 119 120 static PetscErrorCode MatMissingDiagonal_ConstantDiagonal(Mat mat, PetscBool *missing, PetscInt *dd) 121 { 122 PetscFunctionBegin; 123 *missing = PETSC_FALSE; 124 PetscFunctionReturn(PETSC_SUCCESS); 125 } 126 127 static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat) 128 { 129 PetscFunctionBegin; 130 PetscCall(PetscFree(mat->data)); 131 mat->structural_symmetry_eternal = PETSC_FALSE; 132 mat->symmetry_eternal = PETSC_FALSE; 133 PetscFunctionReturn(PETSC_SUCCESS); 134 } 135 136 static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer) 137 { 138 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 139 PetscBool iascii; 140 141 PetscFunctionBegin; 142 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 143 if (iascii) { 144 PetscViewerFormat format; 145 146 PetscCall(PetscViewerGetFormat(viewer, &format)); 147 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS); 148 #if defined(PETSC_USE_COMPLEX) 149 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag))); 150 #else 151 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)(ctx->diag))); 152 #endif 153 } 154 PetscFunctionReturn(PETSC_SUCCESS); 155 } 156 157 static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J, MatAssemblyType mt) 158 { 159 PetscFunctionBegin; 160 PetscFunctionReturn(PETSC_SUCCESS); 161 } 162 163 static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y) 164 { 165 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 166 167 PetscFunctionBegin; 168 PetscCall(VecAXPBY(y, ctx->diag, 0.0, x)); 169 PetscFunctionReturn(PETSC_SUCCESS); 170 } 171 172 static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y) 173 { 174 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 175 176 PetscFunctionBegin; 177 PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x)); 178 PetscFunctionReturn(PETSC_SUCCESS); 179 } 180 181 static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x) 182 { 183 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 184 185 PetscFunctionBegin; 186 PetscCall(VecSet(x, ctx->diag)); 187 PetscFunctionReturn(PETSC_SUCCESS); 188 } 189 190 static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a) 191 { 192 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 193 194 PetscFunctionBegin; 195 ctx->diag += a; 196 PetscFunctionReturn(PETSC_SUCCESS); 197 } 198 199 static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a) 200 { 201 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 202 203 PetscFunctionBegin; 204 ctx->diag *= a; 205 PetscFunctionReturn(PETSC_SUCCESS); 206 } 207 208 static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y) 209 { 210 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 211 212 PetscFunctionBegin; 213 ctx->diag = 0.0; 214 PetscFunctionReturn(PETSC_SUCCESS); 215 } 216 217 static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x) 218 { 219 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data; 220 221 PetscFunctionBegin; 222 if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 223 else matin->factorerrortype = MAT_FACTOR_NOERROR; 224 PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b)); 225 PetscFunctionReturn(PETSC_SUCCESS); 226 } 227 228 PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y) 229 { 230 PetscFunctionBegin; 231 PetscCall(MatSolve_ConstantDiagonal(matin, x, y)); 232 PetscFunctionReturn(PETSC_SUCCESS); 233 } 234 235 PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info) 236 { 237 PetscFunctionBegin; 238 info->block_size = 1.0; 239 info->nz_allocated = 1.0; 240 info->nz_used = 1.0; 241 info->nz_unneeded = 0.0; 242 info->assemblies = A->num_ass; 243 info->mallocs = 0.0; 244 info->memory = 0; /* REVIEW ME */ 245 if (A->factortype) { 246 info->fill_ratio_given = 1.0; 247 info->fill_ratio_needed = 1.0; 248 info->factor_mallocs = 0.0; 249 } else { 250 info->fill_ratio_given = 0; 251 info->fill_ratio_needed = 0; 252 info->factor_mallocs = 0; 253 } 254 PetscFunctionReturn(PETSC_SUCCESS); 255 } 256 257 /*@ 258 MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal 259 260 Collective 261 262 Input Parameters: 263 + comm - MPI communicator 264 . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given) 265 This value should be the same as the local size used in creating the 266 y vector for the matrix-vector product y = Ax. 267 . n - This value should be the same as the local size used in creating the 268 x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have 269 calculated if `N` is given) For square matrices n is almost always `m`. 270 . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given) 271 . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given) 272 - diag - the diagonal value 273 274 Output Parameter: 275 . J - the diagonal matrix 276 277 Level: advanced 278 279 Notes: 280 Only supports square matrices with the same number of local rows and columns 281 282 .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()` 283 @*/ 284 PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J) 285 { 286 PetscFunctionBegin; 287 PetscCall(MatCreate(comm, J)); 288 PetscCall(MatSetSizes(*J, m, n, M, N)); 289 PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL)); 290 PetscCall(MatShift(*J, diag)); 291 PetscCall(MatSetUp(*J)); 292 PetscFunctionReturn(PETSC_SUCCESS); 293 } 294 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->assemblyend = MatAssemblyEnd_ConstantDiagonal; 333 A->ops->destroy = MatDestroy_ConstantDiagonal; 334 A->ops->getinfo = MatGetInfo_ConstantDiagonal; 335 A->ops->equal = MatEqual_ConstantDiagonal; 336 A->ops->axpy = MatAXPY_ConstantDiagonal; 337 338 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL)); 339 PetscFunctionReturn(PETSC_SUCCESS); 340 } 341 342 static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info) 343 { 344 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data; 345 346 PetscFunctionBegin; 347 if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 348 else fact->factorerrortype = MAT_FACTOR_NOERROR; 349 fctx->diag = 1.0 / actx->diag; 350 fact->ops->solve = MatMult_ConstantDiagonal; 351 PetscFunctionReturn(PETSC_SUCCESS); 352 } 353 354 static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info) 355 { 356 PetscFunctionBegin; 357 fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal; 358 PetscFunctionReturn(PETSC_SUCCESS); 359 } 360 361 static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info) 362 { 363 PetscFunctionBegin; 364 fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal; 365 PetscFunctionReturn(PETSC_SUCCESS); 366 } 367 368 PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B) 369 { 370 PetscInt n = A->rmap->n, N = A->rmap->N; 371 372 PetscFunctionBegin; 373 PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B)); 374 375 (*B)->factortype = ftype; 376 (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 377 (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 378 (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 379 (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 380 381 (*B)->ops->shift = NULL; 382 (*B)->ops->scale = NULL; 383 (*B)->ops->mult = NULL; 384 (*B)->ops->sor = NULL; 385 (*B)->ops->zeroentries = NULL; 386 387 PetscCall(PetscFree((*B)->solvertype)); 388 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype)); 389 PetscFunctionReturn(PETSC_SUCCESS); 390 } 391