From: raphael couturier <couturie@extinction>
Date: Tue, 12 Aug 2014 07:20:50 +0000 (+0200)
Subject: ajout d'ex15 qui fonctionne en version 3.5.1
X-Git-Url: https://bilbo.iut-bm.univ-fcomte.fr/and/gitweb/GMRES2stage.git/commitdiff_plain/d085467b95f24aa53166c5284f70be0fd24da9fc

ajout d'ex15 qui fonctionne en version 3.5.1
---

diff --git a/code/ex15.c b/code/ex15.c
new file mode 100644
index 0000000..20144dd
--- /dev/null
+++ b/code/ex15.c
@@ -0,0 +1,405 @@
+
+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\
+Input parameters include:\n\
+  -user_defined_pc : Activate a user-defined preconditioner\n\n";
+
+/*T
+   Concepts: KSP^basic parallel example
+   Concepts: PC^setting a user-defined shell preconditioner
+   Concepts: error handling^Using the macro __FUNCT__ to define routine names;
+   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>
+
+/* Define context for user-provided preconditioner */
+typedef struct {
+  Vec diag;
+} SampleShellPC;
+
+
+
+/*
+   User-defined routines.  Note that immediately before each routine below,
+   we define the macro __FUNCT__ to be a string containing the routine name.
+   If defined, this macro is used in the PETSc error handlers to provide a
+   complete traceback of routine names.  All PETSc library routines use this
+   macro, and users can optionally employ it as well in their application
+   codes.  Note that users can get a traceback of PETSc errors regardless of
+   whether they define __FUNCT__ in application codes; this macro merely
+   provides the added traceback detail of the application routine names.
+*/
+#undef __FUNCT__
+#define __FUNCT__ "main"
+
+
+
+
+
+int KrylovMinimize(Mat A, Vec b, Vec x) {
+  
+
+  //Variables
+
+  PetscScalar  gamma, alpha, oldgamma, beta, t2;
+  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;
+  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;
+  PetscScalar norm_ksp;
+
+
+  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();
+  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);
+
+  return 0;
+
+}
+
+
+
+
+
+
+
+
+int main(int argc,char **args)
+{
+  Vec            x,b,u;   /* approx solution, RHS, exact solution */
+  Mat            A;         /* linear system matrix */
+  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;
+  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);
+  ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
+
+  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+         Compute the matrix and right-hand-side vector that define
+         the linear system, Ax = b.
+     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+  /*
+     Create parallel matrix, specifying only its global dimensions.
+     When using MatCreate(), the matrix format can be specified at
+     runtime. Also, the parallel partioning of the matrix is
+     determined by PETSc at runtime.
+  */
+  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
+  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);
+  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
+  ierr = MatSetUp(A);CHKERRQ(ierr);
+
+  /*
+     Currently, all PETSc parallel matrix formats are partitioned by
+     contiguous chunks of rows across the processors.  Determine which
+     rows of the matrix are locally owned.
+  */
+  ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
+
+  /*
+     Set matrix elements for the 2-D, five-point stencil in parallel.
+      - Each processor needs to insert only elements that it owns
+        locally (but any non-local elements will be sent to the
+        appropriate processor during matrix assembly).
+      - Always specify global rows and columns of matrix entries.
+   */
+  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);}
+    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);}
+    v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);CHKERRQ(ierr);
+  }
+
+  /*
+     Assemble matrix, using the 2-step process:
+       MatAssemblyBegin(), MatAssemblyEnd()
+     Computations can be done while messages are in transition
+     by placing code between these two statements.
+  */
+  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+
+  /*
+     Create parallel vectors.
+      - When using VecCreate() VecSetSizes() and VecSetFromOptions(),
+        we specify only the vector's global
+        dimension; the parallel partitioning is determined at runtime.
+      - Note: We form 1 vector from scratch and then duplicate as needed.
+  */
+  ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
+  ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr);
+  ierr = VecSetFromOptions(u);CHKERRQ(ierr);
+  ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
+  ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
+
+  /*
+     Set exact solution; then compute right-hand-side vector.
+  */
+  ierr = VecSet(u,one);CHKERRQ(ierr);
+  ierr = MatMult(A,u,b);CHKERRQ(ierr);
+
+  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+                Create the linear solver and set various options
+     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+  /*
+     Create linear solver context
+  */
+  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
+
+  /*
+     Set operators. Here the matrix that defines the linear system
+     also serves as the preconditioning matrix.
+  */
+  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
+
+  /*
+     Set linear solver defaults for this problem (optional).
+     - By extracting the KSP and PC contexts from the KSP context,
+       we can then directly call any KSP and PC routines
+       to set various options.
+  */
+  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp,1e-9,1e-9,PETSC_DEFAULT,5000000);CHKERRQ(ierr);
+
+  /*
+    Set runtime options, e.g.,
+        -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
+    These options will override those specified above as long as
+    KSPSetFromOptions() is called _after_ any other customization
+    routines.
+  */
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+
+  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+                      Solve the linear system
+     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+  PetscScalar T1,T2;
+	T1 = MPI_Wtime();
+  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+	T2 = MPI_Wtime();
+  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+                      Check solution and clean up
+     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+  /*
+     Check the error
+  */
+  Vec sol;
+  VecDuplicate(b,&sol);
+  MatMult(A,x,sol);
+  VecAXPY(sol,-1,b);
+  VecNorm(sol, NORM_2, &norm);
+  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n\n\n", T2-T1); 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);
+
+  }
+
+
+
+
+
+  /*
+     Free work space.  All PETSc objects should be destroyed when they
+     are no longer needed.
+  */
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = VecDestroy(&u);CHKERRQ(ierr);  ierr = VecDestroy(&x);CHKERRQ(ierr);
+  ierr = VecDestroy(&b);CHKERRQ(ierr);  ierr = MatDestroy(&A);CHKERRQ(ierr);
+
+  ierr = PetscFinalize();
+  return 0;
+
+}
+