dolfin team mailing list archive

Thread
Date

patch till dolfin

To: dolfin-dev@xxxxxxxxxxxxxxxxxxxxxx
From: Erik Svensson <eriksv@xxxxxxxxxxxxxxxx>
Date: Thu, 01 Jul 2004 12:22:25 +0200
Cc: hoffman@xxxxxxxxxxxx
Reply-to: Discussion of DOLFIN development <dolfin-dev@xxxxxxxxxxxxxxxxxxxxxx>
User-agent: Mozilla/5.0 (X11; U; SunOS sun4u; en-US; rv:1.6) Gecko/20040116

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;

Follow ups

Re: patch till dolfin
From: Anders Logg, 2004-07-01