From: raphael couturier <couturie@extinction>
Date: Wed, 3 Sep 2014 18:28:33 +0000 (+0200)
Subject: mise à jour des codes
X-Git-Url: https://bilbo.iut-bm.univ-fcomte.fr/and/gitweb/GMRES2stage.git/commitdiff_plain/67de434b39a9f2df1984bae1717a1bd52916f54d?ds=inline

mise à jour des codes
---

diff --git a/code/ex10.c b/code/ex10.c
new file mode 100644
index 0000000..265df15
--- /dev/null
+++ b/code/ex10.c
@@ -0,0 +1,657 @@
+
+///home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/crashbasis/crashbasis.bin  -ksp_type gmres  -pc_type none 
+//16s
+//temps simi pour lui et nous
+
+
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/parabolic_fem/parabolic_fem.bin  -ksp_type gmres  -pc_type ilu
+//converge avec nous mais lentement
+//2152s pour lui
+//724s pour nous
+
+
+
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/epb3/epb3.bin  -ksp_type fgmres  -pc_type sor
+//fonctionne avec ilu asm 
+//9s
+//temps simi pour lui et nous
+
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/atmosmodj/atmosmodj.bin  -ksp_type fgmres  -pc_type none
+//aussi sor
+
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/bfwa398/bfwa398.bin  -ksp_type gmres   -pc_type none
+
+
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/torso3/torso3.bin  -ksp_type gmres  -pc_type none
+
+
+static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
+This version first preloads and solves a small system, then loads \n\
+another (larger) system and solves it as well.  This example illustrates\n\
+preloading of instructions with the smaller system so that more accurate\n\
+performance monitoring can be done with the larger one (that actually\n\
+is the system of interest).  See the 'Performance Hints' chapter of the\n\
+users manual for a discussion of preloading.  Input parameters include\n\
+  -f0 <input_file> : first file to load (small system)\n\
+  -f1 <input_file> : second file to load (larger system)\n\n\
+  -nearnulldim <0> : number of vectors in the near-null space immediately following matrix\n\n\
+  -trans  : solve transpose system instead\n\n";
+/*
+  This code can be used to test PETSc interface to other packages.\n\
+  Examples of command line options:       \n\
+   ./ex10 -f0 <datafile> -ksp_type preonly  \n\
+        -help -ksp_view                  \n\
+        -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
+        -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu or superlu_dist or mumps \n\
+        -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package mumps \n\
+   mpiexec -n <np> ./ex10 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
+ \n\n";
+*/
+/*T
+   Concepts: KSP^solving a linear system
+   Processors: n
+T*/
+
+/*
+  Include "petscksp.h" so that we can use KSP solvers.  Note that this file
+  automatically includes:
+     petscsys.h       - base PETSc routines   petscvec.h - vectors
+     petscmat.h - matrices
+     petscis.h     - index sets            petscksp.h - Krylov subspace methods
+     petscviewer.h - viewers               petscpc.h  - preconditioners
+*/
+#include <petscksp.h>
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+
+
+
+
+
+
+
+
+
+int KrylovMinimize(Mat A, Vec b, Vec x) {
+  
+
+  //Variables
+
+  PetscScalar  gamma, alpha, oldgamma, beta;
+  PetscReal norm=20, Eprecision=1e-7, cgprec=1e-40;     
+  PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;
+  PetscErrorCode ierr;
+  PetscScalar T1, T2;
+  KSP ksp;
+  PetscInt total=0;  
+  PetscInt size;
+  PetscInt Istart,Iend;
+  PetscInt i,its;
+  Vec       x_old, residu;
+  Mat S, AS;
+  PetscScalar *array;
+  PetscInt *ind_row;
+  Vec Alpha, p, ss, vect, r, q, Ax;
+
+
+  PetscInt first=1;
+
+  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
+  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
+
+
+
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+
+
+
+
+  VecGetSize(b,&size);
+
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);  
+
+  PetscInt aa,bb;
+  MatGetOwnershipRange(A,&aa,&bb);
+
+  //  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);  
+  //PetscSynchronizedFlush(PETSC_COMM_WORLD);
+
+
+  ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
+  ierr = MatSetSizes(S, bb-aa,  PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
+  ierr = MatSetType(S, MATDENSE); CHKERRQ(ierr);
+  ierr = MatSetUp(S); CHKERRQ(ierr);
+
+  ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
+
+  ierr = VecCreate(PETSC_COMM_WORLD, &Alpha); CHKERRQ(ierr);
+  ierr = VecSetSizes(Alpha, PETSC_DECIDE, ColS); CHKERRQ(ierr);
+  ierr = VecSetFromOptions(Alpha); CHKERRQ(ierr);
+  ierr = VecDuplicate(Alpha, &vect); CHKERRQ(ierr);
+  ierr = VecDuplicate(Alpha, &p); CHKERRQ(ierr);
+  ierr = VecDuplicate(Alpha, &ss); CHKERRQ(ierr);
+  ierr = VecDuplicate(b, &r); CHKERRQ(ierr);
+  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);
+
+
+
+  //indexes of row (these indexes are global)
+  ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
+  for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
+
+  //Initializations
+  //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, 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++){
+
+
+      //Solve
+      ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
+
+      ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
+      total += its; Iiter ++;
+
+
+
+      //Build S'
+      ierr = VecGetArray(x, &array);
+      ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
+      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);
+      //ierr = VecCopy(x, x_old); CHKERRQ(ierr);
+
+
+    }
+
+
+    //minimization step
+    if( norm>Eprecision) {
+
+      MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
+      MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
+
+
+      //Build AS
+      if(first) {
+        MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
+        first=0;
+      }
+      else
+        MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
+
+
+
+
+      //Minimization with CGLS	
+      MatMult(AS, Alpha, r);  VecAYPX(r, -1, b); //r_0 = b-AS*x_0
+
+
+      MatMultTranspose(AS, r, p); //p_0 = AS'*r_0 
+
+
+
+
+      ierr = VecCopy(p, ss); CHKERRQ(ierr); //p_0 = s_0
+      VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
+      while(gamma>cgprec && iter<iterations){
+        MatMult(AS, p, q); //q = AS*p
+        VecNorm(q, NORM_2, &alpha); alpha = alpha * alpha; //alpha = norm2(q)^2
+        alpha = gamma / alpha; 
+        VecAXPY(Alpha, alpha, p); //Alpha += alpha*p
+        VecAXPY(r, -alpha, q); //r -= alpha*q
+        MatMultTranspose(AS, r, ss); // ss = AS'*r
+        oldgamma = gamma;
+        VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
+        beta = gamma / oldgamma;
+        VecAYPX(p, beta, ss); //p = s+beta*p;
+        iter ++;
+      }	
+
+      iter = 0; giter ++;
+      //Minimizer
+      MatMult(S, Alpha, x); //x = S*Alpha;
+    }
+
+  }
+  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);
+
+  return 0;
+
+}
+
+
+
+
+
+
+
+
+
+int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
+  
+
+  //Variables
+
+  PetscScalar  alpha, beta;
+  PetscReal norm=20, Eprecision=1e-7, tol=1e-40;     
+  PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;
+  PetscErrorCode ierr;
+  PetscScalar T1, T2;
+  KSP ksp;
+  PetscInt total=0;  
+  PetscInt size;
+  PetscInt Istart,Iend;
+  PetscInt i,its;
+  Vec       x_old, residu;
+  Mat S, AS;
+  PetscScalar *array;
+  PetscInt *ind_row;
+  Vec  Ax;
+  PetscScalar norm_ksp;
+  Vec u,v,d,uu,vt,zero_long,zero_short,x_lsqr;
+
+  PetscInt first=1;
+
+  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
+  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+
+
+
+
+  VecGetSize(b,&size);
+
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);  
+
+  PetscInt aa,bb;
+  MatGetOwnershipRange(A,&aa,&bb);
+
+  //  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);  
+  //PetscSynchronizedFlush(PETSC_COMM_WORLD);
+
+
+  ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
+  ierr = MatSetSizes(S, bb-aa,  PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
+  ierr = MatSetType(S, MATDENSE); CHKERRQ(ierr);
+  ierr = MatSetUp(S); CHKERRQ(ierr);
+
+  ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
+
+
+
+  ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
+
+  ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
+  ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
+
+
+  //long vector
+  ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
+
+
+  ierr = VecDuplicate(b,&uu);CHKERRQ(ierr);
+  ierr = VecDuplicate(b,&zero_long);CHKERRQ(ierr);
+  ierr = VecSet(zero_long,0);CHKERRQ(ierr);
+
+  //small vector
+  ierr = VecCreate(PETSC_COMM_WORLD, &v); CHKERRQ(ierr);
+  ierr = VecSetSizes(v, PETSC_DECIDE, ColS); CHKERRQ(ierr);
+  ierr = VecSetFromOptions(v); CHKERRQ(ierr);
+  ierr = VecDuplicate(v,&zero_short);CHKERRQ(ierr);
+  ierr = VecSet(zero_short,0);CHKERRQ(ierr);
+  ierr = VecDuplicate(v,&d);CHKERRQ(ierr);
+  ierr = VecDuplicate(v,&vt);CHKERRQ(ierr);
+  ierr = VecDuplicate(v,&x_lsqr);CHKERRQ(ierr);
+
+
+  //indexes of row (these indexes are global)
+  ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
+  for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
+
+  //Initializations
+  //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, 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++){
+
+
+      //Solve
+      ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
+
+      ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
+      total += its; Iiter ++;
+
+
+
+      //Build S'
+      ierr = VecGetArray(x, &array);
+      ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
+      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);
+      //ierr = VecCopy(x, x_old); CHKERRQ(ierr);
+
+
+    }
+
+
+    //minimization step
+    if( norm>Eprecision) {
+
+      MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
+      MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
+
+
+
+
+      //Build AS
+      if(first) {
+        MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
+        
+        first=0;
+      }
+      else
+        MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
+
+
+
+
+
+      //LSQR
+      //LSQR
+      //LSQR
+
+
+
+      PetscScalar n2b,tolb,normr,c,s,phibar,normar,norma,thet,rhot,rho,phi;
+      PetscInt stag;
+      tolb = tol * n2b;
+      VecNorm(b, NORM_2, &n2b); //n2b = norm(b);
+      ierr = VecCopy(b, u); CHKERRQ(ierr);   //u=b
+      VecNorm(u, NORM_2, &beta); // beta=norm(u)
+      normr=beta;
+      if (beta != 0) {
+        VecAYPX(u,1/beta,zero_long); //  u = u / beta;
+      }
+      c=1;
+      s=0;
+      phibar=beta;
+      MatMultTranspose(AS, u, v); //v=A'*u
+      ierr = VecSet(x_lsqr,0);CHKERRQ(ierr);
+      VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
+      if (alpha != 0) {
+        VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
+      }
+      ierr = VecSet(d,0);CHKERRQ(ierr);
+      normar = alpha * beta;
+      norma=0;
+      //stag=0;
+      for(i=0;i<iterations;i++) {
+        MatMult(AS, v, uu); //uu=A*v
+        VecAYPX(u, -alpha, uu); //u = uu-alpha*u;
+        VecNorm(u, NORM_2, &beta); // beta=norm(u)
+        VecAYPX(u,1/beta,zero_long); //  u = u / beta;
+        norma=sqrt(norma*norma+alpha*alpha+beta*beta); // norma = norm([norma alpha beta]);
+        thet = - s * alpha;
+        rhot = c * alpha;
+        rho = sqrt(rhot*rhot + beta*beta);
+        c = rhot / rho;
+        s = - beta / rho;
+        phi = c * phibar;
+        if (phi == 0) {             // stagnation of the method
+          stag = 1;
+        }
+        phibar = s * phibar;
+        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;
+        MatMultTranspose(AS, u, vt); //vt=A'*u;
+        VecAYPX(v,-beta,vt);  // v = vt - beta * v;
+        VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
+        VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
+        normar = alpha * abs( s * phi);
+      }	
+      
+
+
+      iter = 0; giter ++;
+      //Minimizer
+      MatMult(S, x_lsqr, x); //x = S*x_small;
+    }
+
+  }
+  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);
+
+  return 0;
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+int main(int argc,char **args)
+{
+  KSP            ksp;             /* linear solver context */
+  Mat            A,A2;               /* matrix */
+  Vec            x,b,u;           /* approx solution, RHS, exact solution */
+  PetscViewer    fd;              /* viewer */
+  char           file[4][PETSC_MAX_PATH_LEN];     /* input file name */
+  PetscBool      table     =PETSC_FALSE,flg,trans=PETSC_FALSE,initialguess = PETSC_FALSE;
+  PetscBool      outputSoln=PETSC_FALSE;
+  PetscErrorCode ierr;
+  PetscInt       its,num_numfac,m,n,M,nearnulldim = 0;
+  PetscReal      norm;
+  PetscBool      preload=PETSC_TRUE,isSymmetric,cknorm=PETSC_FALSE,initialguessfile = PETSC_FALSE;
+  PetscMPIInt    rank;
+  char           initialguessfilename[PETSC_MAX_PATH_LEN];
+
+  PetscInitialize(&argc,&args,(char*)0,help);
+  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
+
+  ierr = PetscOptionsGetString(NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
+  if (flg) {
+    ierr    = PetscStrcpy(file[1],file[0]);CHKERRQ(ierr);
+    preload = PETSC_FALSE;
+  } else {
+    ierr = PetscOptionsGetString(NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
+    if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
+    ierr = PetscOptionsGetString(NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
+    if (!flg) preload = PETSC_FALSE;   /* don't bother with second system */
+  }
+
+
+  PetscPreLoadBegin(preload,"Load system");
+
+
+
+  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);CHKERRQ(ierr);
+
+  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
+  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
+  ierr = MatLoad(A,fd);CHKERRQ(ierr);
+
+  /*  char           ordering[256] = MATORDERINGRCM;
+  IS             rowperm       = NULL,colperm = NULL;
+  ierr = MatGetOrdering(A2,ordering,&rowperm,&colperm);CHKERRQ(ierr);
+  ierr = MatPermute(A2,rowperm,colperm,&A);CHKERRQ(ierr);
+   */
+
+  ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
+
+  ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
+
+ ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
+      ierr = VecSetSizes(b,m,PETSC_DECIDE);CHKERRQ(ierr);
+      ierr = VecSetFromOptions(b);CHKERRQ(ierr);
+      ierr = VecSet(b,1);CHKERRQ(ierr);
+
+
+
+
+  ierr = MatGetVecs(A,&x,NULL);CHKERRQ(ierr);
+  ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
+
+    ierr = VecSet(x,1.0);CHKERRQ(ierr);
+
+
+  ierr       = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
+  ierr       = KSPSetInitialGuessNonzero(ksp,initialguess);CHKERRQ(ierr);
+  num_numfac = 1;
+
+
+
+
+
+    ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
+
+   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+
+      ierr = KSPSetTolerances(ksp,1e-10,1e-10,PETSC_DEFAULT,5000000);CHKERRQ(ierr);
+  PC             pc;
+      KSPGetPC(ksp, &pc);
+  PCType         type;
+        PCGetType(pc, &type);
+
+         PetscPrintf(PETSC_COMM_WORLD, "PC TYPE %s  \n", type);
+
+
+       PetscScalar T1,T2;
+	T1 = MPI_Wtime();
+        ierr = KSPSetUp(ksp);CHKERRQ(ierr);
+        ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+        T2 = MPI_Wtime();
+
+        
+    ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
+      ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
+
+
+  
+      ierr = MatMult(A,x,u);CHKERRQ(ierr);
+   ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr);
+    ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr);
+
+ {
+
+    Vec x2;
+    Vec sol;
+    VecDuplicate(b,&x2);
+    VecDuplicate(b,&sol);
+    
+    KrylovMinimize(A, b, x2);
+
+
+
+    MatMult(A,x2,sol);
+    VecAXPY(sol,-1,b);
+    VecNorm(sol, NORM_2, &norm);
+    ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Error Krylov Minimization %g\n",norm);
+
+  }
+
+{
+
+    Vec x2;
+    Vec sol;
+    VecDuplicate(b,&x2);
+    VecDuplicate(b,&sol);
+    
+    KrylovMinimizeLSQR(A, b, x2);
+
+
+
+    MatMult(A,x2,sol);
+    VecAXPY(sol,-1,b);
+    VecNorm(sol, NORM_2, &norm);
+    ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Error KrylovLSQR Minimization %g\n",norm);
+
+  }
+
+
+  /*
+     Free work space.  All PETSc objects should be destroyed when they
+     are no longer needed.
+  */
+  ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = VecDestroy(&b);CHKERRQ(ierr);
+  ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr);
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  PetscPreLoadEnd();
+  /* -----------------------------------------------------------
+                      End of linear solver loop
+     ----------------------------------------------------------- */
+
+  ierr = PetscFinalize();
+  return 0;
+}
+
diff --git a/code/ex15.c b/code/ex15.c
index 1fffcbe..820ec8d 100644
--- a/code/ex15.c
+++ b/code/ex15.c
@@ -120,13 +120,14 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, 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++){
 
@@ -310,7 +311,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 16); CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 
 
