Skip to content
Snippets Groups Projects
Commit 1ee4aec2 authored by Oliver Sander's avatar Oliver Sander
Browse files

Removed most of the template metaprogramming code. This simplifies the code, simplifies 

some compiler warnings, and makes FieldVectors of size larger than 500 possible.

I left in the TMP for the norms of FieldVectors, because you can't remove it without
introducing a temporary, and I didn't want to do that without measuring the effects.

Also:  FieldVector::print() is deprecated now.  Please use operator<< instead.

[[Imported from SVN: r4695]]
parent 79c548ef
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
common/fmatrix.hh
+ 38
− 278
View file @ 1ee4aec2
......@@ -30,201 +30,6 @@ namespace Dune {
/** @brief Error thrown if operations of a FieldMatrix fail. */
class FMatrixError : public Exception {};
// template meta program for assignment from scalar
template<int I>
struct fmmeta_assignscalar {
template<class T, class K>
static void assignscalar (T* x, const K& k)
{
fmmeta_assignscalar<I-1>::assignscalar(x,k);
x[I] = k;
}
};
template<>
struct fmmeta_assignscalar<0> {
template<class T, class K>
static void assignscalar (T* x, const K& k)
{
x[0] = k;
}
};
// template meta program for operator+=
template<int I>
struct fmmeta_plusequal {
template<class T>
static void plusequal (T& x, const T& y)
{
x[I] += y[I];
fmmeta_plusequal<I-1>::plusequal(x,y);
}
};
template<>
struct fmmeta_plusequal<0> {
template<class T>
static void plusequal (T& x, const T& y)
{
x[0] += y[0];
}
};
// template meta program for operator-=
template<int I>
struct fmmeta_minusequal {
template<class T>
static void minusequal (T& x, const T& y)
{
x[I] -= y[I];
fmmeta_minusequal<I-1>::minusequal(x,y);
}
};
template<>
struct fmmeta_minusequal<0> {
template<class T>
static void minusequal (T& x, const T& y)
{
x[0] -= y[0];
}
};
// template meta program for operator*=
template<int I>
struct fmmeta_multequal {
template<class T, class K>
static void multequal (T& x, const K& k)
{
x[I] *= k;
fmmeta_multequal<I-1>::multequal(x,k);
}
};
template<>
struct fmmeta_multequal<0> {
template<class T, class K>
static void multequal (T& x, const K& k)
{
x[0] *= k;
}
};
// template meta program for operator/=
template<int I>
struct fmmeta_divequal {
template<class T, class K>
static void divequal (T& x, const K& k)
{
x[I] /= k;
fmmeta_divequal<I-1>::divequal(x,k);
}
};
template<>
struct fmmeta_divequal<0> {
template<class T, class K>
static void divequal (T& x, const K& k)
{
x[0] /= k;
}
};
// template meta program for dot
template<int I>
struct fmmeta_dot {
template<class X, class Y, class K>
static K dot (const X& x, const Y& y)
{
return x[I]*y[I] + fmmeta_dot<I-1>::template dot<X,Y,K>(x,y);
}
};
template<>
struct fmmeta_dot<0> {
template<class X, class Y, class K>
static K dot (const X& x, const Y& y)
{
return x[0]*y[0];
}
};
// template meta program for umv(x,y)
template<int I>
struct fmmeta_umv {
template<class Mat, class X, class Y, int c>
static void umv (const Mat& A, const X& x, Y& y)
{
typedef typename Mat::row_type R;
typedef typename Mat::field_type K;
y[I] += fmmeta_dot<c>::template dot<R,X,K>(A[I],x);
fmmeta_umv<I-1>::template umv<Mat,X,Y,c>(A,x,y);
}
};
template<>
struct fmmeta_umv<0> {
template<class Mat, class X, class Y, int c>
static void umv (const Mat& A, const X& x, Y& y)
{
typedef typename Mat::row_type R;
typedef typename Mat::field_type K;
y[0] += fmmeta_dot<c>::template dot<R,X,K>(A[0],x);
}
};
// template meta program for mmv(x,y)
template<int I>
struct fmmeta_mmv {
template<class Mat, class X, class Y, int c>
static void mmv (const Mat& A, const X& x, Y& y)
{
typedef typename Mat::row_type R;
typedef typename Mat::field_type K;
y[I] -= fmmeta_dot<c>::template dot<R,X,K>(A[I],x);
fmmeta_mmv<I-1>::template mmv<Mat,X,Y,c>(A,x,y);
}
};
template<>
struct fmmeta_mmv<0> {
template<class Mat, class X, class Y, int c>
static void mmv (const Mat& A, const X& x, Y& y)
{
typedef typename Mat::row_type R;
typedef typename Mat::field_type K;
y[0] -= fmmeta_dot<c>::template dot<R,X,K>(A[0],x);
}
};
template<class K, int n, int m, class X, class Y>
inline void fm_mmv (const FieldMatrix<K,n,m>& A, const X& x, Y& y)
{
for (int i=0; i<n; i++)
for (int j=0; j<m; j++)
y[i] -= A[i][j]*x[j];
}
template<class K>
inline void fm_mmv (const FieldMatrix<K,1,1>& A, const FieldVector<K,1>& x, FieldVector<K,1>& y)
{
y[0] -= A[0][0]*x[0];
}
// template meta program for usmv(x,y)
template<int I>
struct fmmeta_usmv {
template<class Mat, class K, class X, class Y, int c>
static void usmv (const Mat& A, const K& alpha, const X& x, Y& y)
{
typedef typename Mat::row_type R;
y[I] += alpha*fmmeta_dot<c>::template dot<R,X,K>(A[I],x);
fmmeta_usmv<I-1>::template usmv<Mat,K,X,Y,c>(A,alpha,x,y);
}
};
template<>
struct fmmeta_usmv<0> {
template<class Mat, class K, class X, class Y, int c>
static void usmv (const Mat& A, const K& alpha, const X& x, Y& y)
{
typedef typename Mat::row_type R;
y[0] += alpha*fmmeta_dot<c>::template dot<R,X,K>(A[0],x);
}
};
// conjugate komplex does nothing for non-complex types
template<class K>
inline K fm_ck (const K& k)
......@@ -425,73 +230,6 @@ namespace Dune {
A[1][1] = temp*detinv;
}
//! left multiplication with matrix
template<class K, int n, int m>
void fm_leftmultiply (const FieldMatrix<K,n,n>& M, FieldMatrix<K,n,m>& A)
{
FieldMatrix<K,n,m> C(A);
for (int i=0; i<n; i++)
for (int j=0; j<m; j++)
{
A[i][j] = 0;
for (int k=0; k<n; k++)
A[i][j] += M[i][k]*C[k][j];
}
}
//! left multiplication with matrix, n=1
template<class K>
void fm_leftmultiply (const FieldMatrix<K,1,1>& M, FieldMatrix<K,1,1>& A)
{
A[0][0] *= M[0][0];
}
//! left multiplication with matrix, n=2
template<class K>
void fm_leftmultiply (const FieldMatrix<K,2,2>& M, FieldMatrix<K,2,2>& A)
{
FieldMatrix<K,2,2> C(A);
A[0][0] = M[0][0]*C[0][0] + M[0][1]*C[1][0];
A[0][1] = M[0][0]*C[0][1] + M[0][1]*C[1][1];
A[1][0] = M[1][0]*C[0][0] + M[1][1]*C[1][0];
A[1][1] = M[1][0]*C[0][1] + M[1][1]*C[1][1];
}
//! right multiplication with matrix
template<class K, int n, int m>
void fm_rightmultiply (const FieldMatrix<K,m,m>& M, FieldMatrix<K,n,m>& A)
{
FieldMatrix<K,n,m> C(A);
for (int i=0; i<n; i++)
for (int j=0; j<m; j++)
{
A[i][j] = 0;
for (int k=0; k<m; k++)
A[i][j] += C[i][k]*M[k][j];
}
}
//! right multiplication with matrix, n=1
template<class K>
void fm_rightmultiply (const FieldMatrix<K,1,1>& M, FieldMatrix<K,1,1>& A)
{
A[0][0] *= M[0][0];
}
//! right multiplication with matrix, n=2
template<class K>
void fm_rightmultiply (const FieldMatrix<K,2,2>& M, FieldMatrix<K,2,2>& A)
{
FieldMatrix<K,2,2> C(A);
A[0][0] = C[0][0]*M[0][0] + C[0][1]*M[1][0];
A[0][1] = C[0][0]*M[0][1] + C[0][1]*M[1][1];
A[1][0] = C[1][0]*M[0][0] + C[1][1]*M[1][0];
A[1][1] = C[1][0]*M[0][1] + C[1][1]*M[1][1];
}
/**
@brief A dense n x m matrix.
......@@ -645,7 +383,8 @@ namespace Dune {
//===== assignment from scalar
FieldMatrix& operator= (const K& k)
{
fmmeta_assignscalar<n-1>::assignscalar(p,k);
for (int i=0; i<n; i++)
p[i] = k;
return *this;
}
......@@ -654,28 +393,32 @@ namespace Dune {
//! vector space addition
FieldMatrix& operator+= (const FieldMatrix& y)
{
fmmeta_plusequal<n-1>::plusequal(*this,y);
for (int i=0; i<n; i++)
p[i] += y.p[i];
return *this;
}
//! vector space subtraction
FieldMatrix& operator-= (const FieldMatrix& y)
{
fmmeta_minusequal<n-1>::minusequal(*this,y);
for (int i=0; i<n; i++)
p[i] -= y.p[i];
return *this;
}
//! vector space multiplication with scalar
FieldMatrix& operator*= (const K& k)
{
fmmeta_multequal<n-1>::multequal(*this,k);
for (int i=0; i<n; i++)
p[i] *= k;
return *this;
}
//! vector space division by scalar
FieldMatrix& operator/= (const K& k)
{
fmmeta_divequal<n-1>::divequal(*this,k);
for (int i=0; i<n; i++)
p[i] /= k;
return *this;
}
......@@ -689,10 +432,9 @@ namespace Dune {
if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
#endif
fmmeta_umv<n-1>::template umv<FieldMatrix,X,Y,m-1>(*this,x,y);
// for (int i=0; i<n; i++)
// for (int j=0; j<m; j++)
// y[i] += (*this)[i][j] * x[j];
for (size_type i=0; i<n; i++)
for (size_type j=0; j<m; j++)
y[i] += (*this)[i][j] * x[j];
}
//! y += A^T x
......@@ -731,8 +473,9 @@ namespace Dune {
if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
#endif
fmmeta_mmv<n-1>::template mmv<FieldMatrix,X,Y,m-1>(*this,x,y);
//fm_mmv(*this,x,y);
for (size_type i=0; i<n; i++)
for (size_type j=0; j<m; j++)
y[i] -= (*this)[i][j] * x[j];
}
//! y -= A^T x
......@@ -771,7 +514,9 @@ namespace Dune {
if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
#endif
fmmeta_usmv<n-1>::template usmv<FieldMatrix,K,X,Y,m-1>(*this,alpha,x,y);
for (size_type i=0; i<n; i++)
for (size_type j=0; j<m; j++)
y[i] += alpha * (*this)[i][j] * x[j];
}
//! y += alpha A^T x
......@@ -863,14 +608,29 @@ namespace Dune {
//! Multiplies M from the left to this matrix
FieldMatrix& leftmultiply (const FieldMatrix<K,n,n>& M)
{
fm_leftmultiply(M,*this);
FieldMatrix<K,n,m> C(*this);
for (int i=0; i<n; i++)
for (int j=0; j<m; j++) {
(*this)[i][j] = 0;
for (int k=0; k<n; k++)
(*this)[i][j] += M[i][k]*C[k][j];
}
return *this;
}
//! Multiplies M from the right to this matrix
FieldMatrix& rightmultiply (const FieldMatrix<K,n,n>& M)
{
fm_rightmultiply(M,*this);
FieldMatrix<K,n,m> C(*this);
for (int i=0; i<n; i++)
for (int j=0; j<m; j++) {
(*this)[i][j] = 0;
for (int k=0; k<m; k++)
(*this)[i][j] += C[i][k]*M[k][j];
}
return *this;
}
......@@ -942,7 +702,7 @@ namespace Dune {
}
private:
// the data, very simply a built in array with rowwise ordering
// the data, very simply a built in array with row-wise ordering
row_type p[(n > 0) ? n : 1];
};
......@@ -1251,7 +1011,7 @@ namespace Dune {
//! infinity norm (row sum norm, how to generalize for blocks?)
double infinity_norm () const
{
return fvmeta_abs(a[0]);
return std::abs(a[0]);
}
//! simplified infinity norm (uses Manhattan norm for complex values)
......
common/fvector.hh
+ 37
− 189
View file @ 1ee4aec2
......@@ -32,165 +32,10 @@ namespace Dune {
// forward declaration of template
template<class K, int n> class FieldVector;
// template meta program for assignment from scalar
template<int I>
struct fvmeta_assignscalar {
template<class K, int n>
static K& assignscalar (FieldVector<K,n>& x, const K& k)
{
x[I] = fvmeta_assignscalar<I-1>::assignscalar(x,k);
return x[I];
}
};
template<>
struct fvmeta_assignscalar<0> {
template<class K, int n>
static K& assignscalar (FieldVector<K,n>& x, const K& k)
{
x[0] = k;
return x[0];
}
};
// template meta program for operator+=
template<int I>
struct fvmeta_plusequal {
template<class K, int n>
static void plusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
{
fvmeta_plusequal<I-1>::plusequal(x,y);
x[I] += y[I];
}
template<class K, int n>
static void plusequal (FieldVector<K,n>& x, const K& y)
{
fvmeta_plusequal<I-1>::plusequal(x,y);
x[I] += y;
}
};
template<>
struct fvmeta_plusequal<0> {
template<class K, int n>
static void plusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
{
x[0] += y[0];
}
template<class K, int n>
static void plusequal (FieldVector<K,n>& x, const K& y)
{
x[0] += y;
}
};
// template meta program for operator-=
template<int I>
struct fvmeta_minusequal {
template<class K, int n>
static void minusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
{
fvmeta_minusequal<I-1>::minusequal(x,y);
x[I] -= y[I];
}
template<class K, int n>
static void minusequal (FieldVector<K,n>& x, const K& y)
{
fvmeta_minusequal<I-1>::minusequal(x,y);
x[I] -= y;
}
};
template<>
struct fvmeta_minusequal<0> {
template<class K, int n>
static void minusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
{
x[0] -= y[0];
}
template<class K, int n>
static void minusequal (FieldVector<K,n>& x, const K& y)
{
x[0] -= y;
}
};
// template meta program for operator*=
template<int I>
struct fvmeta_multequal {
template<class K, int n>
static void multequal (FieldVector<K,n>& x, const K& k)
{
fvmeta_multequal<I-1>::multequal(x,k);
x[I] *= k;
}
};
template<>
struct fvmeta_multequal<0> {
template<class K, int n>
static void multequal (FieldVector<K,n>& x, const K& k)
{
x[0] *= k;
}
};
// template meta program for operator/=
template<int I>
struct fvmeta_divequal {
template<class K, int n>
static void divequal (FieldVector<K,n>& x, const K& k)
{
fvmeta_divequal<I-1>::divequal(x,k);
x[I] /= k;
}
};
template<>
struct fvmeta_divequal<0> {
template<class K, int n>
static void divequal (FieldVector<K,n>& x, const K& k)
{
x[0] /= k;
}
};
#endif
// template meta program for operator==
template<int I>
struct fvmeta_equality {
template<class K, int n>
static bool equality (const FieldVector<K,n>& x, const FieldVector<K,n>& y)
{
return x[I] == y[I] && fvmeta_equality<I-1>::equality(x,y);
}
};
template<>
struct fvmeta_equality<0> {
template<class K, int n>
static bool equality (const FieldVector<K,n>& x, const FieldVector<K,n>& y)
{
return x[0] == y[0];
}
};
#ifndef DUNE_EXPRESSIONTEMPLATES
// template meta program for axpy
template<int I>
struct fvmeta_axpy {
template<class K, int n>
static void axpy (FieldVector<K,n>& x, const K& a, const FieldVector<K,n>& y)
{
fvmeta_axpy<I-1>::axpy(x,a,y);
x[I] += a*y[I];
}
};
template<>
struct fvmeta_axpy<0> {
template<class K, int n>
static void axpy (FieldVector<K,n>& x, const K& a, const FieldVector<K,n>& y)
{
x[0] += a*y[0];
}
};
// template meta program for dot
template<int I>
struct fvmeta_dot {
......@@ -209,23 +54,10 @@ namespace Dune {
}
};
// some abs functions needed for the norms
template<class K>
inline double fvmeta_abs (const K& k)
{
return (k >= 0) ? k : -k;
}
template<class K>
inline double fvmeta_abs (const std::complex<K>& c)
{
return sqrt(c.real()*c.real() + c.imag()*c.imag());
}
template<class K>
inline double fvmeta_absreal (const K& k)
{
return fvmeta_abs(k);
return std::abs(k);
}
template<class K>
......@@ -252,7 +84,7 @@ namespace Dune {
template<class K, int n>
static double one_norm (const FieldVector<K,n>& x)
{
return fvmeta_abs(x[I]) + fvmeta_one_norm<I-1>::one_norm(x);
return std::abs(x[I]) + fvmeta_one_norm<I-1>::one_norm(x);
}
};
template<>
......@@ -260,7 +92,7 @@ namespace Dune {
template<class K, int n>
static double one_norm (const FieldVector<K,n>& x)
{
return fvmeta_abs(x[0]);
return std::abs(x[0]);
}
};
......@@ -306,7 +138,7 @@ namespace Dune {
template<class K, int n>
static double infinity_norm (const FieldVector<K,n>& x)
{
return std::max(fvmeta_abs(x[I]),fvmeta_infinity_norm<I-1>::infinity_norm(x));
return std::max(std::abs(x[I]),fvmeta_infinity_norm<I-1>::infinity_norm(x));
}
};
template<>
......@@ -314,7 +146,7 @@ namespace Dune {
template<class K, int n>
static double infinity_norm (const FieldVector<K,n>& x)
{
return fvmeta_abs(x[0]);
return std::abs(x[0]);
}
};
......@@ -634,7 +466,9 @@ namespace Dune {
//! Assignment operator for scalar
FieldVector& operator= (const K& k)
{
fvmeta_assignscalar<SIZE-1>::assignscalar(*this,k);
//fvmeta_assignscalar<SIZE-1>::assignscalar(*this,k);
for (size_type i=0; i<SIZE; i++)
p[i] = k;
return *this;
}
......@@ -787,14 +621,16 @@ namespace Dune {
//! vector space addition
FieldVector& operator+= (const FieldVector& y)
{
fvmeta_plusequal<SIZE-1>::plusequal(*this,y);
for (size_type i=0; i<SIZE; i++)
p[i] += y.p[i];
return *this;
}
//! vector space subtraction
FieldVector& operator-= (const FieldVector& y)
{
fvmeta_minusequal<SIZE-1>::minusequal(*this,y);
for (size_type i=0; i<SIZE; i++)
p[i] -= y.p[i];
return *this;
}
......@@ -815,28 +651,32 @@ namespace Dune {
//! vector space add scalar to all comps
FieldVector& operator+= (const K& k)
{
fvmeta_plusequal<SIZE-1>::plusequal(*this,k);
for (size_type i=0; i<SIZE; i++)
p[i] += k;
return *this;
}
//! vector space subtract scalar from all comps
FieldVector& operator-= (const K& k)
{
fvmeta_minusequal<SIZE-1>::minusequal(*this,k);
for (size_type i=0; i<SIZE; i++)
p[i] -= k;
return *this;
}
//! vector space multiplication with scalar
FieldVector& operator*= (const K& k)
{
fvmeta_multequal<SIZE-1>::multequal(*this,k);
for (size_type i=0; i<SIZE; i++)
p[i] *= k;
return *this;
}
//! vector space division by scalar
FieldVector& operator/= (const K& k)
{
fvmeta_divequal<SIZE-1>::divequal(*this,k);
for (size_type i=0; i<SIZE; i++)
p[i] /= k;
return *this;
}
......@@ -845,14 +685,19 @@ namespace Dune {
//! Binary vector comparison
bool operator== (const FieldVector& y) const
{
return fvmeta_equality<SIZE-1>::equality(*this,y);
for (size_type i=0; i<SIZE; i++)
if (p[i]!=y.p[i])
return false;
return true;
}
//! vector space axpy operation
FieldVector& axpy (const K& a, const FieldVector& y)
{
#ifndef DUNE_EXPRESSIONTEMPLATES
fvmeta_axpy<SIZE-1>::axpy(*this,a,y);
for (size_type i=0; i<SIZE; i++)
p[i] += a*y.p[i];
#else
*this += a*y;
#endif
......@@ -922,8 +767,9 @@ namespace Dune {
return SIZE;
}
//! Send vector to output stream
void print (std::ostream& s) const
/** \brief Send vector to output stream
\deprecated Use operator << instead */
void print (std::ostream& s) const DUNE_DEPRECATED
{
for (size_type i=0; i<SIZE; i++)
if (i>0)
......@@ -941,7 +787,8 @@ namespace Dune {
template<typename K, int n>
std::ostream& operator<< (std::ostream& s, const FieldVector<K,n>& v)
{
v.print(s);
for (typename FieldVector<K,n>::size_type i=0; i<n; i++)
s << ((i>0) ? " " : "") << v[i];
return s;
}
......@@ -1213,7 +1060,7 @@ namespace Dune {
//! one norm (sum over absolute values of entries)
double one_norm () const
{
return fvmeta_abs(p);
return std::abs(p);
}
//! simplified one norm (uses Manhattan norm for complex values)
......@@ -1237,7 +1084,7 @@ namespace Dune {
//! infinity norm (maximum of absolute values of entries)
double infinity_norm () const
{
return fvmeta_abs(p);
return std::abs(p);
}
//! simplified infinity norm (uses Manhattan norm for complex values)
......@@ -1261,8 +1108,9 @@ namespace Dune {
return 1;
}
//! Send vector to output stream
void print (std::ostream& s) const
/** \brief Send vector to output stream
\deprecated Use operator << instead */
void print (std::ostream& s) const DUNE_DEPRECATED
{
s << p;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment