1 2 /* 3 This provides a simple shell for Fortran (and C programmers) to 4 create a very simple matrix class for use with KSP without coding 5 much of anything. 6 */ 7 8 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 9 #include <petsc/private/vecimpl.h> 10 11 struct _MatShellOps { 12 /* 3 */ PetscErrorCode (*mult)(Mat,Vec,Vec); 13 /* 5 */ PetscErrorCode (*multtranspose)(Mat,Vec,Vec); 14 /* 17 */ PetscErrorCode (*getdiagonal)(Mat,Vec); 15 /* 43 */ PetscErrorCode (*copy)(Mat,Mat,MatStructure); 16 /* 60 */ PetscErrorCode (*destroy)(Mat); 17 }; 18 19 typedef struct { 20 struct _MatShellOps ops[1]; 21 22 PetscScalar vscale,vshift; 23 Vec dshift; 24 Vec left,right; 25 Vec left_work,right_work; 26 Vec left_add_work,right_add_work; 27 Mat axpy; 28 PetscScalar axpy_vscale; 29 PetscBool managescalingshifts; /* The user will manage the scaling and shifts for the MATSHELL, not the default */ 30 /* support for ZeroRows/Columns operations */ 31 IS zrows; 32 IS zcols; 33 Vec zvals; 34 Vec zvals_w; 35 VecScatter zvals_sct_r; 36 VecScatter zvals_sct_c; 37 void *ctx; 38 } Mat_Shell; 39 40 41 /* 42 Store and scale values on zeroed rows 43 xx = [x_1, 0], 0 on zeroed columns 44 */ 45 static PetscErrorCode MatShellPreZeroRight(Mat A,Vec x,Vec *xx) 46 { 47 Mat_Shell *shell = (Mat_Shell*)A->data; 48 PetscErrorCode ierr; 49 50 PetscFunctionBegin; 51 *xx = x; 52 if (shell->zrows) { 53 ierr = VecSet(shell->zvals_w,0.0);CHKERRQ(ierr); 54 ierr = VecScatterBegin(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 55 ierr = VecScatterEnd(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 56 ierr = VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);CHKERRQ(ierr); 57 } 58 if (shell->zcols) { 59 if (!shell->right_work) { 60 ierr = MatCreateVecs(A,&shell->right_work,NULL);CHKERRQ(ierr); 61 } 62 ierr = VecCopy(x,shell->right_work);CHKERRQ(ierr); 63 ierr = VecISSet(shell->right_work,shell->zcols,0.0);CHKERRQ(ierr); 64 *xx = shell->right_work; 65 } 66 PetscFunctionReturn(0); 67 } 68 69 /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */ 70 static PetscErrorCode MatShellPostZeroLeft(Mat A,Vec x) 71 { 72 Mat_Shell *shell = (Mat_Shell*)A->data; 73 PetscErrorCode ierr; 74 75 PetscFunctionBegin; 76 if (shell->zrows) { 77 ierr = VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 78 ierr = VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 79 } 80 PetscFunctionReturn(0); 81 } 82 83 /* 84 Store and scale values on zeroed rows 85 xx = [x_1, 0], 0 on zeroed rows 86 */ 87 static PetscErrorCode MatShellPreZeroLeft(Mat A,Vec x,Vec *xx) 88 { 89 Mat_Shell *shell = (Mat_Shell*)A->data; 90 PetscErrorCode ierr; 91 92 PetscFunctionBegin; 93 *xx = NULL; 94 if (!shell->zrows) { 95 *xx = x; 96 } else { 97 if (!shell->left_work) { 98 ierr = MatCreateVecs(A,NULL,&shell->left_work);CHKERRQ(ierr); 99 } 100 ierr = VecCopy(x,shell->left_work);CHKERRQ(ierr); 101 ierr = VecSet(shell->zvals_w,0.0);CHKERRQ(ierr); 102 ierr = VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 103 ierr = VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 104 ierr = VecScatterBegin(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 105 ierr = VecScatterEnd(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 106 ierr = VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);CHKERRQ(ierr); 107 *xx = shell->left_work; 108 } 109 PetscFunctionReturn(0); 110 } 111 112 /* Zero zero-columns contributions, sum contributions from properly scaled values stored in MatShellPreZeroLeft */ 113 static PetscErrorCode MatShellPostZeroRight(Mat A,Vec x) 114 { 115 Mat_Shell *shell = (Mat_Shell*)A->data; 116 PetscErrorCode ierr; 117 118 PetscFunctionBegin; 119 if (shell->zcols) { 120 ierr = VecISSet(x,shell->zcols,0.0);CHKERRQ(ierr); 121 } 122 if (shell->zrows) { 123 ierr = VecScatterBegin(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 124 ierr = VecScatterEnd(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 125 } 126 PetscFunctionReturn(0); 127 } 128 129 /* 130 xx = diag(left)*x 131 */ 132 static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx) 133 { 134 Mat_Shell *shell = (Mat_Shell*)A->data; 135 PetscErrorCode ierr; 136 137 PetscFunctionBegin; 138 *xx = NULL; 139 if (!shell->left) { 140 *xx = x; 141 } else { 142 if (!shell->left_work) {ierr = VecDuplicate(shell->left,&shell->left_work);CHKERRQ(ierr);} 143 ierr = VecPointwiseMult(shell->left_work,x,shell->left);CHKERRQ(ierr); 144 *xx = shell->left_work; 145 } 146 PetscFunctionReturn(0); 147 } 148 149 /* 150 xx = diag(right)*x 151 */ 152 static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx) 153 { 154 Mat_Shell *shell = (Mat_Shell*)A->data; 155 PetscErrorCode ierr; 156 157 PetscFunctionBegin; 158 *xx = NULL; 159 if (!shell->right) { 160 *xx = x; 161 } else { 162 if (!shell->right_work) {ierr = VecDuplicate(shell->right,&shell->right_work);CHKERRQ(ierr);} 163 ierr = VecPointwiseMult(shell->right_work,x,shell->right);CHKERRQ(ierr); 164 *xx = shell->right_work; 165 } 166 PetscFunctionReturn(0); 167 } 168 169 /* 170 x = diag(left)*x 171 */ 172 static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x) 173 { 174 Mat_Shell *shell = (Mat_Shell*)A->data; 175 PetscErrorCode ierr; 176 177 PetscFunctionBegin; 178 if (shell->left) {ierr = VecPointwiseMult(x,x,shell->left);CHKERRQ(ierr);} 179 PetscFunctionReturn(0); 180 } 181 182 /* 183 x = diag(right)*x 184 */ 185 static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x) 186 { 187 Mat_Shell *shell = (Mat_Shell*)A->data; 188 PetscErrorCode ierr; 189 190 PetscFunctionBegin; 191 if (shell->right) {ierr = VecPointwiseMult(x,x,shell->right);CHKERRQ(ierr);} 192 PetscFunctionReturn(0); 193 } 194 195 /* 196 Y = vscale*Y + diag(dshift)*X + vshift*X 197 198 On input Y already contains A*x 199 */ 200 static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y) 201 { 202 Mat_Shell *shell = (Mat_Shell*)A->data; 203 PetscErrorCode ierr; 204 205 PetscFunctionBegin; 206 if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */ 207 PetscInt i,m; 208 const PetscScalar *x,*d; 209 PetscScalar *y; 210 ierr = VecGetLocalSize(X,&m);CHKERRQ(ierr); 211 ierr = VecGetArrayRead(shell->dshift,&d);CHKERRQ(ierr); 212 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 213 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 214 for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i]; 215 ierr = VecRestoreArrayRead(shell->dshift,&d);CHKERRQ(ierr); 216 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 217 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 218 } else { 219 ierr = VecScale(Y,shell->vscale);CHKERRQ(ierr); 220 } 221 if (shell->vshift != 0.0) {ierr = VecAXPY(Y,shell->vshift,X);CHKERRQ(ierr);} /* if test is for non-square matrices */ 222 PetscFunctionReturn(0); 223 } 224 225 /*@ 226 MatShellGetContext - Returns the user-provided context associated with a shell matrix. 227 228 Not Collective 229 230 Input Parameter: 231 . mat - the matrix, should have been created with MatCreateShell() 232 233 Output Parameter: 234 . ctx - the user provided context 235 236 Level: advanced 237 238 Fortran Notes: 239 To use this from Fortran you must write a Fortran interface definition for this 240 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 241 242 .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext() 243 @*/ 244 PetscErrorCode MatShellGetContext(Mat mat,void *ctx) 245 { 246 PetscErrorCode ierr; 247 PetscBool flg; 248 249 PetscFunctionBegin; 250 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 251 PetscValidPointer(ctx,2); 252 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 253 if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx; 254 else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix"); 255 PetscFunctionReturn(0); 256 } 257 258 static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat,PetscInt nr,PetscInt rows[],PetscInt nc,PetscInt cols[],PetscScalar diag,PetscBool rc) 259 { 260 PetscErrorCode ierr; 261 Mat_Shell *shell = (Mat_Shell*)mat->data; 262 Vec x = NULL,b = NULL; 263 IS is1, is2; 264 const PetscInt *ridxs; 265 PetscInt *idxs,*gidxs; 266 PetscInt cum,rst,cst,i; 267 268 PetscFunctionBegin; 269 if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation"); 270 if (!shell->zvals) { 271 ierr = MatCreateVecs(mat,NULL,&shell->zvals);CHKERRQ(ierr); 272 } 273 if (!shell->zvals_w) { 274 ierr = VecDuplicate(shell->zvals,&shell->zvals_w);CHKERRQ(ierr); 275 } 276 ierr = MatGetOwnershipRange(mat,&rst,NULL);CHKERRQ(ierr); 277 ierr = MatGetOwnershipRangeColumn(mat,&cst,NULL);CHKERRQ(ierr); 278 279 /* Expand/create index set of zeroed rows */ 280 ierr = PetscMalloc1(nr,&idxs);CHKERRQ(ierr); 281 for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst; 282 ierr = ISCreateGeneral(PETSC_COMM_SELF,nr,idxs,PETSC_OWN_POINTER,&is1);CHKERRQ(ierr); 283 ierr = ISSort(is1);CHKERRQ(ierr); 284 ierr = VecISSet(shell->zvals,is1,diag);CHKERRQ(ierr); 285 if (shell->zrows) { 286 ierr = ISSum(shell->zrows,is1,&is2);CHKERRQ(ierr); 287 ierr = ISDestroy(&shell->zrows);CHKERRQ(ierr); 288 ierr = ISDestroy(&is1);CHKERRQ(ierr); 289 shell->zrows = is2; 290 } else shell->zrows = is1; 291 292 /* Create scatters for diagonal values communications */ 293 ierr = VecScatterDestroy(&shell->zvals_sct_c);CHKERRQ(ierr); 294 ierr = VecScatterDestroy(&shell->zvals_sct_r);CHKERRQ(ierr); 295 296 /* row scatter: from/to left vector */ 297 ierr = MatCreateVecs(mat,&x,&b);CHKERRQ(ierr); 298 ierr = VecScatterCreate(b,shell->zrows,shell->zvals_w,shell->zrows,&shell->zvals_sct_r);CHKERRQ(ierr); 299 300 /* col scatter: from right vector to left vector */ 301 ierr = ISGetIndices(shell->zrows,&ridxs);CHKERRQ(ierr); 302 ierr = ISGetLocalSize(shell->zrows,&nr);CHKERRQ(ierr); 303 ierr = PetscMalloc1(nr,&gidxs);CHKERRQ(ierr); 304 for (i = 0, cum = 0; i < nr; i++) { 305 if (ridxs[i] >= mat->cmap->N) continue; 306 gidxs[cum] = ridxs[i]; 307 cum++; 308 } 309 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)mat),cum,gidxs,PETSC_OWN_POINTER,&is1);CHKERRQ(ierr); 310 ierr = VecScatterCreate(x,is1,shell->zvals_w,is1,&shell->zvals_sct_c);CHKERRQ(ierr); 311 ierr = ISDestroy(&is1);CHKERRQ(ierr); 312 ierr = VecDestroy(&x);CHKERRQ(ierr); 313 ierr = VecDestroy(&b);CHKERRQ(ierr); 314 315 /* Expand/create index set of zeroed columns */ 316 if (rc) { 317 ierr = PetscMalloc1(nc,&idxs);CHKERRQ(ierr); 318 for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst; 319 ierr = ISCreateGeneral(PETSC_COMM_SELF,nc,idxs,PETSC_OWN_POINTER,&is1);CHKERRQ(ierr); 320 ierr = ISSort(is1);CHKERRQ(ierr); 321 if (shell->zcols) { 322 ierr = ISSum(shell->zcols,is1,&is2);CHKERRQ(ierr); 323 ierr = ISDestroy(&shell->zcols);CHKERRQ(ierr); 324 ierr = ISDestroy(&is1);CHKERRQ(ierr); 325 shell->zcols = is2; 326 } else shell->zcols = is1; 327 } 328 PetscFunctionReturn(0); 329 } 330 331 static PetscErrorCode MatZeroRows_Shell(Mat mat,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 332 { 333 PetscInt nr, *lrows; 334 PetscErrorCode ierr; 335 336 PetscFunctionBegin; 337 if (x && b) { 338 Vec xt; 339 PetscScalar *vals; 340 PetscInt *gcols,i,st,nl,nc; 341 342 ierr = PetscMalloc1(n,&gcols);CHKERRQ(ierr); 343 for (i = 0, nc = 0; i < n; i++) if (rows[i] < mat->cmap->N) gcols[nc++] = rows[i]; 344 345 ierr = MatCreateVecs(mat,&xt,NULL);CHKERRQ(ierr); 346 ierr = VecCopy(x,xt);CHKERRQ(ierr); 347 ierr = PetscCalloc1(nc,&vals);CHKERRQ(ierr); 348 ierr = VecSetValues(xt,nc,gcols,vals,INSERT_VALUES);CHKERRQ(ierr); /* xt = [x1, 0] */ 349 ierr = PetscFree(vals);CHKERRQ(ierr); 350 ierr = VecAssemblyBegin(xt);CHKERRQ(ierr); 351 ierr = VecAssemblyEnd(xt);CHKERRQ(ierr); 352 ierr = VecAYPX(xt,-1.0,x);CHKERRQ(ierr); /* xt = [0, x2] */ 353 354 ierr = VecGetOwnershipRange(xt,&st,NULL);CHKERRQ(ierr); 355 ierr = VecGetLocalSize(xt,&nl);CHKERRQ(ierr); 356 ierr = VecGetArray(xt,&vals);CHKERRQ(ierr); 357 for (i = 0; i < nl; i++) { 358 PetscInt g = i + st; 359 if (g > mat->rmap->N) continue; 360 if (PetscAbsScalar(vals[i]) == 0.0) continue; 361 ierr = VecSetValue(b,g,diag*vals[i],INSERT_VALUES);CHKERRQ(ierr); 362 } 363 ierr = VecRestoreArray(xt,&vals);CHKERRQ(ierr); 364 ierr = VecAssemblyBegin(b);CHKERRQ(ierr); 365 ierr = VecAssemblyEnd(b);CHKERRQ(ierr); /* b = [b1, x2 * diag] */ 366 ierr = VecDestroy(&xt);CHKERRQ(ierr); 367 ierr = PetscFree(gcols);CHKERRQ(ierr); 368 } 369 ierr = PetscLayoutMapLocal(mat->rmap,n,rows,&nr,&lrows,NULL);CHKERRQ(ierr); 370 ierr = MatZeroRowsColumns_Local_Shell(mat,nr,lrows,0,NULL,diag,PETSC_FALSE);CHKERRQ(ierr); 371 ierr = PetscFree(lrows);CHKERRQ(ierr); 372 PetscFunctionReturn(0); 373 } 374 375 static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat,PetscInt n,const PetscInt rowscols[],PetscScalar diag,Vec x,Vec b) 376 { 377 PetscInt *lrows, *lcols; 378 PetscInt nr, nc; 379 PetscBool congruent; 380 PetscErrorCode ierr; 381 382 PetscFunctionBegin; 383 if (x && b) { 384 Vec xt, bt; 385 PetscScalar *vals; 386 PetscInt *grows,*gcols,i,st,nl; 387 388 ierr = PetscMalloc2(n,&grows,n,&gcols);CHKERRQ(ierr); 389 for (i = 0, nr = 0; i < n; i++) if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i]; 390 for (i = 0, nc = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i]; 391 ierr = PetscCalloc1(n,&vals);CHKERRQ(ierr); 392 393 ierr = MatCreateVecs(mat,&xt,&bt);CHKERRQ(ierr); 394 ierr = VecCopy(x,xt);CHKERRQ(ierr); 395 ierr = VecSetValues(xt,nc,gcols,vals,INSERT_VALUES);CHKERRQ(ierr); /* xt = [x1, 0] */ 396 ierr = VecAssemblyBegin(xt);CHKERRQ(ierr); 397 ierr = VecAssemblyEnd(xt);CHKERRQ(ierr); 398 ierr = VecAXPY(xt,-1.0,x);CHKERRQ(ierr); /* xt = [0, -x2] */ 399 ierr = MatMult(mat,xt,bt);CHKERRQ(ierr); /* bt = [-A12*x2,-A22*x2] */ 400 ierr = VecSetValues(bt,nr,grows,vals,INSERT_VALUES);CHKERRQ(ierr); /* bt = [-A12*x2,0] */ 401 ierr = VecAssemblyBegin(bt);CHKERRQ(ierr); 402 ierr = VecAssemblyEnd(bt);CHKERRQ(ierr); 403 ierr = VecAXPY(b,1.0,bt);CHKERRQ(ierr); /* b = [b1 - A12*x2, b2] */ 404 ierr = VecSetValues(bt,nr,grows,vals,INSERT_VALUES);CHKERRQ(ierr); /* b = [b1 - A12*x2, 0] */ 405 ierr = VecAssemblyBegin(bt);CHKERRQ(ierr); 406 ierr = VecAssemblyEnd(bt);CHKERRQ(ierr); 407 ierr = PetscFree(vals);CHKERRQ(ierr); 408 409 ierr = VecGetOwnershipRange(xt,&st,NULL);CHKERRQ(ierr); 410 ierr = VecGetLocalSize(xt,&nl);CHKERRQ(ierr); 411 ierr = VecGetArray(xt,&vals);CHKERRQ(ierr); 412 for (i = 0; i < nl; i++) { 413 PetscInt g = i + st; 414 if (g > mat->rmap->N) continue; 415 if (PetscAbsScalar(vals[i]) == 0.0) continue; 416 ierr = VecSetValue(b,g,-diag*vals[i],INSERT_VALUES);CHKERRQ(ierr); 417 } 418 ierr = VecRestoreArray(xt,&vals);CHKERRQ(ierr); 419 ierr = VecAssemblyBegin(b);CHKERRQ(ierr); 420 ierr = VecAssemblyEnd(b);CHKERRQ(ierr); /* b = [b1 - A12*x2, x2 * diag] */ 421 ierr = VecDestroy(&xt);CHKERRQ(ierr); 422 ierr = VecDestroy(&bt);CHKERRQ(ierr); 423 ierr = PetscFree2(grows,gcols);CHKERRQ(ierr); 424 } 425 ierr = PetscLayoutMapLocal(mat->rmap,n,rowscols,&nr,&lrows,NULL);CHKERRQ(ierr); 426 ierr = MatHasCongruentLayouts(mat,&congruent);CHKERRQ(ierr); 427 if (congruent) { 428 nc = nr; 429 lcols = lrows; 430 } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */ 431 PetscInt i,nt,*t; 432 433 ierr = PetscMalloc1(n,&t);CHKERRQ(ierr); 434 for (i = 0, nt = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i]; 435 ierr = PetscLayoutMapLocal(mat->cmap,nt,t,&nc,&lcols,NULL);CHKERRQ(ierr); 436 ierr = PetscFree(t);CHKERRQ(ierr); 437 } 438 ierr = MatZeroRowsColumns_Local_Shell(mat,nr,lrows,nc,lcols,diag,PETSC_TRUE);CHKERRQ(ierr); 439 if (!congruent) { 440 ierr = PetscFree(lcols);CHKERRQ(ierr); 441 } 442 ierr = PetscFree(lrows);CHKERRQ(ierr); 443 PetscFunctionReturn(0); 444 } 445 446 PetscErrorCode MatDestroy_Shell(Mat mat) 447 { 448 PetscErrorCode ierr; 449 Mat_Shell *shell = (Mat_Shell*)mat->data; 450 451 PetscFunctionBegin; 452 if (shell->ops->destroy) { 453 ierr = (*shell->ops->destroy)(mat);CHKERRQ(ierr); 454 } 455 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 456 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 457 ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr); 458 ierr = VecDestroy(&shell->left_work);CHKERRQ(ierr); 459 ierr = VecDestroy(&shell->right_work);CHKERRQ(ierr); 460 ierr = VecDestroy(&shell->left_add_work);CHKERRQ(ierr); 461 ierr = VecDestroy(&shell->right_add_work);CHKERRQ(ierr); 462 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 463 ierr = VecDestroy(&shell->zvals_w);CHKERRQ(ierr); 464 ierr = VecDestroy(&shell->zvals);CHKERRQ(ierr); 465 ierr = VecScatterDestroy(&shell->zvals_sct_c);CHKERRQ(ierr); 466 ierr = VecScatterDestroy(&shell->zvals_sct_r);CHKERRQ(ierr); 467 ierr = ISDestroy(&shell->zrows);CHKERRQ(ierr); 468 ierr = ISDestroy(&shell->zcols);CHKERRQ(ierr); 469 ierr = PetscFree(mat->data);CHKERRQ(ierr); 470 PetscFunctionReturn(0); 471 } 472 473 PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str) 474 { 475 Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data; 476 PetscErrorCode ierr; 477 PetscBool matflg; 478 479 PetscFunctionBegin; 480 ierr = PetscObjectTypeCompare((PetscObject)B,MATSHELL,&matflg);CHKERRQ(ierr); 481 if (!matflg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_NOTSAMETYPE,"Output matrix must be a MATSHELL"); 482 483 ierr = PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 484 ierr = PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));CHKERRQ(ierr); 485 486 if (shellA->ops->copy) { 487 ierr = (*shellA->ops->copy)(A,B,str);CHKERRQ(ierr); 488 } 489 shellB->vscale = shellA->vscale; 490 shellB->vshift = shellA->vshift; 491 if (shellA->dshift) { 492 if (!shellB->dshift) { 493 ierr = VecDuplicate(shellA->dshift,&shellB->dshift);CHKERRQ(ierr); 494 } 495 ierr = VecCopy(shellA->dshift,shellB->dshift);CHKERRQ(ierr); 496 } else { 497 ierr = VecDestroy(&shellB->dshift);CHKERRQ(ierr); 498 } 499 if (shellA->left) { 500 if (!shellB->left) { 501 ierr = VecDuplicate(shellA->left,&shellB->left);CHKERRQ(ierr); 502 } 503 ierr = VecCopy(shellA->left,shellB->left);CHKERRQ(ierr); 504 } else { 505 ierr = VecDestroy(&shellB->left);CHKERRQ(ierr); 506 } 507 if (shellA->right) { 508 if (!shellB->right) { 509 ierr = VecDuplicate(shellA->right,&shellB->right);CHKERRQ(ierr); 510 } 511 ierr = VecCopy(shellA->right,shellB->right);CHKERRQ(ierr); 512 } else { 513 ierr = VecDestroy(&shellB->right);CHKERRQ(ierr); 514 } 515 ierr = MatDestroy(&shellB->axpy);CHKERRQ(ierr); 516 if (shellA->axpy) { 517 ierr = PetscObjectReference((PetscObject)shellA->axpy);CHKERRQ(ierr); 518 shellB->axpy = shellA->axpy; 519 shellB->axpy_vscale = shellA->axpy_vscale; 520 } 521 if (shellA->zrows) { 522 ierr = ISDuplicate(shellA->zrows,&shellB->zrows);CHKERRQ(ierr); 523 if (shellA->zcols) { 524 ierr = ISDuplicate(shellA->zcols,&shellB->zcols);CHKERRQ(ierr); 525 } 526 ierr = VecDuplicate(shellA->zvals,&shellB->zvals);CHKERRQ(ierr); 527 ierr = VecCopy(shellA->zvals,shellB->zvals);CHKERRQ(ierr); 528 ierr = VecDuplicate(shellA->zvals_w,&shellB->zvals_w);CHKERRQ(ierr); 529 ierr = PetscObjectReference((PetscObject)shellA->zvals_sct_r);CHKERRQ(ierr); 530 ierr = PetscObjectReference((PetscObject)shellA->zvals_sct_c);CHKERRQ(ierr); 531 shellB->zvals_sct_r = shellA->zvals_sct_r; 532 shellB->zvals_sct_c = shellA->zvals_sct_c; 533 } 534 PetscFunctionReturn(0); 535 } 536 537 PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M) 538 { 539 PetscErrorCode ierr; 540 void *ctx; 541 542 PetscFunctionBegin; 543 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr); 544 ierr = MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);CHKERRQ(ierr); 545 if (op != MAT_DO_NOT_COPY_VALUES) { 546 ierr = MatCopy(mat,*M,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 547 } 548 PetscFunctionReturn(0); 549 } 550 551 PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y) 552 { 553 Mat_Shell *shell = (Mat_Shell*)A->data; 554 PetscErrorCode ierr; 555 Vec xx; 556 PetscObjectState instate,outstate; 557 558 PetscFunctionBegin; 559 if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL"); 560 ierr = MatShellPreZeroRight(A,x,&xx);CHKERRQ(ierr); 561 ierr = MatShellPreScaleRight(A,xx,&xx);CHKERRQ(ierr); 562 ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr); 563 ierr = (*shell->ops->mult)(A,xx,y);CHKERRQ(ierr); 564 ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr); 565 if (instate == outstate) { 566 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 567 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 568 } 569 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 570 ierr = MatShellPostScaleLeft(A,y);CHKERRQ(ierr); 571 ierr = MatShellPostZeroLeft(A,y);CHKERRQ(ierr); 572 573 if (shell->axpy) { 574 if (!shell->left_work) {ierr = MatCreateVecs(A,&shell->left_work,NULL);CHKERRQ(ierr);} 575 ierr = MatMult(shell->axpy,x,shell->left_work);CHKERRQ(ierr); 576 ierr = VecAXPY(y,shell->axpy_vscale,shell->left_work);CHKERRQ(ierr); 577 } 578 PetscFunctionReturn(0); 579 } 580 581 PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z) 582 { 583 Mat_Shell *shell = (Mat_Shell*)A->data; 584 PetscErrorCode ierr; 585 586 PetscFunctionBegin; 587 if (y == z) { 588 if (!shell->right_add_work) {ierr = VecDuplicate(z,&shell->right_add_work);CHKERRQ(ierr);} 589 ierr = MatMult(A,x,shell->right_add_work);CHKERRQ(ierr); 590 ierr = VecAXPY(z,1.0,shell->right_add_work);CHKERRQ(ierr); 591 } else { 592 ierr = MatMult(A,x,z);CHKERRQ(ierr); 593 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 594 } 595 PetscFunctionReturn(0); 596 } 597 598 PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y) 599 { 600 Mat_Shell *shell = (Mat_Shell*)A->data; 601 PetscErrorCode ierr; 602 Vec xx; 603 PetscObjectState instate,outstate; 604 605 PetscFunctionBegin; 606 ierr = MatShellPreZeroLeft(A,x,&xx);CHKERRQ(ierr); 607 ierr = MatShellPreScaleLeft(A,xx,&xx);CHKERRQ(ierr); 608 ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr); 609 if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL"); 610 ierr = (*shell->ops->multtranspose)(A,xx,y);CHKERRQ(ierr); 611 ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr); 612 if (instate == outstate) { 613 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 614 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 615 } 616 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 617 ierr = MatShellPostScaleRight(A,y);CHKERRQ(ierr); 618 ierr = MatShellPostZeroRight(A,y);CHKERRQ(ierr); 619 620 if (shell->axpy) { 621 if (!shell->right_work) {ierr = MatCreateVecs(A,NULL,&shell->right_work);CHKERRQ(ierr);} 622 ierr = MatMultTranspose(shell->axpy,x,shell->right_work);CHKERRQ(ierr); 623 ierr = VecAXPY(y,shell->axpy_vscale,shell->right_work);CHKERRQ(ierr); 624 } 625 PetscFunctionReturn(0); 626 } 627 628 PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z) 629 { 630 Mat_Shell *shell = (Mat_Shell*)A->data; 631 PetscErrorCode ierr; 632 633 PetscFunctionBegin; 634 if (y == z) { 635 if (!shell->left_add_work) {ierr = VecDuplicate(z,&shell->left_add_work);CHKERRQ(ierr);} 636 ierr = MatMultTranspose(A,x,shell->left_add_work);CHKERRQ(ierr); 637 ierr = VecAXPY(z,1.0,shell->left_add_work);CHKERRQ(ierr); 638 } else { 639 ierr = MatMultTranspose(A,x,z);CHKERRQ(ierr); 640 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 641 } 642 PetscFunctionReturn(0); 643 } 644 645 /* 646 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 647 */ 648 PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v) 649 { 650 Mat_Shell *shell = (Mat_Shell*)A->data; 651 PetscErrorCode ierr; 652 653 PetscFunctionBegin; 654 if (shell->ops->getdiagonal) { 655 ierr = (*shell->ops->getdiagonal)(A,v);CHKERRQ(ierr); 656 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)"); 657 ierr = VecScale(v,shell->vscale);CHKERRQ(ierr); 658 if (shell->dshift) { 659 ierr = VecAXPY(v,1.0,shell->dshift);CHKERRQ(ierr); 660 } 661 ierr = VecShift(v,shell->vshift);CHKERRQ(ierr); 662 if (shell->left) {ierr = VecPointwiseMult(v,v,shell->left);CHKERRQ(ierr);} 663 if (shell->right) {ierr = VecPointwiseMult(v,v,shell->right);CHKERRQ(ierr);} 664 if (shell->zrows) { 665 ierr = VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 666 ierr = VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 667 } 668 if (shell->axpy) { 669 if (!shell->left_work) {ierr = VecDuplicate(v,&shell->left_work);CHKERRQ(ierr);} 670 ierr = MatGetDiagonal(shell->axpy,shell->left_work);CHKERRQ(ierr); 671 ierr = VecAXPY(v,shell->axpy_vscale,shell->left_work);CHKERRQ(ierr); 672 } 673 PetscFunctionReturn(0); 674 } 675 676 PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a) 677 { 678 Mat_Shell *shell = (Mat_Shell*)Y->data; 679 PetscErrorCode ierr; 680 681 PetscFunctionBegin; 682 if (shell->left || shell->right) { 683 if (!shell->dshift) { 684 ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);CHKERRQ(ierr); 685 ierr = VecSet(shell->dshift,a);CHKERRQ(ierr); 686 } else { 687 if (shell->left) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 688 if (shell->right) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 689 ierr = VecShift(shell->dshift,a);CHKERRQ(ierr); 690 } 691 if (shell->left) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 692 if (shell->right) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 693 } else shell->vshift += a; 694 if (shell->zrows) { 695 ierr = VecShift(shell->zvals,a);CHKERRQ(ierr); 696 } 697 PetscFunctionReturn(0); 698 } 699 700 PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s) 701 { 702 Mat_Shell *shell = (Mat_Shell*)A->data; 703 PetscErrorCode ierr; 704 705 PetscFunctionBegin; 706 if (!shell->dshift) {ierr = VecDuplicate(D,&shell->dshift);CHKERRQ(ierr);} 707 if (shell->left || shell->right) { 708 if (!shell->right_work) {ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);CHKERRQ(ierr);} 709 if (shell->left && shell->right) { 710 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 711 ierr = VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);CHKERRQ(ierr); 712 } else if (shell->left) { 713 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 714 } else { 715 ierr = VecPointwiseDivide(shell->right_work,D,shell->right);CHKERRQ(ierr); 716 } 717 ierr = VecAXPY(shell->dshift,s,shell->right_work);CHKERRQ(ierr); 718 } else { 719 ierr = VecAXPY(shell->dshift,s,D);CHKERRQ(ierr); 720 } 721 PetscFunctionReturn(0); 722 } 723 724 PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins) 725 { 726 Mat_Shell *shell = (Mat_Shell*)A->data; 727 Vec d; 728 PetscErrorCode ierr; 729 730 PetscFunctionBegin; 731 if (ins == INSERT_VALUES) { 732 if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first"); 733 ierr = VecDuplicate(D,&d);CHKERRQ(ierr); 734 ierr = MatGetDiagonal(A,d);CHKERRQ(ierr); 735 ierr = MatDiagonalSet_Shell_Private(A,d,-1.);CHKERRQ(ierr); 736 ierr = MatDiagonalSet_Shell_Private(A,D,1.);CHKERRQ(ierr); 737 ierr = VecDestroy(&d);CHKERRQ(ierr); 738 if (shell->zrows) { 739 ierr = VecCopy(D,shell->zvals);CHKERRQ(ierr); 740 } 741 } else { 742 ierr = MatDiagonalSet_Shell_Private(A,D,1.);CHKERRQ(ierr); 743 if (shell->zrows) { 744 ierr = VecAXPY(shell->zvals,1.0,D);CHKERRQ(ierr); 745 } 746 } 747 PetscFunctionReturn(0); 748 } 749 750 PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a) 751 { 752 Mat_Shell *shell = (Mat_Shell*)Y->data; 753 PetscErrorCode ierr; 754 755 PetscFunctionBegin; 756 shell->vscale *= a; 757 shell->vshift *= a; 758 if (shell->dshift) { 759 ierr = VecScale(shell->dshift,a);CHKERRQ(ierr); 760 } 761 shell->axpy_vscale *= a; 762 if (shell->zrows) { 763 ierr = VecScale(shell->zvals,a);CHKERRQ(ierr); 764 } 765 PetscFunctionReturn(0); 766 } 767 768 static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right) 769 { 770 Mat_Shell *shell = (Mat_Shell*)Y->data; 771 PetscErrorCode ierr; 772 773 PetscFunctionBegin; 774 if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation"); 775 if (left) { 776 if (!shell->left) { 777 ierr = VecDuplicate(left,&shell->left);CHKERRQ(ierr); 778 ierr = VecCopy(left,shell->left);CHKERRQ(ierr); 779 } else { 780 ierr = VecPointwiseMult(shell->left,shell->left,left);CHKERRQ(ierr); 781 } 782 if (shell->zrows) { 783 ierr = VecPointwiseMult(shell->zvals,shell->zvals,left);CHKERRQ(ierr); 784 } 785 } 786 if (right) { 787 if (!shell->right) { 788 ierr = VecDuplicate(right,&shell->right);CHKERRQ(ierr); 789 ierr = VecCopy(right,shell->right);CHKERRQ(ierr); 790 } else { 791 ierr = VecPointwiseMult(shell->right,shell->right,right);CHKERRQ(ierr); 792 } 793 if (shell->zrows) { 794 if (!shell->left_work) { 795 ierr = MatCreateVecs(Y,NULL,&shell->left_work);CHKERRQ(ierr); 796 } 797 ierr = VecSet(shell->zvals_w,1.0);CHKERRQ(ierr); 798 ierr = VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 799 ierr = VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 800 ierr = VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w);CHKERRQ(ierr); 801 } 802 } 803 PetscFunctionReturn(0); 804 } 805 806 PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t) 807 { 808 Mat_Shell *shell = (Mat_Shell*)Y->data; 809 PetscErrorCode ierr; 810 811 PetscFunctionBegin; 812 if (t == MAT_FINAL_ASSEMBLY) { 813 shell->vshift = 0.0; 814 shell->vscale = 1.0; 815 ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr); 816 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 817 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 818 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 819 ierr = VecScatterDestroy(&shell->zvals_sct_c);CHKERRQ(ierr); 820 ierr = VecScatterDestroy(&shell->zvals_sct_r);CHKERRQ(ierr); 821 ierr = ISDestroy(&shell->zrows);CHKERRQ(ierr); 822 ierr = ISDestroy(&shell->zcols);CHKERRQ(ierr); 823 } 824 PetscFunctionReturn(0); 825 } 826 827 static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d) 828 { 829 PetscFunctionBegin; 830 *missing = PETSC_FALSE; 831 PetscFunctionReturn(0); 832 } 833 834 PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str) 835 { 836 Mat_Shell *shell = (Mat_Shell*)Y->data; 837 PetscErrorCode ierr; 838 839 PetscFunctionBegin; 840 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 841 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 842 shell->axpy = X; 843 shell->axpy_vscale = a; 844 PetscFunctionReturn(0); 845 } 846 847 static struct _MatOps MatOps_Values = {0, 848 0, 849 0, 850 0, 851 /* 4*/ MatMultAdd_Shell, 852 0, 853 MatMultTransposeAdd_Shell, 854 0, 855 0, 856 0, 857 /*10*/ 0, 858 0, 859 0, 860 0, 861 0, 862 /*15*/ 0, 863 0, 864 0, 865 MatDiagonalScale_Shell, 866 0, 867 /*20*/ 0, 868 MatAssemblyEnd_Shell, 869 0, 870 0, 871 /*24*/ MatZeroRows_Shell, 872 0, 873 0, 874 0, 875 0, 876 /*29*/ 0, 877 0, 878 0, 879 0, 880 0, 881 /*34*/ MatDuplicate_Shell, 882 0, 883 0, 884 0, 885 0, 886 /*39*/ MatAXPY_Shell, 887 0, 888 0, 889 0, 890 MatCopy_Shell, 891 /*44*/ 0, 892 MatScale_Shell, 893 MatShift_Shell, 894 MatDiagonalSet_Shell, 895 MatZeroRowsColumns_Shell, 896 /*49*/ 0, 897 0, 898 0, 899 0, 900 0, 901 /*54*/ 0, 902 0, 903 0, 904 0, 905 0, 906 /*59*/ 0, 907 MatDestroy_Shell, 908 0, 909 MatConvertFrom_Shell, 910 0, 911 /*64*/ 0, 912 0, 913 0, 914 0, 915 0, 916 /*69*/ 0, 917 0, 918 MatConvert_Shell, 919 0, 920 0, 921 /*74*/ 0, 922 0, 923 0, 924 0, 925 0, 926 /*79*/ 0, 927 0, 928 0, 929 0, 930 0, 931 /*84*/ 0, 932 0, 933 0, 934 0, 935 0, 936 /*89*/ 0, 937 0, 938 0, 939 0, 940 0, 941 /*94*/ 0, 942 0, 943 0, 944 0, 945 0, 946 /*99*/ 0, 947 0, 948 0, 949 0, 950 0, 951 /*104*/ 0, 952 0, 953 0, 954 0, 955 0, 956 /*109*/ 0, 957 0, 958 0, 959 0, 960 MatMissingDiagonal_Shell, 961 /*114*/ 0, 962 0, 963 0, 964 0, 965 0, 966 /*119*/ 0, 967 0, 968 0, 969 0, 970 0, 971 /*124*/ 0, 972 0, 973 0, 974 0, 975 0, 976 /*129*/ 0, 977 0, 978 0, 979 0, 980 0, 981 /*134*/ 0, 982 0, 983 0, 984 0, 985 0, 986 /*139*/ 0, 987 0, 988 0 989 }; 990 991 /*MC 992 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 993 994 Level: advanced 995 996 .seealso: MatCreateShell() 997 M*/ 998 999 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 1000 { 1001 Mat_Shell *b; 1002 PetscErrorCode ierr; 1003 1004 PetscFunctionBegin; 1005 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1006 1007 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 1008 A->data = (void*)b; 1009 1010 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1011 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1012 1013 b->ctx = 0; 1014 b->vshift = 0.0; 1015 b->vscale = 1.0; 1016 b->managescalingshifts = PETSC_TRUE; 1017 A->assembled = PETSC_TRUE; 1018 A->preallocated = PETSC_FALSE; 1019 1020 ierr = PetscObjectChangeTypeName((PetscObject)A,MATSHELL);CHKERRQ(ierr); 1021 PetscFunctionReturn(0); 1022 } 1023 1024 /*@C 1025 MatCreateShell - Creates a new matrix class for use with a user-defined 1026 private data storage format. 1027 1028 Collective on MPI_Comm 1029 1030 Input Parameters: 1031 + comm - MPI communicator 1032 . m - number of local rows (must be given) 1033 . n - number of local columns (must be given) 1034 . M - number of global rows (may be PETSC_DETERMINE) 1035 . N - number of global columns (may be PETSC_DETERMINE) 1036 - ctx - pointer to data needed by the shell matrix routines 1037 1038 Output Parameter: 1039 . A - the matrix 1040 1041 Level: advanced 1042 1043 Usage: 1044 $ extern int mult(Mat,Vec,Vec); 1045 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 1046 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1047 $ [ Use matrix for operations that have been set ] 1048 $ MatDestroy(mat); 1049 1050 Notes: 1051 The shell matrix type is intended to provide a simple class to use 1052 with KSP (such as, for use with matrix-free methods). You should not 1053 use the shell type if you plan to define a complete matrix class. 1054 1055 Fortran Notes: 1056 To use this from Fortran with a ctx you must write an interface definition for this 1057 function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing 1058 in as the ctx argument. 1059 1060 PETSc requires that matrices and vectors being used for certain 1061 operations are partitioned accordingly. For example, when 1062 creating a shell matrix, A, that supports parallel matrix-vector 1063 products using MatMult(A,x,y) the user should set the number 1064 of local matrix rows to be the number of local elements of the 1065 corresponding result vector, y. Note that this is information is 1066 required for use of the matrix interface routines, even though 1067 the shell matrix may not actually be physically partitioned. 1068 For example, 1069 1070 $ 1071 $ Vec x, y 1072 $ extern int mult(Mat,Vec,Vec); 1073 $ Mat A 1074 $ 1075 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 1076 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 1077 $ VecGetLocalSize(y,&m); 1078 $ VecGetLocalSize(x,&n); 1079 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 1080 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1081 $ MatMult(A,x,y); 1082 $ MatDestroy(&A); 1083 $ VecDestroy(&y); 1084 $ VecDestroy(&x); 1085 $ 1086 1087 1088 MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), MatScale(), MatZeroRows() and MatZeroRowsColumns() internally, so these 1089 operations cannot be overwritten unless MatShellSetManageScalingShifts() is called. 1090 1091 1092 For rectangular matrices do all the scalings and shifts make sense? 1093 1094 Developers Notes: 1095 Regarding shifting and scaling. The general form is 1096 1097 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1098 1099 The order you apply the operations is important. For example if you have a dshift then 1100 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift 1101 you get s*vscale*A + diag(shift) 1102 1103 A is the user provided function. 1104 1105 KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for 1106 for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger 1107 an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat); 1108 each time the MATSHELL matrix has changed. 1109 1110 Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided 1111 with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale(). 1112 1113 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts() 1114 @*/ 1115 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 1116 { 1117 PetscErrorCode ierr; 1118 1119 PetscFunctionBegin; 1120 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1121 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 1122 ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr); 1123 ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr); 1124 ierr = MatSetUp(*A);CHKERRQ(ierr); 1125 PetscFunctionReturn(0); 1126 } 1127 1128 /*@ 1129 MatShellSetContext - sets the context for a shell matrix 1130 1131 Logically Collective on Mat 1132 1133 Input Parameters: 1134 + mat - the shell matrix 1135 - ctx - the context 1136 1137 Level: advanced 1138 1139 Fortran Notes: 1140 To use this from Fortran you must write a Fortran interface definition for this 1141 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 1142 1143 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation() 1144 @*/ 1145 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 1146 { 1147 Mat_Shell *shell = (Mat_Shell*)mat->data; 1148 PetscErrorCode ierr; 1149 PetscBool flg; 1150 1151 PetscFunctionBegin; 1152 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1153 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1154 if (flg) { 1155 shell->ctx = ctx; 1156 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix"); 1157 PetscFunctionReturn(0); 1158 } 1159 1160 /*@ 1161 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately 1162 after MatCreateShell() 1163 1164 Logically Collective on Mat 1165 1166 Input Parameter: 1167 . mat - the shell matrix 1168 1169 Level: advanced 1170 1171 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation() 1172 @*/ 1173 PetscErrorCode MatShellSetManageScalingShifts(Mat A) 1174 { 1175 PetscErrorCode ierr; 1176 Mat_Shell *shell; 1177 PetscBool flg; 1178 1179 PetscFunctionBegin; 1180 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1181 ierr = PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);CHKERRQ(ierr); 1182 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1183 shell = (Mat_Shell*)A->data; 1184 shell->managescalingshifts = PETSC_FALSE; 1185 A->ops->diagonalset = NULL; 1186 A->ops->diagonalscale = NULL; 1187 A->ops->scale = NULL; 1188 A->ops->shift = NULL; 1189 A->ops->axpy = NULL; 1190 PetscFunctionReturn(0); 1191 } 1192 1193 /*@C 1194 MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function. 1195 1196 Logically Collective on Mat 1197 1198 Input Parameters: 1199 + mat - the shell matrix 1200 . f - the function 1201 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 1202 - ctx - an optional context for the function 1203 1204 Output Parameter: 1205 . flg - PETSC_TRUE if the multiply is likely correct 1206 1207 Options Database: 1208 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1209 1210 Level: advanced 1211 1212 Fortran Notes: 1213 Not supported from Fortran 1214 1215 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose() 1216 @*/ 1217 PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 1218 { 1219 PetscErrorCode ierr; 1220 PetscInt m,n; 1221 Mat mf,Dmf,Dmat,Ddiff; 1222 PetscReal Diffnorm,Dmfnorm; 1223 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 1224 1225 PetscFunctionBegin; 1226 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1227 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);CHKERRQ(ierr); 1228 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1229 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);CHKERRQ(ierr); 1230 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 1231 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 1232 1233 ierr = MatComputeOperator(mf,MATAIJ,&Dmf);CHKERRQ(ierr); 1234 ierr = MatComputeOperator(mat,MATAIJ,&Dmat);CHKERRQ(ierr); 1235 1236 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 1237 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 1238 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 1239 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 1240 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 1241 flag = PETSC_FALSE; 1242 if (v) { 1243 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n",(double)(Diffnorm/Dmfnorm));CHKERRQ(ierr); 1244 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 1245 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 1246 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 1247 } 1248 } else if (v) { 1249 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 1250 } 1251 if (flg) *flg = flag; 1252 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 1253 ierr = MatDestroy(&mf);CHKERRQ(ierr); 1254 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 1255 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 1256 PetscFunctionReturn(0); 1257 } 1258 1259 /*@C 1260 MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function. 1261 1262 Logically Collective on Mat 1263 1264 Input Parameters: 1265 + mat - the shell matrix 1266 . f - the function 1267 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 1268 - ctx - an optional context for the function 1269 1270 Output Parameter: 1271 . flg - PETSC_TRUE if the multiply is likely correct 1272 1273 Options Database: 1274 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1275 1276 Level: advanced 1277 1278 Fortran Notes: 1279 Not supported from Fortran 1280 1281 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult() 1282 @*/ 1283 PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 1284 { 1285 PetscErrorCode ierr; 1286 Vec x,y,z; 1287 PetscInt m,n,M,N; 1288 Mat mf,Dmf,Dmat,Ddiff; 1289 PetscReal Diffnorm,Dmfnorm; 1290 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 1291 1292 PetscFunctionBegin; 1293 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1294 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);CHKERRQ(ierr); 1295 ierr = MatCreateVecs(mat,&x,&y);CHKERRQ(ierr); 1296 ierr = VecDuplicate(y,&z);CHKERRQ(ierr); 1297 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1298 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1299 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);CHKERRQ(ierr); 1300 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 1301 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 1302 ierr = MatComputeOperator(mf,MATAIJ,&Dmf);CHKERRQ(ierr); 1303 ierr = MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);CHKERRQ(ierr); 1304 ierr = MatComputeOperatorTranspose(mat,MATAIJ,&Dmat);CHKERRQ(ierr); 1305 1306 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 1307 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 1308 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 1309 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 1310 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 1311 flag = PETSC_FALSE; 1312 if (v) { 1313 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n",(double)(Diffnorm/Dmfnorm));CHKERRQ(ierr); 1314 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 1315 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 1316 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 1317 } 1318 } else if (v) { 1319 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 1320 } 1321 if (flg) *flg = flag; 1322 ierr = MatDestroy(&mf);CHKERRQ(ierr); 1323 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 1324 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 1325 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 1326 ierr = VecDestroy(&x);CHKERRQ(ierr); 1327 ierr = VecDestroy(&y);CHKERRQ(ierr); 1328 ierr = VecDestroy(&z);CHKERRQ(ierr); 1329 PetscFunctionReturn(0); 1330 } 1331 1332 /*@C 1333 MatShellSetOperation - Allows user to set a matrix operation for a shell matrix. 1334 1335 Logically Collective on Mat 1336 1337 Input Parameters: 1338 + mat - the shell matrix 1339 . op - the name of the operation 1340 - f - the function that provides the operation. 1341 1342 Level: advanced 1343 1344 Usage: 1345 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 1346 $ ierr = MatCreateShell(comm,m,n,M,N,ctx,&A); 1347 $ ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 1348 1349 Notes: 1350 See the file include/petscmat.h for a complete list of matrix 1351 operations, which all have the form MATOP_<OPERATION>, where 1352 <OPERATION> is the name (in all capital letters) of the 1353 user interface routine (e.g., MatMult() -> MATOP_MULT). 1354 1355 All user-provided functions (except for MATOP_DESTROY) should have the same calling 1356 sequence as the usual matrix interface routines, since they 1357 are intended to be accessed via the usual matrix interface 1358 routines, e.g., 1359 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 1360 1361 In particular each function MUST return an error code of 0 on success and 1362 nonzero on failure. 1363 1364 Within each user-defined routine, the user should call 1365 MatShellGetContext() to obtain the user-defined context that was 1366 set by MatCreateShell(). 1367 1368 Fortran Notes: 1369 For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not 1370 generate a matrix. See src/mat/examples/tests/ex120f.F 1371 1372 Use MatSetOperation() to set an operation for any matrix type 1373 1374 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts() 1375 @*/ 1376 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void)) 1377 { 1378 PetscBool flg; 1379 Mat_Shell *shell; 1380 PetscErrorCode ierr; 1381 1382 PetscFunctionBegin; 1383 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1384 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1385 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1386 shell = (Mat_Shell*)mat->data; 1387 1388 switch (op) { 1389 case MATOP_DESTROY: 1390 shell->ops->destroy = (PetscErrorCode (*)(Mat))f; 1391 break; 1392 case MATOP_VIEW: 1393 if (!mat->ops->viewnative) { 1394 mat->ops->viewnative = mat->ops->view; 1395 } 1396 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 1397 break; 1398 case MATOP_COPY: 1399 shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f; 1400 break; 1401 case MATOP_DIAGONAL_SET: 1402 case MATOP_DIAGONAL_SCALE: 1403 case MATOP_SHIFT: 1404 case MATOP_SCALE: 1405 case MATOP_AXPY: 1406 case MATOP_ZERO_ROWS: 1407 case MATOP_ZERO_ROWS_COLUMNS: 1408 if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1409 (((void(**)(void))mat->ops)[op]) = f; 1410 break; 1411 case MATOP_GET_DIAGONAL: 1412 if (shell->managescalingshifts) { 1413 shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1414 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1415 } else { 1416 shell->ops->getdiagonal = NULL; 1417 mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1418 } 1419 break; 1420 case MATOP_MULT: 1421 if (shell->managescalingshifts) { 1422 shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1423 mat->ops->mult = MatMult_Shell; 1424 } else { 1425 shell->ops->mult = NULL; 1426 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1427 } 1428 break; 1429 case MATOP_MULT_TRANSPOSE: 1430 if (shell->managescalingshifts) { 1431 shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1432 mat->ops->multtranspose = MatMultTranspose_Shell; 1433 } else { 1434 shell->ops->multtranspose = NULL; 1435 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1436 } 1437 break; 1438 default: 1439 (((void(**)(void))mat->ops)[op]) = f; 1440 break; 1441 } 1442 PetscFunctionReturn(0); 1443 } 1444 1445 /*@C 1446 MatShellGetOperation - Gets a matrix function for a shell matrix. 1447 1448 Not Collective 1449 1450 Input Parameters: 1451 + mat - the shell matrix 1452 - op - the name of the operation 1453 1454 Output Parameter: 1455 . f - the function that provides the operation. 1456 1457 Level: advanced 1458 1459 Notes: 1460 See the file include/petscmat.h for a complete list of matrix 1461 operations, which all have the form MATOP_<OPERATION>, where 1462 <OPERATION> is the name (in all capital letters) of the 1463 user interface routine (e.g., MatMult() -> MATOP_MULT). 1464 1465 All user-provided functions have the same calling 1466 sequence as the usual matrix interface routines, since they 1467 are intended to be accessed via the usual matrix interface 1468 routines, e.g., 1469 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 1470 1471 Within each user-defined routine, the user should call 1472 MatShellGetContext() to obtain the user-defined context that was 1473 set by MatCreateShell(). 1474 1475 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext() 1476 @*/ 1477 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void)) 1478 { 1479 PetscBool flg; 1480 Mat_Shell *shell; 1481 PetscErrorCode ierr; 1482 1483 PetscFunctionBegin; 1484 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1485 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1486 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1487 shell = (Mat_Shell*)mat->data; 1488 1489 switch (op) { 1490 case MATOP_DESTROY: 1491 *f = (void (*)(void))shell->ops->destroy; 1492 break; 1493 case MATOP_VIEW: 1494 *f = (void (*)(void))mat->ops->view; 1495 break; 1496 case MATOP_COPY: 1497 *f = (void (*)(void))shell->ops->copy; 1498 break; 1499 case MATOP_DIAGONAL_SET: 1500 case MATOP_DIAGONAL_SCALE: 1501 case MATOP_SHIFT: 1502 case MATOP_SCALE: 1503 case MATOP_AXPY: 1504 case MATOP_ZERO_ROWS: 1505 case MATOP_ZERO_ROWS_COLUMNS: 1506 *f = (((void (**)(void))mat->ops)[op]); 1507 break; 1508 case MATOP_GET_DIAGONAL: 1509 if (shell->ops->getdiagonal) 1510 *f = (void (*)(void))shell->ops->getdiagonal; 1511 else 1512 *f = (((void (**)(void))mat->ops)[op]); 1513 break; 1514 case MATOP_MULT: 1515 if (shell->ops->mult) 1516 *f = (void (*)(void))shell->ops->mult; 1517 else 1518 *f = (((void (**)(void))mat->ops)[op]); 1519 break; 1520 case MATOP_MULT_TRANSPOSE: 1521 if (shell->ops->multtranspose) 1522 *f = (void (*)(void))shell->ops->multtranspose; 1523 else 1524 *f = (((void (**)(void))mat->ops)[op]); 1525 break; 1526 default: 1527 *f = (((void (**)(void))mat->ops)[op]); 1528 } 1529 PetscFunctionReturn(0); 1530 } 1531