Actual source code: ex25.c

  1: /*$Id: ex25.c,v 1.3 2000/11/15 22:56:05 balay Exp $*/

  3: static char help[] =
  4: "Tests CG, MINRES and SYMMLQ on the symmetric indefinite matrices: afiro and golan\n\
  5: Runtime options: ex25 -fload ~petsc/matrices/indefinite/afiro -pc_type jacobi -pc_jacobi_rowmax\n\
  6: See ~petsc/matrices/indefinite/readme \n\n";

  8: #include <petscksp.h>

 12: int main(int argc,char **args)
 13: {
 14:   Mat         C;
 15:   PetscScalar      v,none = -1.0;
 16:   int         i,j,ierr,Istart,Iend,N,rank,size,its,k;
 17:   double      err_norm,res_norm;
 18:   Vec         x,b,u,u_tmp;
 19:   PetscRandom r;
 20:   PC          pc;
 21:   KSP         ksp;
 22:   PetscViewer      view;
 23:   char        filein[128];     /* input file name */

 25:   PetscInitialize(&argc,&args,(char *)0,help);
 26:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 27:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 28: 
 29:   /* Load the binary data file "filein". Set runtime option: -fload filein */
 30:   PetscPrintf(PETSC_COMM_WORLD,"\n Load dataset ...\n");
 31:   PetscOptionsGetString(PETSC_NULL,"-fload",filein,128,PETSC_NULL);
 32:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filein,FILE_MODE_READ,&view);
 33:   MatCreate(PETSC_COMM_WORLD,&C);
 34:   MatSetType(C,MATMPISBAIJ);
 35:   MatLoad(C,view);
 36:   VecCreate(PETSC_COMM_WORLD,&b);
 37:   VecCreate(PETSC_COMM_WORLD,&u);
 38:   VecLoad(b,view);
 39:   VecLoad(u,view);
 40:   PetscViewerDestroy(&view);
 41:   /* VecView(b,VIEWER_STDOUT_WORLD); */
 42:   /* MatView(C,VIEWER_STDOUT_WORLD); */

 44:   VecDuplicate(u,&u_tmp);

 46:   /* Check accuracy of the data */
 47:   /*
 48:   MatMult(C,u,u_tmp);
 49:   VecAXPY(u_tmp,none,b);
 50:   VecNorm(u_tmp,NORM_2,&res_norm);
 51:   PetscPrintf(PETSC_COMM_WORLD,"Accuracy of the loading data: | b - A*u |_2 : %A \n",res_norm); 
 52:   */

 54:   /* Setup and solve for system */
 55:   VecDuplicate(b,&x);
 56:   for (k=0; k<3; k++){
 57:     if (k == 0){                              /* CG  */
 58:       KSPCreate(PETSC_COMM_WORLD,&ksp);
 59:       KSPSetType(ksp,KSPCG);
 60:       KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);
 61:       PetscPrintf(PETSC_COMM_WORLD,"\n CG: \n");
 62:     } else if (k == 1){                       /* MINRES */
 63:       KSPCreate(PETSC_COMM_WORLD,&ksp);
 64:       KSPSetType(ksp,KSPMINRES);
 65:       KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);
 66:       PetscPrintf(PETSC_COMM_WORLD,"\n MINRES: \n");
 67:     } else {                                 /* SYMMLQ */
 68:       KSPCreate(PETSC_COMM_WORLD,&ksp);
 69:       KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);
 70:       KSPSetType(ksp,KSPSYMMLQ);
 71:       PetscPrintf(PETSC_COMM_WORLD,"\n SYMMLQ: \n");
 72:     }

 74:     KSPGetPC(ksp,&pc);
 75:     PCSetType(pc,PCNONE);
 76:     /* PCSetType(pc,PCJACOBI); */
 77:     KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);

 79:     /*
 80:     Set runtime options, e.g.,
 81:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
 82:                          -pc_type jacobi -pc_jacobi_rowmax
 83:     These options will override those specified above as long as
 84:     KSPSetFromOptions() is called _after_ any other customization routines.
 85:     */
 86:     KSPSetFromOptions(ksp);

 88:     /* Solve linear system; */
 89:     KSPSolve(ksp,b,x);
 90:     KSPGetIterationNumber(ksp,&its);
 91: 
 92:   /* Check error */
 93:     VecCopy(u,u_tmp);
 94:     VecAXPY(u_tmp,none,x);
 95:     VecNorm(u_tmp,NORM_2,&err_norm);
 96:     MatMult(C,x,u_tmp);
 97:     VecAXPY(u_tmp,none,b);
 98:     VecNorm(u_tmp,NORM_2,&res_norm);
 99: 
100:     PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
101:     PetscPrintf(PETSC_COMM_WORLD,"Residual norm: %A;",res_norm);
102:     PetscPrintf(PETSC_COMM_WORLD,"  Error norm: %A.\n",err_norm);

104:     KSPDestroy(&ksp);
105:   }
106: 
107:   /* 
108:        Free work space.  All PETSc objects should be destroyed when they
109:        are no longer needed.
110:   */
111:   VecDestroy(&b);
112:   VecDestroy(&u);
113:   VecDestroy(&x);
114:   VecDestroy(&u_tmp);
115:   MatDestroy(&C);

117:   PetscFinalize();
118:   return 0;
119: }