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