update ex15.c

[GMRES2stage.git] / code / ex15.c

diff --git a/code/ex15.c b/code/ex15.c
index 820ec8d612de26711b0c26cdbe8bf93e0ae33356..04f08f10065869daaa6a6c754d5eb47487bd3fc8 100644 (file)
--- a/code/ex15.c
+++ b/code/ex15.c
@@ -1,3 +1,5 @@
+// /home/couturie/work/petsc-3.5.1_old/arch-linux2-c-debug/bin/mpirun -np 4 -machinefile archi   ./ex15 -m 1000 -n 1000 -ksp_type fgmres -pc_type mg
+
 
 // /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 4    ./ex15 -m 400 -n 400 -ksp_type fgmres
 
@@ -56,7 +58,7 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
   //Variables
 
   PetscScalar  gamma, alpha, oldgamma, beta;
-  PetscReal norm=20, Eprecision=1e-6, cgprec=1e-40;     
+  PetscReal norm=20, Eprecision=1e-3, cgprec=1e-40;     
   PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
   PetscErrorCode ierr;
   PetscScalar T1, T2;
@@ -234,7 +236,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
   //Variables
 
   PetscScalar  alpha, beta;
-  PetscReal norm=20, Eprecision=1e-6, tol=1e-40;     
+  PetscReal norm=20, Eprecision=1e-3, tol=1e-40;     
   PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;
   PetscErrorCode ierr;
   PetscScalar T1, T2;
@@ -575,7 +577,7 @@ int main(int argc,char **args)
        to set various options.
   */
   ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp,1e-7,1e-7,PETSC_DEFAULT,5000000);CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp,4e-6,4e-6,PETSC_DEFAULT,500000);CHKERRQ(ierr);
 
 
   //  PC             pc;
@@ -599,33 +601,64 @@ int main(int argc,char **args)
   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                       Solve the linear system
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-  PetscScalar T1,T2;
-       T1 = MPI_Wtime();
-        ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
-       T2 = MPI_Wtime();
-  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-                      Check solution and clean up
-     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
-  /*
-     Check the error
-  */
-  Vec sol;
-  VecDuplicate(b,&sol);
-  MatMult(A,x,sol);
-  VecAXPY(sol,-1,b);
-  VecNorm(sol, NORM_2, &norm);
-  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
-  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr);
-  ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n\n\n", T2-T1); CHKERRQ(ierr);
-
-
+  /* PetscScalar T1,T2; */
+  /*       T1 = MPI_Wtime(); */
+  /*       ierr = KSPSetUp(ksp); CHKERRQ(ierr); */
+  /*       ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); */
+  /*       T2 = MPI_Wtime(); */
+  /* /\* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+  /*                     Check solution and clean up */
+  /*    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *\/ */
+
+  /* /\* */
+  /*    Check the error */
+  /* *\/ */
+  /* Vec sol; */
+  /* VecDuplicate(b,&sol); */
+  /* MatMult(A,x,sol); */
+  /* VecAXPY(sol,-1,b); */
+  /* VecNorm(sol, NORM_2, &norm); */
+  /* ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); */
+  /* ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr); */
+  /* ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n\n\n", T2-T1); CHKERRQ(ierr); */
+
+
+
+
+  //version to control the error
+ {
 
+    Vec x2;
+    Vec sol;
+    VecDuplicate(b,&x2);
+    VecDuplicate(b,&sol);
+    
+    PetscScalar norm=100;
+    PetscScalar T1,T2;
+    PetscInt total,its;
+    ierr = KSPSetTolerances(ksp,1e-10,1e-10,PETSC_DEFAULT,30);CHKERRQ(ierr);
+    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
+    T1 = MPI_Wtime();
+    while(norm>1e-3) {
+      ierr = KSPSolve(ksp,b,x2);CHKERRQ(ierr);
+      KSPGetResidualNorm(ksp,&norm);
+      ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
+      total += its;
+       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g\n", norm); CHKERRQ(ierr);
+    }
 
+    T2 = MPI_Wtime();
 
+    MatMult(A,x2,sol);
+    VecAXPY(sol,-1,b);
+    VecNorm(sol, NORM_2, &norm);
+     ierr = PetscPrintf(PETSC_COMM_WORLD,"Computed norm of error %g iterations %D\n",(double)norm,total);CHKERRQ(ierr);
+    ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time NORMAL GMRES : %g (s)\n\n\n", T2-T1); CHKERRQ(ierr);
 
+  }
 
 
+ 
  {
 
     Vec x2;
@@ -662,9 +695,8 @@ int main(int argc,char **args)
     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Error Krylov Minimization LSQR %g\n",norm);
 
   }
-
+ 
   
-
   /*
      Free work space.  All PETSc objects should be destroyed when they
      are no longer needed.