相关系数矩阵

      最近一直在做主成分分析的算法实现，其中有一环就是由样本数据来求相关系数矩阵（然后对此实对称矩阵用雅可比方法求出特征值和特征向量），当然我知道我设计的算法肯定不是最有效率的，但我还是想和大家分享一下。
我所处理的数据类似于下面这种格式：

//head file:CMatrix.h
#define COL 1000             //指定最大行数
#define ROW 2000           //指定最大列数
class CMatrix
{
public:
    double GetElement(int row, int col);    //获取data[row][col]的值
    void ChangeElement(int row, int col, double value);    //改变data[row][col]的值为value
    double Cal_Standard_Deviation(int r, int c);    //返回c列变量的标准差
    double Cal_Average(int r, int c);    //返回c列变量的均值
    double* Cal_Mutuality_Matrix(int r, int c);    //返回相关系数矩阵的所有元素值
    double Cal_Exy(int r, int c1, int c2);    //计算E(Xi*Xj)，c1、c2代表两列
private:
    double data[ROW][COL];    //存储矩阵内容
};

//source file:CMatrix.cpp
#include "CMatrix.h"
#include
double CMatrix::GetElement(int row, int col)
{
    return this->data[row][col];
}
void CMatrix::ChangeElement(int row, int col, double value)
{
    this->data[row][col] = value;
}
double CMatrix::Cal_Standard_Deviation(int r, int c)
{
    double result = 0.0;
    double Ex = 0.0;
    double Exx = 0.0;
    Ex = this->Cal_Average(r, c);
    for (int i = 0; i < r; i++)
    {
        Exx += data[i][c - 1] * data[i][c - 1];
    }
    Exx /= r;
    result = sqrt(Exx - Ex * Ex);
    return result;
}
double CMatrix::Cal_Average(int r,int c)
{
    double ave = 0.0; //存储均值
    for(int i = 0; i < r; i ++ )
    {
        ave += data[i][c-1];
    }
    ave /= r;
    return ave;
}
double* CMatrix::Cal_Mutuality_Matrix(int r, int c)
{
    double *Cov = new double[c*c];
    for (int i = 0; i < c; i ++)
    {
        for (int j = 0; j < c; j ++)
        {
            Cov[i * c + j] = ( this->Cal_Exy(r, i + 1, j + 1) - this->Cal_Average(r, i + 1) * this->Cal_Average(r , j + 1) )
                             / ( this->Cal_Standard_Deviation(r, i + 1) * this->Cal_Standard_Deviation(r, j + 1) ) ;
        }
    }

    return Cov;
}
double CMatrix::Cal_Exy(int r, int c1, int c2)
{
    double result = 0.0;    //存储运算结果
    double *xy = new double[r];
    for(int i = 0; i < r; i ++)
    {
        xy[i] = data[i][c1 - 1] * data[i][c2 - 1];
    }
    for (i = 0; i < r; i++)
    {
        result += xy[i];
    }
    result /= r;
    return result;
}