// $Id$
// Version: $Name$
//
-// Copyright (C) 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+// Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
#include <lac/full_matrix.h>
#include <lac/matrix_lib.h>
+#include <lac/householder.h>
#include <multigrid/mg_base.h>
DEAL_II_NAMESPACE_OPEN
const VECTOR &src) const;
private:
- /**
- * Pointer to the coarse grid
- * matrix.
- */
- SmartPointer<const FullMatrix<number> > matrix;
-
/**
* Matrix for QR-factorization.
*/
- mutable FullMatrix<number> work;
-
- /**
- * Auxiliary vector.
- */
- mutable VECTOR aux;
+ Householder<number> householder;
};
/*@}*/
template<typename number, class VECTOR>
MGCoarseGridHouseholder<number, VECTOR>::MGCoarseGridHouseholder(
const FullMatrix<number>* A)
- :
- matrix(A, typeid(*this).name())
-{}
+{
+ if (A != 0) householder.initialize(*A);
+}
MGCoarseGridHouseholder<number, VECTOR>::initialize(
const FullMatrix<number>& A)
{
- matrix = &A;
+ householder.initialize(A);
}
template<typename number, class VECTOR>
void
MGCoarseGridHouseholder<number, VECTOR>::clear()
-{
- matrix = 0;
-}
+{}
VECTOR &dst,
const VECTOR &src) const
{
- work = *matrix;
- aux = src;
- work.least_squares(dst, aux);
+ householder.least_squares(dst, src);
}
#endif // DOXYGEN
// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* the vector space.
*/
number frobenius_norm () const;
-
- /**
- * @deprecated Old name for
- * FullMatrix::frobenius_norm().
- */
- number norm2 () const;
-
+
/**
* Compute the relative norm of
* the skew-symmetric part. The
const FullMatrix<number2> &B,
const number c,
const FullMatrix<number2> &C);
-
- /**
- * @deprecated Simple addition of
- * a scaled matrix,
- * i.e. <tt>*this += a*A</tt>.
- *
- * This function is
- * deprecated. Use <tt>add</tt>
- * instead, since this has the
- * same interface as the other
- * matrix and vector classes in
- * the library.
- */
- template <typename number2>
- void add_scaled (const number a,
- const FullMatrix<number2> &A);
-
+
/**
* Add rectangular block.
*
void invert (const FullMatrix<number2> &M);
- /**
- * @deprecated Use the class Householder
- * to compute a QR-factorization which
- * can be applied to several vectors.
- *
- * QR-factorization of a matrix.
- * The orthogonal transformation
- * Q is applied to the vector y
- * and this matrix.
- *
- * After execution of
- * householder, the upper
- * triangle contains the
- * resulting matrix R, the lower
- * the incomplete factorization
- * matrices.
- */
- template<typename number2>
- void householder (Vector<number2> &y);
-
//@}
///@name Multiplications
//@{
/**
* Forward elimination of lower
* triangle. Inverts the lower
- * triangle of a quadratic matrix
+ * triangle of a rectangular matrix
* for a given right hand side.
*
* If the matrix has more columns
* rows, the upper quadratic part
* of the matrix is considered.
*
- * Note that this function does
- * not fit into this class at
- * all, since it assumes that the
- * elements of this object do not
- * represent a matrix, but rather
- * a decomposition into two
- * factors. Therefore, if this
- * assumption holds, all
- * functions like multiplication
- * by matrices or vectors, norms,
- * etc, have no meaning any
- * more. Conversely, if these
- * functions have a meaning on
- * this object, then the
- * forward() function has no
- * meaning. This bifacial
- * property of this class is
- * probably a design mistake and
- * may once go away by separating
- * the forward() and backward()
- * functions into a class of
- * their own.
+ * @note It is safe to use the
+ * same object for @p dst and @p
+ * src.
*/
template<typename number2>
void forward (Vector<number2> &dst,
* triangle.
*
* See forward()
+ *
+ * @note It is safe to use the
+ * same object for @p dst and @p
+ * src.
*/
template<typename number2>
void backward (Vector<number2> &dst,
const Vector<number2> &src) const;
- /**
- * @deprecated Use the class
- * Householder to solve least
- * squares problems.
- *
- * Least-Squares-Approximation by
- * QR-factorization. The return
- * value is the Euclidean norm of
- * the approximation error.
- */
- template<typename number2>
- double least_squares (Vector<number2> &dst,
- Vector<number2> &src);
//@}
/** @addtogroup Exceptions
// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
-template <typename number>
-template <typename number2>
-void
-FullMatrix<number>::add_scaled (const number s,
- const FullMatrix<number2> &src)
-{
- // this function is
- // deprecated. forward to the other
- // one
- add (s, src);
-}
-
-
-
template <typename number>
template <typename number2>
void
-template <typename number>
-number
-FullMatrix<number>::norm2 () const
-{
- return frobenius_norm();
-}
-
-
-
template <typename number>
number
FullMatrix<number>::relative_symmetry_norm2 () const
}
}
-// QR-transformation cf. Stoer 1 4.8.2 (p. 191)
-
-template <typename number>
-template <typename number2>
-void
-FullMatrix<number>::householder(Vector<number2>& src)
-{
- Assert (!this->empty(), ExcEmptyMatrix());
-
- // m > n, src.n() = m
- Assert (this->n_cols() <= this->n_rows(),
- ExcDimensionMismatch(this->n_cols(), this->n_rows()));
- Assert (src.size() == this->n_rows(),
- ExcDimensionMismatch(src.size(), this->n_rows()));
-
- for (unsigned int j=0 ; j<n() ; ++j)
- {
- number2 sigma = 0;
- unsigned int i;
- for (i=j ; i<m() ; ++i) sigma += this->el(i,j)*this->el(i,j);
- if (std::fabs(sigma) < 1.e-15) return;
- number2 s = this->el(j,j);
- s = (s<0) ? std::sqrt(sigma) : -std::sqrt(sigma);
- number2 dj = s;
-
- number2 beta = 1./(s*this->el(j,j)-sigma);
- this->el(j,j) -= s;
-
- for (unsigned int k=j+1 ; k<n() ; ++k)
- {
- number2 sum = 0.;
- for (i=j ; i<m() ; ++i) sum += this->el(i,j)*this->el(i,k);
- sum *= beta;
-
- for (i=j ; i<m() ; ++i) this->el(i,k) += sum*this->el(i,j);
- }
-
- number2 sum = 0.;
- for (i=j ; i<m() ; ++i) sum += this->el(i,j)*src(i);
- sum *= beta;
-
- for (i=j ; i<m() ; ++i) src(i) += sum*this->el(i,j);
- this->el(j,j) = dj;
- }
-}
-
-
-template <typename number>
-template <typename number2>
-double
-FullMatrix<number>::least_squares (Vector<number2>& dst,
- Vector<number2>& src)
-{
- Assert (!this->empty(), ExcEmptyMatrix());
-
- // m > n, m = src.n, n = dst.n
-
- householder(src);
- backward(dst, src);
-
- number2 sum = 0.;
- for (unsigned int i=n() ; i<m() ; ++i) sum += src(i) * src(i);
- return std::sqrt(sum);
-}
-
template <typename number>
class Householder : private FullMatrix<number>
{
public:
+ /**
+ * Create an empty object.
+ */
+ Householder ();
+
/**
* Create an object holding the
* QR-decomposition of a matrix.
template<typename number2>
Householder (const FullMatrix<number2>&);
+ /**
+ * Compute the QR-decomposition
+ * of another matrix.
+ */
+ template<typename number2>
+ void
+ initialize (const FullMatrix<number2>&);
+
/**
* Solve the least-squares
* problem for the right hand
*/
template<typename number2>
double least_squares (Vector<number2> &dst,
- Vector<number2> &src) const;
+ const Vector<number2> &src) const;
private:
/**
// QR-transformation cf. Stoer 1 4.8.2 (p. 191)
+template <typename number>
+Householder<number>::Householder()
+{}
+
+
template <typename number>
template <typename number2>
-Householder<number>::Householder(const FullMatrix<number2>& M)
- :
- FullMatrix<number>(M),
- diagonal(M.n_rows())
+void
+Householder<number>::initialize(const FullMatrix<number2>& M)
{
-// Assert (!this->empty(), ExcEmptyMatrix());
-
+ this->reinit(M.n_rows(), M.n_cols());
+ this->fill(M);
+ diagonal.resize(M.n_rows());
+ Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
// m > n, src.n() = m
Assert (this->n_cols() <= this->n_rows(),
ExcDimensionMismatch(this->n_cols(), this->n_rows()));
}
+template <typename number>
+template <typename number2>
+Householder<number>::Householder(const FullMatrix<number2>& M)
+{
+ initialize(M);
+}
+
+
template <typename number>
template <typename number2>
double
Householder<number>::least_squares (Vector<number2>& dst,
- Vector<number2>& src) const
+ const Vector<number2>& src) const
{
-// Assert (!this->empty(), ExcEmptyMatrix());
-
+ Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
+ dst = src;
// m > n, m = src.n, n = dst.n
// Multiply Q_n ... Q_2 Q_1 src
// Where Q_i = I-v_i v_i^T
for (unsigned int j=0;j<this->n();++j)
{
- // sum = v_i^T src
- number2 sum = diagonal[j]*src(j);
+ // sum = v_i^T dst
+ number2 sum = diagonal[j]*dst(j);
for (unsigned int i=j+1 ; i<this->m() ; ++i)
- sum += this->el(i,j)*src(i);
- // src -= v * sum
- src(j) -= sum*diagonal[j];
+ sum += this->el(i,j)*dst(i);
+ // dst -= v * sum
+ dst(j) -= sum*diagonal[j];
for (unsigned int i=j+1 ; i<this->m() ; ++i)
- src(i) -= sum*this->el(i,j);
+ dst(i) -= sum*this->el(i,j);
}
-
- this->backward(dst, src);
-
+ // Compute norm of residual
number2 sum = 0.;
for (unsigned int i=this->n() ; i<this->m() ; ++i)
- sum += src(i) * src(i);
+ sum += dst(i) * dst(i);
+ // Compute solution
+ this->backward(dst, dst);
+
return std::sqrt(sum);
}
#include <base/exceptions.h>
#include <base/logstream.h>
#include <base/memory_consumption.h>
+#include <lac/householder.h>
#include <lac/precondition_block.h>
#include <lac/vector.h>
#include <lac/full_matrix.h>
column_cell<this->blocksize; ++column_cell, ++column)
M_cell(row_cell,column_cell)=M(row,column);
}
- M_cell.householder(b_cell);
- M_cell.backward(x_cell,b_cell);
+ Householder<number> house(M_cell);
+ house.least_squares(x_cell,b_cell);
// distribute x_cell to dst
for (row=cell*this->blocksize, row_cell=0;
row_cell<this->blocksize;
}
else
{
- M_cell.householder(b_cell);
- M_cell.backward(x_cell,b_cell);
+ Householder<number> house(M_cell);
+ house.least_squares(x_cell,b_cell);
}
// distribute x_cell to dst
}
else
{
- M_cell.householder(b_cell);
- M_cell.backward(x_cell,b_cell);
+ Householder<number> house(M_cell);
+ house.least_squares(x_cell,b_cell);
}
#include <base/config.h>
#include <base/subscriptor.h>
#include <base/logstream.h>
+#include <lac/householder.h>
#include <lac/solver.h>
#include <lac/solver_control.h>
#include <lac/full_matrix.h>
y.reinit(j);
projected_rhs(0) = beta;
H1.fill(H);
-
- double res = H1.least_squares(y, projected_rhs);
+ Householder<double> house(H1);
+ double res = house.least_squares(y, projected_rhs);
iteration_state = this->control().check(++accumulated_iterations, res);
if (iteration_state != SolverControl::iterate)
break;
// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
const TYPEMAT, const FullMatrix<TYPEMAT2>&,
const TYPEMAT, const FullMatrix<TYPEMAT2>&);
-template void FullMatrix<TYPEMAT>::add_scaled<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (
const FullMatrix<TYPEMAT2>&, double, unsigned, unsigned, unsigned, unsigned);
template void FullMatrix<TYPEMAT>::Tadd<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
template void FullMatrix<TYPEMAT>::backward<TYPEVEC>(
Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
-template void FullMatrix<TYPEMAT>::householder<TYPEVEC>(Vector<TYPEVEC>&);
-template double FullMatrix<TYPEMAT>::least_squares<TYPEVEC>(
- Vector<TYPEVEC>&, Vector<TYPEVEC>&);
template
void FullMatrix<TYPEMAT>::precondition_Jacobi<TYPEVEC> (
Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
template void FullMatrix<TYPEMAT>::backward<TYPEVEC>(
Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
-template void FullMatrix<TYPEMAT>::householder<TYPEVEC>(Vector<TYPEVEC>&);
-template double FullMatrix<TYPEMAT>::least_squares<TYPEVEC>(
- Vector<TYPEVEC>&, Vector<TYPEVEC>&);
template
void FullMatrix<TYPEMAT>::precondition_Jacobi<TYPEVEC> (
Vector<TYPEVEC> &, const Vector<TYPEVEC> &, const TYPEMAT) const;
// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
const TYPEMAT, const FullMatrix<TYPEMAT2>&,
const TYPEMAT, const FullMatrix<TYPEMAT2>&);
-template void FullMatrix<TYPEMAT>::add_scaled<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (
const FullMatrix<TYPEMAT2>&, double, unsigned, unsigned, unsigned, unsigned);
template void FullMatrix<TYPEMAT>::Tadd<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
template void FullMatrix<TYPEMAT>::backward<TYPEVEC>(
Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
-template void FullMatrix<TYPEMAT>::householder<TYPEVEC>(Vector<TYPEVEC>&);
-template double FullMatrix<TYPEMAT>::least_squares<TYPEVEC>(
- Vector<TYPEVEC>&, Vector<TYPEVEC>&);
template
void FullMatrix<TYPEMAT>::precondition_Jacobi<TYPEVEC> (
Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
template void FullMatrix<TYPEMAT>::backward<TYPEVEC>(
Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
-template void FullMatrix<TYPEMAT>::householder<TYPEVEC>(Vector<TYPEVEC>&);
-template double FullMatrix<TYPEMAT>::least_squares<TYPEVEC>(
- Vector<TYPEVEC>&, Vector<TYPEVEC>&);
template
void FullMatrix<TYPEMAT>::precondition_Jacobi<TYPEVEC> (
Vector<TYPEVEC> &, const Vector<TYPEVEC> &, const TYPEMAT) const;
projects / dealii-svn.git / commitdiff