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 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data; 81 82 PetscFunctionBegin; 83 if (a != 1.) PetscObjectStateIncrease((PetscObject)Y); 84 ctx->diag *= a; 85 PetscFunctionReturn(0); 86 } 87 88 static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y) 89 { 90 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data; 91 92 PetscFunctionBegin; 93 if (ctx->diag != 0.0) PetscObjectStateIncrease((PetscObject)Y); 94 ctx->diag = 0.0; 95 PetscFunctionReturn(0); 96 } 97 98 PetscErrorCode MatSOR_ConstantDiagonal(Mat matin,Vec x,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec y) 99 { 100 PetscErrorCode ierr; 101 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)matin->data; 102 103 PetscFunctionBegin; 104 if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 105 else matin->factorerrortype = MAT_FACTOR_NOERROR; 106 ierr = VecAXPBY(y,1.0/ctx->diag,0.0,x);CHKERRQ(ierr); 107 PetscFunctionReturn(0); 108 } 109 110 PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A,MatInfoType flag,MatInfo *info) 111 { 112 PetscFunctionBegin; 113 info->block_size = 1.0; 114 info->nz_allocated = 1.0; 115 info->nz_used = 1.0; 116 info->nz_unneeded = 0.0; 117 info->assemblies = (double)A->num_ass; 118 info->mallocs = 0.0; 119 info->memory = ((PetscObject)A)->mem; 120 if (A->factortype) { 121 info->fill_ratio_given = 1.0; 122 info->fill_ratio_needed = 1.0; 123 info->factor_mallocs = 0.0; 124 } else { 125 info->fill_ratio_given = 0; 126 info->fill_ratio_needed = 0; 127 info->factor_mallocs = 0; 128 } 129 PetscFunctionReturn(0); 130 } 131 132 /*@ 133 MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal 134 135 Collective on MPI_Comm 136 137 Input Parameters: 138 + comm - MPI communicator 139 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 140 This value should be the same as the local size used in creating the 141 y vector for the matrix-vector product y = Ax. 142 . n - This value should be the same as the local size used in creating the 143 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 144 calculated if N is given) For square matrices n is almost always m. 145 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 146 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 147 - diag - the diagonal value 148 149 Output Parameter: 150 . J - the diagonal matrix 151 152 Level: advanced 153 154 Notes: 155 Only supports square matrices with the same number of local rows and columns 156 157 .seealso: MatDestroy(), MATCONSTANTDIAGONAL, MatScale(), MatShift(), MatMult(), MatGetDiagonal(), MatGetFactor(), MatSolve() 158 159 @*/ 160 PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscScalar diag,Mat *J) 161 { 162 PetscErrorCode ierr; 163 164 PetscFunctionBegin; 165 ierr = MatCreate(comm,J);CHKERRQ(ierr); 166 ierr = MatSetSizes(*J,m,n,M,N);CHKERRQ(ierr); 167 ierr = MatSetType(*J,MATCONSTANTDIAGONAL);CHKERRQ(ierr); 168 ierr = MatShift(*J,diag);CHKERRQ(ierr); 169 ierr = MatSetUp(*J);CHKERRQ(ierr); 170 PetscFunctionReturn(0); 171 } 172 173 PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A) 174 { 175 PetscErrorCode ierr; 176 Mat_ConstantDiagonal *ctx; 177 178 PetscFunctionBegin; 179 ierr = PetscNew(&ctx);CHKERRQ(ierr); 180 ctx->diag = 0.0; 181 A->data = (void*)ctx; 182 183 A->assembled = PETSC_TRUE; 184 A->preallocated = PETSC_TRUE; 185 A->ops->mult = MatMult_ConstantDiagonal; 186 A->ops->sor = MatSOR_ConstantDiagonal; 187 A->ops->shift = MatShift_ConstantDiagonal; 188 A->ops->scale = MatScale_ConstantDiagonal; 189 A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal; 190 A->ops->view = MatView_ConstantDiagonal; 191 A->ops->zeroentries = MatZeroEntries_ConstantDiagonal; 192 A->ops->assemblyend = MatAssemblyEnd_ConstantDiagonal; 193 A->ops->destroy = MatDestroy_ConstantDiagonal; 194 A->ops->getinfo = MatGetInfo_ConstantDiagonal; 195 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCONSTANTDIAGONAL);CHKERRQ(ierr); 196 PetscFunctionReturn(0); 197 } 198 199 static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact,Mat A,const MatFactorInfo *info) 200 { 201 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal*)A->data,*fctx = (Mat_ConstantDiagonal*)fact->data; 202 203 PetscFunctionBegin; 204 if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 205 else fact->factorerrortype = MAT_FACTOR_NOERROR; 206 fctx->diag = 1.0/actx->diag; 207 PetscObjectStateIncrease((PetscObject)fact); 208 fact->ops->solve = MatMult_ConstantDiagonal; 209 PetscFunctionReturn(0); 210 } 211 212 static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info) 213 { 214 PetscFunctionBegin; 215 fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal; 216 PetscFunctionReturn(0); 217 } 218 219 static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact,Mat A,IS isrow,const MatFactorInfo *info) 220 { 221 PetscFunctionBegin; 222 fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal; PetscFunctionReturn(0); 223 PetscFunctionReturn(0); 224 } 225 226 PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A,MatFactorType ftype,Mat *B) 227 { 228 PetscInt n = A->rmap->n, N = A->rmap->N; 229 PetscErrorCode ierr; 230 231 PetscFunctionBegin; 232 ierr = MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A),n,n,N,N,0,B);CHKERRQ(ierr); 233 234 (*B)->factortype = ftype; 235 (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 236 (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 237 (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 238 (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 239 240 (*B)->ops->shift = NULL; 241 (*B)->ops->scale = NULL; 242 (*B)->ops->mult = NULL; 243 (*B)->ops->sor = NULL; 244 (*B)->ops->zeroentries = NULL; 245 246 ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr); 247 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250