]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
remove functions in FullMatrix deprecated for about 2 years
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 21 Aug 2007 21:28:55 +0000 (21:28 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 21 Aug 2007 21:28:55 +0000 (21:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@15008 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/multigrid/mg_coarse.h
deal.II/lac/include/lac/full_matrix.h
deal.II/lac/include/lac/full_matrix.templates.h
deal.II/lac/include/lac/householder.h
deal.II/lac/include/lac/precondition_block.templates.h
deal.II/lac/include/lac/solver_gmres.h
deal.II/lac/source/full_matrix.double.cc
deal.II/lac/source/full_matrix.float.cc

index bf0c217ec72716056223c56b1e5bf4b313ad961e..3cc764de656458c78c6092d1eebc706118450cfb 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -16,6 +16,7 @@
 
 #include <lac/full_matrix.h>
 #include <lac/matrix_lib.h>
+#include <lac/householder.h>
 #include <multigrid/mg_base.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -156,21 +157,10 @@ class MGCoarseGridHouseholder : public MGCoarseGridBase<VECTOR>
                     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;
 };
 
 /*@}*/
@@ -274,9 +264,9 @@ MGCoarseGridLACIteration<SOLVER, VECTOR>
 template<typename number, class VECTOR>
 MGCoarseGridHouseholder<number, VECTOR>::MGCoarseGridHouseholder(
   const FullMatrix<number>* A)
-               :
-               matrix(A, typeid(*this).name())
-{}
+{
+  if (A != 0) householder.initialize(*A);
+}
 
 
 
@@ -285,7 +275,7 @@ void
 MGCoarseGridHouseholder<number, VECTOR>::initialize(
   const FullMatrix<number>& A)
 {
-  matrix = &A;
+  householder.initialize(A);
 }
 
 
@@ -293,9 +283,7 @@ MGCoarseGridHouseholder<number, VECTOR>::initialize(
 template<typename number, class VECTOR>
 void
 MGCoarseGridHouseholder<number, VECTOR>::clear()
-{
-  matrix = 0;
-}
+{}
 
 
 
@@ -306,9 +294,7 @@ MGCoarseGridHouseholder<number, VECTOR>::operator() (
   VECTOR       &dst,
   const VECTOR &src) const
 {
-  work = *matrix;
-  aux = src;
-  work.least_squares(dst, aux);
+  householder.least_squares(dst, src);
 }
 
 #endif // DOXYGEN
index 7d6e1cda8f5fbdc4adf48c17a81820bbff083b55..0848a37e6c332281d4873f578b71231cb268abfb 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -485,13 +485,7 @@ class FullMatrix : public Table<2,number>
                                      * 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
@@ -691,23 +685,7 @@ class FullMatrix : public Table<2,number>
              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.
                                      *
@@ -917,26 +895,6 @@ class FullMatrix : public Table<2,number>
     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
 //@{
@@ -1090,7 +1048,7 @@ class FullMatrix : public Table<2,number>
                                     /**
                                      * 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
@@ -1100,28 +1058,9 @@ class FullMatrix : public Table<2,number>
                                      * 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,
@@ -1132,24 +1071,15 @@ class FullMatrix : public Table<2,number>
                                      * 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
index 893aeb3401543649c5789797ff1ffca65a97214d..abaead2eca162d4d911a4175b1183c7877e6586e 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -1609,20 +1609,6 @@ void FullMatrix<number>::Tadd (const FullMatrix<number2> &src,
 
 
 
-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
@@ -1956,15 +1942,6 @@ FullMatrix<number>::frobenius_norm () const
 
 
 
-template <typename number>
-number
-FullMatrix<number>::norm2 () const
-{
-  return frobenius_norm();
-}
-
-
-
 template <typename number>
 number
 FullMatrix<number>::relative_symmetry_norm2 () const
@@ -2321,71 +2298,6 @@ FullMatrix<number>::gauss_jordan ()
     }
 }
 
-// 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>
index f44a4c44268f61326585e8f4e5333a972e761d0b..943548816da3423e2e5254b687086ca57d15385b 100644 (file)
@@ -51,6 +51,11 @@ template<typename number>
 class Householder : private FullMatrix<number>
 {
   public:
+                                    /**
+                                     * Create an empty object.
+                                     */
+    Householder ();
+
                                     /**
                                      * Create an object holding the
                                      * QR-decomposition of a matrix.
@@ -58,6 +63,14 @@ class Householder : private FullMatrix<number>
     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
@@ -77,7 +90,7 @@ class Householder : private FullMatrix<number>
                                      */
     template<typename number2>
     double least_squares (Vector<number2> &dst,
-                         Vector<number2> &src) const;
+                         const Vector<number2> &src) const;
 
   private:
                                     /**
@@ -95,15 +108,20 @@ class Householder : private FullMatrix<number>
 
 // 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()));
@@ -149,35 +167,44 @@ Householder<number>::Householder(const FullMatrix<number2>& M)
 }
 
 
+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);
 }
 
index fd56919f6af3c51f4a6cea432f46fb3b7c0d3fb6..255be19c3e1adc98690a2f4c0c9215fb7061ce1f 100644 (file)
@@ -18,6 +18,7 @@
 #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>
@@ -446,8 +447,8 @@ void PreconditionBlockJacobi<MATRIX,inverse_type>
                   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;
@@ -637,8 +638,8 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::forward (
        }
       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
@@ -751,8 +752,8 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::backward (
        }
       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);
        }
       
       
index ce8ad80c6362fe47f1a4611a802e59a1ebf173a7..9256bb7f08e05416429323d0b01ed17e364e5505 100644 (file)
@@ -18,6 +18,7 @@
 #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>
@@ -930,8 +931,8 @@ SolverFGMRES<VECTOR>::solve (
              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;
index 6a0ca614c144be770eeecbdbb0199aab5f6420d1..29c5e79095a38b2bdfac265f0ee62f58802fc284 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -39,7 +39,6 @@ template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatri
 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>&);
@@ -77,9 +76,6 @@ template void FullMatrix<TYPEMAT>::forward<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> (
@@ -106,9 +102,6 @@ template void FullMatrix<TYPEMAT>::forward<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;
index 14754a8b089ca7093a6f6099355fe8ae38e92ffe..12030d44816a63dab5d5fb5b8435cd0b4f60b91d 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -40,7 +40,6 @@ template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatri
 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>&);
@@ -78,9 +77,6 @@ template void FullMatrix<TYPEMAT>::forward<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> (
@@ -107,9 +103,6 @@ template void FullMatrix<TYPEMAT>::forward<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;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.