elaskjls

[GMRES2stage.git] / code / ex49.c

diff --git a/code/ex49.c b/code/ex49.c
index b93b5c637975fd5e524c97c0ba24475bf0bfd0ec..6e5331d9966195065efcd53aa6fb9ba3458c39a9 100644 (file)
--- a/code/ex49.c
+++ b/code/ex49.c
@@ -1,4 +1,7 @@
-//  /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 4    ./ex49 -mx 900 -my 900 -ksp_type fgmres
+//  /home/couturie/work/petsc-3.5.1_old/arch-linux2-c-debug/bin/mpirun -np 4  -machinefile archi  ./ex49 -mx 900 -my 900 -ksp_type fgmres -pc_type mg
+
+
+
 
 
 static char help[] =  "   Solves the compressible plane strain elasticity equations in 2d on the unit domain using Q1 finite elements. \n\
@@ -88,8 +91,8 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
   //Variables
 
   PetscScalar  gamma, alpha, oldgamma, beta;
-  PetscReal norm=20, Eprecision=1e-6, cgprec=1e-40;     
-  PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
+  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;
   KSP ksp;
@@ -97,7 +100,7 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
   PetscInt size;
   PetscInt Istart,Iend;
   PetscInt i,its;
-  Vec       x_old, residu;
+  Vec       residu;
   Mat S, AS;
   PetscScalar *array;
   PetscInt *ind_row;
@@ -141,7 +144,6 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
   ierr = VecDuplicate(b, &q); CHKERRQ(ierr);
   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
 
-  ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
 
 
@@ -189,7 +191,6 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
 
       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
-      ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 
 
     }
@@ -247,6 +248,13 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
   T2 = MPI_Wtime();
   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = VecDestroy(&r);CHKERRQ(ierr);  
+  ierr = VecDestroy(&vect);CHKERRQ(ierr);  
+  ierr = VecDestroy(&p);CHKERRQ(ierr);  
+  ierr = VecDestroy(&ss);CHKERRQ(ierr);  
+  ierr = VecDestroy(&Ax);CHKERRQ(ierr);  
+  ierr = VecDestroy(&residu);CHKERRQ(ierr);  
 
   return 0;
 
@@ -266,8 +274,8 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
   //Variables
 
   PetscScalar  alpha, beta;
-  PetscReal norm=20, Eprecision=1e-6, tol=1e-40;     
-  PetscInt giter=0, ColS=8, 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=15, Iiter=0;
   PetscErrorCode ierr;
   PetscScalar T1, T2;
   KSP ksp;
@@ -275,7 +283,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
   PetscInt size;
   PetscInt Istart,Iend;
   PetscInt i,its;
-  Vec       x_old, residu;
+  Vec       residu;
   Mat S, AS;
   PetscScalar *array;
   PetscInt *ind_row;
@@ -314,7 +322,6 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
 
-  ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
 
 
@@ -381,12 +388,12 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
 
       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
-      ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 
 
     }
 
 
+
     //minimization step
     if( norm>Eprecision) {
 
@@ -479,6 +486,16 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
   T2 = MPI_Wtime();
   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time LSQR : %g (s)\n", T2-T1); CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations LSQR : %D\n", total); CHKERRQ(ierr);
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = VecDestroy(&Ax);CHKERRQ(ierr);  
+  ierr = VecDestroy(&u);CHKERRQ(ierr);  
+  ierr = VecDestroy(&uu);CHKERRQ(ierr);  
+  ierr = VecDestroy(&zero_long);CHKERRQ(ierr);  
+  ierr = VecDestroy(&d);CHKERRQ(ierr);  
+  ierr = VecDestroy(&residu);CHKERRQ(ierr);  
+  ierr = VecDestroy(&vt);CHKERRQ(ierr);  
+  ierr = VecDestroy(&x_lsqr);CHKERRQ(ierr);  
+
 
   return 0;
 
@@ -1477,12 +1494,16 @@ static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
     KSPGetType(ksp_E,&type);
     PetscPrintf(PETSC_COMM_WORLD, "SOLVER TYPE %s  \n", type);
 
-
+    /*
     T1 = MPI_Wtime();
  ierr = KSPSetUp(ksp_E);CHKERRQ(ierr);
     ierr = KSPSolve(ksp_E,ff,XX);CHKERRQ(ierr);
     T2 = MPI_Wtime();
-    
+     */    
+
+
+
+
     Mat A;
     Vec sol;
     PetscScalar norm;
@@ -1494,20 +1515,65 @@ static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
 
 
     KSPGetOperators(ksp_E,&A,NULL);
-    MatMult(A,XX,sol);
+    /* MatMult(A,XX,sol);
     VecAXPY(sol,-1,ff);
     VecNorm(sol, NORM_2, &norm);
-
+    
 
     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Norm of error : %g\n", (double)norm); CHKERRQ(ierr);
     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
+
+     */
+
+
+        //
+    //
+    //version to control the error
+ {
+
+    Vec x2;
+    Vec sol;
+    VecDuplicate(ff,&x2);
+    VecDuplicate(ff,&sol);
+    
+    PetscScalar norm=100;
+    PetscScalar T1,T2;
+    PetscInt total,its;
+    ierr = KSPSetTolerances(ksp_E,1e-10,1e-10,PETSC_DEFAULT,30);CHKERRQ(ierr);
+    ierr = KSPSetInitialGuessNonzero(ksp_E, PETSC_TRUE); CHKERRQ(ierr);
+    T1 = MPI_Wtime();
+    while(norm>1e-3) {
+      ierr = KSPSolve(ksp_E,ff,x2);CHKERRQ(ierr);
+      KSPGetResidualNorm(ksp_E,&norm);
+      ierr = KSPGetIterationNumber(ksp_E, &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,ff);
+    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);
+
+    ierr = KSPDestroy(&ksp_E);CHKERRQ(ierr);
+    ierr = VecDestroy(&x2);CHKERRQ(ierr);
+    ierr = VecDestroy(&sol);CHKERRQ(ierr);
+  }
+     
+
+
+
+ /*
     PetscInt total;
     ierr = KSPGetIterationNumber(ksp_E, &total); CHKERRQ(ierr);
     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
+  */
 
 
-
-
+  VecCopy(ff,X2);
     KrylovMinimize(A, ff, X2);
     MatMult(A,X2,sol);