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 struct _n_MatShellMatFunctionList { 20 PetscErrorCode (*symbolic)(Mat,Mat,Mat,void**); 21 PetscErrorCode (*numeric)(Mat,Mat,Mat,void*); 22 PetscErrorCode (*destroy)(void*); 23 MatProductType ptype; 24 char *composedname; /* string to identify routine with double dispatch */ 25 char *resultname; /* result matrix type */ 26 27 struct _n_MatShellMatFunctionList *next; 28 }; 29 typedef struct _n_MatShellMatFunctionList *MatShellMatFunctionList; 30 31 typedef struct { 32 struct _MatShellOps ops[1]; 33 34 /* The user will manage the scaling and shifts for the MATSHELL, not the default */ 35 PetscBool managescalingshifts; 36 37 /* support for MatScale, MatShift and MatMultAdd */ 38 PetscScalar vscale,vshift; 39 Vec dshift; 40 Vec left,right; 41 Vec left_work,right_work; 42 Vec left_add_work,right_add_work; 43 44 /* support for MatAXPY */ 45 Mat axpy; 46 PetscScalar axpy_vscale; 47 Vec axpy_left,axpy_right; 48 PetscObjectState axpy_state; 49 50 /* support for ZeroRows/Columns operations */ 51 IS zrows; 52 IS zcols; 53 Vec zvals; 54 Vec zvals_w; 55 VecScatter zvals_sct_r; 56 VecScatter zvals_sct_c; 57 58 /* MatMat operations */ 59 MatShellMatFunctionList matmat; 60 61 /* user defined context */ 62 void *ctx; 63 } Mat_Shell; 64 65 /* 66 Store and scale values on zeroed rows 67 xx = [x_1, 0], 0 on zeroed columns 68 */ 69 static PetscErrorCode MatShellPreZeroRight(Mat A,Vec x,Vec *xx) 70 { 71 Mat_Shell *shell = (Mat_Shell*)A->data; 72 PetscErrorCode ierr; 73 74 PetscFunctionBegin; 75 *xx = x; 76 if (shell->zrows) { 77 ierr = VecSet(shell->zvals_w,0.0);CHKERRQ(ierr); 78 ierr = VecScatterBegin(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 79 ierr = VecScatterEnd(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 80 ierr = VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);CHKERRQ(ierr); 81 } 82 if (shell->zcols) { 83 if (!shell->right_work) { 84 ierr = MatCreateVecs(A,&shell->right_work,NULL);CHKERRQ(ierr); 85 } 86 ierr = VecCopy(x,shell->right_work);CHKERRQ(ierr); 87 ierr = VecISSet(shell->right_work,shell->zcols,0.0);CHKERRQ(ierr); 88 *xx = shell->right_work; 89 } 90 PetscFunctionReturn(0); 91 } 92 93 /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */ 94 static PetscErrorCode MatShellPostZeroLeft(Mat A,Vec x) 95 { 96 Mat_Shell *shell = (Mat_Shell*)A->data; 97 PetscErrorCode ierr; 98 99 PetscFunctionBegin; 100 if (shell->zrows) { 101 ierr = VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 102 ierr = VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 103 } 104 PetscFunctionReturn(0); 105 } 106 107 /* 108 Store and scale values on zeroed rows 109 xx = [x_1, 0], 0 on zeroed rows 110 */ 111 static PetscErrorCode MatShellPreZeroLeft(Mat A,Vec x,Vec *xx) 112 { 113 Mat_Shell *shell = (Mat_Shell*)A->data; 114 PetscErrorCode ierr; 115 116 PetscFunctionBegin; 117 *xx = NULL; 118 if (!shell->zrows) { 119 *xx = x; 120 } else { 121 if (!shell->left_work) { 122 ierr = MatCreateVecs(A,NULL,&shell->left_work);CHKERRQ(ierr); 123 } 124 ierr = VecCopy(x,shell->left_work);CHKERRQ(ierr); 125 ierr = VecSet(shell->zvals_w,0.0);CHKERRQ(ierr); 126 ierr = VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 127 ierr = VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 128 ierr = VecScatterBegin(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 129 ierr = VecScatterEnd(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 130 ierr = VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);CHKERRQ(ierr); 131 *xx = shell->left_work; 132 } 133 PetscFunctionReturn(0); 134 } 135 136 /* Zero zero-columns contributions, sum contributions from properly scaled values stored in MatShellPreZeroLeft */ 137 static PetscErrorCode MatShellPostZeroRight(Mat A,Vec x) 138 { 139 Mat_Shell *shell = (Mat_Shell*)A->data; 140 PetscErrorCode ierr; 141 142 PetscFunctionBegin; 143 if (shell->zcols) { 144 ierr = VecISSet(x,shell->zcols,0.0);CHKERRQ(ierr); 145 } 146 if (shell->zrows) { 147 ierr = VecScatterBegin(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 148 ierr = VecScatterEnd(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 149 } 150 PetscFunctionReturn(0); 151 } 152 153 /* 154 xx = diag(left)*x 155 */ 156 static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx) 157 { 158 Mat_Shell *shell = (Mat_Shell*)A->data; 159 PetscErrorCode ierr; 160 161 PetscFunctionBegin; 162 *xx = NULL; 163 if (!shell->left) { 164 *xx = x; 165 } else { 166 if (!shell->left_work) {ierr = VecDuplicate(shell->left,&shell->left_work);CHKERRQ(ierr);} 167 ierr = VecPointwiseMult(shell->left_work,x,shell->left);CHKERRQ(ierr); 168 *xx = shell->left_work; 169 } 170 PetscFunctionReturn(0); 171 } 172 173 /* 174 xx = diag(right)*x 175 */ 176 static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx) 177 { 178 Mat_Shell *shell = (Mat_Shell*)A->data; 179 PetscErrorCode ierr; 180 181 PetscFunctionBegin; 182 *xx = NULL; 183 if (!shell->right) { 184 *xx = x; 185 } else { 186 if (!shell->right_work) {ierr = VecDuplicate(shell->right,&shell->right_work);CHKERRQ(ierr);} 187 ierr = VecPointwiseMult(shell->right_work,x,shell->right);CHKERRQ(ierr); 188 *xx = shell->right_work; 189 } 190 PetscFunctionReturn(0); 191 } 192 193 /* 194 x = diag(left)*x 195 */ 196 static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x) 197 { 198 Mat_Shell *shell = (Mat_Shell*)A->data; 199 PetscErrorCode ierr; 200 201 PetscFunctionBegin; 202 if (shell->left) {ierr = VecPointwiseMult(x,x,shell->left);CHKERRQ(ierr);} 203 PetscFunctionReturn(0); 204 } 205 206 /* 207 x = diag(right)*x 208 */ 209 static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x) 210 { 211 Mat_Shell *shell = (Mat_Shell*)A->data; 212 PetscErrorCode ierr; 213 214 PetscFunctionBegin; 215 if (shell->right) {ierr = VecPointwiseMult(x,x,shell->right);CHKERRQ(ierr);} 216 PetscFunctionReturn(0); 217 } 218 219 /* 220 Y = vscale*Y + diag(dshift)*X + vshift*X 221 222 On input Y already contains A*x 223 */ 224 static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y) 225 { 226 Mat_Shell *shell = (Mat_Shell*)A->data; 227 PetscErrorCode ierr; 228 229 PetscFunctionBegin; 230 if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */ 231 PetscInt i,m; 232 const PetscScalar *x,*d; 233 PetscScalar *y; 234 ierr = VecGetLocalSize(X,&m);CHKERRQ(ierr); 235 ierr = VecGetArrayRead(shell->dshift,&d);CHKERRQ(ierr); 236 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 237 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 238 for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i]; 239 ierr = VecRestoreArrayRead(shell->dshift,&d);CHKERRQ(ierr); 240 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 241 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 242 } else { 243 ierr = VecScale(Y,shell->vscale);CHKERRQ(ierr); 244 } 245 if (shell->vshift != 0.0) {ierr = VecAXPY(Y,shell->vshift,X);CHKERRQ(ierr);} /* if test is for non-square matrices */ 246 PetscFunctionReturn(0); 247 } 248 249 PetscErrorCode MatShellGetContext_Shell(Mat mat,void *ctx) 250 { 251 PetscFunctionBegin; 252 *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx; 253 PetscFunctionReturn(0); 254 } 255 256 /*@ 257 MatShellGetContext - Returns the user-provided context associated with a shell matrix. 258 259 Not Collective 260 261 Input Parameter: 262 . mat - the matrix, should have been created with MatCreateShell() 263 264 Output Parameter: 265 . ctx - the user provided context 266 267 Level: advanced 268 269 Fortran Notes: 270 To use this from Fortran you must write a Fortran interface definition for this 271 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 272 273 .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext() 274 @*/ 275 PetscErrorCode MatShellGetContext(Mat mat,void *ctx) 276 { 277 PetscErrorCode ierr; 278 279 PetscFunctionBegin; 280 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 281 PetscValidPointer(ctx,2); 282 ierr = PetscUseMethod(mat,"MatShellGetContext_C",(Mat,void*),(mat,ctx));CHKERRQ(ierr); 283 PetscFunctionReturn(0); 284 } 285 286 static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat,PetscInt nr,PetscInt rows[],PetscInt nc,PetscInt cols[],PetscScalar diag,PetscBool rc) 287 { 288 PetscErrorCode ierr; 289 Mat_Shell *shell = (Mat_Shell*)mat->data; 290 Vec x = NULL,b = NULL; 291 IS is1, is2; 292 const PetscInt *ridxs; 293 PetscInt *idxs,*gidxs; 294 PetscInt cum,rst,cst,i; 295 296 PetscFunctionBegin; 297 if (!shell->zvals) { 298 ierr = MatCreateVecs(mat,NULL,&shell->zvals);CHKERRQ(ierr); 299 } 300 if (!shell->zvals_w) { 301 ierr = VecDuplicate(shell->zvals,&shell->zvals_w);CHKERRQ(ierr); 302 } 303 ierr = MatGetOwnershipRange(mat,&rst,NULL);CHKERRQ(ierr); 304 ierr = MatGetOwnershipRangeColumn(mat,&cst,NULL);CHKERRQ(ierr); 305 306 /* Expand/create index set of zeroed rows */ 307 ierr = PetscMalloc1(nr,&idxs);CHKERRQ(ierr); 308 for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst; 309 ierr = ISCreateGeneral(PETSC_COMM_SELF,nr,idxs,PETSC_OWN_POINTER,&is1);CHKERRQ(ierr); 310 ierr = ISSort(is1);CHKERRQ(ierr); 311 ierr = VecISSet(shell->zvals,is1,diag);CHKERRQ(ierr); 312 if (shell->zrows) { 313 ierr = ISSum(shell->zrows,is1,&is2);CHKERRQ(ierr); 314 ierr = ISDestroy(&shell->zrows);CHKERRQ(ierr); 315 ierr = ISDestroy(&is1);CHKERRQ(ierr); 316 shell->zrows = is2; 317 } else shell->zrows = is1; 318 319 /* Create scatters for diagonal values communications */ 320 ierr = VecScatterDestroy(&shell->zvals_sct_c);CHKERRQ(ierr); 321 ierr = VecScatterDestroy(&shell->zvals_sct_r);CHKERRQ(ierr); 322 323 /* row scatter: from/to left vector */ 324 ierr = MatCreateVecs(mat,&x,&b);CHKERRQ(ierr); 325 ierr = VecScatterCreate(b,shell->zrows,shell->zvals_w,shell->zrows,&shell->zvals_sct_r);CHKERRQ(ierr); 326 327 /* col scatter: from right vector to left vector */ 328 ierr = ISGetIndices(shell->zrows,&ridxs);CHKERRQ(ierr); 329 ierr = ISGetLocalSize(shell->zrows,&nr);CHKERRQ(ierr); 330 ierr = PetscMalloc1(nr,&gidxs);CHKERRQ(ierr); 331 for (i = 0, cum = 0; i < nr; i++) { 332 if (ridxs[i] >= mat->cmap->N) continue; 333 gidxs[cum] = ridxs[i]; 334 cum++; 335 } 336 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)mat),cum,gidxs,PETSC_OWN_POINTER,&is1);CHKERRQ(ierr); 337 ierr = VecScatterCreate(x,is1,shell->zvals_w,is1,&shell->zvals_sct_c);CHKERRQ(ierr); 338 ierr = ISDestroy(&is1);CHKERRQ(ierr); 339 ierr = VecDestroy(&x);CHKERRQ(ierr); 340 ierr = VecDestroy(&b);CHKERRQ(ierr); 341 342 /* Expand/create index set of zeroed columns */ 343 if (rc) { 344 ierr = PetscMalloc1(nc,&idxs);CHKERRQ(ierr); 345 for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst; 346 ierr = ISCreateGeneral(PETSC_COMM_SELF,nc,idxs,PETSC_OWN_POINTER,&is1);CHKERRQ(ierr); 347 ierr = ISSort(is1);CHKERRQ(ierr); 348 if (shell->zcols) { 349 ierr = ISSum(shell->zcols,is1,&is2);CHKERRQ(ierr); 350 ierr = ISDestroy(&shell->zcols);CHKERRQ(ierr); 351 ierr = ISDestroy(&is1);CHKERRQ(ierr); 352 shell->zcols = is2; 353 } else shell->zcols = is1; 354 } 355 PetscFunctionReturn(0); 356 } 357 358 static PetscErrorCode MatZeroRows_Shell(Mat mat,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 359 { 360 Mat_Shell *shell = (Mat_Shell*)mat->data; 361 PetscInt nr, *lrows; 362 PetscErrorCode ierr; 363 364 PetscFunctionBegin; 365 if (x && b) { 366 Vec xt; 367 PetscScalar *vals; 368 PetscInt *gcols,i,st,nl,nc; 369 370 ierr = PetscMalloc1(n,&gcols);CHKERRQ(ierr); 371 for (i = 0, nc = 0; i < n; i++) if (rows[i] < mat->cmap->N) gcols[nc++] = rows[i]; 372 373 ierr = MatCreateVecs(mat,&xt,NULL);CHKERRQ(ierr); 374 ierr = VecCopy(x,xt);CHKERRQ(ierr); 375 ierr = PetscCalloc1(nc,&vals);CHKERRQ(ierr); 376 ierr = VecSetValues(xt,nc,gcols,vals,INSERT_VALUES);CHKERRQ(ierr); /* xt = [x1, 0] */ 377 ierr = PetscFree(vals);CHKERRQ(ierr); 378 ierr = VecAssemblyBegin(xt);CHKERRQ(ierr); 379 ierr = VecAssemblyEnd(xt);CHKERRQ(ierr); 380 ierr = VecAYPX(xt,-1.0,x);CHKERRQ(ierr); /* xt = [0, x2] */ 381 382 ierr = VecGetOwnershipRange(xt,&st,NULL);CHKERRQ(ierr); 383 ierr = VecGetLocalSize(xt,&nl);CHKERRQ(ierr); 384 ierr = VecGetArray(xt,&vals);CHKERRQ(ierr); 385 for (i = 0; i < nl; i++) { 386 PetscInt g = i + st; 387 if (g > mat->rmap->N) continue; 388 if (PetscAbsScalar(vals[i]) == 0.0) continue; 389 ierr = VecSetValue(b,g,diag*vals[i],INSERT_VALUES);CHKERRQ(ierr); 390 } 391 ierr = VecRestoreArray(xt,&vals);CHKERRQ(ierr); 392 ierr = VecAssemblyBegin(b);CHKERRQ(ierr); 393 ierr = VecAssemblyEnd(b);CHKERRQ(ierr); /* b = [b1, x2 * diag] */ 394 ierr = VecDestroy(&xt);CHKERRQ(ierr); 395 ierr = PetscFree(gcols);CHKERRQ(ierr); 396 } 397 ierr = PetscLayoutMapLocal(mat->rmap,n,rows,&nr,&lrows,NULL);CHKERRQ(ierr); 398 ierr = MatZeroRowsColumns_Local_Shell(mat,nr,lrows,0,NULL,diag,PETSC_FALSE);CHKERRQ(ierr); 399 if (shell->axpy) { 400 ierr = MatZeroRows(shell->axpy,n,rows,0.0,NULL,NULL);CHKERRQ(ierr); 401 } 402 ierr = PetscFree(lrows);CHKERRQ(ierr); 403 PetscFunctionReturn(0); 404 } 405 406 static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat,PetscInt n,const PetscInt rowscols[],PetscScalar diag,Vec x,Vec b) 407 { 408 Mat_Shell *shell = (Mat_Shell*)mat->data; 409 PetscInt *lrows, *lcols; 410 PetscInt nr, nc; 411 PetscBool congruent; 412 PetscErrorCode ierr; 413 414 PetscFunctionBegin; 415 if (x && b) { 416 Vec xt, bt; 417 PetscScalar *vals; 418 PetscInt *grows,*gcols,i,st,nl; 419 420 ierr = PetscMalloc2(n,&grows,n,&gcols);CHKERRQ(ierr); 421 for (i = 0, nr = 0; i < n; i++) if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i]; 422 for (i = 0, nc = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i]; 423 ierr = PetscCalloc1(n,&vals);CHKERRQ(ierr); 424 425 ierr = MatCreateVecs(mat,&xt,&bt);CHKERRQ(ierr); 426 ierr = VecCopy(x,xt);CHKERRQ(ierr); 427 ierr = VecSetValues(xt,nc,gcols,vals,INSERT_VALUES);CHKERRQ(ierr); /* xt = [x1, 0] */ 428 ierr = VecAssemblyBegin(xt);CHKERRQ(ierr); 429 ierr = VecAssemblyEnd(xt);CHKERRQ(ierr); 430 ierr = VecAXPY(xt,-1.0,x);CHKERRQ(ierr); /* xt = [0, -x2] */ 431 ierr = MatMult(mat,xt,bt);CHKERRQ(ierr); /* bt = [-A12*x2,-A22*x2] */ 432 ierr = VecSetValues(bt,nr,grows,vals,INSERT_VALUES);CHKERRQ(ierr); /* bt = [-A12*x2,0] */ 433 ierr = VecAssemblyBegin(bt);CHKERRQ(ierr); 434 ierr = VecAssemblyEnd(bt);CHKERRQ(ierr); 435 ierr = VecAXPY(b,1.0,bt);CHKERRQ(ierr); /* b = [b1 - A12*x2, b2] */ 436 ierr = VecSetValues(bt,nr,grows,vals,INSERT_VALUES);CHKERRQ(ierr); /* b = [b1 - A12*x2, 0] */ 437 ierr = VecAssemblyBegin(bt);CHKERRQ(ierr); 438 ierr = VecAssemblyEnd(bt);CHKERRQ(ierr); 439 ierr = PetscFree(vals);CHKERRQ(ierr); 440 441 ierr = VecGetOwnershipRange(xt,&st,NULL);CHKERRQ(ierr); 442 ierr = VecGetLocalSize(xt,&nl);CHKERRQ(ierr); 443 ierr = VecGetArray(xt,&vals);CHKERRQ(ierr); 444 for (i = 0; i < nl; i++) { 445 PetscInt g = i + st; 446 if (g > mat->rmap->N) continue; 447 if (PetscAbsScalar(vals[i]) == 0.0) continue; 448 ierr = VecSetValue(b,g,-diag*vals[i],INSERT_VALUES);CHKERRQ(ierr); 449 } 450 ierr = VecRestoreArray(xt,&vals);CHKERRQ(ierr); 451 ierr = VecAssemblyBegin(b);CHKERRQ(ierr); 452 ierr = VecAssemblyEnd(b);CHKERRQ(ierr); /* b = [b1 - A12*x2, x2 * diag] */ 453 ierr = VecDestroy(&xt);CHKERRQ(ierr); 454 ierr = VecDestroy(&bt);CHKERRQ(ierr); 455 ierr = PetscFree2(grows,gcols);CHKERRQ(ierr); 456 } 457 ierr = PetscLayoutMapLocal(mat->rmap,n,rowscols,&nr,&lrows,NULL);CHKERRQ(ierr); 458 ierr = MatHasCongruentLayouts(mat,&congruent);CHKERRQ(ierr); 459 if (congruent) { 460 nc = nr; 461 lcols = lrows; 462 } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */ 463 PetscInt i,nt,*t; 464 465 ierr = PetscMalloc1(n,&t);CHKERRQ(ierr); 466 for (i = 0, nt = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i]; 467 ierr = PetscLayoutMapLocal(mat->cmap,nt,t,&nc,&lcols,NULL);CHKERRQ(ierr); 468 ierr = PetscFree(t);CHKERRQ(ierr); 469 } 470 ierr = MatZeroRowsColumns_Local_Shell(mat,nr,lrows,nc,lcols,diag,PETSC_TRUE);CHKERRQ(ierr); 471 if (!congruent) { 472 ierr = PetscFree(lcols);CHKERRQ(ierr); 473 } 474 ierr = PetscFree(lrows);CHKERRQ(ierr); 475 if (shell->axpy) { 476 ierr = MatZeroRowsColumns(shell->axpy,n,rowscols,0.0,NULL,NULL);CHKERRQ(ierr); 477 } 478 PetscFunctionReturn(0); 479 } 480 481 PetscErrorCode MatDestroy_Shell(Mat mat) 482 { 483 PetscErrorCode ierr; 484 Mat_Shell *shell = (Mat_Shell*)mat->data; 485 MatShellMatFunctionList matmat; 486 487 PetscFunctionBegin; 488 if (shell->ops->destroy) { 489 ierr = (*shell->ops->destroy)(mat);CHKERRQ(ierr); 490 } 491 ierr = PetscMemzero(shell->ops,sizeof(struct _MatShellOps));CHKERRQ(ierr); 492 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 493 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 494 ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr); 495 ierr = VecDestroy(&shell->left_work);CHKERRQ(ierr); 496 ierr = VecDestroy(&shell->right_work);CHKERRQ(ierr); 497 ierr = VecDestroy(&shell->left_add_work);CHKERRQ(ierr); 498 ierr = VecDestroy(&shell->right_add_work);CHKERRQ(ierr); 499 ierr = VecDestroy(&shell->axpy_left);CHKERRQ(ierr); 500 ierr = VecDestroy(&shell->axpy_right);CHKERRQ(ierr); 501 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 502 ierr = VecDestroy(&shell->zvals_w);CHKERRQ(ierr); 503 ierr = VecDestroy(&shell->zvals);CHKERRQ(ierr); 504 ierr = VecScatterDestroy(&shell->zvals_sct_c);CHKERRQ(ierr); 505 ierr = VecScatterDestroy(&shell->zvals_sct_r);CHKERRQ(ierr); 506 ierr = ISDestroy(&shell->zrows);CHKERRQ(ierr); 507 ierr = ISDestroy(&shell->zcols);CHKERRQ(ierr); 508 509 matmat = shell->matmat; 510 while (matmat) { 511 MatShellMatFunctionList next = matmat->next; 512 513 ierr = PetscObjectComposeFunction((PetscObject)mat,matmat->composedname,NULL);CHKERRQ(ierr); 514 ierr = PetscFree(matmat->composedname);CHKERRQ(ierr); 515 ierr = PetscFree(matmat->resultname);CHKERRQ(ierr); 516 ierr = PetscFree(matmat);CHKERRQ(ierr); 517 matmat = next; 518 } 519 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatShellGetContext_C",NULL);CHKERRQ(ierr); 520 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatShellSetContext_C",NULL);CHKERRQ(ierr); 521 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatShellSetVecType_C",NULL);CHKERRQ(ierr); 522 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatShellSetManageScalingShifts_C",NULL);CHKERRQ(ierr); 523 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatShellSetOperation_C",NULL);CHKERRQ(ierr); 524 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatShellGetOperation_C",NULL);CHKERRQ(ierr); 525 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatShellSetMatProductOperation_C",NULL);CHKERRQ(ierr); 526 ierr = PetscFree(mat->data);CHKERRQ(ierr); 527 PetscFunctionReturn(0); 528 } 529 530 typedef struct { 531 PetscErrorCode (*numeric)(Mat,Mat,Mat,void*); 532 PetscErrorCode (*destroy)(void*); 533 void *userdata; 534 Mat B; 535 Mat Bt; 536 Mat axpy; 537 } MatMatDataShell; 538 539 static PetscErrorCode DestroyMatMatDataShell(void *data) 540 { 541 MatMatDataShell *mmdata = (MatMatDataShell *)data; 542 PetscErrorCode ierr; 543 544 PetscFunctionBegin; 545 if (mmdata->destroy) { 546 ierr = (*mmdata->destroy)(mmdata->userdata);CHKERRQ(ierr); 547 } 548 ierr = MatDestroy(&mmdata->B);CHKERRQ(ierr); 549 ierr = MatDestroy(&mmdata->Bt);CHKERRQ(ierr); 550 ierr = MatDestroy(&mmdata->axpy);CHKERRQ(ierr); 551 ierr = PetscFree(mmdata);CHKERRQ(ierr); 552 PetscFunctionReturn(0); 553 } 554 555 static PetscErrorCode MatProductNumeric_Shell_X(Mat D) 556 { 557 PetscErrorCode ierr; 558 Mat_Product *product; 559 Mat A, B; 560 MatMatDataShell *mdata; 561 PetscScalar zero = 0.0; 562 563 PetscFunctionBegin; 564 MatCheckProduct(D,1); 565 product = D->product; 566 if (!product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data empty"); 567 A = product->A; 568 B = product->B; 569 mdata = (MatMatDataShell*)product->data; 570 if (mdata->numeric) { 571 Mat_Shell *shell = (Mat_Shell*)A->data; 572 PetscErrorCode (*stashsym)(Mat) = D->ops->productsymbolic; 573 PetscErrorCode (*stashnum)(Mat) = D->ops->productnumeric; 574 PetscBool useBmdata = PETSC_FALSE, newB = PETSC_TRUE; 575 576 if (shell->managescalingshifts) { 577 if (shell->zcols || shell->zrows) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProduct not supported with zeroed rows/columns"); 578 if (shell->right || shell->left) { 579 useBmdata = PETSC_TRUE; 580 if (!mdata->B) { 581 ierr = MatDuplicate(B,MAT_SHARE_NONZERO_PATTERN,&mdata->B);CHKERRQ(ierr); 582 } else { 583 newB = PETSC_FALSE; 584 } 585 ierr = MatCopy(B,mdata->B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 586 } 587 switch (product->type) { 588 case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */ 589 if (shell->right) { 590 ierr = MatDiagonalScale(mdata->B,shell->right,NULL);CHKERRQ(ierr); 591 } 592 break; 593 case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */ 594 if (shell->left) { 595 ierr = MatDiagonalScale(mdata->B,shell->left,NULL);CHKERRQ(ierr); 596 } 597 break; 598 case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */ 599 if (shell->right) { 600 ierr = MatDiagonalScale(mdata->B,NULL,shell->right);CHKERRQ(ierr); 601 } 602 break; 603 case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */ 604 if (shell->right && shell->left) { 605 PetscBool flg; 606 607 ierr = VecEqual(shell->right,shell->left,&flg);CHKERRQ(ierr); 608 if (!flg) SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name); 609 } 610 if (shell->right) { 611 ierr = MatDiagonalScale(mdata->B,NULL,shell->right);CHKERRQ(ierr); 612 } 613 break; 614 case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */ 615 if (shell->right && shell->left) { 616 PetscBool flg; 617 618 ierr = VecEqual(shell->right,shell->left,&flg);CHKERRQ(ierr); 619 if (!flg) SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name); 620 } 621 if (shell->right) { 622 ierr = MatDiagonalScale(mdata->B,shell->right,NULL);CHKERRQ(ierr); 623 } 624 break; 625 default: SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name); 626 } 627 } 628 /* allow the user to call MatMat operations on D */ 629 D->product = NULL; 630 D->ops->productsymbolic = NULL; 631 D->ops->productnumeric = NULL; 632 633 ierr = (*mdata->numeric)(A,useBmdata ? mdata->B : B,D,mdata->userdata);CHKERRQ(ierr); 634 635 /* clear any leftover user data and restore D pointers */ 636 ierr = MatProductClear(D);CHKERRQ(ierr); 637 D->ops->productsymbolic = stashsym; 638 D->ops->productnumeric = stashnum; 639 D->product = product; 640 641 if (shell->managescalingshifts) { 642 ierr = MatScale(D,shell->vscale);CHKERRQ(ierr); 643 switch (product->type) { 644 case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */ 645 case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */ 646 if (shell->left) { 647 ierr = MatDiagonalScale(D,shell->left,NULL);CHKERRQ(ierr); 648 if (shell->dshift || shell->vshift != zero) { 649 if (!shell->left_work) {ierr = MatCreateVecs(A,NULL,&shell->left_work);CHKERRQ(ierr);} 650 if (shell->dshift) { 651 ierr = VecCopy(shell->dshift,shell->left_work);CHKERRQ(ierr); 652 ierr = VecShift(shell->left_work,shell->vshift);CHKERRQ(ierr); 653 ierr = VecPointwiseMult(shell->left_work,shell->left_work,shell->left);CHKERRQ(ierr); 654 } else { 655 ierr = VecSet(shell->left_work,shell->vshift);CHKERRQ(ierr); 656 } 657 if (product->type == MATPRODUCT_ABt) { 658 MatReuse reuse = mdata->Bt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX; 659 MatStructure str = mdata->Bt ? SUBSET_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN; 660 661 ierr = MatTranspose(mdata->B,reuse,&mdata->Bt);CHKERRQ(ierr); 662 ierr = MatDiagonalScale(mdata->Bt,shell->left_work,NULL);CHKERRQ(ierr); 663 ierr = MatAXPY(D,1.0,mdata->Bt,str);CHKERRQ(ierr); 664 } else { 665 MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN; 666 667 ierr = MatDiagonalScale(mdata->B,shell->left_work,NULL);CHKERRQ(ierr); 668 ierr = MatAXPY(D,1.0,mdata->B,str);CHKERRQ(ierr); 669 } 670 } 671 } 672 break; 673 case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */ 674 if (shell->right) { 675 ierr = MatDiagonalScale(D,shell->right,NULL);CHKERRQ(ierr); 676 if (shell->dshift || shell->vshift != zero) { 677 MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN; 678 679 if (!shell->right_work) {ierr = MatCreateVecs(A,&shell->right_work,NULL);CHKERRQ(ierr);} 680 if (shell->dshift) { 681 ierr = VecCopy(shell->dshift,shell->right_work);CHKERRQ(ierr); 682 ierr = VecShift(shell->right_work,shell->vshift);CHKERRQ(ierr); 683 ierr = VecPointwiseMult(shell->right_work,shell->right_work,shell->right);CHKERRQ(ierr); 684 } else { 685 ierr = VecSet(shell->right_work,shell->vshift);CHKERRQ(ierr); 686 } 687 ierr = MatDiagonalScale(mdata->B,shell->right_work,NULL);CHKERRQ(ierr); 688 ierr = MatAXPY(D,1.0,mdata->B,str);CHKERRQ(ierr); 689 } 690 } 691 break; 692 case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */ 693 case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */ 694 if (shell->dshift || shell->vshift != zero) SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices with diagonal shift",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name); 695 break; 696 default: SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name); 697 } 698 if (shell->axpy && shell->axpy_vscale != zero) { 699 Mat X; 700 PetscObjectState axpy_state; 701 MatStructure str = DIFFERENT_NONZERO_PATTERN; /* not sure it is safe to ever use SUBSET_NONZERO_PATTERN */ 702 703 ierr = MatShellGetContext(shell->axpy,&X);CHKERRQ(ierr); 704 ierr = PetscObjectStateGet((PetscObject)X,&axpy_state);CHKERRQ(ierr); 705 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,...)"); 706 if (!mdata->axpy) { 707 str = DIFFERENT_NONZERO_PATTERN; 708 ierr = MatProductCreate(shell->axpy,B,NULL,&mdata->axpy);CHKERRQ(ierr); 709 ierr = MatProductSetType(mdata->axpy,product->type);CHKERRQ(ierr); 710 ierr = MatProductSetFromOptions(mdata->axpy);CHKERRQ(ierr); 711 ierr = MatProductSymbolic(mdata->axpy);CHKERRQ(ierr); 712 } else { /* May be that shell->axpy has changed */ 713 PetscBool flg; 714 715 ierr = MatProductReplaceMats(shell->axpy,B,NULL,mdata->axpy);CHKERRQ(ierr); 716 ierr = MatHasOperation(mdata->axpy,MATOP_PRODUCTSYMBOLIC,&flg);CHKERRQ(ierr); 717 if (!flg) { 718 str = DIFFERENT_NONZERO_PATTERN; 719 ierr = MatProductSetFromOptions(mdata->axpy);CHKERRQ(ierr); 720 ierr = MatProductSymbolic(mdata->axpy);CHKERRQ(ierr); 721 } 722 } 723 ierr = MatProductNumeric(mdata->axpy);CHKERRQ(ierr); 724 ierr = MatAXPY(D,shell->axpy_vscale,mdata->axpy,str);CHKERRQ(ierr); 725 } 726 } 727 } else SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation"); 728 PetscFunctionReturn(0); 729 } 730 731 static PetscErrorCode MatProductSymbolic_Shell_X(Mat D) 732 { 733 PetscErrorCode ierr; 734 Mat_Product *product; 735 Mat A,B; 736 MatShellMatFunctionList matmat; 737 Mat_Shell *shell; 738 PetscBool flg; 739 char composedname[256]; 740 MatMatDataShell *mdata; 741 742 PetscFunctionBegin; 743 MatCheckProduct(D,1); 744 product = D->product; 745 if (product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data not empty"); 746 A = product->A; 747 B = product->B; 748 shell = (Mat_Shell*)A->data; 749 matmat = shell->matmat; 750 ierr = PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,((PetscObject)B)->type_name);CHKERRQ(ierr); 751 while (matmat) { 752 ierr = PetscStrcmp(composedname,matmat->composedname,&flg);CHKERRQ(ierr); 753 flg = (PetscBool)(flg && (matmat->ptype == product->type)); 754 if (flg) break; 755 matmat = matmat->next; 756 } 757 if (!flg) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Composedname \"%s\" for product type %s not found",composedname,MatProductTypes[product->type]); 758 switch (product->type) { 759 case MATPRODUCT_AB: 760 ierr = MatSetSizes(D,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N);CHKERRQ(ierr); 761 break; 762 case MATPRODUCT_AtB: 763 ierr = MatSetSizes(D,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N);CHKERRQ(ierr); 764 break; 765 case MATPRODUCT_ABt: 766 ierr = MatSetSizes(D,A->rmap->n,B->rmap->n,A->rmap->N,B->rmap->N);CHKERRQ(ierr); 767 break; 768 case MATPRODUCT_RARt: 769 ierr = MatSetSizes(D,B->rmap->n,B->rmap->n,B->rmap->N,B->rmap->N);CHKERRQ(ierr); 770 break; 771 case MATPRODUCT_PtAP: 772 ierr = MatSetSizes(D,B->cmap->n,B->cmap->n,B->cmap->N,B->cmap->N);CHKERRQ(ierr); 773 break; 774 default: SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name); 775 } 776 /* respect users who passed in a matrix for which resultname is the base type */ 777 if (matmat->resultname) { 778 ierr = PetscObjectBaseTypeCompare((PetscObject)D,matmat->resultname,&flg);CHKERRQ(ierr); 779 if (!flg) { 780 ierr = MatSetType(D,matmat->resultname);CHKERRQ(ierr); 781 } 782 } 783 /* If matrix type was not set or different, we need to reset this pointers */ 784 D->ops->productsymbolic = MatProductSymbolic_Shell_X; 785 D->ops->productnumeric = MatProductNumeric_Shell_X; 786 /* attach product data */ 787 ierr = PetscNew(&mdata);CHKERRQ(ierr); 788 mdata->numeric = matmat->numeric; 789 mdata->destroy = matmat->destroy; 790 if (matmat->symbolic) { 791 ierr = (*matmat->symbolic)(A,B,D,&mdata->userdata);CHKERRQ(ierr); 792 } else { /* call general setup if symbolic operation not provided */ 793 ierr = MatSetUp(D);CHKERRQ(ierr); 794 } 795 if (!D->product) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_COR,"Product disappeared after user symbolic phase"); 796 if (D->product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_COR,"Product data not empty after user symbolic phase"); 797 D->product->data = mdata; 798 D->product->destroy = DestroyMatMatDataShell; 799 /* Be sure to reset these pointers if the user did something unexpected */ 800 D->ops->productsymbolic = MatProductSymbolic_Shell_X; 801 D->ops->productnumeric = MatProductNumeric_Shell_X; 802 PetscFunctionReturn(0); 803 } 804 805 static PetscErrorCode MatProductSetFromOptions_Shell_X(Mat D) 806 { 807 PetscErrorCode ierr; 808 Mat_Product *product; 809 Mat A,B; 810 MatShellMatFunctionList matmat; 811 Mat_Shell *shell; 812 PetscBool flg; 813 char composedname[256]; 814 815 PetscFunctionBegin; 816 MatCheckProduct(D,1); 817 product = D->product; 818 A = product->A; 819 B = product->B; 820 ierr = MatIsShell(A,&flg);CHKERRQ(ierr); 821 if (!flg) PetscFunctionReturn(0); 822 shell = (Mat_Shell*)A->data; 823 matmat = shell->matmat; 824 ierr = PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,((PetscObject)B)->type_name);CHKERRQ(ierr); 825 while (matmat) { 826 ierr = PetscStrcmp(composedname,matmat->composedname,&flg);CHKERRQ(ierr); 827 flg = (PetscBool)(flg && (matmat->ptype == product->type)); 828 if (flg) break; 829 matmat = matmat->next; 830 } 831 if (flg) { D->ops->productsymbolic = MatProductSymbolic_Shell_X; } 832 else { ierr = PetscInfo2(D," symbolic product %s not registered for product type %s\n",composedname,MatProductTypes[product->type]);CHKERRQ(ierr); } 833 PetscFunctionReturn(0); 834 } 835 836 static PetscErrorCode MatShellSetMatProductOperation_Private(Mat A,MatProductType ptype,PetscErrorCode (*symbolic)(Mat,Mat,Mat,void**),PetscErrorCode (*numeric)(Mat,Mat,Mat,void*),PetscErrorCode (*destroy)(void*),char *composedname,const char *resultname) 837 { 838 PetscBool flg; 839 PetscErrorCode ierr; 840 Mat_Shell *shell; 841 MatShellMatFunctionList matmat; 842 843 PetscFunctionBegin; 844 if (!numeric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing numeric routine"); 845 if (!composedname) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing composed name"); 846 847 /* add product callback */ 848 shell = (Mat_Shell*)A->data; 849 matmat = shell->matmat; 850 if (!matmat) { 851 ierr = PetscNew(&shell->matmat);CHKERRQ(ierr); 852 matmat = shell->matmat; 853 } else { 854 MatShellMatFunctionList entry = matmat; 855 while (entry) { 856 ierr = PetscStrcmp(composedname,entry->composedname,&flg);CHKERRQ(ierr); 857 flg = (PetscBool)(flg && (entry->ptype == ptype)); 858 if (flg) goto set; 859 matmat = entry; 860 entry = entry->next; 861 } 862 ierr = PetscNew(&matmat->next);CHKERRQ(ierr); 863 matmat = matmat->next; 864 } 865 866 set: 867 matmat->symbolic = symbolic; 868 matmat->numeric = numeric; 869 matmat->destroy = destroy; 870 matmat->ptype = ptype; 871 ierr = PetscFree(matmat->composedname);CHKERRQ(ierr); 872 ierr = PetscFree(matmat->resultname);CHKERRQ(ierr); 873 ierr = PetscStrallocpy(composedname,&matmat->composedname);CHKERRQ(ierr); 874 ierr = PetscStrallocpy(resultname,&matmat->resultname);CHKERRQ(ierr); 875 ierr = PetscInfo3(A,"Composing %s for product type %s with result %s\n",matmat->composedname,MatProductTypes[matmat->ptype],matmat->resultname ? matmat->resultname : "not specified");CHKERRQ(ierr); 876 ierr = PetscObjectComposeFunction((PetscObject)A,matmat->composedname,MatProductSetFromOptions_Shell_X);CHKERRQ(ierr); 877 PetscFunctionReturn(0); 878 } 879 880 /*@C 881 MatShellSetMatProductOperation - Allows user to set a matrix matrix operation for a shell matrix. 882 883 Logically Collective on Mat 884 885 Input Parameters: 886 + A - the shell matrix 887 . ptype - the product type 888 . symbolic - the function for the symbolic phase (can be NULL) 889 . numeric - the function for the numerical phase 890 . destroy - the function for the destruction of the needed data generated during the symbolic phase (can be NULL) 891 . Btype - the matrix type for the matrix to be multiplied against 892 - Ctype - the matrix type for the result (can be NULL) 893 894 Level: advanced 895 896 Usage: 897 $ extern PetscErrorCode usersymbolic(Mat,Mat,Mat,void**); 898 $ extern PetscErrorCode usernumeric(Mat,Mat,Mat,void*); 899 $ extern PetscErrorCode userdestroy(void*); 900 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 901 $ MatShellSetMatProductOperation(A,MATPRODUCT_AB,usersymbolic,usernumeric,userdestroy,MATSEQAIJ,MATDENSE); 902 $ [ create B of type SEQAIJ etc..] 903 $ MatProductCreate(A,B,NULL,&C); 904 $ MatProductSetType(C,MATPRODUCT_AB); 905 $ MatProductSetFromOptions(C); 906 $ MatProductSymbolic(C); -> actually runs the user defined symbolic operation 907 $ MatProductNumeric(C); -> actually runs the user defined numeric operation 908 $ [ use C = A*B ] 909 910 Notes: 911 MATPRODUCT_ABC is not supported yet. Not supported in Fortran. 912 If the symbolic phase is not specified, MatSetUp() is called on the result matrix that must have its type set if Ctype is NULL. 913 Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback. 914 PETSc will take care of calling the user-defined callbacks. 915 It is allowed to specify the same callbacks for different Btype matrix types. 916 The couple (Btype,ptype) uniquely identifies the operation: the last specified callbacks takes precedence. 917 918 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatProductType, MatType, MatSetUp() 919 @*/ 920 PetscErrorCode MatShellSetMatProductOperation(Mat A,MatProductType ptype,PetscErrorCode (*symbolic)(Mat,Mat,Mat,void**),PetscErrorCode (*numeric)(Mat,Mat,Mat,void*),PetscErrorCode (*destroy)(void *),MatType Btype,MatType Ctype) 921 { 922 PetscErrorCode ierr; 923 924 PetscFunctionBegin; 925 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 926 PetscValidLogicalCollectiveEnum(A,ptype,2); 927 if (ptype == MATPRODUCT_ABC) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not for product type %s",MatProductTypes[ptype]); 928 if (!numeric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Missing numeric routine, argument 4"); 929 PetscValidPointer(Btype,6); 930 if (Ctype) PetscValidPointer(Ctype,7); 931 ierr = PetscTryMethod(A,"MatShellSetMatProductOperation_C",(Mat,MatProductType,PetscErrorCode(*)(Mat,Mat,Mat,void**),PetscErrorCode(*)(Mat,Mat,Mat,void*),PetscErrorCode(*)(void*),MatType,MatType),(A,ptype,symbolic,numeric,destroy,Btype,Ctype));CHKERRQ(ierr); 932 PetscFunctionReturn(0); 933 } 934 935 PetscErrorCode MatShellSetMatProductOperation_Shell(Mat A,MatProductType ptype,PetscErrorCode (*symbolic)(Mat,Mat,Mat,void**),PetscErrorCode (*numeric)(Mat,Mat,Mat,void*),PetscErrorCode (*destroy)(void *),MatType Btype,MatType Ctype) 936 { 937 PetscBool flg; 938 PetscErrorCode ierr; 939 char composedname[256]; 940 MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList; 941 PetscMPIInt size; 942 943 PetscFunctionBegin; 944 PetscValidType(A,1); 945 while (Bnames) { /* user passed in the root name */ 946 ierr = PetscStrcmp(Btype,Bnames->rname,&flg);CHKERRQ(ierr); 947 if (flg) break; 948 Bnames = Bnames->next; 949 } 950 while (Cnames) { /* user passed in the root name */ 951 ierr = PetscStrcmp(Ctype,Cnames->rname,&flg);CHKERRQ(ierr); 952 if (flg) break; 953 Cnames = Cnames->next; 954 } 955 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRMPI(ierr); 956 Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype; 957 Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype; 958 ierr = PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,Btype);CHKERRQ(ierr); 959 ierr = MatShellSetMatProductOperation_Private(A,ptype,symbolic,numeric,destroy,composedname,Ctype);CHKERRQ(ierr); 960 PetscFunctionReturn(0); 961 } 962 963 PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str) 964 { 965 Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data; 966 PetscErrorCode ierr; 967 PetscBool matflg; 968 MatShellMatFunctionList matmatA; 969 970 PetscFunctionBegin; 971 ierr = MatIsShell(B,&matflg);CHKERRQ(ierr); 972 if (!matflg) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix %s not derived from MATSHELL",((PetscObject)B)->type_name); 973 974 ierr = PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 975 ierr = PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));CHKERRQ(ierr); 976 977 if (shellA->ops->copy) { 978 ierr = (*shellA->ops->copy)(A,B,str);CHKERRQ(ierr); 979 } 980 shellB->vscale = shellA->vscale; 981 shellB->vshift = shellA->vshift; 982 if (shellA->dshift) { 983 if (!shellB->dshift) { 984 ierr = VecDuplicate(shellA->dshift,&shellB->dshift);CHKERRQ(ierr); 985 } 986 ierr = VecCopy(shellA->dshift,shellB->dshift);CHKERRQ(ierr); 987 } else { 988 ierr = VecDestroy(&shellB->dshift);CHKERRQ(ierr); 989 } 990 if (shellA->left) { 991 if (!shellB->left) { 992 ierr = VecDuplicate(shellA->left,&shellB->left);CHKERRQ(ierr); 993 } 994 ierr = VecCopy(shellA->left,shellB->left);CHKERRQ(ierr); 995 } else { 996 ierr = VecDestroy(&shellB->left);CHKERRQ(ierr); 997 } 998 if (shellA->right) { 999 if (!shellB->right) { 1000 ierr = VecDuplicate(shellA->right,&shellB->right);CHKERRQ(ierr); 1001 } 1002 ierr = VecCopy(shellA->right,shellB->right);CHKERRQ(ierr); 1003 } else { 1004 ierr = VecDestroy(&shellB->right);CHKERRQ(ierr); 1005 } 1006 ierr = MatDestroy(&shellB->axpy);CHKERRQ(ierr); 1007 shellB->axpy_vscale = 0.0; 1008 shellB->axpy_state = 0; 1009 if (shellA->axpy) { 1010 ierr = PetscObjectReference((PetscObject)shellA->axpy);CHKERRQ(ierr); 1011 shellB->axpy = shellA->axpy; 1012 shellB->axpy_vscale = shellA->axpy_vscale; 1013 shellB->axpy_state = shellA->axpy_state; 1014 } 1015 if (shellA->zrows) { 1016 ierr = ISDuplicate(shellA->zrows,&shellB->zrows);CHKERRQ(ierr); 1017 if (shellA->zcols) { 1018 ierr = ISDuplicate(shellA->zcols,&shellB->zcols);CHKERRQ(ierr); 1019 } 1020 ierr = VecDuplicate(shellA->zvals,&shellB->zvals);CHKERRQ(ierr); 1021 ierr = VecCopy(shellA->zvals,shellB->zvals);CHKERRQ(ierr); 1022 ierr = VecDuplicate(shellA->zvals_w,&shellB->zvals_w);CHKERRQ(ierr); 1023 ierr = PetscObjectReference((PetscObject)shellA->zvals_sct_r);CHKERRQ(ierr); 1024 ierr = PetscObjectReference((PetscObject)shellA->zvals_sct_c);CHKERRQ(ierr); 1025 shellB->zvals_sct_r = shellA->zvals_sct_r; 1026 shellB->zvals_sct_c = shellA->zvals_sct_c; 1027 } 1028 1029 matmatA = shellA->matmat; 1030 if (matmatA) { 1031 while (matmatA->next) { 1032 ierr = MatShellSetMatProductOperation_Private(B,matmatA->ptype,matmatA->symbolic,matmatA->numeric,matmatA->destroy,matmatA->composedname,matmatA->resultname);CHKERRQ(ierr); 1033 matmatA = matmatA->next; 1034 } 1035 } 1036 PetscFunctionReturn(0); 1037 } 1038 1039 PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M) 1040 { 1041 PetscErrorCode ierr; 1042 void *ctx; 1043 1044 PetscFunctionBegin; 1045 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr); 1046 ierr = MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);CHKERRQ(ierr); 1047 ierr = PetscObjectChangeTypeName((PetscObject)(*M),((PetscObject)mat)->type_name);CHKERRQ(ierr); 1048 if (op != MAT_DO_NOT_COPY_VALUES) { 1049 ierr = MatCopy(mat,*M,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 1050 } 1051 PetscFunctionReturn(0); 1052 } 1053 1054 PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y) 1055 { 1056 Mat_Shell *shell = (Mat_Shell*)A->data; 1057 PetscErrorCode ierr; 1058 Vec xx; 1059 PetscObjectState instate,outstate; 1060 1061 PetscFunctionBegin; 1062 if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL"); 1063 ierr = MatShellPreZeroRight(A,x,&xx);CHKERRQ(ierr); 1064 ierr = MatShellPreScaleRight(A,xx,&xx);CHKERRQ(ierr); 1065 ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr); 1066 ierr = (*shell->ops->mult)(A,xx,y);CHKERRQ(ierr); 1067 ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr); 1068 if (instate == outstate) { 1069 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 1070 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 1071 } 1072 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 1073 ierr = MatShellPostScaleLeft(A,y);CHKERRQ(ierr); 1074 ierr = MatShellPostZeroLeft(A,y);CHKERRQ(ierr); 1075 1076 if (shell->axpy) { 1077 Mat X; 1078 PetscObjectState axpy_state; 1079 1080 ierr = MatShellGetContext(shell->axpy,&X);CHKERRQ(ierr); 1081 ierr = PetscObjectStateGet((PetscObject)X,&axpy_state);CHKERRQ(ierr); 1082 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,...)"); 1083 1084 ierr = MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);CHKERRQ(ierr); 1085 ierr = VecCopy(x,shell->axpy_right);CHKERRQ(ierr); 1086 ierr = MatMult(shell->axpy,shell->axpy_right,shell->axpy_left);CHKERRQ(ierr); 1087 ierr = VecAXPY(y,shell->axpy_vscale,shell->axpy_left);CHKERRQ(ierr); 1088 } 1089 PetscFunctionReturn(0); 1090 } 1091 1092 PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z) 1093 { 1094 Mat_Shell *shell = (Mat_Shell*)A->data; 1095 PetscErrorCode ierr; 1096 1097 PetscFunctionBegin; 1098 if (y == z) { 1099 if (!shell->right_add_work) {ierr = VecDuplicate(z,&shell->right_add_work);CHKERRQ(ierr);} 1100 ierr = MatMult(A,x,shell->right_add_work);CHKERRQ(ierr); 1101 ierr = VecAXPY(z,1.0,shell->right_add_work);CHKERRQ(ierr); 1102 } else { 1103 ierr = MatMult(A,x,z);CHKERRQ(ierr); 1104 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 1105 } 1106 PetscFunctionReturn(0); 1107 } 1108 1109 PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y) 1110 { 1111 Mat_Shell *shell = (Mat_Shell*)A->data; 1112 PetscErrorCode ierr; 1113 Vec xx; 1114 PetscObjectState instate,outstate; 1115 1116 PetscFunctionBegin; 1117 if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL"); 1118 ierr = MatShellPreZeroLeft(A,x,&xx);CHKERRQ(ierr); 1119 ierr = MatShellPreScaleLeft(A,xx,&xx);CHKERRQ(ierr); 1120 ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr); 1121 ierr = (*shell->ops->multtranspose)(A,xx,y);CHKERRQ(ierr); 1122 ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr); 1123 if (instate == outstate) { 1124 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 1125 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 1126 } 1127 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 1128 ierr = MatShellPostScaleRight(A,y);CHKERRQ(ierr); 1129 ierr = MatShellPostZeroRight(A,y);CHKERRQ(ierr); 1130 1131 if (shell->axpy) { 1132 Mat X; 1133 PetscObjectState axpy_state; 1134 1135 ierr = MatShellGetContext(shell->axpy,&X);CHKERRQ(ierr); 1136 ierr = PetscObjectStateGet((PetscObject)X,&axpy_state);CHKERRQ(ierr); 1137 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,...)"); 1138 ierr = MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);CHKERRQ(ierr); 1139 ierr = VecCopy(x,shell->axpy_left);CHKERRQ(ierr); 1140 ierr = MatMultTranspose(shell->axpy,shell->axpy_left,shell->axpy_right);CHKERRQ(ierr); 1141 ierr = VecAXPY(y,shell->axpy_vscale,shell->axpy_right);CHKERRQ(ierr); 1142 } 1143 PetscFunctionReturn(0); 1144 } 1145 1146 PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z) 1147 { 1148 Mat_Shell *shell = (Mat_Shell*)A->data; 1149 PetscErrorCode ierr; 1150 1151 PetscFunctionBegin; 1152 if (y == z) { 1153 if (!shell->left_add_work) {ierr = VecDuplicate(z,&shell->left_add_work);CHKERRQ(ierr);} 1154 ierr = MatMultTranspose(A,x,shell->left_add_work);CHKERRQ(ierr); 1155 ierr = VecAXPY(z,1.0,shell->left_add_work);CHKERRQ(ierr); 1156 } else { 1157 ierr = MatMultTranspose(A,x,z);CHKERRQ(ierr); 1158 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 1159 } 1160 PetscFunctionReturn(0); 1161 } 1162 1163 /* 1164 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1165 */ 1166 PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v) 1167 { 1168 Mat_Shell *shell = (Mat_Shell*)A->data; 1169 PetscErrorCode ierr; 1170 1171 PetscFunctionBegin; 1172 if (shell->ops->getdiagonal) { 1173 ierr = (*shell->ops->getdiagonal)(A,v);CHKERRQ(ierr); 1174 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)"); 1175 ierr = VecScale(v,shell->vscale);CHKERRQ(ierr); 1176 if (shell->dshift) { 1177 ierr = VecAXPY(v,1.0,shell->dshift);CHKERRQ(ierr); 1178 } 1179 ierr = VecShift(v,shell->vshift);CHKERRQ(ierr); 1180 if (shell->left) {ierr = VecPointwiseMult(v,v,shell->left);CHKERRQ(ierr);} 1181 if (shell->right) {ierr = VecPointwiseMult(v,v,shell->right);CHKERRQ(ierr);} 1182 if (shell->zrows) { 1183 ierr = VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1184 ierr = VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1185 } 1186 if (shell->axpy) { 1187 Mat X; 1188 PetscObjectState axpy_state; 1189 1190 ierr = MatShellGetContext(shell->axpy,&X);CHKERRQ(ierr); 1191 ierr = PetscObjectStateGet((PetscObject)X,&axpy_state);CHKERRQ(ierr); 1192 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,...)"); 1193 ierr = MatCreateVecs(shell->axpy,NULL,shell->axpy_left ? NULL : &shell->axpy_left);CHKERRQ(ierr); 1194 ierr = MatGetDiagonal(shell->axpy,shell->axpy_left);CHKERRQ(ierr); 1195 ierr = VecAXPY(v,shell->axpy_vscale,shell->axpy_left);CHKERRQ(ierr); 1196 } 1197 PetscFunctionReturn(0); 1198 } 1199 1200 PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a) 1201 { 1202 Mat_Shell *shell = (Mat_Shell*)Y->data; 1203 PetscErrorCode ierr; 1204 PetscBool flg; 1205 1206 PetscFunctionBegin; 1207 ierr = MatHasCongruentLayouts(Y,&flg);CHKERRQ(ierr); 1208 if (!flg) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot shift shell matrix if it is not congruent"); 1209 if (shell->left || shell->right) { 1210 if (!shell->dshift) { 1211 ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);CHKERRQ(ierr); 1212 ierr = VecSet(shell->dshift,a);CHKERRQ(ierr); 1213 } else { 1214 if (shell->left) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 1215 if (shell->right) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 1216 ierr = VecShift(shell->dshift,a);CHKERRQ(ierr); 1217 } 1218 if (shell->left) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 1219 if (shell->right) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 1220 } else shell->vshift += a; 1221 if (shell->zrows) { 1222 ierr = VecShift(shell->zvals,a);CHKERRQ(ierr); 1223 } 1224 PetscFunctionReturn(0); 1225 } 1226 1227 PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s) 1228 { 1229 Mat_Shell *shell = (Mat_Shell*)A->data; 1230 PetscErrorCode ierr; 1231 1232 PetscFunctionBegin; 1233 if (!shell->dshift) {ierr = VecDuplicate(D,&shell->dshift);CHKERRQ(ierr);} 1234 if (shell->left || shell->right) { 1235 if (!shell->right_work) {ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);CHKERRQ(ierr);} 1236 if (shell->left && shell->right) { 1237 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 1238 ierr = VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);CHKERRQ(ierr); 1239 } else if (shell->left) { 1240 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 1241 } else { 1242 ierr = VecPointwiseDivide(shell->right_work,D,shell->right);CHKERRQ(ierr); 1243 } 1244 ierr = VecAXPY(shell->dshift,s,shell->right_work);CHKERRQ(ierr); 1245 } else { 1246 ierr = VecAXPY(shell->dshift,s,D);CHKERRQ(ierr); 1247 } 1248 PetscFunctionReturn(0); 1249 } 1250 1251 PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins) 1252 { 1253 Mat_Shell *shell = (Mat_Shell*)A->data; 1254 Vec d; 1255 PetscErrorCode ierr; 1256 PetscBool flg; 1257 1258 PetscFunctionBegin; 1259 ierr = MatHasCongruentLayouts(A,&flg);CHKERRQ(ierr); 1260 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot diagonal set or shift shell matrix if it is not congruent"); 1261 if (ins == INSERT_VALUES) { 1262 if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first"); 1263 ierr = VecDuplicate(D,&d);CHKERRQ(ierr); 1264 ierr = MatGetDiagonal(A,d);CHKERRQ(ierr); 1265 ierr = MatDiagonalSet_Shell_Private(A,d,-1.);CHKERRQ(ierr); 1266 ierr = MatDiagonalSet_Shell_Private(A,D,1.);CHKERRQ(ierr); 1267 ierr = VecDestroy(&d);CHKERRQ(ierr); 1268 if (shell->zrows) { 1269 ierr = VecCopy(D,shell->zvals);CHKERRQ(ierr); 1270 } 1271 } else { 1272 ierr = MatDiagonalSet_Shell_Private(A,D,1.);CHKERRQ(ierr); 1273 if (shell->zrows) { 1274 ierr = VecAXPY(shell->zvals,1.0,D);CHKERRQ(ierr); 1275 } 1276 } 1277 PetscFunctionReturn(0); 1278 } 1279 1280 PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a) 1281 { 1282 Mat_Shell *shell = (Mat_Shell*)Y->data; 1283 PetscErrorCode ierr; 1284 1285 PetscFunctionBegin; 1286 shell->vscale *= a; 1287 shell->vshift *= a; 1288 if (shell->dshift) { 1289 ierr = VecScale(shell->dshift,a);CHKERRQ(ierr); 1290 } 1291 shell->axpy_vscale *= a; 1292 if (shell->zrows) { 1293 ierr = VecScale(shell->zvals,a);CHKERRQ(ierr); 1294 } 1295 PetscFunctionReturn(0); 1296 } 1297 1298 static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right) 1299 { 1300 Mat_Shell *shell = (Mat_Shell*)Y->data; 1301 PetscErrorCode ierr; 1302 1303 PetscFunctionBegin; 1304 if (left) { 1305 if (!shell->left) { 1306 ierr = VecDuplicate(left,&shell->left);CHKERRQ(ierr); 1307 ierr = VecCopy(left,shell->left);CHKERRQ(ierr); 1308 } else { 1309 ierr = VecPointwiseMult(shell->left,shell->left,left);CHKERRQ(ierr); 1310 } 1311 if (shell->zrows) { 1312 ierr = VecPointwiseMult(shell->zvals,shell->zvals,left);CHKERRQ(ierr); 1313 } 1314 } 1315 if (right) { 1316 if (!shell->right) { 1317 ierr = VecDuplicate(right,&shell->right);CHKERRQ(ierr); 1318 ierr = VecCopy(right,shell->right);CHKERRQ(ierr); 1319 } else { 1320 ierr = VecPointwiseMult(shell->right,shell->right,right);CHKERRQ(ierr); 1321 } 1322 if (shell->zrows) { 1323 if (!shell->left_work) { 1324 ierr = MatCreateVecs(Y,NULL,&shell->left_work);CHKERRQ(ierr); 1325 } 1326 ierr = VecSet(shell->zvals_w,1.0);CHKERRQ(ierr); 1327 ierr = VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1328 ierr = VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1329 ierr = VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w);CHKERRQ(ierr); 1330 } 1331 } 1332 if (shell->axpy) { 1333 ierr = MatDiagonalScale(shell->axpy,left,right);CHKERRQ(ierr); 1334 } 1335 PetscFunctionReturn(0); 1336 } 1337 1338 PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t) 1339 { 1340 Mat_Shell *shell = (Mat_Shell*)Y->data; 1341 PetscErrorCode ierr; 1342 1343 PetscFunctionBegin; 1344 if (t == MAT_FINAL_ASSEMBLY) { 1345 shell->vshift = 0.0; 1346 shell->vscale = 1.0; 1347 shell->axpy_vscale = 0.0; 1348 shell->axpy_state = 0; 1349 ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr); 1350 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 1351 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 1352 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 1353 ierr = VecDestroy(&shell->axpy_left);CHKERRQ(ierr); 1354 ierr = VecDestroy(&shell->axpy_right);CHKERRQ(ierr); 1355 ierr = VecScatterDestroy(&shell->zvals_sct_c);CHKERRQ(ierr); 1356 ierr = VecScatterDestroy(&shell->zvals_sct_r);CHKERRQ(ierr); 1357 ierr = ISDestroy(&shell->zrows);CHKERRQ(ierr); 1358 ierr = ISDestroy(&shell->zcols);CHKERRQ(ierr); 1359 } 1360 PetscFunctionReturn(0); 1361 } 1362 1363 static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d) 1364 { 1365 PetscFunctionBegin; 1366 *missing = PETSC_FALSE; 1367 PetscFunctionReturn(0); 1368 } 1369 1370 PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str) 1371 { 1372 Mat_Shell *shell = (Mat_Shell*)Y->data; 1373 PetscErrorCode ierr; 1374 1375 PetscFunctionBegin; 1376 if (X == Y) { 1377 ierr = MatScale(Y,1.0 + a);CHKERRQ(ierr); 1378 PetscFunctionReturn(0); 1379 } 1380 if (!shell->axpy) { 1381 ierr = MatConvertFrom_Shell(X,MATSHELL,MAT_INITIAL_MATRIX,&shell->axpy);CHKERRQ(ierr); 1382 shell->axpy_vscale = a; 1383 ierr = PetscObjectStateGet((PetscObject)X,&shell->axpy_state);CHKERRQ(ierr); 1384 } else { 1385 ierr = MatAXPY(shell->axpy,a/shell->axpy_vscale,X,str);CHKERRQ(ierr); 1386 } 1387 PetscFunctionReturn(0); 1388 } 1389 1390 static struct _MatOps MatOps_Values = {NULL, 1391 NULL, 1392 NULL, 1393 NULL, 1394 /* 4*/ MatMultAdd_Shell, 1395 NULL, 1396 MatMultTransposeAdd_Shell, 1397 NULL, 1398 NULL, 1399 NULL, 1400 /*10*/ NULL, 1401 NULL, 1402 NULL, 1403 NULL, 1404 NULL, 1405 /*15*/ NULL, 1406 NULL, 1407 NULL, 1408 MatDiagonalScale_Shell, 1409 NULL, 1410 /*20*/ NULL, 1411 MatAssemblyEnd_Shell, 1412 NULL, 1413 NULL, 1414 /*24*/ MatZeroRows_Shell, 1415 NULL, 1416 NULL, 1417 NULL, 1418 NULL, 1419 /*29*/ NULL, 1420 NULL, 1421 NULL, 1422 NULL, 1423 NULL, 1424 /*34*/ MatDuplicate_Shell, 1425 NULL, 1426 NULL, 1427 NULL, 1428 NULL, 1429 /*39*/ MatAXPY_Shell, 1430 NULL, 1431 NULL, 1432 NULL, 1433 MatCopy_Shell, 1434 /*44*/ NULL, 1435 MatScale_Shell, 1436 MatShift_Shell, 1437 MatDiagonalSet_Shell, 1438 MatZeroRowsColumns_Shell, 1439 /*49*/ NULL, 1440 NULL, 1441 NULL, 1442 NULL, 1443 NULL, 1444 /*54*/ NULL, 1445 NULL, 1446 NULL, 1447 NULL, 1448 NULL, 1449 /*59*/ NULL, 1450 MatDestroy_Shell, 1451 NULL, 1452 MatConvertFrom_Shell, 1453 NULL, 1454 /*64*/ NULL, 1455 NULL, 1456 NULL, 1457 NULL, 1458 NULL, 1459 /*69*/ NULL, 1460 NULL, 1461 MatConvert_Shell, 1462 NULL, 1463 NULL, 1464 /*74*/ NULL, 1465 NULL, 1466 NULL, 1467 NULL, 1468 NULL, 1469 /*79*/ NULL, 1470 NULL, 1471 NULL, 1472 NULL, 1473 NULL, 1474 /*84*/ NULL, 1475 NULL, 1476 NULL, 1477 NULL, 1478 NULL, 1479 /*89*/ NULL, 1480 NULL, 1481 NULL, 1482 NULL, 1483 NULL, 1484 /*94*/ NULL, 1485 NULL, 1486 NULL, 1487 NULL, 1488 NULL, 1489 /*99*/ NULL, 1490 NULL, 1491 NULL, 1492 NULL, 1493 NULL, 1494 /*104*/ NULL, 1495 NULL, 1496 NULL, 1497 NULL, 1498 NULL, 1499 /*109*/ NULL, 1500 NULL, 1501 NULL, 1502 NULL, 1503 MatMissingDiagonal_Shell, 1504 /*114*/ NULL, 1505 NULL, 1506 NULL, 1507 NULL, 1508 NULL, 1509 /*119*/ NULL, 1510 NULL, 1511 NULL, 1512 NULL, 1513 NULL, 1514 /*124*/ NULL, 1515 NULL, 1516 NULL, 1517 NULL, 1518 NULL, 1519 /*129*/ NULL, 1520 NULL, 1521 NULL, 1522 NULL, 1523 NULL, 1524 /*134*/ NULL, 1525 NULL, 1526 NULL, 1527 NULL, 1528 NULL, 1529 /*139*/ NULL, 1530 NULL, 1531 NULL 1532 }; 1533 1534 PetscErrorCode MatShellSetContext_Shell(Mat mat,void *ctx) 1535 { 1536 Mat_Shell *shell = (Mat_Shell*)mat->data; 1537 1538 PetscFunctionBegin; 1539 shell->ctx = ctx; 1540 PetscFunctionReturn(0); 1541 } 1542 1543 static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype) 1544 { 1545 PetscErrorCode ierr; 1546 1547 PetscFunctionBegin; 1548 ierr = PetscFree(mat->defaultvectype);CHKERRQ(ierr); 1549 ierr = PetscStrallocpy(vtype,(char**)&mat->defaultvectype);CHKERRQ(ierr); 1550 PetscFunctionReturn(0); 1551 } 1552 1553 PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A) 1554 { 1555 Mat_Shell *shell = (Mat_Shell*)A->data; 1556 1557 PetscFunctionBegin; 1558 shell->managescalingshifts = PETSC_FALSE; 1559 A->ops->diagonalset = NULL; 1560 A->ops->diagonalscale = NULL; 1561 A->ops->scale = NULL; 1562 A->ops->shift = NULL; 1563 A->ops->axpy = NULL; 1564 PetscFunctionReturn(0); 1565 } 1566 1567 PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void)) 1568 { 1569 Mat_Shell *shell = (Mat_Shell*)mat->data; 1570 1571 PetscFunctionBegin; 1572 switch (op) { 1573 case MATOP_DESTROY: 1574 shell->ops->destroy = (PetscErrorCode (*)(Mat))f; 1575 break; 1576 case MATOP_VIEW: 1577 if (!mat->ops->viewnative) { 1578 mat->ops->viewnative = mat->ops->view; 1579 } 1580 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 1581 break; 1582 case MATOP_COPY: 1583 shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f; 1584 break; 1585 case MATOP_DIAGONAL_SET: 1586 case MATOP_DIAGONAL_SCALE: 1587 case MATOP_SHIFT: 1588 case MATOP_SCALE: 1589 case MATOP_AXPY: 1590 case MATOP_ZERO_ROWS: 1591 case MATOP_ZERO_ROWS_COLUMNS: 1592 if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1593 (((void(**)(void))mat->ops)[op]) = f; 1594 break; 1595 case MATOP_GET_DIAGONAL: 1596 if (shell->managescalingshifts) { 1597 shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1598 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1599 } else { 1600 shell->ops->getdiagonal = NULL; 1601 mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1602 } 1603 break; 1604 case MATOP_MULT: 1605 if (shell->managescalingshifts) { 1606 shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1607 mat->ops->mult = MatMult_Shell; 1608 } else { 1609 shell->ops->mult = NULL; 1610 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1611 } 1612 break; 1613 case MATOP_MULT_TRANSPOSE: 1614 if (shell->managescalingshifts) { 1615 shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1616 mat->ops->multtranspose = MatMultTranspose_Shell; 1617 } else { 1618 shell->ops->multtranspose = NULL; 1619 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1620 } 1621 break; 1622 default: 1623 (((void(**)(void))mat->ops)[op]) = f; 1624 break; 1625 } 1626 PetscFunctionReturn(0); 1627 } 1628 1629 PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void)) 1630 { 1631 Mat_Shell *shell = (Mat_Shell*)mat->data; 1632 1633 PetscFunctionBegin; 1634 switch (op) { 1635 case MATOP_DESTROY: 1636 *f = (void (*)(void))shell->ops->destroy; 1637 break; 1638 case MATOP_VIEW: 1639 *f = (void (*)(void))mat->ops->view; 1640 break; 1641 case MATOP_COPY: 1642 *f = (void (*)(void))shell->ops->copy; 1643 break; 1644 case MATOP_DIAGONAL_SET: 1645 case MATOP_DIAGONAL_SCALE: 1646 case MATOP_SHIFT: 1647 case MATOP_SCALE: 1648 case MATOP_AXPY: 1649 case MATOP_ZERO_ROWS: 1650 case MATOP_ZERO_ROWS_COLUMNS: 1651 *f = (((void (**)(void))mat->ops)[op]); 1652 break; 1653 case MATOP_GET_DIAGONAL: 1654 if (shell->ops->getdiagonal) 1655 *f = (void (*)(void))shell->ops->getdiagonal; 1656 else 1657 *f = (((void (**)(void))mat->ops)[op]); 1658 break; 1659 case MATOP_MULT: 1660 if (shell->ops->mult) 1661 *f = (void (*)(void))shell->ops->mult; 1662 else 1663 *f = (((void (**)(void))mat->ops)[op]); 1664 break; 1665 case MATOP_MULT_TRANSPOSE: 1666 if (shell->ops->multtranspose) 1667 *f = (void (*)(void))shell->ops->multtranspose; 1668 else 1669 *f = (((void (**)(void))mat->ops)[op]); 1670 break; 1671 default: 1672 *f = (((void (**)(void))mat->ops)[op]); 1673 } 1674 PetscFunctionReturn(0); 1675 } 1676 1677 /*MC 1678 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 1679 1680 Level: advanced 1681 1682 .seealso: MatCreateShell() 1683 M*/ 1684 1685 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 1686 { 1687 Mat_Shell *b; 1688 PetscErrorCode ierr; 1689 1690 PetscFunctionBegin; 1691 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1692 1693 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 1694 A->data = (void*)b; 1695 1696 b->ctx = NULL; 1697 b->vshift = 0.0; 1698 b->vscale = 1.0; 1699 b->managescalingshifts = PETSC_TRUE; 1700 A->assembled = PETSC_TRUE; 1701 A->preallocated = PETSC_FALSE; 1702 1703 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell);CHKERRQ(ierr); 1704 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell);CHKERRQ(ierr); 1705 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellSetVecType_C",MatShellSetVecType_Shell);CHKERRQ(ierr); 1706 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell);CHKERRQ(ierr); 1707 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell);CHKERRQ(ierr); 1708 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell);CHKERRQ(ierr); 1709 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellSetMatProductOperation_C",MatShellSetMatProductOperation_Shell);CHKERRQ(ierr); 1710 ierr = PetscObjectChangeTypeName((PetscObject)A,MATSHELL);CHKERRQ(ierr); 1711 PetscFunctionReturn(0); 1712 } 1713 1714 /*@C 1715 MatCreateShell - Creates a new matrix class for use with a user-defined 1716 private data storage format. 1717 1718 Collective 1719 1720 Input Parameters: 1721 + comm - MPI communicator 1722 . m - number of local rows (must be given) 1723 . n - number of local columns (must be given) 1724 . M - number of global rows (may be PETSC_DETERMINE) 1725 . N - number of global columns (may be PETSC_DETERMINE) 1726 - ctx - pointer to data needed by the shell matrix routines 1727 1728 Output Parameter: 1729 . A - the matrix 1730 1731 Level: advanced 1732 1733 Usage: 1734 $ extern PetscErrorCode mult(Mat,Vec,Vec); 1735 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 1736 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1737 $ [ Use matrix for operations that have been set ] 1738 $ MatDestroy(mat); 1739 1740 Notes: 1741 The shell matrix type is intended to provide a simple class to use 1742 with KSP (such as, for use with matrix-free methods). You should not 1743 use the shell type if you plan to define a complete matrix class. 1744 1745 Fortran Notes: 1746 To use this from Fortran with a ctx you must write an interface definition for this 1747 function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing 1748 in as the ctx argument. 1749 1750 PETSc requires that matrices and vectors being used for certain 1751 operations are partitioned accordingly. For example, when 1752 creating a shell matrix, A, that supports parallel matrix-vector 1753 products using MatMult(A,x,y) the user should set the number 1754 of local matrix rows to be the number of local elements of the 1755 corresponding result vector, y. Note that this is information is 1756 required for use of the matrix interface routines, even though 1757 the shell matrix may not actually be physically partitioned. 1758 For example, 1759 1760 $ 1761 $ Vec x, y 1762 $ extern PetscErrorCode mult(Mat,Vec,Vec); 1763 $ Mat A 1764 $ 1765 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 1766 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 1767 $ VecGetLocalSize(y,&m); 1768 $ VecGetLocalSize(x,&n); 1769 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 1770 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1771 $ MatMult(A,x,y); 1772 $ MatDestroy(&A); 1773 $ VecDestroy(&y); 1774 $ VecDestroy(&x); 1775 $ 1776 1777 MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), MatScale(), MatZeroRows() and MatZeroRowsColumns() internally, so these 1778 operations cannot be overwritten unless MatShellSetManageScalingShifts() is called. 1779 1780 For rectangular matrices do all the scalings and shifts make sense? 1781 1782 Developers Notes: 1783 Regarding shifting and scaling. The general form is 1784 1785 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1786 1787 The order you apply the operations is important. For example if you have a dshift then 1788 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift 1789 you get s*vscale*A + diag(shift) 1790 1791 A is the user provided function. 1792 1793 KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for 1794 for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger 1795 an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat); 1796 each time the MATSHELL matrix has changed. 1797 1798 Matrix product operations (i.e. MatMat, MatTransposeMat etc) can be specified using MatShellSetMatProductOperation() 1799 1800 Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided 1801 with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale(). 1802 1803 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts(), MatShellSetMatProductOperation() 1804 @*/ 1805 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 1806 { 1807 PetscErrorCode ierr; 1808 1809 PetscFunctionBegin; 1810 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1811 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 1812 ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr); 1813 ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr); 1814 ierr = MatSetUp(*A);CHKERRQ(ierr); 1815 PetscFunctionReturn(0); 1816 } 1817 1818 /*@ 1819 MatShellSetContext - sets the context for a shell matrix 1820 1821 Logically Collective on Mat 1822 1823 Input Parameters: 1824 + mat - the shell matrix 1825 - ctx - the context 1826 1827 Level: advanced 1828 1829 Fortran Notes: 1830 To use this from Fortran you must write a Fortran interface definition for this 1831 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 1832 1833 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation() 1834 @*/ 1835 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 1836 { 1837 PetscErrorCode ierr; 1838 1839 PetscFunctionBegin; 1840 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1841 ierr = PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));CHKERRQ(ierr); 1842 PetscFunctionReturn(0); 1843 } 1844 1845 /*@C 1846 MatShellSetVecType - Sets the type of Vec returned by MatCreateVecs() 1847 1848 Logically collective 1849 1850 Input Parameters: 1851 + mat - the shell matrix 1852 - vtype - type to use for creating vectors 1853 1854 Notes: 1855 1856 Level: advanced 1857 1858 .seealso: MatCreateVecs() 1859 @*/ 1860 PetscErrorCode MatShellSetVecType(Mat mat,VecType vtype) 1861 { 1862 PetscErrorCode ierr; 1863 1864 PetscFunctionBegin; 1865 ierr = PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype));CHKERRQ(ierr); 1866 PetscFunctionReturn(0); 1867 } 1868 1869 /*@ 1870 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately 1871 after MatCreateShell() 1872 1873 Logically Collective on Mat 1874 1875 Input Parameter: 1876 . mat - the shell matrix 1877 1878 Level: advanced 1879 1880 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation() 1881 @*/ 1882 PetscErrorCode MatShellSetManageScalingShifts(Mat A) 1883 { 1884 PetscErrorCode ierr; 1885 1886 PetscFunctionBegin; 1887 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1888 ierr = PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));CHKERRQ(ierr); 1889 PetscFunctionReturn(0); 1890 } 1891 1892 /*@C 1893 MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function. 1894 1895 Logically Collective on Mat 1896 1897 Input Parameters: 1898 + mat - the shell matrix 1899 . f - the function 1900 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 1901 - ctx - an optional context for the function 1902 1903 Output Parameter: 1904 . flg - PETSC_TRUE if the multiply is likely correct 1905 1906 Options Database: 1907 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1908 1909 Level: advanced 1910 1911 Fortran Notes: 1912 Not supported from Fortran 1913 1914 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose() 1915 @*/ 1916 PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 1917 { 1918 PetscErrorCode ierr; 1919 PetscInt m,n; 1920 Mat mf,Dmf,Dmat,Ddiff; 1921 PetscReal Diffnorm,Dmfnorm; 1922 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 1923 1924 PetscFunctionBegin; 1925 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1926 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);CHKERRQ(ierr); 1927 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1928 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);CHKERRQ(ierr); 1929 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 1930 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 1931 1932 ierr = MatComputeOperator(mf,MATAIJ,&Dmf);CHKERRQ(ierr); 1933 ierr = MatComputeOperator(mat,MATAIJ,&Dmat);CHKERRQ(ierr); 1934 1935 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 1936 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 1937 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 1938 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 1939 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 1940 flag = PETSC_FALSE; 1941 if (v) { 1942 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); 1943 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 1944 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 1945 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 1946 } 1947 } else if (v) { 1948 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 1949 } 1950 if (flg) *flg = flag; 1951 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 1952 ierr = MatDestroy(&mf);CHKERRQ(ierr); 1953 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 1954 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 1955 PetscFunctionReturn(0); 1956 } 1957 1958 /*@C 1959 MatShellTestMultTranspose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function. 1960 1961 Logically Collective on Mat 1962 1963 Input Parameters: 1964 + mat - the shell matrix 1965 . f - the function 1966 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 1967 - ctx - an optional context for the function 1968 1969 Output Parameter: 1970 . flg - PETSC_TRUE if the multiply is likely correct 1971 1972 Options Database: 1973 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1974 1975 Level: advanced 1976 1977 Fortran Notes: 1978 Not supported from Fortran 1979 1980 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult() 1981 @*/ 1982 PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 1983 { 1984 PetscErrorCode ierr; 1985 Vec x,y,z; 1986 PetscInt m,n,M,N; 1987 Mat mf,Dmf,Dmat,Ddiff; 1988 PetscReal Diffnorm,Dmfnorm; 1989 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 1990 1991 PetscFunctionBegin; 1992 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1993 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);CHKERRQ(ierr); 1994 ierr = MatCreateVecs(mat,&x,&y);CHKERRQ(ierr); 1995 ierr = VecDuplicate(y,&z);CHKERRQ(ierr); 1996 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1997 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1998 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);CHKERRQ(ierr); 1999 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 2000 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 2001 ierr = MatComputeOperator(mf,MATAIJ,&Dmf);CHKERRQ(ierr); 2002 ierr = MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);CHKERRQ(ierr); 2003 ierr = MatComputeOperatorTranspose(mat,MATAIJ,&Dmat);CHKERRQ(ierr); 2004 2005 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 2006 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 2007 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 2008 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 2009 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 2010 flag = PETSC_FALSE; 2011 if (v) { 2012 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); 2013 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 2014 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 2015 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 2016 } 2017 } else if (v) { 2018 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 2019 } 2020 if (flg) *flg = flag; 2021 ierr = MatDestroy(&mf);CHKERRQ(ierr); 2022 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 2023 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 2024 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 2025 ierr = VecDestroy(&x);CHKERRQ(ierr); 2026 ierr = VecDestroy(&y);CHKERRQ(ierr); 2027 ierr = VecDestroy(&z);CHKERRQ(ierr); 2028 PetscFunctionReturn(0); 2029 } 2030 2031 /*@C 2032 MatShellSetOperation - Allows user to set a matrix operation for a shell matrix. 2033 2034 Logically Collective on Mat 2035 2036 Input Parameters: 2037 + mat - the shell matrix 2038 . op - the name of the operation 2039 - g - the function that provides the operation. 2040 2041 Level: advanced 2042 2043 Usage: 2044 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 2045 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 2046 $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 2047 2048 Notes: 2049 See the file include/petscmat.h for a complete list of matrix 2050 operations, which all have the form MATOP_<OPERATION>, where 2051 <OPERATION> is the name (in all capital letters) of the 2052 user interface routine (e.g., MatMult() -> MATOP_MULT). 2053 2054 All user-provided functions (except for MATOP_DESTROY) should have the same calling 2055 sequence as the usual matrix interface routines, since they 2056 are intended to be accessed via the usual matrix interface 2057 routines, e.g., 2058 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 2059 2060 In particular each function MUST return an error code of 0 on success and 2061 nonzero on failure. 2062 2063 Within each user-defined routine, the user should call 2064 MatShellGetContext() to obtain the user-defined context that was 2065 set by MatCreateShell(). 2066 2067 Use MatSetOperation() to set an operation for any matrix type. For matrix product operations (i.e. MatMat, MatTransposeMat etc) use MatShellSetMatProductOperation() 2068 2069 Fortran Notes: 2070 For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not 2071 generate a matrix. See src/mat/tests/ex120f.F 2072 2073 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts(), MatShellSetMatProductOperation() 2074 @*/ 2075 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void)) 2076 { 2077 PetscErrorCode ierr; 2078 2079 PetscFunctionBegin; 2080 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2081 ierr = PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));CHKERRQ(ierr); 2082 PetscFunctionReturn(0); 2083 } 2084 2085 /*@C 2086 MatShellGetOperation - Gets a matrix function for a shell matrix. 2087 2088 Not Collective 2089 2090 Input Parameters: 2091 + mat - the shell matrix 2092 - op - the name of the operation 2093 2094 Output Parameter: 2095 . g - the function that provides the operation. 2096 2097 Level: advanced 2098 2099 Notes: 2100 See the file include/petscmat.h for a complete list of matrix 2101 operations, which all have the form MATOP_<OPERATION>, where 2102 <OPERATION> is the name (in all capital letters) of the 2103 user interface routine (e.g., MatMult() -> MATOP_MULT). 2104 2105 All user-provided functions have the same calling 2106 sequence as the usual matrix interface routines, since they 2107 are intended to be accessed via the usual matrix interface 2108 routines, e.g., 2109 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 2110 2111 Within each user-defined routine, the user should call 2112 MatShellGetContext() to obtain the user-defined context that was 2113 set by MatCreateShell(). 2114 2115 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext() 2116 @*/ 2117 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void)) 2118 { 2119 PetscErrorCode ierr; 2120 2121 PetscFunctionBegin; 2122 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2123 ierr = PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g));CHKERRQ(ierr); 2124 PetscFunctionReturn(0); 2125 } 2126 2127 /*@ 2128 MatIsShell - Inquires if a matrix is derived from MATSHELL 2129 2130 Input Parameter: 2131 . mat - the matrix 2132 2133 Output Parameter: 2134 . flg - the boolean value 2135 2136 Level: developer 2137 2138 Notes: in the future, we should allow the object type name to be changed still using the MatShell data structure for other matrices (i.e. MATTRANSPOSEMAT, MATSCHURCOMPLEMENT etc) 2139 2140 .seealso: MatCreateShell() 2141 @*/ 2142 PetscErrorCode MatIsShell(Mat mat, PetscBool *flg) 2143 { 2144 PetscFunctionBegin; 2145 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2146 PetscValidPointer(flg,2); 2147 *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell); 2148 PetscFunctionReturn(0); 2149 } 2150