wd4.txt

这是本人编写的有关矩阵数值计算的几个C语言程序,方便大家学习<矩阵计算方法>!

/*算法3.1----p204*/
#include"stdio.h"
#include"math.h"
#define TEM 100
float H[TEM][TEM];

void data_in(int m,int n,float arr_In[TEM][TEM]);
void data_out(int m,int n,float arr_Out[TEM][TEM]);
void Housd_T(int k,int m,float arr_H[TEM][TEM]);
void arr_multip(int m,int n,int s,float arr_ms[TEM][TEM],
float arr_sn[TEM][TEM],float arr_mul[TEM][TEM]);
float max_N(int k,int n,int *out_q,float arr[TEM]);

main(void)
{
 int i,j,m,n,k,L,q,t;
 float temp;
 float a[TEM][TEM], p[TEM],gama[TEM],zeta[TEM][TEM];
 float arr_Temp[TEM][TEM],arr_Q[TEM][TEM];

 for(i=1;i<=TEM;i++)                                     /*给矩阵Q赋初值*/
    for(j=1;j<=TEM;j++)
       if(i==j)
          arr_Q[i][j]=1;
       else
          arr_Q[i][j]=0;
Loop:
 printf("\n\n请输入矩阵的行数m：\n");
 scanf("%d",&m);
 printf("\n\n请输入矩阵的列数n：\n");
 scanf("%d",&n);
 if(m<n)
   {
    printf("输入有误！");
    goto Loop;
   }
 data_in(m,n,a);
 k=1;
 for(j=1;j<=n;j++)
    {
     p[j]=j;
     temp=0;
     for(i=1;i<=m;i++)
	temp=(float)(temp+a[i][j]*a[i][j]);
     gama[j]=temp;
    }

Tag_0:
 temp=max_N(k,n,&q,gama);
 L=q;
 gama[L]=temp;
 if(gama[L]==0)
   goto Tag_1;
 temp=p[L];                                                  /*交换pl&pk*/
 p[L]=p[k];
 p[k]=temp;
 temp=gama[L];
 gama[L]=gama[k];
 gama[k]=temp;

 for(i=1;i<=m;i++)
    {
     temp=a[i][L];
     a[i][L]=a[i][k];
     a[i][k]=temp;
    }
 

 Housd_T(k,m,a);					   /*调用Householder变换，将分解后的R存入A中，得AP=Hr...H1R*/

 for(i=1;i<=k-1;i++)                                       /*将diag(I,Hk)存入H中*/
     for(j=1;j<=m;j++)
        {
	 if(i==j)
           H[i][j]=1;
         else
           H[i][j]=0;
        }
 for(i=k;i<=m;i++)
    for(j=1;j<=k-1;j++)
       H[i][j]=0;

 arr_multip(m,m,m,H,arr_Q,arr_Temp);                       /*调用矩阵相乘求Q,即将Householder阵相乘*/
 for(i=1;i<=m;i++)
    for(j=1;j<=m;j++)
       arr_Q[i][j]=arr_Temp[i][j];

 for(i=1;i<=m;i++)                                         /*将diag(I,Hk)与A相乘*/
    for(j=1;j<=n;j++)
       {
        zeta[i][j]=0;
        for(t=1;t<=m;t++)
	   zeta[i][j]=(float)(zeta[i][j]+H[i][t]*a[t][j]);
       }
 for(i=1;i<=m;i++)
    for(j=1;j<=n;j++)
       a[i][j]=zeta[i][j];
 for(j=k+1;j<=n;j++)
    gama[j]=(float)(gama[j]-a[k][j]*a[k][j]);
 k=k+1;
 goto Tag_0;

Tag_1:                                                                       /*调用输出函数输出Q及R*/
 printf("\n\n输出QR分解中的R阵各元素: \n");
 data_out(m,n,a); 

 for(i=1;i<=m;i++)                                                            /*将求得值转置得Q值*/
    for(j=1;j<=m;j++)
       arr_Temp[i][j]=arr_Q[j][i]; 
                                              
 printf("\n\n输出QR分解中的Q阵各元素: \n");
 data_out(m,m,arr_Temp);                                                     
 
 m=1;
 n=2;

}

