# 哪位高手看一下，这是输出什么？

skycncomp 2004-05-08 06:39:53

#include "aband.c"
#include "stdio.h"
main()
{ int i;
static double b[8][5]={ {3.0,-4.0,1.0,0.0,0.0},
{-2.0,-5.0,6.0,1.0,0.0},{1.0,3.0,-1.0,2.0,-3.0},
{2.0,5.0,-5.0,6.0,-1.0},{-3.0,1.0,-1.0,2.0,-5.0},
{6.0,1.0,-3.0,2.0,-9.0},{-4.0,1.0,-1.0,2.0,0.0},
{5.0,1.0,-7.0,0.0,0.0}};
static double d[8][3]={ {13.0,29.0,-13.0},
{-6.0,17.0,-21.0},{-31.0,-6.0,4.0},{64.0,3.0,16.0},
{-20.0,1.0,-5.0},{-22.0,-41.0,56.0},{-29.0,10.0,-21.0},
{7.0,-24.0,20.0}};
if (aband(b,d,8,2,5,3)>0)
for (i=0;i<=7;i++)
printf("x(%d)=%13.7e, %13.7e, %13.7e\n",
i,d[i][0],d[i][1],d[i][2]);
}

#include "math.h"
#include "stdio.h"
int aband(b,d,n,l,il,m)
int n,l,il,m;
double b[],d[];
{ int ls,k,i,j,is,u,v;
double p,t;

if (il!=(2*l+1))
{ printf("fail\n"); return(-2);}
ls=l;

for (k=0;k<=n-2;k++)
{
p=0.0;
for (i=k;i<=ls;i++)
{
t=fabs(b[i*il]);
if (t>p) {p=t; is=i;}
}
if (p+1.0==1.0)
{
printf("fail\n"); return(0);
}

for (j=0;j<=m-1;j++)
{
u=k*m+j; v=is*m+j;
t=d[u]; d[u]=d[v]; d[v]=t;
}

for (j=0;j<=il-1;j++)
{ u=k*il+j; v=is*il+j;
t=b[u]; b[u]=b[v]; b[v]=t;
}
for (j=0;j<=m-1;j++)
{ u=k*m+j; d[u]=d[u]/b[k*il];}
for (j=1;j<=il-1;j++)
{ u=k*il+j; b[u]=b[u]/b[k*il];}
for (i=k+1;i<=ls;i++)
{ t=b[i*il];
for (j=0;j<=m-1;j++)
{ u=i*m+j; v=k*m+j;
d[u]=d[u]-t*d[v];
}
for (j=1;j<=il-1;j++)
{ u=i*il+j; v=k*il+j;
b[u-1]=b[u]-t*b[v];
}
u=i*il+il-1; b[u]=0.0;
}
if (ls!=(n-1)) ls=ls+1;
}
p=b[(n-1)*il];
if (fabs(p)+1.0==1.0)
{ printf("fail\n"); return(0);}
for (j=0;j<=m-1;j++)
{ u=(n-1)*m+j; d[u]=d[u]/p;}
ls=1;
for (i=n-2;i>=0;i--)
{ for (k=0;k<=m-1;k++)
{ u=i*m+k;
for (j=1;j<=ls;j++)
{ v=i*il+j; is=(i+j)*m+k;
d[u]=d[u]-b[v]*d[is];
}
}
if (ls!=(il-1)) ls=ls+1;
}
return(2);
}

...全文
41 21 打赏 收藏 举报

21 条回复

improgrammer 2004-05-11

{
{ 3.0,-4.0, 1.0, 0.0, 0.0},
{-2.0,-5.0, 6.0, 1.0, 0.0},
{ 1.0, 3.0,-1.0, 2.0,-3.0},
{ 2.0, 5.0,-5.0, 6.0,-1.0},
{-3.0, 1.0,-1.0, 2.0,-5.0},
{ 6.0, 1.0,-3.0, 2.0,-9.0},
{-4.0, 1.0,-1.0, 2.0, 0.0},
{ 5.0, 1.0,-7.0, 0.0, 0.0}
};

