c++实现计算特征值和特征向量

今天实现了一下特征值和特征向量，参考https://blog.youkuaiyun.com/webzhuce/article/details/85013301
这里的特征值可以理解为对应特征向量的贡献量。每一个特征向量理解为一个线性的子空间。

bool My_Jacobi(const vector<vector<double>>& matrix, vector<vector<double>>& eigenvectors, vector<double>& eigenvalues, double precision, int max)
{
 int dim = matrix.size();
 vector<vector<double>> mat;
 vector<double> t;
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 mat.push_back(t);
 t.clear();
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 mat.push_back(t);
 t.clear();
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 mat.push_back(t);
 t.clear();
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 mat.push_back(t);
 t.clear();
    // 初始化特征矩阵
    for (int i = 0; i < dim; i++)
    {
        for (int j = 0; j < dim; j++)
        {
   mat[i][j] = matrix[i][j];
            if(i==j)
            {
                eigenvectors[i][j] = 1.0;
            }
            else
            {
                eigenvectors[i][j] = 0.0;
            }
        }
    }
    int nCount = 0;  //current iteration
    while(1)
    {
        // 搜索矩阵绝对值最大元素及下标
  double dbMax = mat[0][1];
        int nRow = 0;
        int nCol = 1;
        for (int i = 0; i < dim; i++) //row
        {
            for (int j = 0; j < dim; j++)
            {
    double d = fabs(mat[i][j]);
                if((i != j) && (d > dbMax))
                {
     dbMax = d;
     nRow = i;
     nCol = j;
                }
            }
        }
  cout << dbMax << "," << nRow << "," << nCol << endl;
  // 阈值条件判断
  if (dbMax < precision)     //precision check
   break;
  if (nCount > max)       //iterations check
   break;
  nCount++;
  double dbApp = mat[nRow][nRow];
  double dbApq = mat[nRow][nCol];
  double dbAqq = mat[nCol][nCol];
  // 计算旋转矩阵
  double dbAngle = 0.5*atan2(-2 * dbApq, dbAqq - dbApp);
  double dbSinTheta = sin(dbAngle);
  double dbCosTheta = cos(dbAngle);
  double dbSin2Theta = sin(2 * dbAngle);
  double dbCos2Theta = cos(2 * dbAngle);
  mat[nRow][nRow] = dbApp*dbCosTheta*dbCosTheta + dbAqq*dbSinTheta*dbSinTheta + 2 * dbApq*dbCosTheta*dbSinTheta;
  mat[nCol][nCol] = dbApp*dbSinTheta*dbSinTheta + dbAqq*dbCosTheta*dbCosTheta - 2 * dbApq*dbCosTheta*dbSinTheta;
  mat[nRow][nCol] = 0.5*(dbAqq - dbApp)*dbSin2Theta + dbApq*dbCos2Theta;
  mat[nCol][nRow] = mat[nRow][nCol];
  for (int i = 0; i < dim; i++) 
  {
   if ((i != nCol) && (i != nRow)) 
   {
    dbMax = mat[i][nRow];
    mat[i][nRow] = mat[i][nCol] * dbSinTheta + dbMax*dbCosTheta;
    mat[i][nCol] = mat[i][nCol] * dbCosTheta - dbMax*dbSinTheta;
   }
  }
  for (int j = 0; j < dim; j++) {
   if ((j != nCol) && (j != nRow)) {
    dbMax = mat[nRow][j];
    mat[nRow][j] = mat[nCol][j] * dbSinTheta + dbMax*dbCosTheta;
    mat[nCol][j] = mat[nCol][j] * dbCosTheta - dbMax*dbSinTheta;
   }
  }
  //compute eigenvector
  for (int i = 0; i < dim; i++) 
  {
   dbMax = eigenvectors[i][nRow];
   eigenvectors[i][nRow] = eigenvectors[i][nCol] * dbSinTheta + dbMax*dbCosTheta;
   eigenvectors[i][nCol] = eigenvectors[i][nCol] * dbCosTheta - dbMax*dbSinTheta;
  }
    }
 // 特征值排序
 std::map<double, int> mapEigen;
 for (int i = 0; i < dim; i++) 
 {
  eigenvalues[i] = mat[i][i];
  mapEigen.insert(make_pair(eigenvalues[i], i));
 }
 double *pdbTmpVec = new double[dim*dim];
 std::map<double, int>::reverse_iterator iter = mapEigen.rbegin();
 for (int j = 0; iter != mapEigen.rend(), j < dim; ++iter, ++j) {
  for (int i = 0; i < dim; i++) {
   pdbTmpVec[i*dim + j] = eigenvectors[i][iter->second];
  }
  eigenvalues[j] = iter->first;
 }
 for (int i = 0; i < dim; i++) 
 {
  double dSumVec = 0;
  for (int j = 0; j < dim; j++)
   dSumVec += pdbTmpVec[j * dim + i];
  if (dSumVec < 0) 
  {
   for (int j = 0; j < dim; j++)
    pdbTmpVec[j * dim + i] *= -1;
  }
 }
 for (int i = 0; i < dim; i++)
 {
  for (int j = 0; j < dim; j++)
  {
   eigenvectors[i][j] = pdbTmpVec[i * dim + j];
  }
 }
 delete[]pdbTmpVec;
 mat.clear();
 t.clear();
    return true;
}

int main()
{
 vector<vector<double>> temp;
 vector<double> t;
 t.push_back(4.0);
 t.push_back(-30.0);
 t.push_back(60.0);
 t.push_back(-35.0);
 temp.push_back(t);
 t.clear();
 t.push_back(-30.0);
 t.push_back(300.0);
 t.push_back(-675.0);
 t.push_back(420.0);
 temp.push_back(t);
 t.clear();
 t.push_back(60.0);
 t.push_back(-675.0);
 t.push_back(1620.0);
 t.push_back(-1050.0);
 temp.push_back(t);
 t.clear();
 t.push_back(-35.0);
 t.push_back(420.0);
 t.push_back(-1050.0);
 t.push_back(700.0);
 temp.push_back(t);
 t.clear();
 vector<std::vector<double>> eigenvectors;
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 eigenvectors.push_back(t);
 t.clear();
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 eigenvectors.push_back(t);
 t.clear();
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 eigenvectors.push_back(t);
 t.clear();
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 t.push_back(0.0);
 eigenvectors.push_back(t);
 t.clear();
 vector<double> eigenvalues;
 eigenvalues.resize(4, 0.0);
 int n = 4;
 double eps = 1e-10; 
 int T = 10000;
 bool re = My_Jacobi(temp, eigenvectors, eigenvalues, eps, T);
    if (re) 
 {
        cout << "Matrix:" << endl;
        for (int i = 0; i < n; i++) 
  {
   for (int j = 0; j < n; j++)
   {
    cout << setw(15) << temp[i][j];
   }
            cout << endl;
        }
        cout << "eigenvectors:" << endl;
        for (int i = 0; i < n; i++) 
  {
   for (int j = 0; j < n; j++)
   {
    cout << setw(15) << eigenvectors[i][j];
   }
            cout << endl;
        }
        cout << "eigenvalues:" << endl;
        for (int i = 0; i < n; i++) 
  {
   cout << setw(15) << eigenvalues[i];
            cout << endl;
        }
    }
 else
 {
  cout << "false" << endl;
 }
    system("pause");
}