帮忙用c 语言把矩阵进行QR分解,写出源程序

帮忙用c 语言把矩阵进行QR分解,写出源程序,第1张

void QR(double a[N][N],double q[N][N],double r1[N][N],int n) /*QR分解*/

{

int i,j,k,r,m

double temp,sum,dr,cr,hr

double ur[N],pr[N],wr[N]

double q1[N][N],emp[N][N]

for(i=1i<n+1i++)//将a放入temp中

for(j=1j<n+1j++)

{

emp[i][j]=a[i][j]

}

for(i=1i<n+1i++)//定义单位矩阵

for(j=1j<n+1j++)

{

if(i==j)q[i][j]=1

else q[i][j]=0

}

for(r=1r<nr++)

{

temp=0

for(k=r+1k<n+1k++)

temp+=fabs(a[k][r])

if(temp>=ZERO)

{

sum=0

for(k=rk<n+1k++)

sum+=a[k][r]*a[k][r]

dr=sqrt(sum)

if(a[r][r]>ZERO)m=-1

else m=1

cr=m*dr

hr=cr*(cr-a[r][r])

for(i=1i<n+1i++)//定义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=1i<n+1i++)//定义wr

{

sum=0

for(j=1j<n+1j++)

sum+=q[i][j]*ur[j]

wr[i]=sum

}

for(i=1i<n+1i++)//定义qr

for(j=1j<n+1j++)

{

q1[i][j]=q[i][j]-wr[i]*ur[j]/hr

}

for(i=1i<n+1i++)//定义qr+1

for(j=1j<n+1j++)

{

q[i][j]=q1[i][j]

}

for(i=1i<n+1i++)//定义pr

{

sum=0

for(j=1j<n+1j++)

sum+=a[j][i]*ur[j]

pr[i]=sum/hr

}

for(i=1i<n+1i++)

for(j=1j<n+1j++)

{

a[i][j]=a[i][j]-ur[i]*pr[j]

}

}

}

for(i=1i<n+1i++)

for(j=1j<n+1j++)

{

if(fabs(a[i][j])<ZERO)a[i][j]=0

}

for(i=1i<n+1i++)

for(j=1j<n+1j++)

{

r1[i][j]=a[i][j]

}

for(i=1i<n+1i++)//将a取出

for(j=1j<n+1j++)

{

a[i][j]=emp[i][j]

}}

这是我编的一个子函数,直接调用就可以了

给你一段matlab代码,你可以照着改成 c++

    n = size(V,1)

    k = size(V,2)

    U = zeros(n,k)

    U(:,1) = V(:,1)/sqrt(V(:,1)'*V(:,1))

    for i = 2:k

      U(:,i) = V(:,i)

      for j = 1:i-1

        U(:,i) = U(:,i) - ( U(:,i)'*U(:,j) )/( U(:,j)'*U(:,j) )*U(:,j)

      end

      U(:,i) = U(:,i)/sqrt(U(:,i)'*U(:,i))

    end

楼主的问题是自己写程序完成矩阵的QR分解,既然是迭代实现QR分解,就与矩阵论中说的计算特征值和特征向量的方法有些区别了。大体的步骤应该是首先将矩阵化成双对角矩阵,然后追赶计算特征值和特征向量,程序代码可以参考 徐士良编的 常用数值算法 c语言描述


欢迎分享,转载请注明来源:内存溢出

原文地址:https://54852.com/yw/12196521.html

(0)
打赏 微信扫一扫微信扫一扫 支付宝扫一扫支付宝扫一扫
上一篇 2023-05-21
下一篇2023-05-21

发表评论

登录后才能评论

评论列表(0条)

    保存