Thanks for you convenient library. I just found an amazing thing: Eigen performs very well in multiplying complex matrices.

This is my *.cpp

Code: Select all

#include <cstdio>
#include <ctime>
#include <iostream>
#include <Eigen/Core>
using namespace std;
using namespace Eigen;

unsigned int i, res = 1, size, timer;
MatrixXcd X, Y, Z, ZZ;
MatrixXd A, B, C, D, P, Q, R, S;

void program()
{   
   X = MatrixXcd::Random(size, size);
   Y = MatrixXcd::Random(size, size);
   Z.resize(size, size);
   printf("%d\t", size);
   
   timer = clock();
   for(i = 0; i < res; ++i)
   ZZ = X*Y;
   printf("%d\t", clock() - timer);
   
   timer = clock();
   for(i = 0; i < res; ++i)
   {
      A = X.real();
      B = X.imag();
      C = Y.real();
      D = Y.imag();
      P = A*C;
      Q = B*D;
      Z.real() = P - Q;
      Z.imag() = (A + B)*(C + D) - P - Q;
   }
   printf("%d\t", clock() - timer);
   cout << (Z - ZZ).norm() << endl;
   return;
}

int main()
{
   srand(time(NULL));
   for(size = 1; ; size <<= 1)
   program();
   return 0;
}


And it output:

Code: Select all

1       0       0       0
2       0       0       2.02064e-016
4       0       0       7.35654e-016
8       0       0       3.55112e-015
16      0       0       1.17275e-014
32      0       16      4.25013e-014
64      62      63      1.6616e-013
128     454     437     6.55415e-013
256     3578    3484    2.56773e-012
512     30704   32437   9.30642e-012
...


I think it indicates that Eigen's complex multiplication takes ONLY 3 REAL MULTS!!! (instead of 4)

However, when it comes to matrix chain multiplication, it seems that Eigen just multiplies from left to right. Doesn't a MatrixBase stores a "matrix expression?" Is it possible to optimize matrix chain multiplication without any function other than "binary operator* ?"

Note that Eigen vectorizes operations on complexes, including matrix products involving real and complex objects.
Regarding A*B*C we currently follow the standard C++ rules to order the operations, so you can still use parentheses if you know in advance what's the best order. In the near future we will be able to perform such optimization automatically.

Thanks for you reply.

I'm going to settle for my ugly function expression for now.

(Methods provided by Eigen is far more elegant!!!)