Actual source code: ex25.c

petsc-3.12.2 2019-11-22
Report Typos and Errors
  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;
 11:   PetscInt       its,k;
 12:   PetscReal      err_norm,res_norm;
 13:   Vec            x,b,u,u_tmp;
 14:   PC             pc;
 15:   KSP            ksp;
 16:   PetscViewer    view;
 17:   char           filein[128];     /* input file name */

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

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

 38:   VecDuplicate(u,&u_tmp);

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

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

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

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

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

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

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

 96:     KSPDestroy(&ksp);
 97:   }

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

109:   PetscFinalize();
110:   return ierr;
111: }

113: /*TEST

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

119: TEST*/