/*************************************************************************/
void data_out(int m,int n,float arr_Out[TEM][TEM])            /*数据输出函数*/
{
 int i,j;
 
 for(i=1;i<=m;i++)
    {
     for(j=1;j<=n;j++)
     printf("%.3f  ",arr_Out[i][j]);
     printf("\n");
    }
 
 return;
}

/*************************************************************************/
void data_in(int m,int n,float arr_In[TEM][TEM])              /*数据输入函数*/
{
 int i,j;
 printf("\n\n请输入mxn阵各元素: \n");
 for(i=1;i<=m;i++)
    {
     for(j=1;j<=n;j++)
	{
	 printf("Input a[%d][%d]> ",i,j);
	 scanf("%f",&arr_In[i][j]);
	}
    }
 return;
}

/*************************************************************************/
float max_N(int k,int n,int *out_q,float arr[TEM])            /*求最大值函数*/
{
 int i;
 float temp;
 temp=arr[k];
 *out_q=k;
 for(i=k;i<=n;i++)
    if(temp<arr[i])
      {
       temp=arr[i];
       *out_q=i;
      }
 return (temp);
}

/*************************************************************************/
void Housd_T(int k,int m,float arr_H[TEM][TEM])             /*构造一个Householder变换*/
{
 int i,j,d;
 float y,beta,alfa,temp_H;
 float s[TEM],tem[TEM];
 for(i=k;i<=m;i++)
    s[i]=arr_H[i][k];
 y=max_N(k,m,&d,s);

 if(y!=0)
   {
    for(i=k;i<=m;i++)
       s[i]=(float)(s[i]/y);
    alfa=0;
    for(i=k;i<=m;i++)
       alfa=(float)(alfa+s[i]*s[i]);
    alfa=(float)pow((float)alfa,0.5);

    if(s[k]<0)
      alfa=-alfa;
    s[k]=(float)(s[k]+alfa);
    beta=(float)(1/(alfa*s[k]));
    alfa=(float)(y*alfa);
   }
 else
   beta=0;                         /* H=I-beta*(s[k],s[k+1],...,s[m])'(s[k],s[k+1],...,s[m])*/
 for(i=k;i<=m;i++)
    {
     for(j=k;j<=m;j++)
	if(i==j)
          H[i][j]=1-beta*s[i]*s[j];
        else
	  H[i][j]=-beta*s[i]*s[j];
    }


 printf("\n\n(k=%d)检验Householder变换的正确性:\n",k);
 for(i=k;i<=m;i++)
    {
     temp_H=0;
     for(j=k;j<=m;j++)
	temp_H=temp_H+H[i][j]*arr_H[j][k];
     tem[i]=temp_H;
     printf("%.4f ",tem[i]);
    }

 return;
}
/*************************************************************************/

void arr_multip(int m,int n,int s,float arr_ms[TEM][TEM],
float arr_sn[TEM][TEM],float arr_mul[TEM][TEM])
{                                                                          /*两个矩阵相乘*/
 int i,j,t;
 float temp;
 for(i=1;i<=m;i++)
     for(j=1;j<=n;j++)
        {
         temp=0;
         for(t=1;t<=s;t++)
            temp=temp+arr_ms[i][t]*arr_sn[t][j];
         arr_mul[i][j]=temp;
        }
 return;
}


A:
2 2 -1 
2 -1 2
-1 2 2
1 0 -4
4 3 2
R:
5.385 0.743 1.114 
0.000 -5.045 -2.215
0.000 0.000 -3.443
0.000 0.000 0.000
0.000 0.000 0.000

Q:
-0.186 -0.424 -0.368 -0.683 -0.428
0.371 -0.342 0.630 0.131 0.575 
0.371 0.253 -0.623 0.394 -0.504 
-0.743 -0.308 -0.0043 0.552 -0.217
0.371 -0.738 -0.276 0.236 0.430