Eig

bool matrix::Eig(

  matrix&  eigen_vectors,     // 本征向量矩阵

  vector&  eigen_values       // 本征向量向量

   );

#property script_show_inputs

//--- 输入参数

input int  InpSize1   =512;

input int  InpSize2   =256;

input int  InpSize3   =1024;

//+------------------------------------------------------------------+

//| 填充测试平方可逆矩阵                                                |

//+------------------------------------------------------------------+

template<typename T>

void MatrixFill(matrix<T> &matrix_a)

 {

  ulong size_m=matrix_a.Rows();

  ulong size_k=matrix_a.Cols();

  T     value=0.0;

//--- 填充矩阵

  for(ulong i=0; i<size_m; i++)

   {

    for(ulong j=0; j<size_k; j++)

     {

      if(i==j)

        matrix_a[i][j]=T(1.0+i);

      else

       {

        value+=1.0;

        matrix_a[i][j]=value;

       }

     }

   }

 }

//+------------------------------------------------------------------+

//| 脚本程序 start 函数                                                |

//+------------------------------------------------------------------+

int OnStart()

 {

  int errors=0;



  errors+=TestEigen<double>(InpSize1);

  errors+=TestEigen<double>(InpSize2);

  errors+=TestEigen<double>(InpSize3);



  errors+=TestEigen<float>(InpSize1);

  errors+=TestEigen<float>(InpSize2);

  errors+=TestEigen<float>(InpSize3);

//---

  Print("Test ", errors?"failed":"passed");

  return(errors);

 }



/*

结果



Eigen solver of double matrix 512 x 512  passed  vectors=489  time=4268.506 ms

Eigen solver of double matrix 256 x 256  passed  vectors=251  time=417.610 ms

Eigen solver of double matrix 1024 x 1024  passed  vectors=916  time=43708.280 ms

Eigen solver of float matrix 512 x 512  passed  vectors=1  time=2508.357 ms

Eigen solver of float matrix 256 x 256  passed  vectors=1  time=188.859 ms

Eigen solver of float matrix 1024 x 1024  passed  vectors=1  time=27209.666 ms

测试 通过

*/



//+------------------------------------------------------------------+

//|  测试 Eig 方法                                                    |

//+------------------------------------------------------------------+

template<typename T>

int TestEigen(const int size_m)

 {

  int vectors=0;

  matrix<T> matrix_a(size_m, size_m);

  matrix<T> matrix_v(size_m, size_m);

  vector<T> vector_e(size_m);

//--- 填充平方矩阵

  MatrixFill(matrix_a);

//--- 将测量单位微秒

  ulong t1=GetMicrosecondCount();

//--- 本征解算器

  matrix_a.Eig(matrix_v, vector_e);

//--- 测量

  ulong t2=GetMicrosecondCount();

//--- 验算是否 A * v = lambda * v

  for(ulong n=0; n<vector_e.Size(); n++)

   {

    vector<T> eigen_vector=matrix_v.Col(n);

    vector<T> vector_c1   =eigen_vector*vector_e[n];

    vector<T> vector_c2   =matrix_a.MatMul(eigen_vector);



    //--- 偏差太多，将精度检查降低到第 10 位

    ulong errors=vector_c1.CompareByDigits(vector_c2, sizeof(T)==sizeof(double) ? 10 : 5);

    if(int(errors)<size_m/10)

      vectors++;

   }

  double elapsed_time=double(t2-t1)/1000.0;

  printf("Eigen solver of %s matrix %d x %d  %s  vectors=%d  time=%.3f ms",

         typename(T), size_m, size_m, vectors>0?"passed":"failed", vectors, elapsed_time);

  return(vectors==0);

 }