elaskjls

[GMRES2stage.git] / code / ex29.c

diff --git a/code/ex29.c b/code/ex29.c
index 5e21c65b386f6c654ad382198ecf76891fa9250f..427f70a6af61ed6968bf45a2d8583e9f0e78b130 100644 (file)
--- a/code/ex29.c
+++ b/code/ex29.c
@@ -1,6 +1,6 @@
 
 
-//  /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 3 ex29 -da_grid_x 900 -da_grid_y 900
+//  /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 3 ex29 -da_grid_x 600 -da_grid_y 600 -ksp_type fgmres 
 
 
 
@@ -64,8 +64,8 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
   //Variables
 
   PetscScalar  gamma, alpha, oldgamma, beta;
-  PetscReal norm=20, Eprecision=1e-5, cgprec=1e-40;     
-  PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
+  PetscReal norm=20, Eprecision=1e-6, cgprec=1e-40;     
+  PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
   PetscErrorCode ierr;
   PetscScalar T1, T2;
   KSP ksp;
@@ -128,13 +128,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++){
 
@@ -154,11 +155,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);
@@ -239,8 +242,8 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
   //Variables
 
   PetscScalar  alpha, beta;
-  PetscReal norm=20, Eprecision=1e-5, tol=1e-40;     
-  PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
+  PetscReal norm=20, Eprecision=1e-6, tol=1e-40;     
+  PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
   PetscErrorCode ierr;
   PetscScalar T1, T2;
   KSP ksp;
@@ -316,7 +319,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);
 
 
@@ -324,6 +327,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++){
 
@@ -342,12 +346,14 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
       VecRestoreArray(x, &array);
 
 
-
+      
+      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);
@@ -428,13 +434,6 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
         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;
@@ -483,6 +482,11 @@ int main(int argc,char **argv)
 
   PetscInitialize(&argc,&argv,(char*)0,help);
 
+
+  PetscMPIInt size;
+  MPI_Comm_size(PETSC_COMM_WORLD,&size);
+  PetscPrintf(PETSC_COMM_WORLD,"Number of processors = %d\n",size);
+
   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
   ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);CHKERRQ(ierr);
   ierr = DMDASetUniformCoordinates(da,0,1,0,1,0,0);CHKERRQ(ierr);
@@ -502,10 +506,22 @@ int main(int argc,char **argv)
   ierr = KSPSetComputeOperators(ksp,ComputeMatrix,&user);CHKERRQ(ierr);
   ierr = KSPSetDM(ksp,da);CHKERRQ(ierr);
   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
-  ierr = KSPSetUp(ksp);CHKERRQ(ierr);
 
-  ierr = KSPSetTolerances(ksp, 1e-5, 1e-5, PETSC_DEFAULT, 50000000); 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);
+
+
+  ierr = KSPSetTolerances(ksp, 1e-6, 1e-6, PETSC_DEFAULT, 50000000); CHKERRQ(ierr);
+  
   T1 = MPI_Wtime();
+    ierr = KSPSetUp(ksp);CHKERRQ(ierr);
   ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);
   T2 = MPI_Wtime();