求解广义逆矩阵
double MatrixL[3*m*n*l];
int line=l,queue=m*n*3;
double *const FKMatrixL = new double[(line+queue)*(line+queue)]; //矩阵
double *const _MatrixLdp = new double[line*line]; //矩阵
double *const _MatrixLc = new double[line*line]; //矩阵
double *const _MatrixLct = new double[line*line]; //矩阵
double *const _MatrixLctc = new double[line*line]; //矩阵
double *const MatrixLctc = new double[line*line]; //矩阵
double *const MatrixLctcct = new double[line*line]; //矩阵
double *const _MatrixLdq = new double[queue*queue]; //矩阵
double *const _MatrixLd = new double[queue*queue]; //矩阵
double *const _MatrixLdt = new double[queue*queue]; //矩阵
double *const _MatrixLddt = new double[queue*queue]; //矩阵
double *const MatrixLddt = new double[queue*queue]; //矩阵
double *const MatrixLdtddt = new double[queue*queue]; //矩阵
double *const Ajia = new double[line*queue]; //矩阵
double *const CMatrixL = new double[line*line]; //C矩阵
double *const DMatrixL = new double[queue*queue]; //D矩阵
void DisplayMatrix(double *const lMatrix, const int &line,const int &queue)
{
for (int i=0; i<line; i++)
{
for (int j=0; j<queue; j++)
{
cout << lMatrix[i*queue+j]<<" " ;
}
cout << endl;
}
cout << endl;
}
//打印Czong矩阵
void DisplayCzongMatrix(double *const pl,const int &queue,const int line)
{
for (int i=0; i<queue; i++)
{
for (int j=0; j<line; j++)
{
cout << pl[i*line+j]<<" " ;
}
cout << endl;
}
cout << endl;
}
//打印初等变换后D、Dt、DDt矩阵
//求分块矩阵fkMatrix，并对其进行初等变换，存在_pMatrix中
void FenkuaiCDMatrix(double *const pMatrix, double *const _pMatrixcp,double *const _pMatrixc,double *const _pMatrixct,double *const _pMatrixctc,double *const _pMatrixdq,double *const _pMatrixd,double *const _pMatrixdt,double *const _pMatrixddt, const int &line,const int &queue)
{
double *FKMatrix = new double[(line+queue)*(line+queue)];
for (int i=0; i<line; i++)
{
for (int j=0; j<queue; j++)
FKMatrix[i*(line+queue)+j] = pMatrix[i*queue+j];
}
for (int i=0; i<line; i++)
{
for (int j=queue; j<(queue+line); j++)
FKMatrix[i*(line+queue)+j] = 0.0;
FKMatrix[queue+i*(queue+line+1)] =1.0;
}
for (int i=line; i<line+queue; i++)
{
for (int j=0; j<queue+line; j++)
FKMatrix[i*(line+queue)+j] = 0.0;
FKMatrix[i*(queue+line+1)-line] = 1.0;
}
//fkMatrix Initialization over!
for (int i=0; i<queue; i++)//i小于line和queue中小的那个数
{
double base = FKMatrix[i*(line+queue)+i];
if (fabs(base) < 1E-300)
{
cout << endl << "求行变换过程中被零除，无法求解!" << endl;
exit(0);
}
for (int j=0; j<queue; j++)//j小于line和queue中小的那个数
{
if (j == i) continue;
double times = FKMatrix[i*(line+queue)+j]/base;
for (int k=0; k<(line+queue); k++)//col
{
FKMatrix[k*(line+queue)+j] = FKMatrix[k*(line+queue)+j] - times*FKMatrix[k*(line+queue)+i];
}
}
for (int k=0; k<(line+queue); k++)
{
FKMatrix[k*(line+queue)+i] /= base;
}
}
for (int i=queue; i<line; i++)//i等于line和queue中小的那个数
{
for (int j=0; j<queue; j++)//j小于line和queue中小的那个数
{
FKMatrix[i*(line+queue)+queue+j]=-FKMatrix[i*(line+queue)+j];
FKMatrix[i*(line+queue)+j]=0;
}
}
//求P矩阵
for (int i=0; i<line; i++)
{
for (int j=0; j<line; j++)
_pMatrixcp[i*line+j] = FKMatrix[i*(line+queue)+queue+j];
}
//求C矩阵
double *tMatrixc = new double[2*line*line];
for ( int i=0; i<line; i++){
for (int j=0; j<line; j++)
tMatrixc[i*line*2+j] = _pMatrixcp[i*line+j];
}
for (int i=0; i<line; i++)
{
for (int j=line; j<line*2; j++)
tMatrixc[i*line*2+j] = 0.0;
tMatrixc[i*line*2+line+i] = 1.0;
}
//Initialization over!
for (int i=0; i<line; i++)//Process Cols
{
double base = tMatrixc[i*line*2+i];
if (fabs(base) < 1E-300)
{
cout << endl << "求逆矩阵过程中被零除，无法求解!" << endl;
exit(0);
}
for (int j=0; j<line; j++)//row
{
if (j == i) continue;
double times = tMatrixc[j*line*2+i]/base;
for (int k=0; k<line*2; k++)//col
{
tMatrixc[j*line*2+k] = tMatrixc[j*line*2+k] - times*tMatrixc[i*line*2+k];
}
}
for (int k=0; k<line*2; k++){
tMatrixc[i*line*2+k] /= base;
}
}
for (int i=0; i<line; i++)
{
for (int j=0; j<queue; j++)
_pMatrixc[i*queue+j] = tMatrixc[i*(line+line)+line+j];
}
//求Ct矩阵
for (int i=0; i<queue; i++)
{
for (int j=0; j<line; j++)
_pMatrixct[i*line+j] = _pMatrixc[j*queue+i] ;
}
//求CtC矩阵
for (int i=0; i<queue; i++)
{ double sum=0;
for (int j=0; j<queue; j++)
{
for (int k=0; k<line; k++)
{
sum=sum+(_pMatrixct[i*line+k]*_pMatrixct[j*line+k]);
}
_pMatrixctc[i*queue+j]=sum;sum=0 ;
}
}
//求Q矩阵
for (int i=0; i<queue; i++)
{
for (int j=0; j<queue; j++)
_pMatrixdq[i*queue+j] = FKMatrix[(i+line)*(line+queue)+j];
}
//求D矩阵
double *tMatrixd = new double[2*queue*queue];
for ( int i=0; i<queue; i++)
{
for (int j=0; j<queue; j++)
tMatrixd[i*queue*2+j] = _pMatrixdq[i*queue+j];
}
for (int i=0; i<queue; i++){
for (int j=queue; j<queue*2; j++)
tMatrixd[i*queue*2+j] = 0.0;
tMatrixd[i*queue*2+queue+i] = 1.0;
}
//Initialization over!
for (int i=0; i<queue; i++)//Process Cols
{
double base = tMatrixd[i*queue*2+i];
if (fabs(base) < 1E-300)
{
cout << endl << "求逆矩阵过程中被零除，无法求解!" << endl;
exit(0);
}
for (int j=0; j<queue; j++)//row
{
if (j == i) continue;
double times = tMatrixd[j*queue*2+i]/base;
for (int k=0; k<queue*2; k++)//col
{
tMatrixd[j*queue*2+k] = tMatrixd[j*queue*2+k] - times*tMatrixd[i*queue*2+k];
}
}
for (int k=0; k<queue*2; k++){
tMatrixd[i*queue*2+k] /= base;
}
}
for (int i=0; i<queue; i++)
{
for (int j=0; j<queue; j++)
_pMatrixd[i*queue+j] = tMatrixd[i*queue*2+j+queue];
}
//求Dt矩阵
for (int i=0; i<queue; i++)
{
for (int j=0; j<queue; j++)
_pMatrixdt[i*queue+j] = _pMatrixd[j*queue+i] ;
}
//求DDt积矩阵
for (int i=0; i<queue; i++)
{ double sum=0;
for (int j=0; j<queue; j++)
{
for (int k=0; k<queue; k++)
{
sum=sum+(_pMatrixd[i*queue+k]*_pMatrixd[j*queue+k]);
}
_pMatrixddt[i*queue+j]=sum;
sum=0 ;
}
}
delete[] FKMatrix;
}
//求_pMatrix的逆矩阵，并存结果于矩阵pMatrix中
void ContraryMatrix(double *const _pMatrixctc, double *const pMatrixctc, const int &queue)
{
double *tMatrix = new double[2*queue*queue];
for (int i=0; i<queue; i++){
for (int j=0; j<queue; j++)
tMatrix[i*queue*2+j] = _pMatrixctc[i*queue+j];
}
for (int i=0; i<queue; i++){
for (int j=queue; j<queue*2; j++)
tMatrix[i*queue*2+j] = 0.0;
tMatrix[i*queue*2+queue+i] = 1.0;
}
//Initialization over!
for (int i=0; i<queue; i++)//Process Cols
{
double base = tMatrix[i*queue*2+i];
if (fabs(base) < 1E-300){
cout << endl << "求逆矩阵过程中被零除，无法求解!" << endl;
exit(0);
}
for (int j=0; j<queue; j++)//row
{
if (j == i) continue;
double times = tMatrix[j*queue*2+i]/base;
for (int k=0; k<queue*2; k++)//col
{
tMatrix[j*queue*2+k] = tMatrix[j*queue*2+k] - times*tMatrix[i*queue*2+k];
}
}
for (int k=0; k<queue*2; k++){
tMatrix[i*queue*2+k] /= base;
}
}
for (int i=0; i<queue; i++)
{
for (int j=0; j<queue; j++)
pMatrixctc[i*queue+j] = tMatrix[i*queue*2+j+queue];
}
delete[] tMatrix;
}
//求D总的矩阵积
void DianjiD(double *const _MatrixLdt,double *const MatrixLddt,double *const MatrixLdtddt, const int queue)
{
for (int i=0; i<queue; i++)
{ double sum=0;
for (int j=0; j<queue; j++)
{
for (int k=0; k<queue; k++)
{
sum=sum+(_MatrixLdt[i*queue+k]*MatrixLddt[j*queue+k]);
}
MatrixLdtddt[i*queue+j]=sum;
sum=0 ;
}
}
}
//求C总的矩阵积
void DianjiC(double *const _MatrixLct,double *const MatrixLctc,double *const MatrixLctcct, const int queue,const int line)
{
for (int i=0; i<queue; i++)
{ double sum=0;
for (int j=0; j<line; j++)
{
for (int k=0; k<queue; k++)
{
sum=sum+(MatrixLctc[i*queue+k]*_MatrixLct[j+k*line]);
}
MatrixLctcct[i*line+j]=sum;
sum=0 ;
}
}
}
//求D总和C总的矩阵积
void DianjiAjia(double *const MatrixLdtddt,double *const MatrixLctcct,double *const Ajia, const int queue,const int line)
{
for (int i=0; i<queue; i++)
{ double sum=0;
for (int j=0; j<line; j++)
{
for (int k=0; k<queue; k++)
{
sum=sum+(MatrixLdtddt[i*queue+k]*MatrixLctcct[j+k*line]);
}
Ajia[i*line+j]=sum;
sum=0 ;
}
}
}
DisplayMatrix(MatrixL, line, queue);
FenkuaiCDMatrix(MatrixL,_MatrixLdp, _MatrixLc,_MatrixLct,_MatrixLctc,_MatrixLdq,_MatrixLd,_MatrixLdt,_MatrixLddt,line,queue);//求分块矩阵
ContraryMatrix(_MatrixLctc, MatrixLctc, queue);//求逆
ContraryMatrix(_MatrixLddt, MatrixLddt, queue);//求逆
DianjiD(_MatrixLdt,MatrixLddt,MatrixLdtddt,queue);
DianjiC(_MatrixLct,MatrixLctc,MatrixLctcct,queue,line);
DianjiAjia(MatrixLdtddt,MatrixLctcct,Ajia,queue,line);
cout << "A的广义逆矩阵为：" << endl;
DisplayCzongMatrix(Ajia, queue,line);
delete[] _MatrixLc;
delete[] _MatrixLd;
网上下的一段广义逆矩阵程序，问题如下：
当取90行，3列时矩阵可求逆，同样元素数变成6行，45列的矩阵求逆时，则会在上标红处报错，提示内存访问冲突中断，这种问题一般是内存不够，可是其矩阵定义为const，具体该怎么改呢？
请高手指教具体怎么改？
[解决办法]
求广义逆矩阵最重要的是计算SVD分解
推荐http://my.oschina.net/zmjerry/blog/3774
