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