archi

[GMRES2stage.git] / code / ex15.c

diff --git a/code/ex15.c b/code/ex15.c
index 0e0a8c467b785a3ebbaae3f4135057fa2abf1390..04f08f10065869daaa6a6c754d5eb47487bd3fc8 100644 (file)
--- a/code/ex15.c
+++ b/code/ex15.c
@@ -1,5 +1,7 @@
+// /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
+
+// /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-8, 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;
@@ -120,13 +122,14 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, PETSC_DEFAULT, 16); CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 
 
 
   //GMRES WITH MINIMIZATION
   T1 = MPI_Wtime();
+  ierr = KSPSetUp(ksp); CHKERRQ(ierr);
   while(giter<Emaxiter && norm>Eprecision ){
     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 
@@ -146,11 +149,13 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
 
 
-      //Error
+      KSPGetResidualNorm(ksp,&norm);
+
+      /*      //Error
       ierr = VecCopy(x, residu); CHKERRQ(ierr);
       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
-
+       */
 
 
       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
@@ -231,8 +236,8 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
   //Variables
 
   PetscScalar  alpha, beta;
-  PetscReal norm=20, Eprecision=1e-8, tol=1e-40;     
-  PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
+  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;
   KSP ksp;
@@ -308,7 +313,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, PETSC_DEFAULT, 16); CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 
 
@@ -316,6 +321,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   //GMRES WITH MINIMIZATION
   T1 = MPI_Wtime();
+  ierr = KSPSetUp(ksp); CHKERRQ(ierr);
   while(giter<Emaxiter && norm>Eprecision ){
     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 
@@ -335,11 +341,14 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
 
 
+      KSPGetResidualNorm(ksp,&norm);
+
+      /*
       //Error
       ierr = VecCopy(x, residu); CHKERRQ(ierr);
       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
-
+       */
 
 
       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
@@ -419,14 +428,14 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
         VecAYPX(d,-thet,v);       //d = (v - thet * d);
         VecAYPX(d,1/rho,zero_short);     //d=d/ rho;
 
-
+        /*
         if (normar/(norma*normr) <= tol) { // check for convergence in min{|b-A*x|}
           break;
         }
         if (normr <= tolb) {           // check for convergence in A*x=b
           break;
         }
-
+         */
 
         VecAXPY(x_lsqr,phi,d);     // x_lsqr=x_lsqr+phi*d
         normr = abs(s) * normr;
@@ -475,6 +484,10 @@ int main(int argc,char **args)
   ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
 
+  PetscMPIInt size;
+  MPI_Comm_size(PETSC_COMM_WORLD,&size);
+  PetscPrintf(PETSC_COMM_WORLD,"Number of processors = %d\n",size);
+
   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
          Compute the matrix and right-hand-side vector that define
          the linear system, Ax = b.
@@ -506,7 +519,8 @@ int main(int argc,char **args)
    */
   for (Ii=Istart; Ii<Iend; Ii++) {
     v = -1.0; i = Ii/n; j = Ii - i*n;
-    if (i>0)   {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
+    PetscScalar v2=-1.;
+    if (i>0)   {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v2,INSERT_VALUES);CHKERRQ(ierr);}
     if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
     if (j>0)   {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
     if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
@@ -563,7 +577,17 @@ int main(int argc,char **args)
        to set various options.
   */
   ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp,1e-9,1e-9,PETSC_DEFAULT,5000000);CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp,4e-6,4e-6,PETSC_DEFAULT,500000);CHKERRQ(ierr);
+
+
+  //  PC             pc;
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+  KSPGetPC(ksp, &pc);
+  PCType         type;
+  PCGetType(pc, &type);
+
+  PetscPrintf(PETSC_COMM_WORLD, "PC TYPE %s  \n", type);
+
 
   /*
     Set runtime options, e.g.,
@@ -577,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;
@@ -630,7 +685,7 @@ int main(int argc,char **args)
     VecDuplicate(b,&x2);
     VecDuplicate(b,&sol);
     
-    KrylovMinimize(A, b, x2);
+    KrylovMinimizeLSQR(A, b, x2);
 
 
 
@@ -640,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.