jacobian

commit 88f6f81d4b6441ea80fab2f13d35d8c8df5b0f3f
parent 12e0a6fff7be7c5916ae861fe2b303fdb16870df
Author: David Freifeld <freifeld.david@gmail.com>
Date:   Wed,  1 Jul 2020 12:20:07 -0700

Strassen now works but has incorrect quadrants

Diffstat:
M
tests/strassen.cpp
 | 
21
+++++++--------------

1 file changed, 7 insertions(+), 14 deletions(-)
diff --git a/tests/strassen.cpp b/tests/strassen.cpp
@@ -13,12 +13,11 @@
 // Not remotely extensible, just initial test.
 Eigen::MatrixXd strassen_mul(Eigen::MatrixXd a, Eigen::MatrixXd b)
 {
-  int block_len = 1250;
+  int block_len = a.rows()/2;
   Eigen::MatrixXd result (a.rows(), a.cols());
+
   // Nigh unreadable code follows.
-  std::cout << "STRASSEN END INIT\n";
-  
-  Eigen::MatrixXd m1 = (a.block(0,0, block_len, block_len) + a.block(a.rows()-block_len,a.cols()-block_len, block_len, block_len)) * (b.block(0,0, block_len, block_len) + b.block(b.rows()-block_len,b.cols()-block_len, block_len, block_len));
+  Eigen::MatrixXd m1 = (a.block(0,0, block_len, block_len)) + a.block(a.rows()-block_len,a.cols()-block_len, block_len, block_len) * (b.block(0,0, block_len, block_len) + b.block(b.rows()-block_len,b.cols()-block_len, block_len, block_len));
   Eigen::MatrixXd m2 = (a.block(a.rows()-block_len, 0, block_len, block_len) + a.block(a.rows()-block_len,a.cols()-block_len, block_len, block_len)) * (b.block(0,0, block_len, block_len));
   Eigen::MatrixXd m3 = a.block(0,0, block_len, block_len) * (b.block(0,b.cols()-block_len, block_len, block_len) - b.block(b.rows()-block_len,b.cols()-block_len, block_len, block_len));
   Eigen::MatrixXd m4 = a.block(a.rows()-block_len,a.cols()-block_len, block_len, block_len) * (b.block(b.rows()-block_len,0, block_len, block_len) - b.block(0,0, block_len, block_len));
@@ -26,33 +25,27 @@ Eigen::MatrixXd strassen_mul(Eigen::MatrixXd a, Eigen::MatrixXd b)
   Eigen::MatrixXd m6 = (a.block(a.rows()-block_len,0, block_len, block_len) - a.block(0,0, block_len, block_len)) * (b.block(0,0, block_len, block_len) + b.block(0,b.cols()-block_len, block_len, block_len));
   Eigen::MatrixXd m7 = (a.block(0,a.cols()-block_len, block_len, block_len) - a.block(a.rows()-block_len,a.cols()-block_len, block_len, block_len)) * (b.block(a.rows()-block_len,0, block_len, block_len) + b.block(b.rows()-block_len,b.cols()-block_len, block_len, block_len));
 
-  std::cout << "STRASSEN END m\n";
-  
   result.block(0,0, block_len, block_len) =  m1 + m4 - m5 + m7;
   result.block(0,result.cols()-block_len, block_len, block_len) =  m3 + m5;
   result.block(result.rows()-block_len,0, block_len, block_len) =  m2 + m4;
   result.block(result.rows()-block_len,result.cols()-block_len, block_len, block_len) =  m1 - m2 + m3 + m6;
 
-  std::cout << "STRASSEN END CALC\n";
   return result;
 }
 
 int main()
 {
-  Eigen::MatrixXd a = Eigen::MatrixXd::Random(5000, 5000);
-  Eigen::MatrixXd b = Eigen::MatrixXd::Random(5000, 5000);
-
-  std::cout << "END INIT\n";
+  Eigen::MatrixXd a = Eigen::MatrixXd::Random(16, 16);
+  Eigen::MatrixXd b = Eigen::MatrixXd::Random(16, 16);
   
   auto eigen_begin = std::chrono::high_resolution_clock::now();
   Eigen::MatrixXd product = a * b;
   auto eigen_end = std::chrono::high_resolution_clock::now();
   
-  std::cout << "END EIGEN\n";
- 
   auto strassen_begin = std::chrono::high_resolution_clock::now();
   Eigen::MatrixXd sproduct = strassen_mul(a, b);
   auto strassen_end = std::chrono::high_resolution_clock::now();
   
-  std::cout << "EIGEN: " << std::chrono::duration_cast<std::chrono::nanoseconds>(eigen_end - eigen_begin).count() / pow(10,9) << " STRASSEN: " << std::chrono::duration_cast<std::chrono::nanoseconds>(strassen_end - strassen_begin).count() / pow(10,9);
+  std::cout << "EIGEN: " << std::chrono::duration_cast<std::chrono::nanoseconds>(eigen_end - eigen_begin).count() / pow(10,9) << " STRASSEN: " << std::chrono::duration_cast<std::chrono::nanoseconds>(strassen_end - strassen_begin).count() / pow(10,9) << "\n";
+  std::cout << "A:\n" << a << "\nB:\n" << b << "\nEigen:\n" << product << "\nStrassen:\n" << sproduct << "\n";
 }