From bc9c5a57e4b0b66714c6ac200b4aff7c26b38608 Mon Sep 17 00:00:00 2001
From: raphael couturier <couturie@extinction>
Date: Wed, 13 Aug 2014 14:35:10 +0200
Subject: [PATCH] ajout ex29, lsqr ne semble pas bon sur cet exemple...

---
 code/ex15.c | 267 +++++++++++++++++++-
 code/ex29.c | 700 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 960 insertions(+), 7 deletions(-)
 create mode 100644 code/ex29.c

diff --git a/code/ex15.c b/code/ex15.c
index 20144dd..0e0a8c4 100644
--- a/code/ex15.c
+++ b/code/ex15.c
@@ -1,4 +1,8 @@
 
+// /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 4    ./ex15 -m 400 -n 400
+
+
+
 static char help[] = "Solves a linear system in parallel with KSP.  Also\n\
 illustrates setting a user-defined shell preconditioner and using the\n\
 macro __FUNCT__ to define routine names for use in error handling.\n\
@@ -51,11 +55,11 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
   //Variables
 
-  PetscScalar  gamma, alpha, oldgamma, beta, t2;
+  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;
   PetscErrorCode ierr;
-  PetscScalar T1, T2, t1;
+  PetscScalar T1, T2;
   KSP ksp;
   PetscInt total=0;  
   PetscInt size;
@@ -66,7 +70,6 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
   PetscScalar *array;
   PetscInt *ind_row;
   Vec Alpha, p, ss, vect, r, q, Ax;
-  PetscScalar norm_ksp;
 
 
   PetscInt first=1;
@@ -207,7 +210,6 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
   }
   T2 = MPI_Wtime();
-  t1 = T2 - T1;
   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);
 
@@ -222,6 +224,242 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
 
 
+
+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;
+  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, MATMPIDENSE); 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, 16); CHKERRQ(ierr);
+  ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
+
+
+
+
+  //GMRES WITH MINIMIZATION
+  T1 = MPI_Wtime();
+  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);
+
+
+
+      //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)
 {
   Vec            x,b,u;   /* approx solution, RHS, exact solution */
@@ -229,11 +467,9 @@ int main(int argc,char **args)
   KSP            ksp;      /* linear solver context */
   PC             pc;        /* preconditioner context */
   PetscReal      norm;      /* norm of solution error */
-  SampleShellPC  *shell;    /* user-defined preconditioner context */
-  PetscScalar    v,one = 1.0,none = -1.0;
+  PetscScalar    v,one = 1.0;
   PetscInt       i,j,Ii,J,Istart,Iend,m = 8,n = 7,its;
   PetscErrorCode ierr;
-  PetscBool      user_defined_pc = PETSC_FALSE;
 
   PetscInitialize(&argc,&args,(char*)0,help);
   ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
@@ -387,6 +623,23 @@ int main(int argc,char **args)
   }
 
 
+ {
+
+    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 LSQR %g\n",norm);
+
+  }
 
 
 
