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