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