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