@@ -318,6 +319,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++){
 
@@ -515,7 +517,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);}
@@ -574,6 +577,16 @@ int main(int argc,char **args)
   ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
   ierr = KSPSetTolerances(ksp,1e-7,1e-7,PETSC_DEFAULT,5000000);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.,
         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
@@ -588,7 +601,7 @@ int main(int argc,char **args)
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
   PetscScalar T1,T2;
 	T1 = MPI_Wtime();
-  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+        ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 	T2 = MPI_Wtime();
   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                       Check solution and clean up
diff --git a/code/ex29.c b/code/ex29.c
index ec89b02..427f70a 100644
--- a/code/ex29.c
+++ b/code/ex29.c
@@ -128,13 +128,14 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, 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++){
 
@@ -318,7 +319,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 16); CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 
 
@@ -326,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++){
 
@@ -504,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);
+
+
+
+  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();
 
diff --git a/code/ex45.c b/code/ex45.c
index dd3f3a6..d3e5fb4 100644
--- a/code/ex45.c
+++ b/code/ex45.c
@@ -45,7 +45,7 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
   PetscScalar  gamma, alpha, oldgamma, beta;
   PetscReal norm=20, Eprecision=1e-8, cgprec=1e-40;     
-  PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
+  PetscInt giter=0, ColS=6, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
   PetscErrorCode ierr;
   PetscScalar T1, T2;
   KSP ksp;
