$args0 = $args;
$args0 =~ s/\*[^,]*,/\*,/g;
$args0 =~ s/\*[^,]*$/\*/g;
+ # First, do the general template None of these functions is
+ # implemented, but they allow us to link for instance with
+ # long double lapack support
+ my $numbers = 1;
+ my $argst = $args0;
+ my $typet = $type;
+ while ($argst =~ s/double/number$numbers/)
+ {
+ $numbers++;
+ }
+ $typet =~ s/double/number1/g;
+ $templates .= "\n\n//--------------------------------------------------------------------------------//\ntemplate<typename number1";
+ for (my $i=2;$i<$numbers;++$i)
+ {
+ $templates .= ", typename number$i";
+ }
+ $templates .= ">\ninline $typet\n$name ($argst)\n";
+ $templates .= "{\n Assert (false, ExcNotImplemented());\n}";
+
+ # The specialization for double. Note that we always have two
+ # versions, one implemented and calling LAPACK, the other not
+ # implemented
$templates .= "\n\n#ifdef HAVE_D$capname\_";
$templates .= "\ninline $type\n$name ($args)\n{\n d$name\_ ($args2);\n}\n";
$templates .= "#else\ninline $type\n$name ($args0)\n";
-namespace internal
-{
- // a namespace into which we import
- // the global definitions of the
- // LAPACK functions getrf for
- // various data types and add
- // general templates for all other
- // types. this allows to call the
- // getrf function for all data
- // types but generated an exception
- // if something is called for a
- // data type not supported by
- // LAPACK.
- namespace gemm_switch
- {
- using ::dealii::gemm;
-
- template <typename T, typename T2>
- void
- gemm (const char*, const char*, const int*, const int*, const int*,
- const T*, const T2*, const int*, const T*, const int*,
- const T*, T2*, const int*)
- {
- Assert (false, LAPACKSupport::ExcMissing("dgemm for this data type"));
- }
-
- }
-}
-
-
-
template <typename number>
template <typename number2>
void FullMatrix<number>::mmult (FullMatrix<number2> &dst,
// Use the BLAS function gemm for
// calculating the matrix-matrix
// product.
- internal::gemm_switch::gemm(notrans, notrans, &m, &n, &k, &alpha, &src(0,0), &m,
- &this->values[0], &k, &beta, &dst(0,0), &m);
+ gemm(notrans, notrans, &m, &n, &k, &alpha, &src(0,0), &m,
+ &this->values[0], &k, &beta, &dst(0,0), &m);
return;
}
// Use the BLAS function gemm for
// calculating the matrix-matrix
// product.
- internal::gemm_switch::gemm(notrans, trans, &m, &n, &k, &alpha, &src(0,0), &m,
- &this->values[0], &n, &beta, &dst(0,0), &m);
+ gemm(notrans, trans, &m, &n, &k, &alpha, &src(0,0), &m,
+ &this->values[0], &n, &beta, &dst(0,0), &m);
return;
}
// Use the BLAS function gemm for
// calculating the matrix-matrix
// product.
- internal::gemm_switch::gemm(trans, notrans, &m, &n, &k, &alpha, &src(0,0), &k,
- &this->values[0], &k, &beta, &dst(0,0), &m);
+ gemm(trans, notrans, &m, &n, &k, &alpha, &src(0,0), &k,
+ &this->values[0], &k, &beta, &dst(0,0), &m);
return;
}
// Use the BLAS function gemm for
// calculating the matrix-matrix
// product.
- internal::gemm_switch::gemm(trans, trans, &m, &n, &k, &alpha, &src(0,0), &k,
- &this->values[0], &n, &beta, &dst(0,0), &m);
+ gemm(trans, trans, &m, &n, &k, &alpha, &src(0,0), &k,
+ &this->values[0], &n, &beta, &dst(0,0), &m);
return;
}
}
-namespace internal
-{
- // a namespace into which we import
- // the global definitions of the
- // LAPACK functions getrf for
- // various data types and add
- // general templates for all other
- // types. this allows to call the
- // getrf function for all data
- // types but generated an exception
- // if something is called for a
- // data type not supported by
- // LAPACK.
- namespace getrf_switch
- {
- using ::dealii::getrf;
-
- template <typename T>
- void
- getrf (const int*, const int*, T*, const int*, int*, int*)
- {
- Assert (false, LAPACKSupport::ExcMissing("dgetrf for this data type"));
- }
-
-
- using ::dealii::getri;
-
- template <typename T>
- void
- getri (const int*, T*, const int*, int*, T*, const int*, int*)
- {
- Assert (false, LAPACKSupport::ExcMissing("dgetri for this data type"));
- }
- }
-}
-
-
-
template <typename number>
void
FullMatrix<number>::gauss_jordan ()
// Use the LAPACK function getrf for
// calculating the LU factorization.
- internal::getrf_switch::getrf(&nn, &nn, &this->values[0], &nn, &ipiv[0], &info);
+ getrf(&nn, &nn, &this->values[0], &nn, &ipiv[0], &info);
Assert(info >= 0, ExcInternalError());
Assert(info == 0, LACExceptions::ExcSingular());
// Use the LAPACK function getri for
// calculating the actual inverse using
// the LU factorization.
- internal::getrf_switch::getri(&nn, &this->values[0], &nn, &ipiv[0], &inv_work[0], &nn, &info);
+ getri(&nn, &this->values[0], &nn, &ipiv[0], &inv_work[0], &nn, &info);
Assert(info >= 0, ExcInternalError());
Assert(info == 0, LACExceptions::ExcSingular());
const number factor = 1.,
const bool transpose = false);
+
/**
* Matrix-vector-multiplication.
*
*
* Source and destination must
* not be the same vector.
+ *
+ * @note This template only
+ * exists for compile-time
+ * compatibility with
+ * FullMatrix. Implementation is
+ * only available for <tt>VECTOR=Vector<number></tt>
*/
- void vmult (Vector<number> &w,
- const Vector<number> &v,
- const bool adding=false) const;
+ template <class VECTOR>
+ void vmult(VECTOR& dst, const VECTOR& src, const bool adding = false) const;
/**
* Adding Matrix-vector-multiplication.
* <i>w += A*v</i>
*
* Source and destination must
* not be the same vector.
+ *
+ * @note This template only
+ * exists for compile-time
+ * compatibility with
+ * FullMatrix. Implementation is
+ * only available for <tt>VECTOR=Vector<number></tt>
*/
- void vmult_add (Vector<number> &w,
- const Vector<number> &v) const;
+ template <class VECTOR>
+ void vmult_add (VECTOR& w, const VECTOR& v) const;
/**
* Transpose
*
* Source and destination must
* not be the same vector.
- */
- void Tvmult (Vector<number> &w,
- const Vector<number> &v,
+ *
+ * @note This template only
+ * exists for compile-time
+ * compatibility with
+ * FullMatrix. Implementation is
+ * only available for <tt>VECTOR=Vector<number></tt>
+ */
+ template <class VECTOR>
+ void Tvmult (VECTOR& w, const VECTOR& v,
const bool adding=false) const;
/**
*
* Source and destination must
* not be the same vector.
- */
+ *
+ * @note This template only
+ * exists for compile-time
+ * compatibility with
+ * FullMatrix. Implementation is
+ * only available for <tt>VECTOR=Vector<number></tt>
+ */
+ template <class VECTOR>
+ void Tvmult_add (VECTOR& w, const VECTOR& v) const;
+
+ void vmult (Vector<number> &w,
+ const Vector<number> &v,
+ const bool adding=false) const;
+ void vmult_add (Vector<number> &w,
+ const Vector<number> &v) const;
+ void Tvmult (Vector<number> &w,
+ const Vector<number> &v,
+ const bool adding=false) const;
void Tvmult_add (Vector<number> &w,
const Vector<number> &v) const;
-
/**
* Compute the LU factorization
* of the matrix using LAPACK
template <typename number>
template <class MATRIX>
-void
+inline void
LAPACKFullMatrix<number>::copy_from (const MATRIX& M)
{
this->reinit (M.m(), M.n());
template <typename number>
template <class MATRIX>
-void
+inline void
LAPACKFullMatrix<number>::fill (
const MATRIX& M,
const unsigned int dst_offset_i,
template <typename number>
-std::complex<number>
+template <class VECTOR>
+inline void
+LAPACKFullMatrix<number>::vmult(VECTOR&, const VECTOR&, bool) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+
+template <typename number>
+template <class VECTOR>
+inline void
+LAPACKFullMatrix<number>::vmult_add(VECTOR&, const VECTOR&) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+
+template <typename number>
+template <class VECTOR>
+inline void
+LAPACKFullMatrix<number>::Tvmult(VECTOR&, const VECTOR&, bool) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+
+template <typename number>
+template <class VECTOR>
+inline void
+LAPACKFullMatrix<number>::Tvmult_add(VECTOR&, const VECTOR&) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+
+template <typename number>
+inline std::complex<number>
LAPACKFullMatrix<number>::eigenvalue (const unsigned int i) const
{
Assert (state & LAPACKSupport::eigenvalues, ExcInvalidState());
template <typename number>
-number
+inline number
LAPACKFullMatrix<number>::singular_value (const unsigned int i) const
{
Assert (state == LAPACKSupport::svd || state == LAPACKSupport::inverse_svd, LAPACKSupport::ExcState(state));
case PreconditionBlockBase<inverse_type>::householder:
this->inverse_householder(0).initialize(M_cell);
break;
+ case PreconditionBlockBase<inverse_type>::svd:
+ this->inverse_svd(0) = M_cell;
+ this->inverse_svd(0).compute_inverse_svd(0.);
+ break;
default:
Assert(false, ExcNotImplemented());
case PreconditionBlockBase<inverse_type>::householder:
this->inverse_householder(cell).initialize(M_cell);
break;
+ case PreconditionBlockBase<inverse_type>::svd:
+ this->inverse_svd(cell) = M_cell;
+ this->inverse_svd(cell).compute_inverse_svd(0.);
+ break;
default:
Assert(false, ExcNotImplemented());
#include <base/smartpointer.h>
#include <base/memory_consumption.h>
#include <lac/householder.h>
+#include <lac/lapack_full_matrix.h>
#include <vector>
* Use QR decomposition of
* the Householder class.
*/
- householder
+ householder,
+ /**
+ * Use the singular value
+ * decomposition of LAPACKFullMatrix.
+ */
+ svd
};
/**
*/
Householder<number>& inverse_householder (unsigned int i);
+ /**
+ * Access to the inverse diagonal
+ * blocks if Inversion is #householder.
+ */
+ LAPACKFullMatrix<number>& inverse_svd (unsigned int i);
+
/**
* Access to the inverse diagonal
* blocks.
* #gauss_jordan is used. Using
* <tt>number=float</tt> saves
* memory in comparison with
- * <tt>number=double</tt>.
+ * <tt>number=double</tt>, but
+ * may introduce numerical instability.
*/
std::vector<FullMatrix<number> > var_inverse_full;
* Storage of the inverse
* matrices of the diagonal
* blocks matrices as
- * <tt>Householder<number></tt>
+ * <tt>Householder</tt>
* matrices if Inversion
* #householder is used. Using
* <tt>number=float</tt> saves
* memory in comparison with
- * <tt>number=double</tt>.
+ * <tt>number=double</tt>, but
+ * may introduce numerical instability.
*/
std::vector<Householder<number> > var_inverse_householder;
+ /**
+ * Storage of the inverse
+ * matrices of the diagonal
+ * blocks matrices as
+ * <tt>LAPACKFullMatrix</tt>
+ * matrices if Inversion
+ * #svd is used. Using
+ * <tt>number=float</tt> saves
+ * memory in comparison with
+ * <tt>number=double</tt>, but
+ * may introduce numerical instability.
+ */
+ std::vector<LAPACKFullMatrix<number> > var_inverse_svd;
+
/**
* Storage of the original diagonal blocks.
*
{
if (var_inverse_full.size()!=0)
var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end());
- if (var_inverse_full.size()!=0)
+ if (var_inverse_householder.size()!=0)
var_inverse_householder.erase(var_inverse_householder.begin(), var_inverse_householder.end());
+ if (var_inverse_svd.size()!=0)
+ var_inverse_svd.erase(var_inverse_svd.begin(), var_inverse_svd.end());
if (var_diagonal.size()!=0)
var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
var_same_diagonal = false;
case householder:
var_inverse_householder.resize(1);
break;
+ case svd:
+ var_inverse_svd.resize(1);
+ var_inverse_svd[0].reinit(b,b);
+ break;
default:
Assert(false, ExcNotImplemented());
}
case householder:
var_inverse_householder.resize(n);
break;
+ case svd:
+ {
+ std::vector<LAPACKFullMatrix<number> >
+ tmp(n, LAPACKFullMatrix<number>(b));
+ var_inverse_svd.swap (tmp);
+ break;
+ }
default:
Assert(false, ExcNotImplemented());
}
AssertIndexRange (ii, var_inverse_householder.size());
var_inverse_householder[ii].vmult(dst, src);
break;
+ case svd:
+ AssertIndexRange (ii, var_inverse_svd.size());
+ var_inverse_svd[ii].vmult(dst, src);
+ break;
default:
Assert(false, ExcNotImplemented());
}
AssertIndexRange (ii, var_inverse_householder.size());
var_inverse_householder[ii].Tvmult(dst, src);
break;
+ case svd:
+ AssertIndexRange (ii, var_inverse_svd.size());
+ var_inverse_svd[ii].Tvmult(dst, src);
+ break;
default:
Assert(false, ExcNotImplemented());
}
}
+template <typename number>
+inline
+LAPACKFullMatrix<number>&
+PreconditionBlockBase<number>::inverse_svd(unsigned int i)
+{
+ if (same_diagonal())
+ return var_inverse_svd[0];
+
+ AssertIndexRange (i, var_inverse_svd.size());
+ return var_inverse_svd[i];
+}
+
+
template <typename number>
inline
FullMatrix<number>&
case PreconditionBlockBase<inverse_type>::householder:
this->inverse_householder(block).initialize(M_cell);
break;
+ case PreconditionBlockBase<inverse_type>::svd:
+ this->inverse_svd(block).reinit(bs, bs);
+ this->inverse_svd(block) = M_cell;
+ this->inverse_svd(block).compute_inverse_svd(0.);
+ break;
default:
Assert(false, ExcNotImplemented());
}
template LAPACKFullMatrix<float> &
LAPACKFullMatrix<float>::operator = (const FullMatrix<float>& M);
+template class LAPACKFullMatrix<long double>;
+template LAPACKFullMatrix<long double> &
+LAPACKFullMatrix<long double>::operator = (const FullMatrix<long double>& M);
+template LAPACKFullMatrix<long double> &
+LAPACKFullMatrix<long double>::operator = (const FullMatrix<double>& M);
+template LAPACKFullMatrix<long double> &
+LAPACKFullMatrix<long double>::operator = (const FullMatrix<float>& M);
+
template class PreconditionLU<double>;
template class PreconditionLU<float>;
projects / dealii.git / commitdiff