+3.000 +5.000 +0.000
-1.000 -3.000 +3.000
+0.000 +2.000 -1.000
-5.000 -0.000 +0.000
+7.000 -0.000 +2.000
+1.000 +1.000 -3.000
+2.000 -1.000 +0.000
+0.000 +4.000 -5.000

{
{ 13.0, 29.0,-13.0},
{ -6.0, 17.0,-21.0},
{-31.0, -6.0, 4.0},
{ 64.0, 3.0, 16.0},
{-20.0, 1.0, -5.0},
{-22.0,-41.0, 56.0},
{-29.0, 10.0,-21.0},
{ 7.0,-24.0, 20.0}
};

• 打赏
• 举报

improgrammer 2004-05-11

• 打赏
• 举报

skycncomp 2004-05-11

• 打赏
• 举报

improgrammer 2004-05-10
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

//判断p是否充分接近零
int near_zero(double p)
{
//——只要浮点数p与1.0对齐后没有有效数字，就判p充分接近零
return p+1.0==1.0;
}

//打印两矩阵
void my_print(double b[],double d[],const int N,const int IL,const int M)
{
int i,j;
for (i=0;i<N;i++)
{
printf("row[%d]",i);
for(j=0;j<IL;++j)
{
printf(" %+7.3f",b[i*IL+j]);
}
printf(" |");
for(j=0;j<M;++j)
{
printf(" %+7.3f",d[i*M+j]);
}
printf("\n");
}
printf("\n");
}

//向量点积
double vec_dm(const double *v1,const double *v2,int d1,int d2,int n)
{
double r=0.0;
while((n--) > 0)
{
r += (*v1)*(*v2);
v1+=d1; v2+=d2;
}
return r;
}

//v1 <--> v2
void vec_swap(double *v1,double *v2,int d1,int d2,int n)
{
double t;
while((n--) > 0)
{
t=*v1;*v1=*v2;*v2=t;
v1+=d1;v2+=d2;
}
}

// vec *= scale
void vec_scale(double *vec,int d,int n,double scale)
{
while((n--) > 0)
{
*vec *= scale;
vec+=d;
}
}

// v1 += v2 * scale
void vec_scale_add(double *v1,const double *v2,int d1,int d2,int n,double scale)
{
while((n--) > 0)
{
(*v1) += (*v2) * scale;
v1+=d1;v2+=d2;
}
}

//求向量中绝对值最大元素之索引
int vec_max_index(const double *v,const int d,const int n)
{
int i,r=0;
double t,p=0.0;
for(i=0;i<n;++i)
{
t=fabs(*v);
if(t>p)
{
r=i;
p=t;
}
v+=d;
}
return r;
}

//向量左移一个元素，最左元素被丢弃，最右元素写为0
void vec_lshift(double *v,int d,int n)
{
while((n--)>1)
{
*v = *(v+d);
v+=d;
}
*v=0.0;
}