diff --git a/code/ex29.c b/code/ex29.c
new file mode 100644
index 0000000..5e21c65
--- /dev/null
+++ b/code/ex29.c
@@ -0,0 +1,700 @@
+
+
+//  /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 3 ex29 -da_grid_x 900 -da_grid_y 900
+
+
+
+/*T
+   Concepts: KSP^solving a system of linear equations
+   Concepts: KSP^Laplacian, 2d
+   Processors: n
+T*/
+
+/*
+Added at the request of Marc Garbey.
+
+Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation
+
+   -div \rho grad u = f,  0 < x,y < 1,
+
+with forcing function
+
+   f = e^{-x^2/\nu} e^{-y^2/\nu}
+
+with Dirichlet boundary conditions
+
+   u = f(x,y) for x = 0, x = 1, y = 0, y = 1
+
+or pure Neumman boundary conditions
+
+This uses multigrid to solve the linear system
+*/
+
+static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";
+
+#include <petscdm.h>
+#include <petscdmda.h>
+#include <petscksp.h>
+
+extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
+extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
+
+typedef enum {DIRICHLET, NEUMANN} BCType;
+
+typedef struct {
+  PetscReal rho;
+  PetscReal nu;
+  BCType    bcType;
+} UserContext;
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+
+
+
+
+
+
+
+
+
+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;
+  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, MATMPIDENSE); 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, 16); CHKERRQ(ierr);
+  ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
+
+
+
+  //GMRES WITH MINIMIZATION
+  T1 = MPI_Wtime();
+  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);
+
+
+
+      //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-5, 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;
+  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, MATMPIDENSE); 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, 16); CHKERRQ(ierr);
+  ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
+
+
+
+
+  //GMRES WITH MINIMIZATION
+  T1 = MPI_Wtime();
+  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);
+
+
+
+      //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 **argv)
+{
+  KSP            ksp;
+  DM             da;
+  UserContext    user;
+  const char     *bcTypes[2] = {"dirichlet","neumann"};
+  PetscErrorCode ierr;
+  PetscInt       bc;
+  Vec            b,x;
+  Mat A;
+  PetscScalar T1,T2,norm;
+  Vec sol;
+  PetscInt total;
+
+  PetscInitialize(&argc,&argv,(char*)0,help);
+
+  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);
+  ierr = DMDASetFieldName(da,0,"Pressure");CHKERRQ(ierr);
+
+  ierr        = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMqq");
+  user.rho    = 1.0;
+  ierr        = PetscOptionsReal("-rho", "The conductivity", "ex29.c", user.rho, &user.rho, NULL);CHKERRQ(ierr);
+  user.nu     = 0.1;
+  ierr        = PetscOptionsReal("-nu", "The width of the Gaussian source", "ex29.c", user.nu, &user.nu, NULL);CHKERRQ(ierr);
+  bc          = (PetscInt)DIRICHLET;
+  ierr        = PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,NULL);CHKERRQ(ierr);
+  user.bcType = (BCType)bc;
+  ierr        = PetscOptionsEnd();
+
+  ierr = KSPSetComputeRHS(ksp,ComputeRHS,&user);CHKERRQ(ierr);
+  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);
+  T1 = MPI_Wtime();
+  ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);
+  T2 = MPI_Wtime();
+
+
+
+  ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr);
+  ierr = KSPGetRhs(ksp,&b);CHKERRQ(ierr);
+  KSPGetOperators(ksp,&A,NULL);
+
+  VecDuplicate(x,&sol);
+  MatMult(A,x,sol);
+  VecAXPY(sol,-1,b);
+  VecNorm(sol, NORM_2, &norm);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "Error %g\n",norm);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
+
+  ierr = KSPGetIterationNumber(ksp, &total); CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
+
+ {
+
+    Vec x;
+    VecDuplicate(b,&x);
+    
+    KrylovMinimize(A, b, x);
+
+
+
+    MatMult(A,x,sol);
+
+
+    VecAXPY(sol,-1,b);
+    VecNorm(sol, NORM_2, &norm);
+    ierr = PetscPrintf(PETSC_COMM_WORLD, "Error2 %g\n",norm);
+
+  }
+
+
+ {
+
+    Vec x;
+    VecDuplicate(b,&x);
+    
+    KrylovMinimizeLSQR(A, b, x);
+
+
+
+    MatMult(A,x,sol);
+
+
+    VecAXPY(sol,-1,b);
+    VecNorm(sol, NORM_2, &norm);
+    ierr = PetscPrintf(PETSC_COMM_WORLD, "Error LSQR %g\n",norm);
+
+  }
+
+
+
+
+  ierr = DMDestroy(&da);CHKERRQ(ierr);
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = PetscFinalize();
+
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ComputeRHS"
+PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
+{
+  UserContext    *user = (UserContext*)ctx;
+  PetscErrorCode ierr;
+  PetscInt       i,j,mx,my,xm,ym,xs,ys;
+  PetscScalar    Hx,Hy;
+  PetscScalar    **array;
+  DM             da;
+
+  PetscFunctionBeginUser;
+  ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
+  ierr = DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
+  Hx   = 1.0 / (PetscReal)(mx-1);
+  Hy   = 1.0 / (PetscReal)(my-1);
+  ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
+  ierr = DMDAVecGetArray(da, b, &array);CHKERRQ(ierr);
+  for (j=ys; j<ys+ym; j++) {
+    for (i=xs; i<xs+xm; i++) {
+      array[j][i] = PetscExpScalar(-((PetscReal)i*Hx)*((PetscReal)i*Hx)/user->nu)*PetscExpScalar(-((PetscReal)j*Hy)*((PetscReal)j*Hy)/user->nu)*Hx*Hy;
+    }
+  }
+  ierr = DMDAVecRestoreArray(da, b, &array);CHKERRQ(ierr);
+  ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
+
+  /* force right hand side to be consistent for singular matrix */
+  /* note this is really a hack, normally the model would provide you with a consistent right handside */
+  if (user->bcType == NEUMANN) {
+    MatNullSpace nullspace;
+
+    ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
+    ierr = MatNullSpaceRemove(nullspace,b);CHKERRQ(ierr);
+    ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "ComputeRho"
+PetscErrorCode ComputeRho(PetscInt i, PetscInt j, PetscInt mx, PetscInt my, PetscReal centerRho, PetscReal *rho)
+{
+  PetscFunctionBeginUser;
+  if ((i > mx/3.0) && (i < 2.0*mx/3.0) && (j > my/3.0) && (j < 2.0*my/3.0)) {
+    *rho = centerRho;
+  } else {
+    *rho = 1.0;
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ComputeMatrix"
+PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx)
+{
+  UserContext    *user = (UserContext*)ctx;
+  PetscReal      centerRho;
+  PetscErrorCode ierr;
+  PetscInt       i,j,mx,my,xm,ym,xs,ys;
+  PetscScalar    v[5];
+  PetscReal      Hx,Hy,HydHx,HxdHy,rho;
+  MatStencil     row, col[5];
+  DM             da;
+
+  PetscFunctionBeginUser;
+  ierr      = KSPGetDM(ksp,&da);CHKERRQ(ierr);
+  centerRho = user->rho;
+  ierr      = DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
+  Hx        = 1.0 / (PetscReal)(mx-1);
+  Hy        = 1.0 / (PetscReal)(my-1);
+  HxdHy     = Hx/Hy;
+  HydHx     = Hy/Hx;
+  ierr      = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
+  for (j=ys; j<ys+ym; j++) {
+    for (i=xs; i<xs+xm; i++) {
+      row.i = i; row.j = j;
+      ierr  = ComputeRho(i, j, mx, my, centerRho, &rho);CHKERRQ(ierr);
+      if (i==0 || j==0 || i==mx-1 || j==my-1) {
+        if (user->bcType == DIRICHLET) {
+          v[0] = 2.0*rho*(HxdHy + HydHx);
+          ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
+        } else if (user->bcType == NEUMANN) {
+          PetscInt numx = 0, numy = 0, num = 0;
+          if (j!=0) {
+            v[num] = -rho*HxdHy;              col[num].i = i;   col[num].j = j-1;
+            numy++; num++;
+          }
+          if (i!=0) {
+            v[num] = -rho*HydHx;              col[num].i = i-1; col[num].j = j;
+            numx++; num++;
+          }
+          if (i!=mx-1) {
+            v[num] = -rho*HydHx;              col[num].i = i+1; col[num].j = j;
+            numx++; num++;
+          }
+          if (j!=my-1) {
+            v[num] = -rho*HxdHy;              col[num].i = i;   col[num].j = j+1;
+            numy++; num++;
+          }
+          v[num] = numx*rho*HydHx + numy*rho*HxdHy; col[num].i = i;   col[num].j = j;
+          num++;
+          ierr = MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);
+        }
+      } else {
+        v[0] = -rho*HxdHy;              col[0].i = i;   col[0].j = j-1;
+        v[1] = -rho*HydHx;              col[1].i = i-1; col[1].j = j;
+        v[2] = 2.0*rho*(HxdHy + HydHx); col[2].i = i;   col[2].j = j;
+        v[3] = -rho*HydHx;              col[3].i = i+1; col[3].j = j;
+        v[4] = -rho*HxdHy;              col[4].i = i;   col[4].j = j+1;
+        ierr = MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);CHKERRQ(ierr);
+      }
+    }
+  }
+  ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  if (user->bcType == NEUMANN) {
+    MatNullSpace nullspace;
+
+    ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
+    ierr = MatSetNullSpace(jac,nullspace);CHKERRQ(ierr);
+    ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
-- 
2.39.5