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