1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 PetscScalar diag; 6 } Mat_ConstantDiagonal; 7 8 9 /* ----------------------------------------------------------------------------------------*/ 10 static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat) 11 { 12 PetscErrorCode ierr; 13 14 PetscFunctionBegin; 15 ierr = PetscFree(mat->data);CHKERRQ(ierr); 16 PetscFunctionReturn(0); 17 } 18 19 static PetscErrorCode MatView_ConstantDiagonal(Mat J,PetscViewer viewer) 20 { 21 PetscErrorCode ierr; 22 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data; 23 PetscBool iascii; 24 25 PetscFunctionBegin; 26 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 27 if (iascii) { 28 PetscViewerFormat format; 29 30 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 31 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(0); 32 #if !defined(PETSC_USE_COMPLEX) 33 ierr = PetscViewerASCIIPrintf(viewer,"Diagonal value: %g\n",ctx->diag);CHKERRQ(ierr); 34 #else 35 ierr = PetscViewerASCIIPrintf(viewer,"Diagonal value: %g + i %g\n",PetscRealPart(ctx->diag),PetscImaginaryPart(ctx->diag));CHKERRQ(ierr); 36 #endif 37 } 38 PetscFunctionReturn(0); 39 } 40 41 static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J,MatAssemblyType mt) 42 { 43 PetscFunctionBegin; 44 PetscFunctionReturn(0); 45 } 46 47 static PetscErrorCode MatMult_ConstantDiagonal(Mat J,Vec x,Vec y) 48 { 49 PetscErrorCode ierr; 50 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data; 51 52 PetscFunctionBegin; 53 ierr = VecAXPBY(y,ctx->diag,0.0,x);CHKERRQ(ierr); 54 PetscFunctionReturn(0); 55 } 56 57 PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J,Vec x) 58 { 59 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data; 60 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 ierr = VecSet(x,ctx->diag);CHKERRQ(ierr); 64 PetscFunctionReturn(0); 65 } 66 67 static PetscErrorCode MatShift_ConstantDiagonal(Mat Y,PetscScalar a) 68 { 69 PetscErrorCode ierr; 70 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data; 71 72 PetscFunctionBegin; 73 if (a != 0.) {ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr);} 74 ctx->diag += a; 75 PetscFunctionReturn(0); 76 } 77 78 static PetscErrorCode MatScale_ConstantDiagonal(Mat Y,PetscScalar a) 79 { 80 PetscErrorCode ierr; 81 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data; 82 83 PetscFunctionBegin; 84 if (a != 1.) {ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr);} 85 ctx->diag *= a; 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y) 90 { 91 PetscErrorCode ierr; 92 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data; 93 94 PetscFunctionBegin; 95 if (ctx->diag != 0.0) {ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr);} 96 ctx->diag = 0.0; 97 PetscFunctionReturn(0); 98 } 99 100 PetscErrorCode MatSOR_ConstantDiagonal(Mat matin,Vec x,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec y) 101 { 102 PetscErrorCode ierr; 103 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)matin->data; 104 105 PetscFunctionBegin; 106 if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 107 else matin->factorerrortype = MAT_FACTOR_NOERROR; 108 ierr = VecAXPBY(y,1.0/ctx->diag,0.0,x);CHKERRQ(ierr); 109 PetscFunctionReturn(0); 110 } 111 112 PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A,MatInfoType flag,MatInfo *info) 113 { 114 PetscFunctionBegin; 115 info->block_size = 1.0; 116 info->nz_allocated = 1.0; 117 info->nz_used = 1.0; 118 info->nz_unneeded = 0.0; 119 info->assemblies = A->num_ass; 120 info->mallocs = 0.0; 121 info->memory = ((PetscObject)A)->mem; 122 if (A->factortype) { 123 info->fill_ratio_given = 1.0; 124 info->fill_ratio_needed = 1.0; 125 info->factor_mallocs = 0.0; 126 } else { 127 info->fill_ratio_given = 0; 128 info->fill_ratio_needed = 0; 129 info->factor_mallocs = 0; 130 } 131 PetscFunctionReturn(0); 132 } 133 134 /*@ 135 MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal 136 137 Collective 138 139 Input Parameters: 140 + comm - MPI communicator 141 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 142 This value should be the same as the local size used in creating the 143 y vector for the matrix-vector product y = Ax. 144 . n - This value should be the same as the local size used in creating the 145 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 146 calculated if N is given) For square matrices n is almost always m. 147 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 148 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 149 - diag - the diagonal value 150 151 Output Parameter: 152 . J - the diagonal matrix 153 154 Level: advanced 155 156 Notes: 157 Only supports square matrices with the same number of local rows and columns 158 159 .seealso: MatDestroy(), MATCONSTANTDIAGONAL, MatScale(), MatShift(), MatMult(), MatGetDiagonal(), MatGetFactor(), MatSolve() 160 161 @*/ 162 PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscScalar diag,Mat *J) 163 { 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 ierr = MatCreate(comm,J);CHKERRQ(ierr); 168 ierr = MatSetSizes(*J,m,n,M,N);CHKERRQ(ierr); 169 ierr = MatSetType(*J,MATCONSTANTDIAGONAL);CHKERRQ(ierr); 170 ierr = MatShift(*J,diag);CHKERRQ(ierr); 171 ierr = MatSetUp(*J);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175 PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A) 176 { 177 PetscErrorCode ierr; 178 Mat_ConstantDiagonal *ctx; 179 180 PetscFunctionBegin; 181 ierr = PetscNew(&ctx);CHKERRQ(ierr); 182 ctx->diag = 0.0; 183 A->data = (void*)ctx; 184 185 A->assembled = PETSC_TRUE; 186 A->preallocated = PETSC_TRUE; 187 A->ops->mult = MatMult_ConstantDiagonal; 188 A->ops->sor = MatSOR_ConstantDiagonal; 189 A->ops->shift = MatShift_ConstantDiagonal; 190 A->ops->scale = MatScale_ConstantDiagonal; 191 A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal; 192 A->ops->view = MatView_ConstantDiagonal; 193 A->ops->zeroentries = MatZeroEntries_ConstantDiagonal; 194 A->ops->assemblyend = MatAssemblyEnd_ConstantDiagonal; 195 A->ops->destroy = MatDestroy_ConstantDiagonal; 196 A->ops->getinfo = MatGetInfo_ConstantDiagonal; 197 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCONSTANTDIAGONAL);CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200 201 static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact,Mat A,const MatFactorInfo *info) 202 { 203 PetscErrorCode ierr; 204 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal*)A->data,*fctx = (Mat_ConstantDiagonal*)fact->data; 205 206 PetscFunctionBegin; 207 if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 208 else fact->factorerrortype = MAT_FACTOR_NOERROR; 209 fctx->diag = 1.0/actx->diag; 210 ierr = PetscObjectStateIncrease((PetscObject)fact);CHKERRQ(ierr); 211 fact->ops->solve = MatMult_ConstantDiagonal; 212 PetscFunctionReturn(0); 213 } 214 215 static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info) 216 { 217 PetscFunctionBegin; 218 fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal; 219 PetscFunctionReturn(0); 220 } 221 222 static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact,Mat A,IS isrow,const MatFactorInfo *info) 223 { 224 PetscFunctionBegin; 225 fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal; PetscFunctionReturn(0); 226 PetscFunctionReturn(0); 227 } 228 229 PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A,MatFactorType ftype,Mat *B) 230 { 231 PetscInt n = A->rmap->n, N = A->rmap->N; 232 PetscErrorCode ierr; 233 234 PetscFunctionBegin; 235 ierr = MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A),n,n,N,N,0,B);CHKERRQ(ierr); 236 237 (*B)->factortype = ftype; 238 (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 239 (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 240 (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 241 (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 242 243 (*B)->ops->shift = NULL; 244 (*B)->ops->scale = NULL; 245 (*B)->ops->mult = NULL; 246 (*B)->ops->sor = NULL; 247 (*B)->ops->zeroentries = NULL; 248 249 ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr); 250 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr); 251 PetscFunctionReturn(0); 252 } 253