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