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) PetscCall(PetscFree(*cols)); 40 if (vals) PetscCall(PetscFree(*vals)); 41 PetscFunctionReturn(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 117 } 118 119 static PetscErrorCode MatMissingDiagonal_ConstantDiagonal(Mat mat,PetscBool *missing,PetscInt *dd) 120 { 121 PetscFunctionBegin; 122 *missing = PETSC_FALSE; 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat) 127 { 128 PetscFunctionBegin; 129 PetscCall(PetscFree(mat->data)); 130 PetscFunctionReturn(0); 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(0); 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(0); 152 } 153 154 static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J,MatAssemblyType mt) 155 { 156 PetscFunctionBegin; 157 PetscFunctionReturn(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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 = ((PetscObject)A)->mem; 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(0); 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: `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()` 264 265 @*/ 266 PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscScalar diag,Mat *J) 267 { 268 PetscFunctionBegin; 269 PetscCall(MatCreate(comm,J)); 270 PetscCall(MatSetSizes(*J,m,n,M,N)); 271 PetscCall(MatSetType(*J,MATCONSTANTDIAGONAL)); 272 PetscCall(MatShift(*J,diag)); 273 PetscCall(MatSetUp(*J)); 274 PetscFunctionReturn(0); 275 } 276 277 PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A) 278 { 279 Mat_ConstantDiagonal *ctx; 280 281 PetscFunctionBegin; 282 PetscCall(PetscNew(&ctx)); 283 ctx->diag = 0.0; 284 A->data = (void*)ctx; 285 286 A->assembled = PETSC_TRUE; 287 A->preallocated = PETSC_TRUE; 288 289 A->ops->mult = MatMult_ConstantDiagonal; 290 A->ops->multadd = MatMultAdd_ConstantDiagonal; 291 A->ops->multtranspose = MatMultTranspose_ConstantDiagonal; 292 A->ops->multtransposeadd = MatMultTransposeAdd_ConstantDiagonal; 293 A->ops->norm = MatNorm_ConstantDiagonal; 294 A->ops->createsubmatrices= MatCreateSubMatrices_ConstantDiagonal; 295 A->ops->duplicate = MatDuplicate_ConstantDiagonal; 296 A->ops->missingdiagonal = MatMissingDiagonal_ConstantDiagonal; 297 A->ops->getrow = MatGetRow_ConstantDiagonal; 298 A->ops->restorerow = MatRestoreRow_ConstantDiagonal; 299 A->ops->sor = MatSOR_ConstantDiagonal; 300 A->ops->shift = MatShift_ConstantDiagonal; 301 A->ops->scale = MatScale_ConstantDiagonal; 302 A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal; 303 A->ops->view = MatView_ConstantDiagonal; 304 A->ops->zeroentries = MatZeroEntries_ConstantDiagonal; 305 A->ops->assemblyend = MatAssemblyEnd_ConstantDiagonal; 306 A->ops->destroy = MatDestroy_ConstantDiagonal; 307 A->ops->getinfo = MatGetInfo_ConstantDiagonal; 308 A->ops->axpy = MatAXPY_ConstantDiagonal; 309 310 PetscCall(PetscObjectChangeTypeName((PetscObject)A,MATCONSTANTDIAGONAL)); 311 PetscFunctionReturn(0); 312 } 313 314 static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact,Mat A,const MatFactorInfo *info) 315 { 316 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal*)A->data,*fctx = (Mat_ConstantDiagonal*)fact->data; 317 318 PetscFunctionBegin; 319 if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 320 else fact->factorerrortype = MAT_FACTOR_NOERROR; 321 fctx->diag = 1.0/actx->diag; 322 fact->ops->solve = MatMult_ConstantDiagonal; 323 PetscFunctionReturn(0); 324 } 325 326 static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info) 327 { 328 PetscFunctionBegin; 329 fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal; 330 PetscFunctionReturn(0); 331 } 332 333 static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact,Mat A,IS isrow,const MatFactorInfo *info) 334 { 335 PetscFunctionBegin; 336 fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal; 337 PetscFunctionReturn(0); 338 } 339 340 PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A,MatFactorType ftype,Mat *B) 341 { 342 PetscInt n = A->rmap->n, N = A->rmap->N; 343 344 PetscFunctionBegin; 345 PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A),n,n,N,N,0,B)); 346 347 (*B)->factortype = ftype; 348 (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 349 (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 350 (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 351 (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 352 353 (*B)->ops->shift = NULL; 354 (*B)->ops->scale = NULL; 355 (*B)->ops->mult = NULL; 356 (*B)->ops->sor = NULL; 357 (*B)->ops->zeroentries = NULL; 358 359 PetscCall(PetscFree((*B)->solvertype)); 360 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype)); 361 PetscFunctionReturn(0); 362 } 363