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