Ax=B，改写成Ly=B，Ux=y的方谨空旁程组。就相当于将A=LU分解成了两个矩阵。称为矩阵A的三角分解，或LU分解。如果L为单位下三角阵，则叫Doolittle分解，若祥橡U为单位上三角阵，则叫Crout分解。只要A的各顺序主子式不为零，则A可唯一分解成一个单位下三角阵L与一亏码个上三角阵U的乘积。
•设Ax=b，A=LU，则Ax=LUx=b
于是令Ux=y，则Ly=b
这样原来方程能化为两个简单方程组
下面是LU分解的Fortran子程序 希望可以有所帮助
!求解au=b,u
!n表示为方程维数
subroutine lu(a,b,n,u)
implicit real(8) (a-h,o-z)
real(8)::a(n,n),b(n),u(n),a_bak(n,n),b1(n),aL(n,n),aU(n,n),y(n)
!exchange rows
do i=1,n
tmpMax=0.d0
do ic=i,n
if(tmpMax<dabs(a(ic,i))) then
tmpMax=dabs(a(ic,i))
i_rec=ic
endif
enddo
if(i_rec.ne.i) then
do jc=i,n
tmp=a(i,jc)
a(i,jc)=a(i_rec,jc)
a(i_rec,jc)=tmp
enddo
tmp=b(i)
b(i)=b(i_rec)
b(i_rec)=tmp
endif
!decomposition
do j=i,n
tmp=0.d0
do k=1,i-1
tmp=tmp+aL(i,k)*aU(k,j)
enddo
aU(i,j)=a(i,j)-tmp
tmp=0.d0
do k=1,i-1
tmp=tmp+aL(j,k)*aU(k,i)
enddo
aL(j,i)=(a(j,i)-tmp)/aU(i,i)
enddo
enddo
!find answer
do i=1,n
tmp=0.d0
do j=1,i-1
tmp=tmp+aL(i,j)*y(j)
enddo
y(i)=b(i)-tmp
enddo
do i=n,1,-1
tmp=0.d0
do j=i+1,n
tmp=tmp+aU(i,j)*u(j)
enddo
u(i)=(y(i)-tmp)/aU(i,i)
enddo
end
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]
}}
这是我编的一个子函数，返斗戚直接调用就可以了
给，下凯空迟面是Cholesky分解法的C++经典算法：//-------------------------------------------------------------------
// Cholesky分解法
//-------------------------------------------------------------------
template <class T>
int cholesky(Matrix<T>& mat, double epsilon=EPSILON) {
size_t i, j, k
for (i=0 i<mat.Rows() ++i) {
// 计算第 i 轮主元
for (k=0 k<i ++k) mat[i][i] -= mat[k][i]*mat[k][i]
mat[i][i] = sqrt(mat[i][i])
// 计算第 i 轮主元结束
if (fabs(mat[i][i])<epsilon) break
//计算第 i 列
// for (j=i+1 j<mat.Rows() ++j) {
// for (k=0 k<i ++k) mat[j][i] -= mat[j][k]*mat[i][k]
// }
// 计算第 i 列结束
/盯李/ 计算第 i 行
for (j=i+1 j<mat.Cols() ++j) {
for (k=0 k<i ++k) mat[i][j] -= mat[k][i]*mat[k][j]
mat[i][j] /= mat[i][i]
}
// 计算第 i 行结束
}
return (i==mat.Rows())
}
下面这是另一种实现方法，输入输出语句自己根据需要写一下吧：
#include<malloc.h>
#include<math.h>
void cholesky(double **a,double *b,int n,double *x)
{
int i,j,m,k
double **L
L=(double **)malloc(n*sizeof(double))
for(i=0i<ni++)
L[i]=(double *)malloc(n*sizeof(double))
for(i=0i<ni++)
for(j=0j<nj++)
L[i][j]=0
for(k=0k<nk++)
{
L[k][k]=a[k][k]
for(m=0m<亏判km++)
L[k][k]-=L[k][m]*L[k][m]
L[k][k]=sqrt(L[k][k])
for(i=k+1i<ni++)
{
L[i][k]=a[i][k]
for(m=0m<km++)
L[i][k]-=L[i][m]*L[k][m]
L[i][k]/=L[k][k]
}
}
for(i=0i<ni++)
{
x[i]=b[i]
for(m=0m<im++)
x[i]-=L[i][m]*x[m]
x[i]/=L[i][i]
}
for(i=n-1i>=0i--)
{
for(m=i+1m<nm++)
x[i]-=L[m][i]*x[m]
x[i]/=L[i][i]
}
}
void main()
{
int i,j,n=3
double **a,*b,*x
a=(double **)malloc(n*sizeof(double))
for(i=0i<ni++)
a[i]=(double *)malloc(n*sizeof(double))
for(i=0i<ni++)
for(j=0j<nj++)
a[i][j]=1
for(i=0i<ni++)
a[i][i]=100
b=(double *)malloc(n*sizeof(double))
x=(double *)malloc(n*sizeof(double))
for(i=0i<ni++)
b[i]=3
cholesky(a,b,n,x)
}
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）