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) break; 860 matmat = entry; 861 entry = entry->next; 862 } 863 if (!flg) { 864 ierr = PetscNew(&matmat->next);CHKERRQ(ierr); 865 matmat = matmat->next; 866 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen"); 867 } 868 869 matmat->symbolic = symbolic; 870 matmat->numeric = numeric; 871 matmat->destroy = destroy; 872 matmat->ptype = ptype; 873 ierr = PetscFree(matmat->composedname);CHKERRQ(ierr); 874 ierr = PetscFree(matmat->resultname);CHKERRQ(ierr); 875 ierr = PetscStrallocpy(composedname,&matmat->composedname);CHKERRQ(ierr); 876 ierr = PetscStrallocpy(resultname,&matmat->resultname);CHKERRQ(ierr); 877 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); 878 ierr = PetscObjectComposeFunction((PetscObject)A,matmat->composedname,MatProductSetFromOptions_Shell_X);CHKERRQ(ierr); 879 PetscFunctionReturn(0); 880 } 881 882 /*@C 883 MatShellSetMatProductOperation - Allows user to set a matrix matrix operation for a shell matrix. 884 885 Logically Collective on Mat 886 887 Input Parameters: 888 + A - the shell matrix 889 . ptype - the product type 890 . symbolic - the function for the symbolic phase (can be NULL) 891 . numeric - the function for the numerical phase 892 . destroy - the function for the destruction of the needed data generated during the symbolic phase (can be NULL) 893 . Btype - the matrix type for the matrix to be multiplied against 894 - Ctype - the matrix type for the result (can be NULL) 895 896 Level: advanced 897 898 Usage: 899 $ extern PetscErrorCode usersymbolic(Mat,Mat,Mat,void**); 900 $ extern PetscErrorCode usernumeric(Mat,Mat,Mat,void*); 901 $ extern PetscErrorCode userdestroy(void*); 902 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 903 $ MatShellSetMatProductOperation(A,MATPRODUCT_AB,usersymbolic,usernumeric,userdestroy,MATSEQAIJ,MATDENSE); 904 $ [ create B of type SEQAIJ etc..] 905 $ MatProductCreate(A,B,NULL,&C); 906 $ MatProductSetType(C,MATPRODUCT_AB); 907 $ MatProductSetFromOptions(C); 908 $ MatProductSymbolic(C); -> actually runs the user defined symbolic operation 909 $ MatProductNumeric(C); -> actually runs the user defined numeric operation 910 $ [ use C = A*B ] 911 912 Notes: 913 MATPRODUCT_ABC is not supported yet. Not supported in Fortran. 914 If the symbolic phase is not specified, MatSetUp() is called on the result matrix that must have its type set if Ctype is NULL. 915 Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback. 916 PETSc will take care of calling the user-defined callbacks. 917 It is allowed to specify the same callbacks for different Btype matrix types. 918 The couple (Btype,ptype) uniquely identifies the operation: the last specified callbacks takes precedence. 919 920 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatProductType, MatType, MatSetUp() 921 @*/ 922 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) 923 { 924 PetscErrorCode ierr; 925 926 PetscFunctionBegin; 927 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 928 PetscValidLogicalCollectiveEnum(A,ptype,2); 929 if (ptype == MATPRODUCT_ABC) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not for product type %s",MatProductTypes[ptype]); 930 if (!numeric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Missing numeric routine, argument 4"); 931 PetscValidPointer(Btype,6); 932 if (Ctype) PetscValidPointer(Ctype,7); 933 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); 934 PetscFunctionReturn(0); 935 } 936 937 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) 938 { 939 PetscBool flg; 940 PetscErrorCode ierr; 941 char composedname[256]; 942 MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList; 943 PetscMPIInt size; 944 945 PetscFunctionBegin; 946 PetscValidType(A,1); 947 while (Bnames) { /* user passed in the root name */ 948 ierr = PetscStrcmp(Btype,Bnames->rname,&flg);CHKERRQ(ierr); 949 if (flg) break; 950 Bnames = Bnames->next; 951 } 952 while (Cnames) { /* user passed in the root name */ 953 ierr = PetscStrcmp(Ctype,Cnames->rname,&flg);CHKERRQ(ierr); 954 if (flg) break; 955 Cnames = Cnames->next; 956 } 957 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 958 Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype; 959 Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype; 960 ierr = PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,Btype);CHKERRQ(ierr); 961 ierr = MatShellSetMatProductOperation_Private(A,ptype,symbolic,numeric,destroy,composedname,Ctype);CHKERRQ(ierr); 962 PetscFunctionReturn(0); 963 } 964 965 PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str) 966 { 967 Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data; 968 PetscErrorCode ierr; 969 PetscBool matflg; 970 MatShellMatFunctionList matmatA; 971 972 PetscFunctionBegin; 973 ierr = MatIsShell(B,&matflg);CHKERRQ(ierr); 974 if (!matflg) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix %s not derived from MATSHELL",((PetscObject)B)->type_name); 975 976 ierr = PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 977 ierr = PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));CHKERRQ(ierr); 978 979 if (shellA->ops->copy) { 980 ierr = (*shellA->ops->copy)(A,B,str);CHKERRQ(ierr); 981 } 982 shellB->vscale = shellA->vscale; 983 shellB->vshift = shellA->vshift; 984 if (shellA->dshift) { 985 if (!shellB->dshift) { 986 ierr = VecDuplicate(shellA->dshift,&shellB->dshift);CHKERRQ(ierr); 987 } 988 ierr = VecCopy(shellA->dshift,shellB->dshift);CHKERRQ(ierr); 989 } else { 990 ierr = VecDestroy(&shellB->dshift);CHKERRQ(ierr); 991 } 992 if (shellA->left) { 993 if (!shellB->left) { 994 ierr = VecDuplicate(shellA->left,&shellB->left);CHKERRQ(ierr); 995 } 996 ierr = VecCopy(shellA->left,shellB->left);CHKERRQ(ierr); 997 } else { 998 ierr = VecDestroy(&shellB->left);CHKERRQ(ierr); 999 } 1000 if (shellA->right) { 1001 if (!shellB->right) { 1002 ierr = VecDuplicate(shellA->right,&shellB->right);CHKERRQ(ierr); 1003 } 1004 ierr = VecCopy(shellA->right,shellB->right);CHKERRQ(ierr); 1005 } else { 1006 ierr = VecDestroy(&shellB->right);CHKERRQ(ierr); 1007 } 1008 ierr = MatDestroy(&shellB->axpy);CHKERRQ(ierr); 1009 shellB->axpy_vscale = 0.0; 1010 shellB->axpy_state = 0; 1011 if (shellA->axpy) { 1012 ierr = PetscObjectReference((PetscObject)shellA->axpy);CHKERRQ(ierr); 1013 shellB->axpy = shellA->axpy; 1014 shellB->axpy_vscale = shellA->axpy_vscale; 1015 shellB->axpy_state = shellA->axpy_state; 1016 } 1017 if (shellA->zrows) { 1018 ierr = ISDuplicate(shellA->zrows,&shellB->zrows);CHKERRQ(ierr); 1019 if (shellA->zcols) { 1020 ierr = ISDuplicate(shellA->zcols,&shellB->zcols);CHKERRQ(ierr); 1021 } 1022 ierr = VecDuplicate(shellA->zvals,&shellB->zvals);CHKERRQ(ierr); 1023 ierr = VecCopy(shellA->zvals,shellB->zvals);CHKERRQ(ierr); 1024 ierr = VecDuplicate(shellA->zvals_w,&shellB->zvals_w);CHKERRQ(ierr); 1025 ierr = PetscObjectReference((PetscObject)shellA->zvals_sct_r);CHKERRQ(ierr); 1026 ierr = PetscObjectReference((PetscObject)shellA->zvals_sct_c);CHKERRQ(ierr); 1027 shellB->zvals_sct_r = shellA->zvals_sct_r; 1028 shellB->zvals_sct_c = shellA->zvals_sct_c; 1029 } 1030 1031 matmatA = shellA->matmat; 1032 if (matmatA) { 1033 while (matmatA->next) { 1034 ierr = MatShellSetMatProductOperation_Private(B,matmatA->ptype,matmatA->symbolic,matmatA->numeric,matmatA->destroy,matmatA->composedname,matmatA->resultname);CHKERRQ(ierr); 1035 matmatA = matmatA->next; 1036 } 1037 } 1038 PetscFunctionReturn(0); 1039 } 1040 1041 PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M) 1042 { 1043 PetscErrorCode ierr; 1044 void *ctx; 1045 1046 PetscFunctionBegin; 1047 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr); 1048 ierr = MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);CHKERRQ(ierr); 1049 ierr = PetscObjectChangeTypeName((PetscObject)(*M),((PetscObject)mat)->type_name);CHKERRQ(ierr); 1050 if (op != MAT_DO_NOT_COPY_VALUES) { 1051 ierr = MatCopy(mat,*M,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 1052 } 1053 PetscFunctionReturn(0); 1054 } 1055 1056 PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y) 1057 { 1058 Mat_Shell *shell = (Mat_Shell*)A->data; 1059 PetscErrorCode ierr; 1060 Vec xx; 1061 PetscObjectState instate,outstate; 1062 1063 PetscFunctionBegin; 1064 if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL"); 1065 ierr = MatShellPreZeroRight(A,x,&xx);CHKERRQ(ierr); 1066 ierr = MatShellPreScaleRight(A,xx,&xx);CHKERRQ(ierr); 1067 ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr); 1068 ierr = (*shell->ops->mult)(A,xx,y);CHKERRQ(ierr); 1069 ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr); 1070 if (instate == outstate) { 1071 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 1072 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 1073 } 1074 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 1075 ierr = MatShellPostScaleLeft(A,y);CHKERRQ(ierr); 1076 ierr = MatShellPostZeroLeft(A,y);CHKERRQ(ierr); 1077 1078 if (shell->axpy) { 1079 Mat X; 1080 PetscObjectState axpy_state; 1081 1082 ierr = MatShellGetContext(shell->axpy,(void *)&X);CHKERRQ(ierr); 1083 ierr = PetscObjectStateGet((PetscObject)X,&axpy_state);CHKERRQ(ierr); 1084 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,...)"); 1085 1086 ierr = MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);CHKERRQ(ierr); 1087 ierr = VecCopy(x,shell->axpy_right);CHKERRQ(ierr); 1088 ierr = MatMult(shell->axpy,shell->axpy_right,shell->axpy_left);CHKERRQ(ierr); 1089 ierr = VecAXPY(y,shell->axpy_vscale,shell->axpy_left);CHKERRQ(ierr); 1090 } 1091 PetscFunctionReturn(0); 1092 } 1093 1094 PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z) 1095 { 1096 Mat_Shell *shell = (Mat_Shell*)A->data; 1097 PetscErrorCode ierr; 1098 1099 PetscFunctionBegin; 1100 if (y == z) { 1101 if (!shell->right_add_work) {ierr = VecDuplicate(z,&shell->right_add_work);CHKERRQ(ierr);} 1102 ierr = MatMult(A,x,shell->right_add_work);CHKERRQ(ierr); 1103 ierr = VecAXPY(z,1.0,shell->right_add_work);CHKERRQ(ierr); 1104 } else { 1105 ierr = MatMult(A,x,z);CHKERRQ(ierr); 1106 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 1107 } 1108 PetscFunctionReturn(0); 1109 } 1110 1111 PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y) 1112 { 1113 Mat_Shell *shell = (Mat_Shell*)A->data; 1114 PetscErrorCode ierr; 1115 Vec xx; 1116 PetscObjectState instate,outstate; 1117 1118 PetscFunctionBegin; 1119 if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL"); 1120 ierr = MatShellPreZeroLeft(A,x,&xx);CHKERRQ(ierr); 1121 ierr = MatShellPreScaleLeft(A,xx,&xx);CHKERRQ(ierr); 1122 ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr); 1123 ierr = (*shell->ops->multtranspose)(A,xx,y);CHKERRQ(ierr); 1124 ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr); 1125 if (instate == outstate) { 1126 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 1127 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 1128 } 1129 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 1130 ierr = MatShellPostScaleRight(A,y);CHKERRQ(ierr); 1131 ierr = MatShellPostZeroRight(A,y);CHKERRQ(ierr); 1132 1133 if (shell->axpy) { 1134 Mat X; 1135 PetscObjectState axpy_state; 1136 1137 ierr = MatShellGetContext(shell->axpy,(void *)&X);CHKERRQ(ierr); 1138 ierr = PetscObjectStateGet((PetscObject)X,&axpy_state);CHKERRQ(ierr); 1139 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,...)"); 1140 ierr = MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);CHKERRQ(ierr); 1141 ierr = VecCopy(x,shell->axpy_left);CHKERRQ(ierr); 1142 ierr = MatMultTranspose(shell->axpy,shell->axpy_left,shell->axpy_right);CHKERRQ(ierr); 1143 ierr = VecAXPY(y,shell->axpy_vscale,shell->axpy_right);CHKERRQ(ierr); 1144 } 1145 PetscFunctionReturn(0); 1146 } 1147 1148 PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z) 1149 { 1150 Mat_Shell *shell = (Mat_Shell*)A->data; 1151 PetscErrorCode ierr; 1152 1153 PetscFunctionBegin; 1154 if (y == z) { 1155 if (!shell->left_add_work) {ierr = VecDuplicate(z,&shell->left_add_work);CHKERRQ(ierr);} 1156 ierr = MatMultTranspose(A,x,shell->left_add_work);CHKERRQ(ierr); 1157 ierr = VecAXPY(z,1.0,shell->left_add_work);CHKERRQ(ierr); 1158 } else { 1159 ierr = MatMultTranspose(A,x,z);CHKERRQ(ierr); 1160 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 1161 } 1162 PetscFunctionReturn(0); 1163 } 1164 1165 /* 1166 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1167 */ 1168 PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v) 1169 { 1170 Mat_Shell *shell = (Mat_Shell*)A->data; 1171 PetscErrorCode ierr; 1172 1173 PetscFunctionBegin; 1174 if (shell->ops->getdiagonal) { 1175 ierr = (*shell->ops->getdiagonal)(A,v);CHKERRQ(ierr); 1176 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)"); 1177 ierr = VecScale(v,shell->vscale);CHKERRQ(ierr); 1178 if (shell->dshift) { 1179 ierr = VecAXPY(v,1.0,shell->dshift);CHKERRQ(ierr); 1180 } 1181 ierr = VecShift(v,shell->vshift);CHKERRQ(ierr); 1182 if (shell->left) {ierr = VecPointwiseMult(v,v,shell->left);CHKERRQ(ierr);} 1183 if (shell->right) {ierr = VecPointwiseMult(v,v,shell->right);CHKERRQ(ierr);} 1184 if (shell->zrows) { 1185 ierr = VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1186 ierr = VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1187 } 1188 if (shell->axpy) { 1189 Mat X; 1190 PetscObjectState axpy_state; 1191 1192 ierr = MatShellGetContext(shell->axpy,(void *)&X);CHKERRQ(ierr); 1193 ierr = PetscObjectStateGet((PetscObject)X,&axpy_state);CHKERRQ(ierr); 1194 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,...)"); 1195 ierr = MatCreateVecs(shell->axpy,NULL,shell->axpy_left ? NULL : &shell->axpy_left);CHKERRQ(ierr); 1196 ierr = MatGetDiagonal(shell->axpy,shell->axpy_left);CHKERRQ(ierr); 1197 ierr = VecAXPY(v,shell->axpy_vscale,shell->axpy_left);CHKERRQ(ierr); 1198 } 1199 PetscFunctionReturn(0); 1200 } 1201 1202 PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a) 1203 { 1204 Mat_Shell *shell = (Mat_Shell*)Y->data; 1205 PetscErrorCode ierr; 1206 PetscBool flg; 1207 1208 PetscFunctionBegin; 1209 ierr = MatHasCongruentLayouts(Y,&flg);CHKERRQ(ierr); 1210 if (!flg) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot shift shell matrix if it is not congruent"); 1211 if (shell->left || shell->right) { 1212 if (!shell->dshift) { 1213 ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);CHKERRQ(ierr); 1214 ierr = VecSet(shell->dshift,a);CHKERRQ(ierr); 1215 } else { 1216 if (shell->left) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 1217 if (shell->right) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 1218 ierr = VecShift(shell->dshift,a);CHKERRQ(ierr); 1219 } 1220 if (shell->left) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 1221 if (shell->right) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 1222 } else shell->vshift += a; 1223 if (shell->zrows) { 1224 ierr = VecShift(shell->zvals,a);CHKERRQ(ierr); 1225 } 1226 PetscFunctionReturn(0); 1227 } 1228 1229 PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s) 1230 { 1231 Mat_Shell *shell = (Mat_Shell*)A->data; 1232 PetscErrorCode ierr; 1233 1234 PetscFunctionBegin; 1235 if (!shell->dshift) {ierr = VecDuplicate(D,&shell->dshift);CHKERRQ(ierr);} 1236 if (shell->left || shell->right) { 1237 if (!shell->right_work) {ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);CHKERRQ(ierr);} 1238 if (shell->left && shell->right) { 1239 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 1240 ierr = VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);CHKERRQ(ierr); 1241 } else if (shell->left) { 1242 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 1243 } else { 1244 ierr = VecPointwiseDivide(shell->right_work,D,shell->right);CHKERRQ(ierr); 1245 } 1246 ierr = VecAXPY(shell->dshift,s,shell->right_work);CHKERRQ(ierr); 1247 } else { 1248 ierr = VecAXPY(shell->dshift,s,D);CHKERRQ(ierr); 1249 } 1250 PetscFunctionReturn(0); 1251 } 1252 1253 PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins) 1254 { 1255 Mat_Shell *shell = (Mat_Shell*)A->data; 1256 Vec d; 1257 PetscErrorCode ierr; 1258 PetscBool flg; 1259 1260 PetscFunctionBegin; 1261 ierr = MatHasCongruentLayouts(A,&flg);CHKERRQ(ierr); 1262 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot diagonal set or shift shell matrix if it is not congruent"); 1263 if (ins == INSERT_VALUES) { 1264 if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first"); 1265 ierr = VecDuplicate(D,&d);CHKERRQ(ierr); 1266 ierr = MatGetDiagonal(A,d);CHKERRQ(ierr); 1267 ierr = MatDiagonalSet_Shell_Private(A,d,-1.);CHKERRQ(ierr); 1268 ierr = MatDiagonalSet_Shell_Private(A,D,1.);CHKERRQ(ierr); 1269 ierr = VecDestroy(&d);CHKERRQ(ierr); 1270 if (shell->zrows) { 1271 ierr = VecCopy(D,shell->zvals);CHKERRQ(ierr); 1272 } 1273 } else { 1274 ierr = MatDiagonalSet_Shell_Private(A,D,1.);CHKERRQ(ierr); 1275 if (shell->zrows) { 1276 ierr = VecAXPY(shell->zvals,1.0,D);CHKERRQ(ierr); 1277 } 1278 } 1279 PetscFunctionReturn(0); 1280 } 1281 1282 PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a) 1283 { 1284 Mat_Shell *shell = (Mat_Shell*)Y->data; 1285 PetscErrorCode ierr; 1286 1287 PetscFunctionBegin; 1288 shell->vscale *= a; 1289 shell->vshift *= a; 1290 if (shell->dshift) { 1291 ierr = VecScale(shell->dshift,a);CHKERRQ(ierr); 1292 } 1293 shell->axpy_vscale *= a; 1294 if (shell->zrows) { 1295 ierr = VecScale(shell->zvals,a);CHKERRQ(ierr); 1296 } 1297 PetscFunctionReturn(0); 1298 } 1299 1300 static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right) 1301 { 1302 Mat_Shell *shell = (Mat_Shell*)Y->data; 1303 PetscErrorCode ierr; 1304 1305 PetscFunctionBegin; 1306 if (left) { 1307 if (!shell->left) { 1308 ierr = VecDuplicate(left,&shell->left);CHKERRQ(ierr); 1309 ierr = VecCopy(left,shell->left);CHKERRQ(ierr); 1310 } else { 1311 ierr = VecPointwiseMult(shell->left,shell->left,left);CHKERRQ(ierr); 1312 } 1313 if (shell->zrows) { 1314 ierr = VecPointwiseMult(shell->zvals,shell->zvals,left);CHKERRQ(ierr); 1315 } 1316 } 1317 if (right) { 1318 if (!shell->right) { 1319 ierr = VecDuplicate(right,&shell->right);CHKERRQ(ierr); 1320 ierr = VecCopy(right,shell->right);CHKERRQ(ierr); 1321 } else { 1322 ierr = VecPointwiseMult(shell->right,shell->right,right);CHKERRQ(ierr); 1323 } 1324 if (shell->zrows) { 1325 if (!shell->left_work) { 1326 ierr = MatCreateVecs(Y,NULL,&shell->left_work);CHKERRQ(ierr); 1327 } 1328 ierr = VecSet(shell->zvals_w,1.0);CHKERRQ(ierr); 1329 ierr = VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1330 ierr = VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1331 ierr = VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w);CHKERRQ(ierr); 1332 } 1333 } 1334 if (shell->axpy) { 1335 ierr = MatDiagonalScale(shell->axpy,left,right);CHKERRQ(ierr); 1336 } 1337 PetscFunctionReturn(0); 1338 } 1339 1340 PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t) 1341 { 1342 Mat_Shell *shell = (Mat_Shell*)Y->data; 1343 PetscErrorCode ierr; 1344 1345 PetscFunctionBegin; 1346 if (t == MAT_FINAL_ASSEMBLY) { 1347 shell->vshift = 0.0; 1348 shell->vscale = 1.0; 1349 shell->axpy_vscale = 0.0; 1350 shell->axpy_state = 0; 1351 ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr); 1352 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 1353 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 1354 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 1355 ierr = VecDestroy(&shell->axpy_left);CHKERRQ(ierr); 1356 ierr = VecDestroy(&shell->axpy_right);CHKERRQ(ierr); 1357 ierr = VecScatterDestroy(&shell->zvals_sct_c);CHKERRQ(ierr); 1358 ierr = VecScatterDestroy(&shell->zvals_sct_r);CHKERRQ(ierr); 1359 ierr = ISDestroy(&shell->zrows);CHKERRQ(ierr); 1360 ierr = ISDestroy(&shell->zcols);CHKERRQ(ierr); 1361 } 1362 PetscFunctionReturn(0); 1363 } 1364 1365 static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d) 1366 { 1367 PetscFunctionBegin; 1368 *missing = PETSC_FALSE; 1369 PetscFunctionReturn(0); 1370 } 1371 1372 PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str) 1373 { 1374 Mat_Shell *shell = (Mat_Shell*)Y->data; 1375 PetscErrorCode ierr; 1376 1377 PetscFunctionBegin; 1378 if (X == Y) { 1379 ierr = MatScale(Y,1.0 + a);CHKERRQ(ierr); 1380 PetscFunctionReturn(0); 1381 } 1382 if (!shell->axpy) { 1383 ierr = MatConvertFrom_Shell(X,MATSHELL,MAT_INITIAL_MATRIX,&shell->axpy);CHKERRQ(ierr); 1384 shell->axpy_vscale = a; 1385 ierr = PetscObjectStateGet((PetscObject)X,&shell->axpy_state);CHKERRQ(ierr); 1386 } else { 1387 ierr = MatAXPY(shell->axpy,a/shell->axpy_vscale,X,str);CHKERRQ(ierr); 1388 } 1389 PetscFunctionReturn(0); 1390 } 1391 1392 static struct _MatOps MatOps_Values = {NULL, 1393 NULL, 1394 NULL, 1395 NULL, 1396 /* 4*/ MatMultAdd_Shell, 1397 NULL, 1398 MatMultTransposeAdd_Shell, 1399 NULL, 1400 NULL, 1401 NULL, 1402 /*10*/ NULL, 1403 NULL, 1404 NULL, 1405 NULL, 1406 NULL, 1407 /*15*/ NULL, 1408 NULL, 1409 NULL, 1410 MatDiagonalScale_Shell, 1411 NULL, 1412 /*20*/ NULL, 1413 MatAssemblyEnd_Shell, 1414 NULL, 1415 NULL, 1416 /*24*/ MatZeroRows_Shell, 1417 NULL, 1418 NULL, 1419 NULL, 1420 NULL, 1421 /*29*/ NULL, 1422 NULL, 1423 NULL, 1424 NULL, 1425 NULL, 1426 /*34*/ MatDuplicate_Shell, 1427 NULL, 1428 NULL, 1429 NULL, 1430 NULL, 1431 /*39*/ MatAXPY_Shell, 1432 NULL, 1433 NULL, 1434 NULL, 1435 MatCopy_Shell, 1436 /*44*/ NULL, 1437 MatScale_Shell, 1438 MatShift_Shell, 1439 MatDiagonalSet_Shell, 1440 MatZeroRowsColumns_Shell, 1441 /*49*/ NULL, 1442 NULL, 1443 NULL, 1444 NULL, 1445 NULL, 1446 /*54*/ NULL, 1447 NULL, 1448 NULL, 1449 NULL, 1450 NULL, 1451 /*59*/ NULL, 1452 MatDestroy_Shell, 1453 NULL, 1454 MatConvertFrom_Shell, 1455 NULL, 1456 /*64*/ NULL, 1457 NULL, 1458 NULL, 1459 NULL, 1460 NULL, 1461 /*69*/ NULL, 1462 NULL, 1463 MatConvert_Shell, 1464 NULL, 1465 NULL, 1466 /*74*/ NULL, 1467 NULL, 1468 NULL, 1469 NULL, 1470 NULL, 1471 /*79*/ NULL, 1472 NULL, 1473 NULL, 1474 NULL, 1475 NULL, 1476 /*84*/ NULL, 1477 NULL, 1478 NULL, 1479 NULL, 1480 NULL, 1481 /*89*/ NULL, 1482 NULL, 1483 NULL, 1484 NULL, 1485 NULL, 1486 /*94*/ NULL, 1487 NULL, 1488 NULL, 1489 NULL, 1490 NULL, 1491 /*99*/ NULL, 1492 NULL, 1493 NULL, 1494 NULL, 1495 NULL, 1496 /*104*/ NULL, 1497 NULL, 1498 NULL, 1499 NULL, 1500 NULL, 1501 /*109*/ NULL, 1502 NULL, 1503 NULL, 1504 NULL, 1505 MatMissingDiagonal_Shell, 1506 /*114*/ NULL, 1507 NULL, 1508 NULL, 1509 NULL, 1510 NULL, 1511 /*119*/ NULL, 1512 NULL, 1513 NULL, 1514 NULL, 1515 NULL, 1516 /*124*/ NULL, 1517 NULL, 1518 NULL, 1519 NULL, 1520 NULL, 1521 /*129*/ NULL, 1522 NULL, 1523 NULL, 1524 NULL, 1525 NULL, 1526 /*134*/ NULL, 1527 NULL, 1528 NULL, 1529 NULL, 1530 NULL, 1531 /*139*/ NULL, 1532 NULL, 1533 NULL 1534 }; 1535 1536 PetscErrorCode MatShellSetContext_Shell(Mat mat,void *ctx) 1537 { 1538 Mat_Shell *shell = (Mat_Shell*)mat->data; 1539 1540 PetscFunctionBegin; 1541 shell->ctx = ctx; 1542 PetscFunctionReturn(0); 1543 } 1544 1545 static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype) 1546 { 1547 PetscErrorCode ierr; 1548 1549 PetscFunctionBegin; 1550 ierr = PetscFree(mat->defaultvectype);CHKERRQ(ierr); 1551 ierr = PetscStrallocpy(vtype,(char**)&mat->defaultvectype);CHKERRQ(ierr); 1552 PetscFunctionReturn(0); 1553 } 1554 1555 PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A) 1556 { 1557 Mat_Shell *shell = (Mat_Shell*)A->data; 1558 1559 PetscFunctionBegin; 1560 shell->managescalingshifts = PETSC_FALSE; 1561 A->ops->diagonalset = NULL; 1562 A->ops->diagonalscale = NULL; 1563 A->ops->scale = NULL; 1564 A->ops->shift = NULL; 1565 A->ops->axpy = NULL; 1566 PetscFunctionReturn(0); 1567 } 1568 1569 PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void)) 1570 { 1571 Mat_Shell *shell = (Mat_Shell*)mat->data; 1572 1573 PetscFunctionBegin; 1574 switch (op) { 1575 case MATOP_DESTROY: 1576 shell->ops->destroy = (PetscErrorCode (*)(Mat))f; 1577 break; 1578 case MATOP_VIEW: 1579 if (!mat->ops->viewnative) { 1580 mat->ops->viewnative = mat->ops->view; 1581 } 1582 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 1583 break; 1584 case MATOP_COPY: 1585 shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f; 1586 break; 1587 case MATOP_DIAGONAL_SET: 1588 case MATOP_DIAGONAL_SCALE: 1589 case MATOP_SHIFT: 1590 case MATOP_SCALE: 1591 case MATOP_AXPY: 1592 case MATOP_ZERO_ROWS: 1593 case MATOP_ZERO_ROWS_COLUMNS: 1594 if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1595 (((void(**)(void))mat->ops)[op]) = f; 1596 break; 1597 case MATOP_GET_DIAGONAL: 1598 if (shell->managescalingshifts) { 1599 shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1600 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1601 } else { 1602 shell->ops->getdiagonal = NULL; 1603 mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1604 } 1605 break; 1606 case MATOP_MULT: 1607 if (shell->managescalingshifts) { 1608 shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1609 mat->ops->mult = MatMult_Shell; 1610 } else { 1611 shell->ops->mult = NULL; 1612 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1613 } 1614 break; 1615 case MATOP_MULT_TRANSPOSE: 1616 if (shell->managescalingshifts) { 1617 shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1618 mat->ops->multtranspose = MatMultTranspose_Shell; 1619 } else { 1620 shell->ops->multtranspose = NULL; 1621 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1622 } 1623 break; 1624 default: 1625 (((void(**)(void))mat->ops)[op]) = f; 1626 break; 1627 } 1628 PetscFunctionReturn(0); 1629 } 1630 1631 PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void)) 1632 { 1633 Mat_Shell *shell = (Mat_Shell*)mat->data; 1634 1635 PetscFunctionBegin; 1636 switch (op) { 1637 case MATOP_DESTROY: 1638 *f = (void (*)(void))shell->ops->destroy; 1639 break; 1640 case MATOP_VIEW: 1641 *f = (void (*)(void))mat->ops->view; 1642 break; 1643 case MATOP_COPY: 1644 *f = (void (*)(void))shell->ops->copy; 1645 break; 1646 case MATOP_DIAGONAL_SET: 1647 case MATOP_DIAGONAL_SCALE: 1648 case MATOP_SHIFT: 1649 case MATOP_SCALE: 1650 case MATOP_AXPY: 1651 case MATOP_ZERO_ROWS: 1652 case MATOP_ZERO_ROWS_COLUMNS: 1653 *f = (((void (**)(void))mat->ops)[op]); 1654 break; 1655 case MATOP_GET_DIAGONAL: 1656 if (shell->ops->getdiagonal) 1657 *f = (void (*)(void))shell->ops->getdiagonal; 1658 else 1659 *f = (((void (**)(void))mat->ops)[op]); 1660 break; 1661 case MATOP_MULT: 1662 if (shell->ops->mult) 1663 *f = (void (*)(void))shell->ops->mult; 1664 else 1665 *f = (((void (**)(void))mat->ops)[op]); 1666 break; 1667 case MATOP_MULT_TRANSPOSE: 1668 if (shell->ops->multtranspose) 1669 *f = (void (*)(void))shell->ops->multtranspose; 1670 else 1671 *f = (((void (**)(void))mat->ops)[op]); 1672 break; 1673 default: 1674 *f = (((void (**)(void))mat->ops)[op]); 1675 } 1676 PetscFunctionReturn(0); 1677 } 1678 1679 /*MC 1680 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 1681 1682 Level: advanced 1683 1684 .seealso: MatCreateShell() 1685 M*/ 1686 1687 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 1688 { 1689 Mat_Shell *b; 1690 PetscErrorCode ierr; 1691 1692 PetscFunctionBegin; 1693 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1694 1695 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 1696 A->data = (void*)b; 1697 1698 b->ctx = NULL; 1699 b->vshift = 0.0; 1700 b->vscale = 1.0; 1701 b->managescalingshifts = PETSC_TRUE; 1702 A->assembled = PETSC_TRUE; 1703 A->preallocated = PETSC_FALSE; 1704 1705 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell);CHKERRQ(ierr); 1706 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell);CHKERRQ(ierr); 1707 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellSetVecType_C",MatShellSetVecType_Shell);CHKERRQ(ierr); 1708 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell);CHKERRQ(ierr); 1709 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell);CHKERRQ(ierr); 1710 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell);CHKERRQ(ierr); 1711 ierr = PetscObjectComposeFunction((PetscObject)A,"MatShellSetMatProductOperation_C",MatShellSetMatProductOperation_Shell);CHKERRQ(ierr); 1712 ierr = PetscObjectChangeTypeName((PetscObject)A,MATSHELL);CHKERRQ(ierr); 1713 PetscFunctionReturn(0); 1714 } 1715 1716 /*@C 1717 MatCreateShell - Creates a new matrix class for use with a user-defined 1718 private data storage format. 1719 1720 Collective 1721 1722 Input Parameters: 1723 + comm - MPI communicator 1724 . m - number of local rows (must be given) 1725 . n - number of local columns (must be given) 1726 . M - number of global rows (may be PETSC_DETERMINE) 1727 . N - number of global columns (may be PETSC_DETERMINE) 1728 - ctx - pointer to data needed by the shell matrix routines 1729 1730 Output Parameter: 1731 . A - the matrix 1732 1733 Level: advanced 1734 1735 Usage: 1736 $ extern PetscErrorCode mult(Mat,Vec,Vec); 1737 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 1738 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1739 $ [ Use matrix for operations that have been set ] 1740 $ MatDestroy(mat); 1741 1742 Notes: 1743 The shell matrix type is intended to provide a simple class to use 1744 with KSP (such as, for use with matrix-free methods). You should not 1745 use the shell type if you plan to define a complete matrix class. 1746 1747 Fortran Notes: 1748 To use this from Fortran with a ctx you must write an interface definition for this 1749 function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing 1750 in as the ctx argument. 1751 1752 PETSc requires that matrices and vectors being used for certain 1753 operations are partitioned accordingly. For example, when 1754 creating a shell matrix, A, that supports parallel matrix-vector 1755 products using MatMult(A,x,y) the user should set the number 1756 of local matrix rows to be the number of local elements of the 1757 corresponding result vector, y. Note that this is information is 1758 required for use of the matrix interface routines, even though 1759 the shell matrix may not actually be physically partitioned. 1760 For example, 1761 1762 $ 1763 $ Vec x, y 1764 $ extern PetscErrorCode mult(Mat,Vec,Vec); 1765 $ Mat A 1766 $ 1767 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 1768 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 1769 $ VecGetLocalSize(y,&m); 1770 $ VecGetLocalSize(x,&n); 1771 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 1772 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1773 $ MatMult(A,x,y); 1774 $ MatDestroy(&A); 1775 $ VecDestroy(&y); 1776 $ VecDestroy(&x); 1777 $ 1778 1779 1780 MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), MatScale(), MatZeroRows() and MatZeroRowsColumns() internally, so these 1781 operations cannot be overwritten unless MatShellSetManageScalingShifts() is called. 1782 1783 1784 For rectangular matrices do all the scalings and shifts make sense? 1785 1786 Developers Notes: 1787 Regarding shifting and scaling. The general form is 1788 1789 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1790 1791 The order you apply the operations is important. For example if you have a dshift then 1792 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift 1793 you get s*vscale*A + diag(shift) 1794 1795 A is the user provided function. 1796 1797 KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for 1798 for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger 1799 an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat); 1800 each time the MATSHELL matrix has changed. 1801 1802 Matrix product operations (i.e. MatMat, MatTranposeMat etc) can be specified using MatShellSetMatProductOperation() 1803 1804 Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided 1805 with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale(). 1806 1807 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts(), MatShellSetMatProductOperation() 1808 @*/ 1809 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 1810 { 1811 PetscErrorCode ierr; 1812 1813 PetscFunctionBegin; 1814 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1815 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 1816 ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr); 1817 ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr); 1818 ierr = MatSetUp(*A);CHKERRQ(ierr); 1819 PetscFunctionReturn(0); 1820 } 1821 1822 /*@ 1823 MatShellSetContext - sets the context for a shell matrix 1824 1825 Logically Collective on Mat 1826 1827 Input Parameters: 1828 + mat - the shell matrix 1829 - ctx - the context 1830 1831 Level: advanced 1832 1833 Fortran Notes: 1834 To use this from Fortran you must write a Fortran interface definition for this 1835 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 1836 1837 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation() 1838 @*/ 1839 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 1840 { 1841 PetscErrorCode ierr; 1842 1843 PetscFunctionBegin; 1844 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1845 ierr = PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));CHKERRQ(ierr); 1846 PetscFunctionReturn(0); 1847 } 1848 1849 /*@C 1850 MatShellSetVecType - Sets the type of Vec returned by MatCreateVecs() 1851 1852 Logically collective 1853 1854 Input Parameters: 1855 + mat - the shell matrix 1856 - vtype - type to use for creating vectors 1857 1858 Notes: 1859 1860 Level: advanced 1861 1862 .seealso: MatCreateVecs() 1863 @*/ 1864 PetscErrorCode MatShellSetVecType(Mat mat,VecType vtype) 1865 { 1866 PetscErrorCode ierr; 1867 1868 PetscFunctionBegin; 1869 ierr = PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype));CHKERRQ(ierr); 1870 PetscFunctionReturn(0); 1871 } 1872 1873 /*@ 1874 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately 1875 after MatCreateShell() 1876 1877 Logically Collective on Mat 1878 1879 Input Parameter: 1880 . mat - the shell matrix 1881 1882 Level: advanced 1883 1884 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation() 1885 @*/ 1886 PetscErrorCode MatShellSetManageScalingShifts(Mat A) 1887 { 1888 PetscErrorCode ierr; 1889 1890 PetscFunctionBegin; 1891 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1892 ierr = PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));CHKERRQ(ierr); 1893 PetscFunctionReturn(0); 1894 } 1895 1896 /*@C 1897 MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function. 1898 1899 Logically Collective on Mat 1900 1901 Input Parameters: 1902 + mat - the shell matrix 1903 . f - the function 1904 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 1905 - ctx - an optional context for the function 1906 1907 Output Parameter: 1908 . flg - PETSC_TRUE if the multiply is likely correct 1909 1910 Options Database: 1911 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1912 1913 Level: advanced 1914 1915 Fortran Notes: 1916 Not supported from Fortran 1917 1918 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose() 1919 @*/ 1920 PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 1921 { 1922 PetscErrorCode ierr; 1923 PetscInt m,n; 1924 Mat mf,Dmf,Dmat,Ddiff; 1925 PetscReal Diffnorm,Dmfnorm; 1926 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 1927 1928 PetscFunctionBegin; 1929 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1930 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);CHKERRQ(ierr); 1931 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1932 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);CHKERRQ(ierr); 1933 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 1934 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 1935 1936 ierr = MatComputeOperator(mf,MATAIJ,&Dmf);CHKERRQ(ierr); 1937 ierr = MatComputeOperator(mat,MATAIJ,&Dmat);CHKERRQ(ierr); 1938 1939 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 1940 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 1941 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 1942 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 1943 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 1944 flag = PETSC_FALSE; 1945 if (v) { 1946 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); 1947 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 1948 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 1949 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 1950 } 1951 } else if (v) { 1952 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 1953 } 1954 if (flg) *flg = flag; 1955 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 1956 ierr = MatDestroy(&mf);CHKERRQ(ierr); 1957 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 1958 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 1959 PetscFunctionReturn(0); 1960 } 1961 1962 /*@C 1963 MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function. 1964 1965 Logically Collective on Mat 1966 1967 Input Parameters: 1968 + mat - the shell matrix 1969 . f - the function 1970 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 1971 - ctx - an optional context for the function 1972 1973 Output Parameter: 1974 . flg - PETSC_TRUE if the multiply is likely correct 1975 1976 Options Database: 1977 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1978 1979 Level: advanced 1980 1981 Fortran Notes: 1982 Not supported from Fortran 1983 1984 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult() 1985 @*/ 1986 PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 1987 { 1988 PetscErrorCode ierr; 1989 Vec x,y,z; 1990 PetscInt m,n,M,N; 1991 Mat mf,Dmf,Dmat,Ddiff; 1992 PetscReal Diffnorm,Dmfnorm; 1993 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 1994 1995 PetscFunctionBegin; 1996 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1997 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);CHKERRQ(ierr); 1998 ierr = MatCreateVecs(mat,&x,&y);CHKERRQ(ierr); 1999 ierr = VecDuplicate(y,&z);CHKERRQ(ierr); 2000 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 2001 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 2002 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);CHKERRQ(ierr); 2003 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 2004 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 2005 ierr = MatComputeOperator(mf,MATAIJ,&Dmf);CHKERRQ(ierr); 2006 ierr = MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);CHKERRQ(ierr); 2007 ierr = MatComputeOperatorTranspose(mat,MATAIJ,&Dmat);CHKERRQ(ierr); 2008 2009 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 2010 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 2011 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 2012 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 2013 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 2014 flag = PETSC_FALSE; 2015 if (v) { 2016 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); 2017 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 2018 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 2019 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 2020 } 2021 } else if (v) { 2022 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 2023 } 2024 if (flg) *flg = flag; 2025 ierr = MatDestroy(&mf);CHKERRQ(ierr); 2026 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 2027 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 2028 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 2029 ierr = VecDestroy(&x);CHKERRQ(ierr); 2030 ierr = VecDestroy(&y);CHKERRQ(ierr); 2031 ierr = VecDestroy(&z);CHKERRQ(ierr); 2032 PetscFunctionReturn(0); 2033 } 2034 2035 /*@C 2036 MatShellSetOperation - Allows user to set a matrix operation for a shell matrix. 2037 2038 Logically Collective on Mat 2039 2040 Input Parameters: 2041 + mat - the shell matrix 2042 . op - the name of the operation 2043 - g - the function that provides the operation. 2044 2045 Level: advanced 2046 2047 Usage: 2048 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 2049 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 2050 $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 2051 2052 Notes: 2053 See the file include/petscmat.h for a complete list of matrix 2054 operations, which all have the form MATOP_<OPERATION>, where 2055 <OPERATION> is the name (in all capital letters) of the 2056 user interface routine (e.g., MatMult() -> MATOP_MULT). 2057 2058 All user-provided functions (except for MATOP_DESTROY) should have the same calling 2059 sequence as the usual matrix interface routines, since they 2060 are intended to be accessed via the usual matrix interface 2061 routines, e.g., 2062 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 2063 2064 In particular each function MUST return an error code of 0 on success and 2065 nonzero on failure. 2066 2067 Within each user-defined routine, the user should call 2068 MatShellGetContext() to obtain the user-defined context that was 2069 set by MatCreateShell(). 2070 2071 Use MatSetOperation() to set an operation for any matrix type. For matrix product operations (i.e. MatMat, MatTransposeMat etc) use MatShellSetMatProductOperation() 2072 2073 Fortran Notes: 2074 For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not 2075 generate a matrix. See src/mat/tests/ex120f.F 2076 2077 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts(), MatShellSetMatProductOperation() 2078 @*/ 2079 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void)) 2080 { 2081 PetscErrorCode ierr; 2082 2083 PetscFunctionBegin; 2084 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2085 ierr = PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));CHKERRQ(ierr); 2086 PetscFunctionReturn(0); 2087 } 2088 2089 /*@C 2090 MatShellGetOperation - Gets a matrix function for a shell matrix. 2091 2092 Not Collective 2093 2094 Input Parameters: 2095 + mat - the shell matrix 2096 - op - the name of the operation 2097 2098 Output Parameter: 2099 . g - the function that provides the operation. 2100 2101 Level: advanced 2102 2103 Notes: 2104 See the file include/petscmat.h for a complete list of matrix 2105 operations, which all have the form MATOP_<OPERATION>, where 2106 <OPERATION> is the name (in all capital letters) of the 2107 user interface routine (e.g., MatMult() -> MATOP_MULT). 2108 2109 All user-provided functions have the same calling 2110 sequence as the usual matrix interface routines, since they 2111 are intended to be accessed via the usual matrix interface 2112 routines, e.g., 2113 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 2114 2115 Within each user-defined routine, the user should call 2116 MatShellGetContext() to obtain the user-defined context that was 2117 set by MatCreateShell(). 2118 2119 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext() 2120 @*/ 2121 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void)) 2122 { 2123 PetscErrorCode ierr; 2124 2125 PetscFunctionBegin; 2126 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2127 ierr = PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g));CHKERRQ(ierr); 2128 PetscFunctionReturn(0); 2129 } 2130 2131 /*@ 2132 MatIsShell - Inquires if a matrix is derived from MATSHELL 2133 2134 Input Parameter: 2135 . mat - the matrix 2136 2137 Output Parameter: 2138 . flg - the boolean value 2139 2140 Level: developer 2141 2142 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) 2143 2144 .seealso: MatCreateShell() 2145 @*/ 2146 PetscErrorCode MatIsShell(Mat mat, PetscBool *flg) 2147 { 2148 PetscFunctionBegin; 2149 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2150 PetscValidPointer(flg,2); 2151 *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell); 2152 PetscFunctionReturn(0); 2153 } 2154