181e6777dSBarry Smith 281e6777dSBarry Smith #ifndef lint 3*b9fa9cd0SBarry Smith static char vcid[] = "$Id: snesmfj.c,v 1.12 1995/06/29 23:54:14 bsmith Exp bsmith $"; 481e6777dSBarry Smith #endif 581e6777dSBarry Smith 6b1f0a012SBarry Smith #include "draw.h" /*I "draw.h" I*/ 7b1f0a012SBarry Smith #include "snes.h" /*I "snes.h" I*/ 8*b9fa9cd0SBarry Smith #include "ptscimpl.h" 9*b9fa9cd0SBarry Smith #include "plog.h" 1081e6777dSBarry Smith 1139e2f89bSBarry Smith typedef struct { 1239e2f89bSBarry Smith SNES snes; 1339e2f89bSBarry Smith Vec w; 1439e2f89bSBarry Smith } MFCtx_Private; 1539e2f89bSBarry Smith 16*b9fa9cd0SBarry Smith int SNESMatrixFreeDestroy_Private(void *ptr) 17*b9fa9cd0SBarry Smith { 18*b9fa9cd0SBarry Smith int ierr; 19*b9fa9cd0SBarry Smith MFCtx_Private *ctx = (MFCtx_Private* ) ptr; 20*b9fa9cd0SBarry Smith ierr = VecDestroy(ctx->w); CHKERRQ(ierr); 21*b9fa9cd0SBarry Smith PETSCFREE(ptr); 22*b9fa9cd0SBarry Smith return 0; 23*b9fa9cd0SBarry Smith } 2439e2f89bSBarry Smith /* 2539e2f89bSBarry Smith SNESMatrixFreeMult_Private - Default matrix free form of A*u. 2639e2f89bSBarry Smith 2739e2f89bSBarry Smith */ 2839e2f89bSBarry Smith int SNESMatrixFreeMult_Private(void *ptr,Vec dx,Vec y) 2939e2f89bSBarry Smith { 3039e2f89bSBarry Smith MFCtx_Private *ctx = (MFCtx_Private* ) ptr; 3139e2f89bSBarry Smith SNES snes = ctx->snes; 3239e2f89bSBarry Smith double norm,epsilon = 1.e-8; /* assumes double precision */ 3319a167f6SBarry Smith Scalar h,dot; 3419a167f6SBarry Smith double sum; 3539e2f89bSBarry Smith Scalar mone = -1.0; 3639e2f89bSBarry Smith Vec w = ctx->w,U,F; 3739e2f89bSBarry Smith int ierr; 3839e2f89bSBarry Smith 3978b31e54SBarry Smith ierr = SNESGetSolution(snes,&U); CHKERRQ(ierr); 4078b31e54SBarry Smith ierr = SNESGetFunction(snes,&F); CHKERRQ(ierr); 4139e2f89bSBarry Smith /* determine a "good" step size */ 4239e2f89bSBarry Smith VecDot(U,dx,&dot); VecASum(dx,&sum); VecNorm(dx,&norm); 43edd2f0e1SBarry Smith if (sum == 0.0) {dot = 1.0; norm = 1.0;} 4419a167f6SBarry Smith #if defined(PETSC_COMPLEX) 4519a167f6SBarry Smith else if (abs(dot) < 1.e-16*sum && real(dot) >= 0.0) dot = 1.e-16*sum; 4619a167f6SBarry Smith else if (abs(dot) < 0.0 && real(dot) > 1.e-16*sum) dot = -1.e-16*sum; 4719a167f6SBarry Smith #else 48edd2f0e1SBarry Smith else if (dot < 1.e-16*sum && dot >= 0.0) dot = 1.e-16*sum; 4939e2f89bSBarry Smith else if (dot < 0.0 && dot > 1.e-16*sum) dot = -1.e-16*sum; 5019a167f6SBarry Smith #endif 5139e2f89bSBarry Smith h = epsilon*dot/(norm*norm); 5239e2f89bSBarry Smith 5339e2f89bSBarry Smith /* evaluate function at F(x + dx) */ 5439e2f89bSBarry Smith VecWAXPY(&h,dx,U,w); 5578b31e54SBarry Smith ierr = SNESComputeFunction(snes,w,y); CHKERRQ(ierr); 5639e2f89bSBarry Smith VecAXPY(&mone,F,y); 5739e2f89bSBarry Smith h = -1.0/h; 5839e2f89bSBarry Smith VecScale(&h,y); 5939e2f89bSBarry Smith return 0; 6039e2f89bSBarry Smith } 6181e6777dSBarry Smith /*@ 625392566eSBarry Smith SNESDefaultMatrixFreeMatCreate - Creates a matrix-free matrix 635392566eSBarry Smith for use with SNES solver. You may use this matrix as 645392566eSBarry Smith Jacobian argument for the routine SNESSetJacobian. This is 655392566eSBarry Smith most useful when you are using finite differences for a 665392566eSBarry Smith matrix free Newton method but explictly are forming a 675392566eSBarry Smith preconditioner matrix. 685392566eSBarry Smith 695392566eSBarry Smith Input Parameters: 705392566eSBarry Smith . x - vector where SNES solution is to be stored. 715392566eSBarry Smith 725392566eSBarry Smith Output Parameters: 735392566eSBarry Smith . J - the matrix-free matrix 745392566eSBarry Smith 755392566eSBarry Smith @*/ 765392566eSBarry Smith int SNESDefaultMatrixFreeMatCreate(SNES snes,Vec x, Mat *J) 775392566eSBarry Smith { 785392566eSBarry Smith MPI_Comm comm; 795392566eSBarry Smith MFCtx_Private *mfctx; 805392566eSBarry Smith int n,ierr; 815392566eSBarry Smith 825392566eSBarry Smith mfctx = (MFCtx_Private *) PETSCMALLOC(sizeof(MFCtx_Private));CHKPTRQ(mfctx); 835392566eSBarry Smith mfctx->snes = snes; 845392566eSBarry Smith ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 855392566eSBarry Smith PetscObjectGetComm((PetscObject)x,&comm); 865392566eSBarry Smith VecGetSize(x,&n); 875392566eSBarry Smith ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERRQ(ierr); 885392566eSBarry Smith MatShellSetMult(*J,SNESMatrixFreeMult_Private); 89*b9fa9cd0SBarry Smith MatShellSetDestroy(*J,SNESMatrixFreeDestroy_Private); 90*b9fa9cd0SBarry Smith PLogObjectParent(*J,mfctx->w); 9181e6777dSBarry Smith return 0; 9281e6777dSBarry Smith } 9381e6777dSBarry Smith 94