]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make some small optimization in vmult of FullMatrix.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 28 May 2009 09:55:05 +0000 (09:55 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 28 May 2009 09:55:05 +0000 (09:55 +0000)
git-svn-id: https://svn.dealii.org/trunk@18884 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/full_matrix.templates.h
deal.II/lac/include/lac/householder.h
deal.II/lac/include/lac/sparse_matrix.templates.h

index 936884e238f5094422423b2c9fd8d70f00aadd3b..4b7744d5e09c87ab581a0eb0c491365ed0fb5e28 100644 (file)
@@ -198,27 +198,17 @@ FullMatrix<number>::vmult (Vector<number2>& dst,
   Assert (&src != &dst, ExcSourceEqualsDestination());
 
   const number* e = this->data();
-  const unsigned int size_m = m(),
-                    size_n = n();
-  if (!adding)
-    {
-      for (unsigned int i=0; i<size_m; ++i)
-       {
-         number2 s = 0.;
-         for (unsigned int j=0; j<size_n; ++j)
-           s += number2(src(j)) * number2(*(e++));
-         dst(i) = s;
-       };
-    }
-  else
+                                  // get access to the data in order to
+                                  // avoid copying it when using the ()
+                                  // operator
+  const number2* src_ptr = &(*const_cast<Vector<number2>*>(&src))(0);
+  const unsigned int size_m = m(), size_n = n();
+  for (unsigned int i=0; i<size_m; ++i)
     {
-      for (unsigned int i=0; i<size_m; ++i)
-       {
-         number2 s = 0.;
-         for (unsigned int j=0; j<size_n; ++j)
-           s += number2(src(j)) * number2(*(e++));
-         dst(i) += s;
-       }
+      number2 s = adding ? dst(i) : 0.;
+      for (unsigned int j=0; j<size_n; ++j)
+       s += src_ptr[j] * number2(*(e++));
+      dst(i) = s;
     }
 }
 
@@ -231,35 +221,29 @@ void FullMatrix<number>::Tvmult (Vector<number2>       &dst,
                                 const bool             adding) const
 {
   Assert (!this->empty(), ExcEmptyMatrix());
-  
+
   Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n()));
   Assert(src.size() == m(), ExcDimensionMismatch(src.size(), m()));
 
   Assert (&src != &dst, ExcSourceEqualsDestination());
 
-  const unsigned int size_m = m(),
-                    size_n = n();
+  const number* e = this->data();
+  number2* dst_ptr = &dst(0);
+  const unsigned int size_m = m(), size_n = n();
 
+                                  // zero out data if we are not adding
   if (!adding)
+    for (unsigned int j=0; j<size_n; ++j)
+      dst_ptr[j] = 0.;
+
+                                  // write the loop in a way that we can
+                                  // access the data contiguously
+  for (unsigned int i=0; i<size_m; ++i)
     {
-      for (unsigned int i=0; i<size_n; ++i)
-       {
-         number s = 0.;
-         for (unsigned int j=0; j<size_m; ++j)
-           s += number(src(j)) * (*this)(j,i);
-         dst(i) = s;
-       };
-    }
-  else
-    {
-      for (unsigned int i=0; i<size_n; ++i)
-       {
-         number s = 0.;
-         for (unsigned int j=0; j<size_m; ++j)
-           s += number(src(j)) * (*this)(j,i);
-         dst(i) += s;
-       }
-    }
+      const number2 d = src(i);
+      for (unsigned int j=0; j<size_n; ++j)
+       dst_ptr[j] += d * number2(*(e++));
+    };
 }
 
 
index 3162e71f39e12094708f9259a433a14545f08eee..28833bf57966c39c2a750e62e3cfd8f82ec8c355 100644 (file)
@@ -114,25 +114,28 @@ Householder<number>::Householder()
 {}
 
 
