eignvalue.txt

一般实矩阵的特征值计算

#include<stdio.h>
#include<conio.h>
#include<float.h>
#include<math.h>
#define N 2
#define Eps 1.0e-15

typedef struct{
               double re;
               double im;
              }complex;

void main()
{
 int Matrix_qr();
 int rd_data();
 void Matrix_trans();

 double a[N][N],*a1[N];
// int i,j;
 int i;
 complex u[N];
 FILE *fp;
    char buffer[100];

	if((fp = fopen("transqr.dat", "wb")) == NULL)
	{
		printf("Can't open file\n");
		return;
	}

	gets(buffer);
	
	fputs(buffer, fp);
	if(ferror(fp))
		printf("Fail to write file\n");
	
	fclose(fp);

 for(i=0;i<N;i++)
     a1[i]=a[i];

 if(rd_data("transqr.dat",N,a1)==0)
    {
     Matrix_trans(N,a1,Eps);
     if(Matrix_qr(N,a1,u,Eps)==0)
        {
         printf("\nThe result is:\n");
            for(i=0;i<N;i++)
               {
                printf("%6.2f+i%6.2f",u[i].re,u[i].im);
                printf("\n");
               }
        }
        else
        printf("\nThe precision isn't enough.\n");
    }
    else
    printf("\007The datafile isn't matched.\n");
    getch();
}

int rd_data(tmp,n,x)
char *tmp;
int n;
double *x[];
{
 int i,j;
 FILE *fp;

 if((fp=fopen(tmp,"r"))==NULL)
    return(-1);
 for(j=0;j<n;j++)
    for(i=0;i<n;i++)
       fscanf(fp,"%lf",&x[j][i]);
 fclose(fp);

 return(0);
}

//#include<float.h>
//#include<math.h>

/*typedef struct{
              double re;
              double im;
              }complex; */

int Matrix_qr(n,a,u)
    int n;
    double *a[];
    /*double Eps;*/
    complex *u;
{
 int i,j,k,l,m,it;
 double b,c,xy,w,p,q,r,x,y,s,e,f,g,z;

 m=n-1;
 it=0;

 for(;;)
    {
     if(m==-1)
       return(0);
     l=m+1;
     do{
        l=l-1;
       }while(l>0&&fabs(a[l][l-1])>Eps*(fabs(a[l-1][l-1])+fabs(a[l][l])));
       if(l==m)
         {
          u[l].re=a[l][l];
          u[l].im=0.0;
          m=m-1;
          it=0;
         }
        else
         {
          if(l==m-1)
            {
             b=-(a[m][m]+a[m-1][m-1]);
             c=a[m][m]*a[m-1][m-1]-a[m][m-1]*a[m-1][m];
             w=b*b-4.0*c;
             y=sqrt(fabs(w));
             if(w>Eps)
               {
                xy=1.0;
                if(b<Eps)xy=-1.0;
                u[m].re=(-b-xy*y)/2.0;
                u[m-1].re=c/u[m].re;
                u[m].im=0.0;
                u[m-1].im=0.0;
               }
               else
                  {
                   u[m].re=-b/2.0;
                   u[m-1].re=u[m].re;
                   u[m].im=y/2.0;
                   u[m-1].im=-u[m].im;
                  }
                  m=m-2;
                  it=0;
            }
            else
               {
                if(it>200)
                  return(-1);
                it++;
                for(j=l+2;j<=m;j++)
                   a[j][j-2]=0.0;
                for(j=l+3;j<=m;j++)
                   a[j][j-3]=0.0;
                for(k=l;k<=m-1;k++)
                   {
                    if(k!=l)
                      {
                       p=a[k][k-1];
                       q=a[k+1][k-1];
                       r=0.0;
                       if(k!=m-1)
                          r=a[k+2][k-1];
                      }
                      else
                        {
                         x=a[m][m]+a[m-1][m-1];
                         y=a[m-1][m-1]*a[m][m]-a[m-1][m]*a[m][m-1];
                         p=a[l][l]*(a[l][l]-x)+a[l][l+1]*a[l+1][l]+y;
                         q=a[l+1][l]*(a[l][l]+a[l+1][l+1]-x);
                         r=a[l+1][l]*a[l+2][l+1];
                        }
                        if(fabs(p)+fabs(q)+fabs(r)>Eps)
                          {
                           xy=1.0;
                           if(p<Eps)
                             xy=-1.0;
                           s=xy*sqrt(p*p+q*q+r*r);
                           if(k!=l)
                             a[k][k-1]=-s;
                           e=-q/s;
                           f=-r/s;
                           x=-p/s;
                           y=-x-f*r/(p+s);
                           g=e*r/(p+s);
                           z=-x-e*q/(p+s);
                           for(j=k;j<=m;j++)
                              {
                               p=x*a[k][j]+e*a[k+1][j];
                               q=e*a[k][j]+y*a[k+1][j];
                               r=f*a[k][j]+g*a[k+1][j];
                               if(k!=m-1)
                               {
                                p+=f*a[k+2][j];
                                q+=g*a[k+2][j];
                                r+=z*a[k+2][j];
                                a[k+2][j]=r;
                               }
                               a[k+1][j]=q;
                               a[k][j]=p;
                              }
                           j=k+3;
                           if(j>=m)j=m;
                           for(i=l;i<=j;i++)
                              {
                               p=x*a[i][k]+e*a[i][k+1];
                               q=e*a[i][k]+y*a[i][k+1];
                               r=f*a[i][k]+g*a[i][k+1];
                               if(k!=m-1)
                                 {
                                  p+=f*a[i][k+2];
                                  q+=f*a[i][k+2];
                                  r+=z*a[i][k+2];
                                  a[i][k+2]=r;
                                 }
                                 a[i][k+1]=q;
                                 a[i][k]=p;
                              }
                          }
                   }
               }
         }
    }
}

//#include<float.h>
//#include<math.h>

void Matrix_trans(n,a)
int n;
double *a[];
/*double Eps;*/
{
 int k,j,i;
 double d,t;
 void dswap();

 for(k=1;k<n-1;k++)
    {
     d=0.0;
     for(j=k;j<n;j++)
        {
         if(fabs(a[j][k-1])>fabs(d))
            {
             d=a[j][k-1];
             i=j;
            }
        }
        if(fabs(d)>Eps)
          {
           if(i!=k)
              {
               for(j=k-1;j<n;j++)
                  dswap(&a[i][j],&a[k][j]);
               for(j=0;j<n;j++)
                  dswap(&a[j][i],&a[j][k]);
              }
              for(i=k+1;i<n;i++)
                 {
                  t=a[i][k-1]/d;
                  a[i][k-1]=0.0;
                  for(j=k;j<n;j++)
                     a[i][j]-=t*a[k][j];
                  for(j=0;j<n;j++)
                     a[j][k]+=t*a[j][i];
                 }
          }
    }
}

void dswap(x,y)
double *x,*y;
  {
   double temp;
   temp=*x;
   *x=*y;
   *y=temp;
  }