modif de LSQR qui fonctionne quand meme mieux (un poil moins rapide que CGLS quand...

diff --git a/code/ex15.c b/code/ex15.c
index 0e0a8c467b785a3ebbaae3f4135057fa2abf1390..94c93659b6acb2e0ca63c874f44ed14e79f649c4 100644 (file)
--- a/code/ex15.c
+++ b/code/ex15.c
@@ -1,5 +1,5 @@
 
 

-// /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

 
 
 
 
 
 


@@ -232,7 +232,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   PetscScalar  alpha, beta;
   PetscReal norm=20, Eprecision=1e-8, tol=1e-40;     
 
   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;
+  PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;

   PetscErrorCode ierr;
   PetscScalar T1, T2;
   KSP ksp;
   PetscErrorCode ierr;
   PetscScalar T1, T2;
   KSP ksp;


@@ -419,14 +419,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;
 
         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;
         }
         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;
 
         VecAXPY(x_lsqr,phi,d);     // x_lsqr=x_lsqr+phi*d
         normr = abs(s) * normr;


@@ -630,7 +630,7 @@ int main(int argc,char **args)
     VecDuplicate(b,&x2);
     VecDuplicate(b,&sol);
     
     VecDuplicate(b,&x2);
     VecDuplicate(b,&sol);
     

-    KrylovMinimize(A, b, x2);
+    KrylovMinimizeLSQR(A, b, x2);

 
 
 
 
 
 

diff --git a/code/ex29.c b/code/ex29.c
index 5e21c65b386f6c654ad382198ecf76891fa9250f..ede23559af1f7c1b3857d7b638259ae9ee7f9878 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 900 -da_grid_y 900 -ksp_type fgmres

 
 
 
 
 
 


@@ -240,7 +240,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   PetscScalar  alpha, beta;
   PetscReal norm=20, Eprecision=1e-5, tol=1e-40;     
 
   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;
+  PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;

   PetscErrorCode ierr;
   PetscScalar T1, T2;
   KSP ksp;
   PetscErrorCode ierr;
   PetscScalar T1, T2;
   KSP ksp;


@@ -428,13 +428,6 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
         VecAYPX(d,1/rho,zero_short);     //d=d/ rho;
 
 
         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;
 
         VecAXPY(x_lsqr,phi,d);     // x_lsqr=x_lsqr+phi*d
         normr = abs(s) * normr;

diff --git a/code/ex45.c b/code/ex45.c
index a44a2d09174416a6ae6d2a08056a78d5770b4126..11fd5c28340202ee669d6fea574e220a1d8a771d 100644 (file)
--- a/code/ex45.c
+++ b/code/ex45.c
@@ -408,14 +408,6 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
         VecAYPX(d,1/rho,zero_short);     //d=d/ rho;
 
 
         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;
         MatMultTranspose(AS, u, vt); //vt=A'*u;
         VecAXPY(x_lsqr,phi,d);     // x_lsqr=x_lsqr+phi*d
         normr = abs(s) * normr;
         MatMultTranspose(AS, u, vt); //vt=A'*u;