+
 template <typename number>
 template <typename number2>
 void
 Householder<number>::initialize(const FullMatrix<number2>& M)
 {
-  this->reinit(M.n_rows(), M.n_cols());
+  const unsigned int m = M.n_rows(), n = M.n_cols();
+  this->reinit(m, n);
   this->fill(M);
-  diagonal.resize(M.n_rows());
   Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
+  diagonal.resize(m);
+
                                   // m > n, src.n() = m
   Assert (this->n_cols() <= this->n_rows(),
          ExcDimensionMismatch(this->n_cols(), this->n_rows()));
 
-  for (unsigned int j=0 ; j<this->n() ; ++j)
+  for (unsigned int j=0 ; j<n ; ++j)
   {
     number2 sigma = 0;
     unsigned int i;
                                     // sigma = ||v||^2
-    for (i=j ; i<this->m() ; ++i)
+    for (i=j ; i<m ; ++i)
       sigma += this->el(i,j)*this->el(i,j);
                                     // We are ready if the column is
                                     // empty. Are we?
@@ -148,20 +151,20 @@ Householder<number>::initialize(const FullMatrix<number2>& M)
     diagonal[j] = beta*(this->el(j,j) - s);
     this->el(j,j) = s;
 
-    for (i=j+1 ; i<this->m() ; ++i)
+    for (i=j+1 ; i<m ; ++i)
       this->el(i,j) *= beta;
 
 
                                     // For all subsequent columns do
                                     // the Householder reflexion
-    for (unsigned int k=j+1 ; k<this->n() ; ++k)
+    for (unsigned int k=j+1 ; k<n ; ++k)
     {
       number2 sum = diagonal[j]*this->el(j,k);
-      for (i=j+1 ; i<this->m() ; ++i)
+      for (i=j+1 ; i<m ; ++i)
        sum += this->el(i,j)*this->el(i,k);
 
       this->el(j,k) -= sum*this->diagonal[j];
-      for (i=j+1 ; i<this->m() ; ++i)
+      for (i=j+1 ; i<m ; ++i)
        this->el(i,k) -= sum*this->el(i,j);
     }    
   }
@@ -185,6 +188,8 @@ Householder<number>::least_squares (Vector<number2>& dst,
   Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
   AssertDimension(dst.size(), this->n());
   AssertDimension(src.size(), this->m());
+
+  const unsigned int m = this->m(), n = this->n();
   
   GrowingVectorMemory<Vector<number2> > mem;
   Vector<number2>* aux = mem.alloc();
@@ -194,20 +199,20 @@ Householder<number>::least_squares (Vector<number2>& dst,
 
                                   // 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)
+  for (unsigned int j=0;j<n;++j)
     {
                                       // sum = v_i^T dst
       number2 sum = diagonal[j]* (*aux)(j);
-      for (unsigned int i=j+1 ; i<this->m() ; ++i)
+      for (unsigned int i=j+1 ; i<m ; ++i)
        sum += this->el(i,j)*(*aux)(i);
                                       // dst -= v * sum
       (*aux)(j) -= sum*diagonal[j];
-      for (unsigned int i=j+1 ; i<this->m() ; ++i)
+      for (unsigned int i=j+1 ; i<m ; ++i)
        (*aux)(i) -= sum*this->el(i,j);
     }
                                   // Compute norm of residual
   number2 sum = 0.;
-  for (unsigned int i=this->n() ; i<this->m() ; ++i)
+  for (unsigned int i=n ; i<m ; ++i)
     sum += (*aux)(i) * (*aux)(i);
                                   // Compute solution
   this->backward(dst, *aux);
index df34ec954645afa55e0fb022983f240f694b5adb..5b7ee61e36ff451a9025b09f12b55d25251121f9 100644 (file)
@@ -374,11 +374,11 @@ namespace internal
       else
        for (unsigned int row=begin_row; row<end_row; ++row)
          {
-           typename OutVector::value_type s = 0.;
+           typename OutVector::value_type s = *dst_ptr;
            const number *const val_end_of_row = &values[rowstart[row+1]];
            while (val_ptr != val_end_of_row)
              s += *val_ptr++ * src(*colnum_ptr++);
-           *dst_ptr++ += s;
+           *dst_ptr++ = s;
          }
     }
   }

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.