void
QR(double
**a,int a_row,int
a_colunm,double** q,double
**r1)
{
int n=0;
if (a_row>a_colunm){n=a_colunm;}else{n=a_row;};
int i,j,k,r,m;
double temp,sum,dr,cr,hr;
double *ur=new
double[n];
double *pr=new
double[n];
double *wr=new
double[n];
double **q1=new
double*[n];
double **emp=new
double*[n];
for (int i=0;i<n;i++)
{
q1[i]=new
double[n];
emp[i]=new
double[n];
}
for(i=0;i<a_row;i++)//将a放入temp中
for(j=0;j<a_colunm;j++)
{
emp[i][j]=a[i][j];
};
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(a[k][r]);
if(temp>=0.0)
{
sum=0;
for(k=r;k<n;k++)
sum+=a[k][r]*a[k][r];
dr=sqrt(sum);
if(a[r][r]>0.0)m=-1;
else m=1;
cr=m*dr;
hr=cr*(cr-a[r][r]);
for(i=0;i<n;i++)//定义ur
{
if(i<r)ur[i]=0;
if(i==r)ur[i]=a[r][r]-cr;
if(i>r)ur[i]=a[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<a_row;i++)//定义qr
for(j=0;j<a_colunm;j++)
{
q1[i][j]=q[i][j]-wr[i]*ur[j]/hr;
};
for(i=0;i<a_row;i++)//定义qr+1
for(j=0;j<a_colunm;j++)
{
q[i][j]=q1[i][j];
};
for(i=0;i<a_colunm;i++)//定义pr
{
sum=0;
for(j=0;j<n;j++)
sum+=a[j][i]*ur[j];
pr[i]=sum/hr;
};
for(i=0;i<a_row;i++)
for(j=0;j<a_colunm;j++)
{
a[i][j]=a[i][j]-ur[i]*pr[j];
};
};
};
for(i=0;i<a_row;i++)
for(j=0;j<a_colunm;j++)
{
if(fabs(a[i][j])<0.0)a[i][j]=0;
};
for(i=0;i<a_row;i++)
for(j=0;j<a_colunm;j++)
{
r1[i][j]=a[i][j];
};
for(i=0;i<a_row;i++)//将a取出
for(j=0;j<a_colunm;j++)
{
a[i][j]=emp[i][j];
}}
int
main()
{
double ** A=new
double*[3];
for (int i=0;i<3;i++)
{
A[i]=new
double[4];
}
for (int i=0;i<3;i++)
{
for (int j=0;j<4;j++)
{
A[i][j]=1+i*i+j+9;
}
}
printf("A=\n");
for (int i=0;i<3;i++)
{
for (int j=0;j<4;j++)
{
printf("%
10f",A[i][j]);
}
printf("\n");
}
double **qrtq=new
double*[3];
double **qrtt=new
double*[3];
for (int i=0;i<3;i++)
{
qrtq[i]=new
double[4];
qrtt[i]=new
double[4];
}
QR(A,3,4,qrtq,qrtt);
printf("qrqt=\n");
for (int i=0;i<3;i++)
{
for (int j=0;j<3;j++)
{
printf("%12f",qrtq[i][j]);
}
printf("\n");
}
printf("qrtt=\n");
for (int i=0;i<3;i++)
{
for (int j=0;j<4;j++)
{
printf("%12f",qrtt[i][j]);
}
printf("\n");
}
system("pause");
return 0;
}
void
QR(double
*a,int a_row,int
a_colunm,double* q,double
*r1)
{
int i,j,k,r,m;
double temp,sum,dr,cr,hr;
double *ur=new
double[a_row*a_row];
double *pr=new
double[a_row*a_row];
double *wr=new
double[a_row*a_row];
double *q1=new
double[a_row*a_row];
double *emp=new
double[a_row*a_colunm];
for(i=0;i<a_row;i++)//将a放入temp中
for(j=0;j<a_colunm;j++)
{
emp[i*a_colunm+j]=a[i*a_colunm+j];
};
for(i=0;i<a_row;i++)//定义单位矩阵
for(j=0;j<a_row;j++)
{
if(i==j)q[i*a_row+j]=1;
else q[i*a_row+j]=0;
};
for(r=0;r<a_colunm;r++)
{
temp=0;
for(k=r;k<a_row;k++)
temp+=fabs(a[k*a_colunm+r]);
if(temp>=0.0)
{
sum=0;
for(k=r;k<a_row;k++)
sum+=a[k*a_colunm+r]*a[k*a_colunm+r];
dr=sqrt(sum);
if(a[r*a_colunm+r]>0.0)m=-1;
else m=1;
cr=m*dr;
hr=cr*(cr-a[r*a_colunm+r]);
for(i=0;i<a_row;i++)//定义ur
{
if(i<r)ur[i]=0;
if(i==r)ur[i]=a[r*a_colunm+r]-cr;
if(i>r)ur[i]=a[i*a_colunm+r];
};
for(i=0;i<a_row;i++)//定义wr
{
sum=0;
for(j=0;j<a_row;j++)
sum+=q[i*a_row+j]*ur[j];
wr[i]=sum;
};
for(i=0;i<a_row;i++)//定义qr
for(j=0;j<a_row;j++)
{
q1[i*a_row+j]=q[i*a_row+j]-wr[i]*ur[j]/hr;
};
for(i=0;i<a_row;i++)//定义qr+1
for(j=0;j<a_row;j++)
{
q[i*a_row+j]=q1[i*a_row+j];
};
for(i=0;i<a_row;i++)//定义pr
{
sum=0;
for(j=0;j<a_row;j++)
sum+=a[j*a_colunm+i]*ur[j];
pr[i]=sum/hr;
};
for(i=0;i<a_row;i++)
for(j=0;j<a_colunm;j++)
{
a[i*a_colunm+j]=a[i*a_colunm+j]-ur[i]*pr[j];
};
};
};
for(i=0;i<a_row;i++)
for(j=0;j<a_colunm;j++)
{
if(fabs(a[i*a_colunm+j])<0.0)a[i*a_colunm+j]=0;
};
for(i=0;i<a_row;i++)
for(j=0;j<a_colunm;j++)
{
r1[i*a_colunm+j]=a[i*a_colunm+j];
};
for(i=0;i<a_row;i++)//将a取出
for(j=0;j<a_colunm;j++)
{
a[i*a_colunm+j]=emp[i*a_colunm+j];
}
}
加载中…
喜欢