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