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 MatGetRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[]) 19 { 20 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data; 21 22 PetscFunctionBegin; 23 if (ncols) *ncols = 1; 24 if (cols) { 25 PetscCall(PetscMalloc1(1, cols)); 26 (*cols)[0] = row; 27 } 28 if (vals) { 29 PetscCall(PetscMalloc1(1, vals)); 30 (*vals)[0] = ctx->diag; 31 } 32 PetscFunctionReturn(PETSC_SUCCESS); 33 } 34 35 static PetscErrorCode MatRestoreRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[]) 36 { 37 PetscFunctionBegin; 38 if (ncols) *ncols = 0; 39 if (cols) PetscCall(PetscFree(*cols)); 40 if (vals) PetscCall(PetscFree(*vals)); 41 PetscFunctionReturn(PETSC_SUCCESS); 42 } 43 44 static PetscErrorCode MatMultTranspose_ConstantDiagonal(Mat A, Vec x, Vec y) 45 { 46 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data; 47 48 PetscFunctionBegin; 49 PetscCall(VecAXPBY(y, ctx->diag, 0.0, x)); 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 MatMultTransposeAdd_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, ctx->diag, 1.0, v1)); 73 } else { 74 PetscCall(VecAXPBYPCZ(v3, 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 PetscFunctionReturn(PETSC_SUCCESS); 131 } 132 133 static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer) 134 { 135 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 136 PetscBool iascii; 137 138 PetscFunctionBegin; 139 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 140 if (iascii) { 141 PetscViewerFormat format; 142 143 PetscCall(PetscViewerGetFormat(viewer, &format)); 144 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS); 145 #if defined(PETSC_USE_COMPLEX) 146 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag))); 147 #else 148 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)(ctx->diag))); 149 #endif 150 } 151 PetscFunctionReturn(PETSC_SUCCESS); 152 } 153 154 static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J, MatAssemblyType mt) 155 { 156 PetscFunctionBegin; 157 PetscFunctionReturn(PETSC_SUCCESS); 158 } 159 160 static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y) 161 { 162 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 163 164 PetscFunctionBegin; 165 PetscCall(VecAXPBY(y, ctx->diag, 0.0, x)); 166 PetscFunctionReturn(PETSC_SUCCESS); 167 } 168 169 PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x) 170 { 171 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data; 172 173 PetscFunctionBegin; 174 PetscCall(VecSet(x, ctx->diag)); 175 PetscFunctionReturn(PETSC_SUCCESS); 176 } 177 178 static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a) 179 { 180 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 181 182 PetscFunctionBegin; 183 ctx->diag += a; 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 187 static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a) 188 { 189 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 190 191 PetscFunctionBegin; 192 ctx->diag *= a; 193 PetscFunctionReturn(PETSC_SUCCESS); 194 } 195 196 static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y) 197 { 198 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data; 199 200 PetscFunctionBegin; 201 ctx->diag = 0.0; 202 PetscFunctionReturn(PETSC_SUCCESS); 203 } 204 205 PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y) 206 { 207 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data; 208 209 PetscFunctionBegin; 210 if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 211 else matin->factorerrortype = MAT_FACTOR_NOERROR; 212 PetscCall(VecAXPBY(y, 1.0 / ctx->diag, 0.0, x)); 213 PetscFunctionReturn(PETSC_SUCCESS); 214 } 215 216 PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info) 217 { 218 PetscFunctionBegin; 219 info->block_size = 1.0; 220 info->nz_allocated = 1.0; 221 info->nz_used = 1.0; 222 info->nz_unneeded = 0.0; 223 info->assemblies = A->num_ass; 224 info->mallocs = 0.0; 225 info->memory = 0; /* REVIEW ME */ 226 if (A->factortype) { 227 info->fill_ratio_given = 1.0; 228 info->fill_ratio_needed = 1.0; 229 info->factor_mallocs = 0.0; 230 } else { 231 info->fill_ratio_given = 0; 232 info->fill_ratio_needed = 0; 233 info->factor_mallocs = 0; 234 } 235 PetscFunctionReturn(PETSC_SUCCESS); 236 } 237 238 /*@ 239 MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal 240 241 Collective 242 243 Input Parameters: 244 + comm - MPI communicator 245 . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given) 246 This value should be the same as the local size used in creating the 247 y vector for the matrix-vector product y = Ax. 248 . n - This value should be the same as the local size used in creating the 249 x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have 250 calculated if `N` is given) For square matrices n is almost always `m`. 251 . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given) 252 . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given) 253 - diag - the diagonal value 254 255 Output Parameter: 256 . J - the diagonal matrix 257 258 Level: advanced 259 260 Notes: 261 Only supports square matrices with the same number of local rows and columns 262 263 .seealso: [](chapter_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()` 264 @*/ 265 PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J) 266 { 267 PetscFunctionBegin; 268 PetscCall(MatCreate(comm, J)); 269 PetscCall(MatSetSizes(*J, m, n, M, N)); 270 PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL)); 271 PetscCall(MatShift(*J, diag)); 272 PetscCall(MatSetUp(*J)); 273 PetscFunctionReturn(PETSC_SUCCESS); 274 } 275 276 PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A) 277 { 278 Mat_ConstantDiagonal *ctx; 279 280 PetscFunctionBegin; 281 PetscCall(PetscNew(&ctx)); 282 ctx->diag = 0.0; 283 A->data = (void *)ctx; 284 285 A->assembled = PETSC_TRUE; 286 A->preallocated = PETSC_TRUE; 287 288 A->ops->mult = MatMult_ConstantDiagonal; 289 A->ops->multadd = MatMultAdd_ConstantDiagonal; 290 A->ops->multtranspose = MatMultTranspose_ConstantDiagonal; 291 A->ops->multtransposeadd = MatMultTransposeAdd_ConstantDiagonal; 292 A->ops->norm = MatNorm_ConstantDiagonal; 293 A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal; 294 A->ops->duplicate = MatDuplicate_ConstantDiagonal; 295 A->ops->missingdiagonal = MatMissingDiagonal_ConstantDiagonal; 296 A->ops->getrow = MatGetRow_ConstantDiagonal; 297 A->ops->restorerow = MatRestoreRow_ConstantDiagonal; 298 A->ops->sor = MatSOR_ConstantDiagonal; 299 A->ops->shift = MatShift_ConstantDiagonal; 300 A->ops->scale = MatScale_ConstantDiagonal; 301 A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal; 302 A->ops->view = MatView_ConstantDiagonal; 303 A->ops->zeroentries = MatZeroEntries_ConstantDiagonal; 304 A->ops->assemblyend = MatAssemblyEnd_ConstantDiagonal; 305 A->ops->destroy = MatDestroy_ConstantDiagonal; 306 A->ops->getinfo = MatGetInfo_ConstantDiagonal; 307 A->ops->axpy = MatAXPY_ConstantDiagonal; 308 309 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL)); 310 PetscFunctionReturn(PETSC_SUCCESS); 311 } 312 313 static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info) 314 { 315 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data; 316 317 PetscFunctionBegin; 318 if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 319 else fact->factorerrortype = MAT_FACTOR_NOERROR; 320 fctx->diag = 1.0 / actx->diag; 321 fact->ops->solve = MatMult_ConstantDiagonal; 322 PetscFunctionReturn(PETSC_SUCCESS); 323 } 324 325 static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info) 326 { 327 PetscFunctionBegin; 328 fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal; 329 PetscFunctionReturn(PETSC_SUCCESS); 330 } 331 332 static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info) 333 { 334 PetscFunctionBegin; 335 fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal; 336 PetscFunctionReturn(PETSC_SUCCESS); 337 } 338 339 PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B) 340 { 341 PetscInt n = A->rmap->n, N = A->rmap->N; 342 343 PetscFunctionBegin; 344 PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B)); 345 346 (*B)->factortype = ftype; 347 (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 348 (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 349 (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 350 (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 351 352 (*B)->ops->shift = NULL; 353 (*B)->ops->scale = NULL; 354 (*B)->ops->mult = NULL; 355 (*B)->ops->sor = NULL; 356 (*B)->ops->zeroentries = NULL; 357 358 PetscCall(PetscFree((*B)->solvertype)); 359 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype)); 360 PetscFunctionReturn(PETSC_SUCCESS); 361 } 362