//哪位高手看一下，这是输出什么？
//——应该是解方程
int aband(double b[],double d[],const int N,const int IL,const int M)
{
const int L=(IL+1)/2;
int i,j,is,_ls;
double p;
//第一步：好象是做一个消元法，完成后b数组的第一列为主对角线元素
for (i=0;i<N-1;i++)
{
_ls=min(i+L,N);
//从b[i][0],b[i+1][0]..b[_ls-1][0]找最大元素[绝对值最大]
//——最大元素的索引保存在is中
is=i+vec_max_index(b+i*IL,IL,_ls-i);
p=b[is*IL];
assert(!near_zero(p));//最大元素应不为0
//对调第i,is行
vec_swap(d+i*M,d+is*M,1,1,M);
vec_swap(b+i*IL,b+is*IL,1,1,IL);
//d[N][M]第i行除以b[N][IL]的第i行排头元素
vec_scale(d+i*M,1,M,1/p);
//b[N][M]第i行（排头除外）除以排头元素
vec_scale(b+i*IL+1,1,IL-1,1/p);
//
for (j=i+1;j<_ls;j++)
{
p=b[j*IL];
vec_lshift(b+j*IL,1,IL);
}
}
//d[N-1][0..M-1] /= p;
p=b[(N-1)*IL];
assert(!near_zero(fabs(p)));
vec_scale(d+(N-1)*M,1,M,1/p);
//第二步：某种卷积运算
for (i = N - 2;i >= 0;i--)
{
const int R = min(N-i-1,IL-1);
for (j = 0;j < M;j++)
{
d[i*M+j] -= vec_dm(b+i*IL+1,d+(i+1)*M+j,1,M,R);
}
}
return(2);
}

void work()
{
static double b[8][5]=
{
{ 3.0,-4.0, 1.0, 0.0, 0.0},
{-2.0,-5.0, 6.0, 1.0, 0.0},
{ 1.0, 3.0,-1.0, 2.0,-3.0},
{ 2.0, 5.0,-5.0, 6.0,-1.0},
{-3.0, 1.0,-1.0, 2.0,-5.0},
{ 6.0, 1.0,-3.0, 2.0,-9.0},
{-4.0, 1.0,-1.0, 2.0, 0.0},
{ 5.0, 1.0,-7.0, 0.0, 0.0}
};
static double d[8][3]=
{
{ 13.0, 29.0,-13.0},
{ -6.0, 17.0,-21.0},
{-31.0, -6.0, 4.0},
{ 64.0, 3.0, 16.0},
{-20.0, 1.0, -5.0},
{-22.0,-41.0, 56.0},
{-29.0, 10.0,-21.0},
{ 7.0,-24.0, 20.0}
};
if (aband(b,d,8,5,3)>0)
{
my_print(b,d,8,5,3);
}
}
• 打赏
• 举报

skycncomp 2004-05-10

• 打赏
• 举报

skycncomp 2004-05-10

• 打赏
• 举报

improgrammer 2004-05-10

B*X=D
B是N*N方阵
X是N*M矩阵
D是N*M矩阵

B是对角方阵（宽度是IL），所以b[N][IL]采用了压缩存储法。

• 打赏
• 举报

skycncomp 2004-05-09

• 打赏
• 举报

holyfair 2004-05-09

x(0)=3.0000000e+000, 5.0000000e+000, -1.2767565e-015
x(1)=-1.0000000e+000, -3.0000000e+000, 3.0000000e+000
x(2)=8.3266727e-016, 2.0000000e+000, -1.0000000e+000
x(3)=-5.0000000e+000, -2.6645353e-015, -1.7763568e-01
x(4)=7.0000000e+000, -8.8817842e-016, 2.0000000e+000
x(5)=1.0000000e+000, 1.0000000e+000, -3.0000000e+000
x(6)=2.0000000e+000, -1.0000000e+000, -8.8817842e-016
x(7)=1.5199238e-016, 4.0000000e+000, -5.0000000e+000

• 打赏
• 举报

skycncomp 2004-05-09

• 打赏
• 举报

aband.c原文件在那里
• 打赏
• 举报

liangenbo 2004-05-09

• 打赏
• 举报

skycncomp 2004-05-09

• 打赏
• 举报

improgrammer 2004-05-09

• 打赏
• 举报

skycncomp 2004-05-09

• 打赏
• 举报

julyclyde 2004-05-09

• 打赏
• 举报

cancer001 2004-05-08

• 打赏
• 举报

wangnewton 2004-05-08

ps：一个注释都没有

• 打赏
• 举报

sten 2004-05-08

• 打赏
• 举报

tianyaht 2004-05-08

• 打赏
• 举报

C语言

6.6w+

C语言相关问题讨论

2004-05-08 06:39