dolfin team mailing list archive

Thread
Date
Re: patch till dolfin

To: Discussion of DOLFIN development <dolfin-dev@xxxxxxxxxxxxxxxxxxxxxx>
From: Anders Logg <logg@xxxxxxxxxxxxxxxx>
Date: Thu, 1 Jul 2004 13:21:23 +0200
In-reply-to: <40E3E5E1.5030309@math.chalmers.se>
Mail-followup-to: Discussion of DOLFIN development <dolfin-dev@xxxxxxxxxxxxxxxxxxxxxx>
Reply-to: Discussion of DOLFIN development <dolfin-dev@xxxxxxxxxxxxxxxxxxxxxx>
User-agent: Mutt/1.5.5.1+cvs20040105i
Great, I will apply the patch when I get back from my vacation.

What about GMRES?

/Anders

On Thu, Jul 01, 2004 at 12:22:25PM +0200, Erik Svensson wrote:

> Här kommer en patch.
> 
> Gjorda tillägg i: FGMRES, GaussSeidle och LinearSolver.
> 
> /Erik
> 

> ? Makefile
> ? config.log
> ? config.status
> ? dolfin-erik-2004-07-01.patch
> ? src/Makefile
> ? src/config/Makefile
> ? src/config/dolfin-config
> ? src/config/var.tmp
> ? src/demo/Makefile
> ? src/demo/bench/dolfin-bench
> ? src/demo/elementfunction/dolfin-elementfunction
> ? src/demo/fem/dolfin-fem
> ? src/demo/la/generic_matrix/dolfin-la-generic_matrix
> ? src/demo/la/test/dolfin-test
> ? src/demo/mesh/dolfin-mesh
> ? src/demo/ode/dolfin-tanganyika
> ? src/demo/quadrature/dolfin-quadrature
> ? src/demo/solvers/Makefile
> ? src/demo/solvers/convdiff/dolfin-convdiff
> ? src/demo/solvers/elasticity/dolfin-elasticity
> ? src/demo/solvers/elasticity-stationary/dolfin-elasticity-stationary
> ? src/demo/solvers/euler/dolfin-euler
> ? src/demo/solvers/heat/dolfin-heat
> ? src/demo/solvers/navierstokes/Makefile
> ? src/demo/solvers/navierstokes/benchmark/dolfin-navierstokes-benchmark
> ? src/demo/solvers/ode/Makefile
> ? src/demo/solvers/ode/bench/dolfin-ode-bench
> ? src/demo/solvers/ode/bistable/dolfin-ode-bistable
> ? src/demo/solvers/ode/lorenz/dolfin-lorenz
> ? src/demo/solvers/ode/mechanical/Makefile
> ? src/demo/solvers/ode/mechanical/spring/dolfin-ode-spring
> ? src/demo/solvers/ode/modeling/Makefile
> ? src/demo/solvers/ode/modeling/lattice/dolfin-ode-modeling-lattice
> ? src/demo/solvers/ode/modeling/simple/dolfin-ode-modeling-simple
> ? src/demo/solvers/ode/newton/dolfin-ode-newton
> ? src/demo/solvers/ode/performance/dolfin-ode-perf-test
> ? src/demo/solvers/ode/stiff/stabilization/dolfin-ode-stiff-stabilization
> ? src/demo/solvers/ode/stiff/testproblems/dolfin-ode-stiff-testproblems
> ? src/demo/solvers/ode/test/dolfin-ode-test
> ? src/demo/solvers/poisson/dolfin-poisson
> ? src/demo/solvers/poisson-multigrid/dolfin-poisson-multigrid
> ? src/demo/solvers/template/dolfin-template
> ? src/demo/solvers/wave/dolfin-wave
> ? src/demo/solvers/wave-vector/dolfin-wavevector
> ? src/demo/test/dolfin-test
> ? src/greeting/Makefile
> ? src/greeting/var.tmp
> ? src/kernel/Makefile
> ? src/kernel/common/.deps
> ? src/kernel/common/Makefile
> ? src/kernel/common/dolfin/Makefile
> ? src/kernel/element/.deps
> ? src/kernel/element/Makefile
> ? src/kernel/element/dolfin/Makefile
> ? src/kernel/fem/.deps
> ? src/kernel/fem/dolfin/Makefile
> ? src/kernel/form/.deps
> ? src/kernel/form/Makefile
> ? src/kernel/form/dolfin/Makefile
> ? src/kernel/function/.deps
> ? src/kernel/function/Makefile
> ? src/kernel/function/dolfin/Makefile
> ? src/kernel/io/.deps
> ? src/kernel/io/Makefile
> ? src/kernel/io/dolfin/Makefile
> ? src/kernel/la/.deps
> ? src/kernel/la/Makefile
> ? src/kernel/la/dolfin/Makefile
> ? src/kernel/log/.deps
> ? src/kernel/log/Makefile
> ? src/kernel/log/dolfin/Makefile
> ? src/kernel/main/.deps
> ? src/kernel/main/Makefile
> ? src/kernel/main/dolfin_modules.h
> ? src/kernel/main/dolfin/Makefile
> ? src/kernel/map/.deps
> ? src/kernel/map/Makefile
> ? src/kernel/map/dolfin/Makefile
> ? src/kernel/math/.deps
> ? src/kernel/math/Makefile
> ? src/kernel/math/dolfin/Makefile
> ? src/kernel/mesh/.deps
> ? src/kernel/mesh/Makefile
> ? src/kernel/mesh/dolfin/Makefile
> ? src/kernel/ode/.deps
> ? src/kernel/ode/Makefile
> ? src/kernel/ode/dolfin/Makefile
> ? src/kernel/quadrature/.deps
> ? src/kernel/quadrature/Makefile
> ? src/kernel/quadrature/dolfin/Makefile
> ? src/kernel/settings/.deps
> ? src/kernel/settings/Makefile
> ? src/kernel/settings/dolfin/Makefile
> ? src/modules/Makefile
> ? src/modules/modules.list
> ? src/modules/convdiff/.deps
> ? src/modules/convdiff/Makefile
> ? src/modules/elasticity/.deps
> ? src/modules/elasticity/Makefile
> ? src/modules/elasticity-stationary/.deps
> ? src/modules/elasticity-stationary/Makefile
> ? src/modules/euler/.deps
> ? src/modules/euler/Makefile
> ? src/modules/heat/.deps
> ? src/modules/heat/Makefile
> ? src/modules/navierstokes/.deps
> ? src/modules/navierstokes/Makefile
> ? src/modules/odesolver/.deps
> ? src/modules/odesolver/Makefile
> ? src/modules/poisson/.deps
> ? src/modules/poisson/Makefile
> ? src/modules/poisson-multigrid/.deps
> ? src/modules/poisson-multigrid/Makefile
> ? src/modules/template/.deps
> ? src/modules/template/Makefile
> ? src/modules/wave/.deps
> ? src/modules/wave/Makefile
> ? src/modules/wave-vector/.deps
> ? src/modules/wave-vector/Makefile
> ? src/post/Makefile
> ? src/pre/Makefile
> ? src/utils/Makefile
> ? src/utils/inp2dx/.deps
> ? src/utils/inp2dx/Makefile
> ? src/utils/inp2dx/inp2dx
> ? src/utils/inp2dx/inp2dxanim
> Index: src/demo/la/Makefile
> ===================================================================
> RCS file: /cvs/dolfin/src/demo/la/Makefile,v
> retrieving revision 1.128
> diff -r1.128 Makefile
> 19c19
> < prefix = /usr/local
> > prefix = /home/erik/local
> 65c65
> < CXXFLAGS = -g -O2
> > CXXFLAGS = -DDEBUG=1 -g -O2 -Wall -Werror -pedantic -ansi -std=c++98
> Index: src/demo/solvers/ode/stiff/Makefile
> ===================================================================
> RCS file: /cvs/dolfin/src/demo/solvers/ode/stiff/Makefile,v
> retrieving revision 1.83
> diff -r1.83 Makefile
> 19c19
> < prefix = /usr/local
> > prefix = /home/erik/local
> 65c65
> < CXXFLAGS = -g -O2
> > CXXFLAGS = -DDEBUG=1 -g -O2 -Wall -Werror -pedantic -ansi -std=c++98
> Index: src/kernel/fem/Makefile
> ===================================================================
> RCS file: /cvs/dolfin/src/kernel/fem/Makefile,v
> retrieving revision 1.131
> diff -r1.131 Makefile
> 19c19
> < prefix = /usr/local
> > prefix = /home/erik/local
> 65c65
> < CXXFLAGS = -g -O2
> > CXXFLAGS = -DDEBUG=1 -g -O2 -Wall -Werror -pedantic -ansi -std=c++98
> Index: src/kernel/la/FGMRES.cpp
> ===================================================================
> RCS file: /cvs/dolfin/src/kernel/la/FGMRES.cpp,v
> retrieving revision 1.1
> diff -r1.1 FGMRES.cpp
> 4a5,6
> > #include <dolfin/dolfin_math.h>
> > #include <dolfin/dolfin_settings.h>
> 5a8,9
> > #include <dolfin/Matrix.h>
> > #include <dolfin/Vector.h>
> 10,11c14,16
> < FGMRES::FGMRES(const Matrix& A, real tol, unsigned int maxiter)
> <   : A(A), tol(tol), maxiter(maxiter)
> > FGMRES::FGMRES(const Matrix& A, unsigned int restarts, unsigned int maxiter, 
> > 	   real tol, Preconditioner& pc)
> >   : Preconditioner(), A(A), pc(pc)
> 13c18,21
> <   dolfin_error("Not implemented.");
> >   this -> restarts = restarts;
> >   this -> maxiter = maxiter;
> >   this -> tol = tol;
> >   solve2convergence = true;
> 16c24,26
> < FGMRES::~FGMRES()
> > FGMRES::FGMRES(const Matrix& A, unsigned int maxiter, 
> > 	   real tol, Preconditioner& pc)
> >   : Preconditioner(), A(A), pc(pc)
> 18c28,31
> <   dolfin_error("Not implemented.");
> >   restarts = 1;
> >   this -> maxiter = maxiter;
> >   this -> tol = tol;
> >   solve2convergence = false;
> 21c34
> < void FGMRES::solve(Vector& x, const Vector& b)
> > FGMRES::~FGMRES()
> 23c36
> <   dolfin_error("Not implemented.");
> >  
> 27c40,41
> < 		  real tol, unsigned int maxiter)
> > 		      unsigned int restarts, unsigned int maxiter, real tol,
> > 		      Preconditioner& pc)
> 30,31c44,45
> <   FGMRES fgmres(A, tol, maxiter);
> < 
> >   FGMRES fgmres(A, restarts, maxiter, tol, pc);
> >   
> 35a50,227
> > void FGMRES::solve(Vector& x, const Vector& b)
> > {
> >   // Check compatibility in the matrix and vectors.
> >   check(A, x, b);
> > 
> >   bnorm = b.norm();
> >   // Check if b=0 => x=0.
> >   if ( bnorm < DOLFIN_EPS ){
> >     x = 0.0;
> >     return;
> >   }
> >   
> >   // Compute residual
> >   Vector r(b.size());
> >   real rnorm = residual(A,x,b,r);
> > 
> >   // Restarted GMRES
> >   for ( unsigned int i = 0; i < restarts; i++ ){
> >     
> >     // Krylov iterations
> >     unsigned int noiter = iterator(x,b,r);
> >     
> >     // Check stopping criterion
> >     if ( solve2convergence ) { 
> >       
> >       // Evaluate residual
> >       rnorm = residual(A,x,b,r);
> > 
> >       // Check residual
> >       if ( rnorm < tol*bnorm ){
> > 	dolfin_info("Restarted FGMRES converged after %i iterations (residual %1.5e)",
> > 		    i*maxiter+noiter, rnorm);
> > 	break;
> >       }
> >       if ( i == restarts-1 )
> > 	dolfin_error2("Restarted FGMRES did not converge (%i iterations, residual %1.5e)",
> > 		     i*maxiter+noiter, rnorm);
> >     }  
> >   }
> > }
> > //-----------------------------------------------------------------------------
> > unsigned int FGMRES::iterator(Vector& x, const Vector& b, Vector& r)
> > {
> >   // Number of unknowns
> >   unsigned int n = x.size();
> >   
> >   // Initialization of temporary variables 
> >   Matrix mat_h(maxiter+1,maxiter+1, Matrix::dense);
> >   Matrix mat_v(1,n, Matrix::dense);
> >   Matrix mat_z(1,n, Matrix::dense);
> > 
> >   Vector vec_s(maxiter+1);  
> >   Vector vec_c(maxiter+1); 
> >   Vector vec_g(maxiter+1); 
> >   Vector vec_z(n);
> >   Vector vec_v(n);
> > 
> >   real tmp1, tmp2, nu;
> >   real htmp = 0.0;
> >   
> >   // Compute start residual = b - Ax.
> >   real rnorm = residual(A,x,b,r);
> > 
> >   // Set start values for v. 
> >   for (unsigned int i = 0; i < n; i++)
> >     mat_v[0][i] = vec_v(i) = r(i)/rnorm;
> >   
> >   vec_g = 0.0;
> >   vec_g(0) = rnorm;
> > 
> >   unsigned int k = 0;
> >   unsigned int k_end = 0;
> >   
> >   for (k = 0; k < maxiter; k++) {
> > 
> >     // Use the modified Gram-Schmidt process with reorthogonalization to
> >     // find orthogonal Krylov vectors.
> >     if (k != 0) mat_v.addrow(vec_v);
> >     
> >     pc.solve(vec_z,vec_v);
> >     A.mult(vec_z,vec_v);
> >          
> >     if (k == 0) 
> >       for (unsigned int i=0;i<n;i++) mat_z[0][i] = vec_z(i); 
> >     else
> >       mat_z.addrow(vec_z);
> >     
> >     // Modified Gram-Schmit orthogonalization. 
> >     for (unsigned int j = 0; j < k+1; j++) {
> >       htmp = 0.0;
> >       for (unsigned int i = 0; i < n; i++) htmp += vec_v(i) * mat_v[j][i];
> >       mat_h[j][k] = htmp;
> >       for (unsigned int i = 0; i < n; i++) vec_v(i) -= htmp * mat_v[j][i];
> >     }
> >     mat_h[k+1][k] = vec_v.norm();
> >     
> >     // Test for non orthogonality and reorthogonalise if necessary.
> >     if ( reorthog(mat_v,vec_v,k) ){ 
> >       for (unsigned int j = 0; j < k+1; j++) {
> > 	for (unsigned int i = 0; i < n; i++) htmp = vec_v(i)*mat_v[j][i];
> > 	mat_h(j,k) += htmp;
> > 	for (unsigned int i = 0; i < n; i++) vec_v(i) -= htmp*mat_v[j][i];
> >       }	  
> >       mat_h[k+1][k] = vec_v.norm();
> >     }
> >     
> >     // Normalize
> >     vec_v *= 1.0/mat_h[k+1][k]; 
> > 	 
> >     // If k > 0, solve least squares problem using QR factorization
> >     // and Givens rotations.  
> >     if (k>0) {
> >       for (unsigned int j = 0; j < k; j++) {
> > 	tmp1 = mat_h[j][k];
> > 	tmp2 = mat_h[j+1][k];
> > 	mat_h[j][k]   = vec_c(j)*tmp1 - vec_s(j)*tmp2;
> > 	mat_h[j+1][k] = vec_s(j)*tmp1 + vec_c(j)*tmp2; 
> >       }
> >     }
> >     
> >     nu = sqrt( sqr(mat_h[k][k]) + sqr(mat_h[k+1][k]) );
> >     
> >     vec_c(k) =   mat_h[k][k]   / nu;
> >     vec_s(k) = - mat_h[k+1][k] / nu;
> >     
> >     mat_h[k][k] = vec_c(k)*mat_h[k][k] - vec_s(k)*mat_h[k+1][k];
> >     mat_h[k+1][k] = 0.0;
> >     
> >     tmp1 = vec_c(k)*vec_g(k) - vec_s(k)*vec_g(k+1);
> >     vec_g(k+1) = vec_s(k)*vec_g(k) + vec_c(k)*vec_g(k+1);
> >     vec_g(k) = tmp1;
> >     
> >     // If residual rho = Ax-b less than tolerance (normalized with ||b||)
> >     // we are done.  
> > 
> >     if ( (fabs(vec_g(k+1)) < tol*bnorm) || (k == maxiter-1) ){
> >       k_end = k;
> >       break;
> >     }
> >   }
> >   k = k_end;
> > 
> >   // Postprocess to obtain solution by solving system Ry=w.
> >   Vector vec_w(vec_g);
> >   Vector vec_y(k_end+1); 
> >   
> >   // solve triangular system ry = w
> >   vec_y(k) = vec_w(k)/mat_h[k][k];
> >   for (unsigned int i = 1; i < k+1; i++) {
> >     vec_y(k-i) = vec_w(k-i);
> >     for (unsigned int j = 0; j < i; j++) vec_y(k-i) -= mat_h[k-i][k-j]*vec_y(k-j);
> >     vec_y(k-i) /= mat_h[k-i][k-i];
> >   }
> > 
> >   vec_z = 0.0;
> >   mat_z.multt(vec_y,vec_z);
> >   x += vec_z;
> >   
> >   return k_end;
> > }
> > //-----------------------------------------------------------------------------
> > bool FGMRES::reorthog(Matrix& v, Vector &x, int k)
> > {
> >   // reorthogonalize if ||Av_k||+delta*||v_k+1||=||Av_k|| to working precision 
> >   // with  delta \approx 10^-3
> >   Vector Av(v.size(1));
> >   Vector w(v.size(1));
> >  
> >   w = v(k,all);
> > 
> >   A.mult(w,Av);
> >   
> >   real delta = 1.0e-3;
> >   real Avnorm = Av.norm();
> > 
> >   return ((Avnorm + delta*x.norm()) - Avnorm) < DOLFIN_EPS;
> > }
> > //-----------------------------------------------------------------------------
> Index: src/kernel/la/GaussSeidel.cpp
> ===================================================================
> RCS file: /cvs/dolfin/src/kernel/la/GaussSeidel.cpp,v
> retrieving revision 1.1
> diff -r1.1 GaussSeidel.cpp
> 12c12
> <   : Preconditioner(), A(A), tol(tol), maxiter(maxiter)
> >   : Preconditioner(), A(A)
> 14c14,15
> <   // Do nothing
> >   this -> tol = tol;
> >   this -> maxiter = maxiter;
> 60a62
> > 
> Index: src/kernel/la/LinearSolver.cpp
> ===================================================================
> RCS file: /cvs/dolfin/src/kernel/la/LinearSolver.cpp,v
> retrieving revision 1.2
> diff -r1.2 LinearSolver.cpp
> 33a34,50
> > real LinearSolver::residual(const Matrix& A, Vector& x, const Vector& b, 
> > 			    Vector& r) const
> > {
> >   int n = b.size();
> > 
> >   r.init(n);
> >   A.mult(x,r);
> > 
> >   real sum = 0.0;
> >   for (int i = 0; i < n; i++) {
> >     r(i) = b(i) - r(i);
> >     sum += r(i)*r(i);
> >   }
> >   
> >   return sqrt(sum);
> > }
> > //-----------------------------------------------------------------------------
> Index: src/kernel/la/dolfin/FGMRES.h
> ===================================================================
> RCS file: /cvs/dolfin/src/kernel/la/dolfin/FGMRES.h,v
> retrieving revision 1.1
> diff -r1.1 FGMRES.h
> 16,17d15
> < 
> <   /// This is just a template. Write documentation here.
> 19,24c17,33
> <   class FGMRES : public Preconditioner, public LinearSolver
> <   {
> <   public:
> < 
> <     /// Create GMRES preconditioner/solver for a given matrix
> <     FGMRES(const Matrix& A, real tol, unsigned int maxiter);
> >   /// This is a flexible GMRES with right preconditioning,
> >   /// due to Y. Saad, SIAM J. Sci. Comput., 14 (1993), 461-469.
> >     
> >     class FGMRES : public Preconditioner, public LinearSolver
> >     {
> >     public:
> > 	
> >     /// Creates a FGMRES solver for a given matrix. Return an error 
> >     /// message if convergence is not obtained.
> >     FGMRES(const Matrix& A, unsigned int restarts, unsigned int maxiter, 
> > 	   real tol, Preconditioner& pc);
> > 
> >     /// Creates a FGMRES solver/preconditioner for a given matrix.
> >     /// Does not return an error message if convergence is not 
> >     /// obtained. Recommended for preconditioning.
> >     FGMRES(const Matrix& A, unsigned int restarts, real tol, 
> > 	   Preconditioner& pc);
> 28c37
> < 
> > 	  
> 31c40
> < 
> > 	    
> 34,38c43,54
> < 		      real tol, unsigned int maxiter);
> < 
> <   private:
> <     
> <     // The matrix
> > 		      unsigned int restarts, unsigned int maxiter, real tol,
> > 		      Preconditioner& pc);
> >     	    
> >     private:
> > 
> >     // Main FGMRES iterator.
> >     unsigned int iterator(Vector& x, const Vector& b, Vector& r);
> > 
> >     // Reorthogonalize.
> >     bool reorthog(Matrix& v, Vector &x, int k);
> > 	      
> >     // The matrix.
> 40,42d55
> <     
> <     // Tolerance
> <     real tol;
> 44c57,60
> <     // Maximum number of iterations
> >     // Maximum number of restarts.
> >     unsigned int restarts;
> > 
> >     // Number of iterations per restart.
> 47c63,64
> <   };
> >     // Tolerance
> >     real tol;  
> 48a66,75
> >     // The preconditioner.
> >     Preconditioner& pc;
> > 
> >     // Switch, solve to convergence. 
> >     bool solve2convergence; 
> > 
> >     // Norm of right hand side.
> >     real bnorm;
> >       
> >    };
> 51a79
> > 
> Index: src/kernel/la/dolfin/GaussSeidel.h
> ===================================================================
> RCS file: /cvs/dolfin/src/kernel/la/dolfin/GaussSeidel.h,v
> retrieving revision 1.1
> diff -r1.1 GaussSeidel.h
> 29a30,32
> >     /// Create Gauss-Seidel preconditioner/solver for a given matrix
> >     GaussSeidel(const Matrix& A, unsigned int maxiter);
> > 
> Index: src/kernel/la/dolfin/LinearSolver.h
> ===================================================================
> RCS file: /cvs/dolfin/src/kernel/la/dolfin/LinearSolver.h,v
> retrieving revision 1.1
> diff -r1.1 LinearSolver.h
> 28a29,31
> >     /// Compute l2-norm of residual and the residual vector
> >     real residual(const Matrix& A, Vector& x, const Vector& b, 
> > 		  Vector& r) const;
References

patch till dolfin
From: Erik Svensson, 2004-07-01