主题：有程序能获得LU分解的L和U吗？
timwan
[专家分：0] 发布于 2010-03-22 11:46:00
使用过IMSL的库函数LFTCG，但是输出只有1个矩阵，并没有输出L和U。
我按照IMSL说明文档也没能从这1个矩阵凑出L和U，找不到L*U等于原来矩阵A的。
比如
A(1,1:3)=(/1,2,1/)
A(2,1:3)=(/2,2,3/)
A(3,1:3)=(/-1,-3,0/)
调用LFTCG进行LU分解后输出：
B(1,1:3)=(/2,2,3/)
B(2,1:3)=(/-0.5,-2,1.5/)
B(3,1:3)=(/0.5,0.5,0.25/)
从B里我怎么也找不到L和U的信息。
是什么原因呢？谢谢！
回复列表 （共8个回复）
沙发
trimtrim [专家分：1640] 发布于 2010-03-22 17:14:00
看手册，觉得应该可以的。
板凳
anney169 [专家分：130] 发布于 2010-03-24 10:36:00
找 徐士良编著的《FORTRAN常用算法程序集》这本书，里面有源代码。
若没找到给我发邮件annie16934@yahoo.com.cn
3 楼
doctorlive [专家分：800] 发布于 2010-03-25 17:47:00
这是QU分解
请输入系数矩阵的行数、列数
3 3
请输入系数矩阵(按行输入)
1 2 1
2 2 3
-1 -3 0
3
2.000000 2.000000 3.000000
0.0000000E+00 -2.000000 1.500000
0.0000000E+00 0.0000000E+00 0.2500000
输出Q
1.000000 0.0000000E+00 0.0000000E+00
0.0000000E+00 -1.000000 0.0000000E+00
0.0000000E+00 0.0000000E+00 1.000000
输出U
2.000000 2.000000 3.000000
0.0000000E+00 2.000000 -1.500000
0.0000000E+00 0.0000000E+00 0.2500000
Press any key to continue
4 楼
timwan [专家分：0] 发布于 2010-03-26 15:01:00
[quote]找 徐士良编著的《FORTRAN常用算法程序集》这本书，里面有源代码。
若没找到给我发邮件annie16934@yahoo.com.cn[/quote]
谢谢，我去找找看。
5 楼
timwan [专家分：0] 发布于 2010-03-26 15:03:00
[quote]这是QU分解
请输入系数矩阵的行数、列数
3 3
请输入系数矩阵(按行输入)
1 2 1
2 2 3
-1 -3 0
3
2.000000 2.000000 3.000000
0.0000000E+00 -2.000000 1.500000
0.0000000E+00 0.0000000E+00 0.2500000
输出Q
1.000000 0.0000000E+00 0.0000000E+00
0.0000000E+00 -1.000000 0.0000000E+00
0.0000000E+00 0.0000000E+00 1.000000
输出U
2.000000 2.000000 3.000000
0.0000000E+00 2.000000 -1.500000
0.0000000E+00 0.0000000E+00 0.2500000
Press any key to continue[/quote]
您好，QU是什么意思呢？
您的这个U正好是我程序输出矩阵的上三角，但是Q能拿来作L用么？
6 楼
doctorlive [专家分：800] 发布于 2010-03-30 19:26:00
这是一种满秩分解
我在网上搜了个C语言写的LU分解程序，自己未亲测，你试试看。
一个N阶矩阵的LU分解
键盘输入一个矩阵，对他进行LU分解。源程序可编译运行通过：
#define N 4
#include<stdio.h>
main()
{
int i,j,r,k;
float a[N][N],l[N][N],u[N][N];
float de=0;
printf("请输入%d*%d矩阵A:\n");
for(i=1;i<N;i++)
for(j=1;j<N;j++)
scanf("%f",&a[i][j]);
getchar();
for(j=1;j<N;j++)
{
u[1][j]=a[1][j];
l[j][1]=a[j][1]/u[1][1];
printf("l[%d][1] is %f\n",j,l[j][1]);
}
getchar();
for(r=2;r<N;r++)
{
for(j=r;j<N;j++)
{
for(k=1;k<r;k++)
{
de+=l[r][k]*u[k][j];
}
u[r][j]=a[r][j]-de;
de=0;
printf("u[%d][%d] is %f\n",r,j,u[r][j]);
}
getchar();
for(i=r+1;i<N;i++)
{
for(k=1;k<r;k++)
{
de+=l[i][k]*u[k][r];
}
l[i][r]=a[i][r]-de;
l[i][r]=l[i][r]/u[r][r];
de=0;
printf("l[%d][%d] is %f\n",i,r,l[i][r]);
}
}
printf("The L matrix is:\n");
for(i=1;i<N;i++)
{
l[i][i]=1;
}
for(i=1;i<N;i++)
{
for(j=1;j<N;j++)
printf("%f ",l[i][j]);
printf("\n");
}
printf("The u matrix is:\n");
for(i=1;i<N;i++)
{
for(j=1;j<N;j++)
printf("%f ",u[i][j]);
printf("\n");
}
getchar();
}
7 楼
timwan [专家分：0] 发布于 2010-03-31 00:26:00
[quote]
这是一种满秩分解
我在网上搜了个C语言写的LU分解程序，自己未亲测，你试试看。
一个N阶矩阵的LU分解
键盘输入一个矩阵，对他进行LU分解。源程序可编译运行通过：
#define N 4
#include<stdio.h>
main()
{
int i,j,r,k;
float a[N][N],l[N][N],u[N][N];
float de=0;
printf("请输入%d*%d矩阵A:\n");
for(i=1;i<N;i++)
for(j=1;j<N;j++)
scanf("%f",&a[i][j]);
getchar();
for(j=1;j<N;j++)
{
u[1][j]=a[1][j];
l[j][1]=a[j][1]/u[1][1];
printf("l[%d][1] is %f\n",j,l[j][1]);
}
getchar();
for(r=2;r<N;r++)
{
for(j=r;j<N;j++)
{
for(k=1;k<r;k++)
{
de+=l[r][k]*u[k][j];
}
u[r][j]=a[r][j]-de;
de=0;
printf("u[%d][%d] is %f\n",r,j,u[r][j]);
}
getchar();
for(i=r+1;i<N;i++)
{
for(k=1;k<r;k++)
{
de+=l[i][k]*u[k][r];
}
l[i][r]=a[i][r]-de;
l[i][r]=l[i][r]/u[r][r];
de=0;
printf("l[%d][%d] is %f\n",i,r,l[i][r]);
}
}
printf("The L matrix is:\n");
for(i=1;i<N;i++)
{
l[i][i]=1;
}
for(i=1;i<N;i++)
{
for(j=1;j<N;j++)
printf("%f ",l[i][j]);
printf("\n");
}
printf("The u matrix is:\n");
for(i=1;i<N;i++)
{
for(j=1;j<N;j++)
printf("%f ",u[i][j]);
printf("\n");
}
getchar();
}
[/quote]
谢谢您！我去试下~
8 楼
他寻欢欢 [专家分：90] 发布于 2010-05-09 16:22:00
[quote][quote]这是QU分解
请输入系数矩阵的行数、列数
3 3
请输入系数矩阵(按行输入)
1 2 1
2 2 3
-1 -3 0
3
2.000000 2.000000 3.000000
0.0000000E+00 -2.000000 1.500000
0.0000000E+00 0.0000000E+00 0.2500000
输出Q
1.000000 0.0000000E+00 0.0000000E+00
0.0000000E+00 -1.000000 0.0000000E+00
0.0000000E+00 0.0000000E+00 1.000000
输出U
2.000000 2.000000 3.000000
0.0000000E+00 2.000000 -1.500000
0.0000000E+00 0.0000000E+00 0.2500000
Press any key to continue[/quote]
您好，QU是什么意思呢？
您的这个U正好是我程序输出矩阵的上三角，但是Q能拿来作L用么？
[/quote]
不能，可惜，
L U分解适用于所有满秩矩阵，Q U分解是对称正定矩阵的正交分解，不同的，还是
前面的建议，研究一下，徐世亮的那本书，学习不要偷懒。
我来回复
您所在位置：