@@ -108,13 +108,14 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, 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++){
 
@@ -221,7 +222,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   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=6, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
   PetscErrorCode ierr;
   PetscScalar T1, T2;
   KSP ksp;
@@ -297,7 +298,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 16); CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 
 
@@ -305,6 +306,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++){
 
@@ -466,8 +468,22 @@ int main(int argc,char **argv)
   ierr = KSPSetFromOptions(ksp);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);
+  KSPGetType(ksp,&type);
+  PetscPrintf(PETSC_COMM_WORLD, "SOLVER TYPE %s  \n", type);
+
+
+
   ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, PETSC_DEFAULT, 50000000); CHKERRQ(ierr);
   T1 = MPI_Wtime();
+  ierr = KSPSetUp(ksp);CHKERRQ(ierr);
   ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);
   T2 = MPI_Wtime();
 
diff --git a/code/ex49.c b/code/ex49.c
index 260effb..b93b5c6 100644
--- a/code/ex49.c
+++ b/code/ex49.c
@@ -152,13 +152,14 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, 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++){
 
@@ -342,7 +343,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 16); CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 
 
@@ -350,6 +351,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++){
 
@@ -1445,9 +1447,9 @@ static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
 
 
   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_E);CHKERRQ(ierr);
-  ierr = KSPSetOptionsPrefix(ksp_E,"elas_");CHKERRQ(ierr);  /* elasticity */
+  // ierr = KSPSetOptionsPrefix(ksp_E,"elas_");CHKERRQ(ierr);  /* elasticity */
 
-  ierr = PetscOptionsGetBool(NULL,"-use_nonsymbc",&use_nonsymbc,&flg);CHKERRQ(ierr);
+  //ierr = PetscOptionsGetBool(NULL,"-use_nonsymbc",&use_nonsymbc,&flg);CHKERRQ(ierr);
   /* solve */
   if (!use_nonsymbc) {
     Mat        AA;
@@ -1461,12 +1463,23 @@ static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
     ierr = KSPSetOperators(ksp_E,AA,AA);CHKERRQ(ierr);
     ierr = KSPSetFromOptions(ksp_E);CHKERRQ(ierr);
 
- ierr = KSPSetFromOptions(ksp_E);CHKERRQ(ierr);
 
    PetscScalar T1,T2;
     ierr = KSPSetTolerances(ksp_E, 1e-7, 1e-7, PETSC_DEFAULT, 50000000); CHKERRQ(ierr);
-    T1 = MPI_Wtime();
 
+
+    PC             pc;
+    KSPGetPC(ksp_E, &pc);
+    PCType         type;
+    PCGetType(pc, &type);
+
+    PetscPrintf(PETSC_COMM_WORLD, "PC TYPE %s  \n", type);
+    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();