C++ 关于QR矩阵分解

void QR(const Matrix& m, Matrix& Q,Matrix& R)
{
int n = 0;
int row = m.rows();
int col = m.cols();
if (row>col){ n =col; }
else{ n =row; }
Matrix mtemp = m;
Matrix q1(row, col);
Q.resize(row,row);
R.resize(row,col);

int i, j, k, r, l;
double temp, sum, dr, cr, hr;
double *ur = new double[n];
double *pr = new double[n];
double *wr = new double[n];

for (i = 0; i<n; i++)//定义单位矩阵
for (j = 0; j<n; j++)
{
    if (i == j)Q[i][j] = 1;
    else Q[i][j] = 0;
}
for (r = 0; r<n; r++)
{
    temp = 0;
    for (k = r; k<n; k++)
        temp += fabs(mtemp[k][r]);
    if (temp >= 0.0)
    {
        sum = 0;
        for (k = r; k<n; k++)
            sum += mtemp[k][r] * mtemp[k][r];
        dr = sqrt(sum);
        if (mtemp[r][r]>0.0) l = -1;
        else l = 1;
        cr = l*dr;
        hr = cr*(cr - mtemp[r][r]);

        for (i = 0; i<n; i++)//定义ur
        {
            if (i<r)ur[i] = 0;
            if (i == r)ur[i] = mtemp[r][r] - cr;
            if (i>r)ur[i] = mtemp[i][r];
        }
        for (i = 0; i<n; i++)//定义wr
        {
            sum = 0;
            for (j = 0; j<n; j++)
                sum += Q[i][j] * ur[j];
            wr[i] = sum;
        }
        for (i = 0; i<row; i++)//定义qr
        for (j = 0; j<col; j++)
        {
            q1[i][j] = Q[i][j] - wr[i] * ur[j] / hr;
        }
        for (i = 0; i<row; i++)//定义qr+1
        for (j = 0; j<col; j++)
        {
            Q[i][j] = q1[i][j];
        };
        for (i = 0; i<col; i++)//定义pr
        {
            sum = 0;
            for (j = 0; j<n; j++)
                sum += mtemp[j][i] * ur[j];
            pr[i] = sum / hr;
        };
        for (i = 0; i<row; i++)
        for (j = 0; j<col; j++)
        {
            mtemp[i][j] = mtemp[i][j] - ur[i] * pr[j];
        };
    };
};
for (i = 0; i<row; i++)
for (j = 0; j<col; j++)
{
    if (fabs(mtemp[i][j])<0.0) mtemp[i][j] = 0;
};
for (i = 0; i<row; i++)
for (j = 0; j<col; j++)
{
    R[i][j] = mtemp[i][j];
}

}

什么问题呢?难道是分享资料不成!

哈,楼主还个代码写的还可以,就是不知道什么问题

有问题吗?还是只是来分享代码呢?