]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Improve our own mmult implementation.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 30 May 2009 12:29:22 +0000 (12:29 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 30 May 2009 12:29:22 +0000 (12:29 +0000)
git-svn-id: https://svn.dealii.org/trunk@18889 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 4b7744d5e09c87ab581a0eb0c491365ed0fb5e28..f56acaf1c44a502c87cf4fa5fbb0ed45134b0b8b 100644 (file)
@@ -626,24 +626,26 @@ void FullMatrix<number>::mmult (FullMatrix<number2>       &dst,
 
 #endif
 
+  const unsigned int m = this->m(), n = src.n(), l = this->n();
+
+                                       // arrange the loops in a way 
+                                       // that we can access contiguous 
+                                       // fields at the innermost loop 
+                                       // (even though this introduces 
+                                       // a lot of writing into the 
+                                       // destination matrix)
   if (!adding)
-    for (unsigned int i=0; i<m(); i++)
-      for (unsigned int j=0; j<src.n(); j++)
-       {
-         number s = 0.;
-         for (unsigned k=0; k<n(); k++)
-           s += (*this)(i,k) * number(src(k,j));
-         dst(i,j) = s;
-       }
-  else
-    for (unsigned int i=0; i<m(); i++)
-      for (unsigned int j=0; j<src.n(); j++)
-       {
-         number s = 0.;
-         for (unsigned k=0; k<n(); k++)
-           s += (*this)(i,k) * number(src(k,j));
-         dst(i,j) += s;
+    dst = (number2)0;
+  for (unsigned int i=0; i<m; i++)
+    {
+      number2 * dst_ptr = &dst(i,0);
+      for (unsigned int k=0; k<l; k++)
+        {
+         const number2 mult = (number2)(*this)(i,k);
+         for (unsigned int j=0; j<n; j++)
+           dst_ptr[j] += mult * (number2)(src(k,j));
        }
+    }
 }
 
 
@@ -716,25 +718,26 @@ void FullMatrix<number>::Tmmult (FullMatrix<number2>       &dst,
 
 #endif
 
+  const unsigned int m = n(), n = src.n(), l = this->m();
 
+                                       // arrange the loops in a way 
+                                       // that we can access contiguous 
+                                       // fields at the innermost loop 
+                                       // (even though this introduces 
+                                       // a lot of writing into the 
+                                       // destination matrix)
   if (!adding)
-    for (unsigned int i=0; i<n(); i++)
-      for (unsigned int j=0; j<src.n(); j++)
-       {
-         number s = 0;
-         for (unsigned int k=0; k<m(); k++)
-           s += (*this)(k,i) * number(src(k,j));
-         dst(i,j) = s;
-       }
-  else
-    for (unsigned int i=0; i<n(); i++)
-      for (unsigned int j=0; j<src.n(); j++)
-       {
-         number s = 0;
-         for (unsigned int k=0; k<m(); k++)
-           s += (*this)(k,i) * number(src(k,j));
-         dst(i,j) += s;
+    dst = (number2)0;
+  for (unsigned int i=0; i<m; i++)
+    {
+      number2 * dst_ptr = &dst(i,0);
+      for (unsigned int k=0; k<l; k++)
+        {
+         const number2 mult = (number2)(*this)(k,i);
+         for (unsigned int j=0; j<n; j++)
+           dst_ptr[j] += mult * (number2)(src(k,j));
        }
+    }
 }
 
 

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.