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