Actual source code: ex25.c

  1: static char help[] = "Tests CG, MINRES and SYMMLQ on the symmetric indefinite matrices: afiro \n\n";

  3: #include <petscksp.h>

  5: int main(int argc,char **args)
  6: {
  7:   Mat            C;
  8:   PetscScalar    none = -1.0;
  9:   PetscMPIInt    rank,size;
 10:   PetscInt       its,k;
 11:   PetscReal      err_norm,res_norm;
 12:   Vec            x,b,u,u_tmp;
 13:   PC             pc;
 14:   KSP            ksp;
 15:   PetscViewer    view;
 16:   char           filein[128];     /* input file name */

 18:   PetscInitialize(&argc,&args,(char*)0,help);
 19:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 20:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 22:   /* Load the binary data file "filein". Set runtime option: -f filein */
 23:   PetscPrintf(PETSC_COMM_WORLD,"\n Load dataset ...\n");
 24:   PetscOptionsGetString(NULL,NULL,"-f",filein,sizeof(filein),NULL);
 25:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filein,FILE_MODE_READ,&view);
 26:   MatCreate(PETSC_COMM_WORLD,&C);
 27:   MatSetType(C,MATMPISBAIJ);
 28:   MatLoad(C,view);
 29:   VecCreate(PETSC_COMM_WORLD,&b);
 30:   VecCreate(PETSC_COMM_WORLD,&u);
 31:   VecLoad(b,view);
 32:   VecLoad(u,view);
 33:   PetscViewerDestroy(&view);
 34:   /* VecView(b,VIEWER_STDOUT_WORLD); */
 35:   /* MatView(C,VIEWER_STDOUT_WORLD); */

 37:   VecDuplicate(u,&u_tmp);

 39:   /* Check accuracy of the data */
 40:   /*
 41:   MatMult(C,u,u_tmp);
 42:   VecAXPY(u_tmp,none,b);
 43:   VecNorm(u_tmp,NORM_2,&res_norm);
 44:   PetscPrintf(PETSC_COMM_WORLD,"Accuracy of the loading data: | b - A*u |_2 : %g \n",(double)res_norm);
 45:   */

 47:   /* Setup and solve for system */
 48:   VecDuplicate(b,&x);
 49:   for (k=0; k<3; k++) {
 50:     if (k == 0) {                              /* CG  */
 51:       KSPCreate(PETSC_COMM_WORLD,&ksp);
 52:       KSPSetType(ksp,KSPCG);
 53:       KSPSetOperators(ksp,C,C);
 54:       PetscPrintf(PETSC_COMM_WORLD,"\n CG: \n");
 55:     } else if (k == 1) {                       /* MINRES */
 56:       KSPCreate(PETSC_COMM_WORLD,&ksp);
 57:       KSPSetType(ksp,KSPMINRES);
 58:       KSPSetOperators(ksp,C,C);
 59:       PetscPrintf(PETSC_COMM_WORLD,"\n MINRES: \n");
 60:     } else {                                 /* SYMMLQ */
 61:       KSPCreate(PETSC_COMM_WORLD,&ksp);
 62:       KSPSetOperators(ksp,C,C);
 63:       KSPSetType(ksp,KSPSYMMLQ);
 64:       PetscPrintf(PETSC_COMM_WORLD,"\n SYMMLQ: \n");
 65:     }

 67:     KSPGetPC(ksp,&pc);
 68:     PCSetType(pc,PCNONE);

 70:     /*
 71:     Set runtime options, e.g.,
 72:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
 73:                          -pc_type jacobi -pc_jacobi_type rowmax
 74:     These options will override those specified above as long as
 75:     KSPSetFromOptions() is called _after_ any other customization routines.
 76:     */
 77:     KSPSetFromOptions(ksp);

 79:     /* Solve linear system; */
 80:     KSPSolve(ksp,b,x);
 81:     KSPGetIterationNumber(ksp,&its);

 83:     /* Check error */
 84:     VecCopy(u,u_tmp);
 85:     VecAXPY(u_tmp,none,x);
 86:     VecNorm(u_tmp,NORM_2,&err_norm);
 87:     MatMult(C,x,u_tmp);
 88:     VecAXPY(u_tmp,none,b);
 89:     VecNorm(u_tmp,NORM_2,&res_norm);

 91:     PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
 92:     PetscPrintf(PETSC_COMM_WORLD,"Residual norm: %g;",(double)res_norm);
 93:     PetscPrintf(PETSC_COMM_WORLD,"  Error norm: %g.\n",(double)err_norm);

 95:     KSPDestroy(&ksp);
 96:   }

 98:   /*
 99:        Free work space.  All PETSc objects should be destroyed when they
100:        are no longer needed.
101:   */
102:   VecDestroy(&b);
103:   VecDestroy(&u);
104:   VecDestroy(&x);
105:   VecDestroy(&u_tmp);
106:   MatDestroy(&C);

108:   PetscFinalize();
109:   return 0;
110: }

112: /*TEST

114:     test:
115:       args: -f ${DATAFILESPATH}/matrices/indefinite/afiro -ksp_rtol 1.e-3
116:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)

118: TEST*/