Actual source code: ex53.c

petsc-3.3-p7 2013-05-11
  1: static const char help[] = "Read linear variational inequality from file and solve it.\n\n";
  2: #include <petscsnes.h>

  4: typedef struct {
  5:   Vec         q,zz,lb,ub;
  6:   Mat         M,Jac;
  7: } AppCtx;

  9: /* 
 10:    User-defined routines
 11: */
 12: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 13: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);

 17: int main(int argc,char **argv)
 18: {
 19:   SNES           snes;         /* nonlinear solver context */
 20:   Vec            r;          /* solution, residual vectors */
 22:   AppCtx         user;         /* user-defined work context */
 23:   PetscViewer    viewer;

 25:   PetscInitialize(&argc,&argv,(char *)0,help);
 26:   PetscViewerBinaryOpen(PETSC_COMM_SELF,"videfinition",FILE_MODE_READ,&viewer);
 27:   MatCreate(PETSC_COMM_SELF,&user.M); MatLoad(user.M,viewer);
 28:   MatDuplicate(user.M,MAT_COPY_VALUES,&user.Jac);
 29:   VecCreate(PETSC_COMM_SELF,&user.q); VecLoad(user.q,viewer);
 30:   VecCreate(PETSC_COMM_SELF,&user.lb); VecLoad(user.lb,viewer);
 31:   VecCreate(PETSC_COMM_SELF,&user.ub);VecLoad(user.ub,viewer);
 32:   VecCreate(PETSC_COMM_SELF,&user.zz);VecLoad(user.zz,viewer);
 33:   VecView(user.zz,PETSC_VIEWER_STDOUT_SELF);
 34:   /*  VecSet(user.zz,0.01);*/
 35:   PetscViewerDestroy(&viewer);
 36:   VecDuplicate(user.q,&r);

 38:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 39:      Create nonlinear solver context
 40:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 41:   SNESCreate(PETSC_COMM_WORLD,&snes);

 43:   SNESSetFunction(snes,r,FormFunction1,&user);

 45:   SNESSetJacobian(snes,user.Jac,user.Jac,FormJacobian1,&user);

 47:   SNESVISetVariableBounds(snes,user.lb,user.ub);
 48:   SNESSetFromOptions(snes);
 49:   SNESSolve(snes,PETSC_NULL,user.zz);
 50:   PetscObjectSetName((PetscObject)user.zz,"x*");
 51:   VecView(user.zz,PETSC_VIEWER_STDOUT_SELF);
 52:   PetscObjectSetName((PetscObject)r,"f(x*)");
 53:   FormFunction1(snes,user.zz,r,&user);
 54:   VecView(r,PETSC_VIEWER_STDOUT_SELF);
 55:   PetscFinalize();
 56:   return 0;
 57: }

 59: /* ------------------------------------------------------------------- */
 62: /* 
 63:    FormFunction1 - Evaluates nonlinear function, F(x).

 65:    Input Parameters:
 66: .  snes - the SNES context
 67: .  x    - input vector
 68: .  ctx  - optional user-defined context

 70:    Output Parameter:
 71: .  f - function vector
 72:  */
 73: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
 74: {
 76:   AppCtx         *user = (AppCtx*)ctx;

 78:   MatMult(user->M,x,f);
 79:   VecAXPY(f,1.0,user->q);
 80:   return 0;
 81: }
 82: /* ------------------------------------------------------------------- */
 85: /*
 86:    FormJacobian1 - Evaluates Jacobian matrix.

 88:    Input Parameters:
 89: .  snes - the SNES context
 90: .  x - input vector
 91: .  ctx - optional user-defined context

 93:    Output Parameters:
 94: .  jac - Jacobian matrix
 95: .  B - optionally different preconditioning matrix
 96: .  flag - flag indicating matrix structure
 97: */
 98: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *ctx)
 99: {
101:   AppCtx         *user = (AppCtx*)ctx;
102:   MatCopy(user->M,*jac,DIFFERENT_NONZERO_PATTERN);
103:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
104:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);

106:   return 0;
107: }