jacobian

commit 67393ec4b94378dac2c8ae450c5309227f285dae
parent 58dab65b73329cbb51485c00b420846b8445785f
Author: David Freifeld <freifeld.david@gmail.com>
Date:   Sat,  4 Jul 2020 10:37:51 -0700

Attempt to get strassen pad working

Diffstat:
M
tests/strassen.cpp
 | 
24
++++++++++++++++++++----
M
utils.cpp
 | 
6
+++++-

2 files changed, 25 insertions(+), 5 deletions(-)
diff --git a/tests/strassen.cpp b/tests/strassen.cpp
@@ -11,9 +11,24 @@
 #include <iostream>
 
 // Not remotely extensible, just initial test.
-Eigen::MatrixXd strassen_mul(Eigen::MatrixXd a, Eigen::MatrixXd b)
+Eigen::MatrixXd strassen_mul(Eigen::MatrixXd x, Eigen::MatrixXd y)
 {
-  int block_len = a.rows()/2;
+  int power = 0;
+  int largest;
+  if (x.rows() >= x.cols() || x.rows() >= y.cols() || x.rows() >= y.rows()) largest = x.rows();
+  else if (x.cols() >= x.rows() || x.cols() >= y.cols() || x.cols() >= y.rows()) largest = x.cols();
+  else if (y.cols() >= x.rows() || y.cols() >= x.cols() || y.cols() >= y.rows()) largest = y.cols();
+  else if (y.rows() >= x.rows() || y.rows() >= x.cols() || y.rows() >= y.cols()) largest = y.rows();
+  for (; pow(2, power) < largest; power++) {
+    std::cout << pow(2, power) << " " << power <<"\n";
+  }
+  std::cout << "POWER: " << power << " LARGEST: " << largest << "\n";
+  Eigen::MatrixXd a ((int)pow(2,power), (int)pow(2,power));
+  Eigen::MatrixXd b ((int)pow(2,power), (int)pow(2,power));
+  a.block(0,0,x.rows(), x.cols()) = x;
+  b.block(0,0,y.rows(), y.cols()) = y;
+  //Eigen::MatrixXd a ()
+  int block_len = largest/2;
   Eigen::MatrixXd result (a.rows(), a.cols());
 
   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));
@@ -23,6 +38,7 @@ Eigen::MatrixXd strassen_mul(Eigen::MatrixXd a, Eigen::MatrixXd b)
   Eigen::MatrixXd m5 = (a.block(0, 0, block_len, block_len) + a.block(0,a.cols()-block_len, block_len, block_len)) * (b.block(b.rows()-block_len,b.cols()-block_len, block_len, block_len));
   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));
+
   
   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;
@@ -37,8 +53,8 @@ int main()
   srand((unsigned int) time(0));
   int sz;
   std::cin >> sz;
-  Eigen::MatrixXd a = Eigen::MatrixXd::Random(sz, sz);
-  Eigen::MatrixXd b = Eigen::MatrixXd::Random(sz, sz);
+  Eigen::MatrixXd a = Eigen::MatrixXd::Random(sz, 3);
+  Eigen::MatrixXd b = Eigen::MatrixXd::Random(3, sz);
   
   auto eigen_begin = std::chrono::high_resolution_clock::now();
   Eigen::MatrixXd product = a * b;
diff --git a/utils.cpp b/utils.cpp
@@ -57,8 +57,12 @@ std::function<double(double)> rectifier(double (*activation)(double))
   return rectified;
 }
 
-Eigen::MatrixXd strassen_mul(Eigen::MatrixXd a, Eigen::MatrixXd b)
+Eigen::MatrixXd strassen_mul(Eigen::MatrixXd x, Eigen::MatrixXd y)
 {
+  int apower;
+  int bpower;
+  for (;a
+  //Eigen::MatrixXd a ()
   int block_len = a.rows()/2;
   Eigen::MatrixXd result (a.rows(), a.cols());