det.cpp

利用C语言实现矩阵的行列式的计算

#include "iostream"
#include "math.h"
using namespace std;
  void hhbg(int n,double *a)
  { int i,j,k,u,v;
    double d,t;
    for (k=1; k<=n-2; k++)
      { d=0.0;
        for (j=k; j<=n-1; j++)
          { u=j*n+k-1; t=a[u];
            if (fabs(t)>fabs(d))
              { d=t; i=j;}
          }
        if (fabs(d)+1.0!=1.0)
          { if (i!=k)
              { for (j=k-1; j<=n-1; j++)
                  { u=i*n+j; v=k*n+j;
                    t=a[u]; a[u]=a[v]; a[v]=t;
                  }
                for (j=0; j<=n-1; j++)
                  { u=j*n+i; v=j*n+k;
                    t=a[u]; a[u]=a[v]; a[v]=t;
                  }
              }
            for (i=k+1; i<=n-1; i++)
              { u=i*n+k-1; t=a[u]/d; a[u]=0.0;
                for (j=k; j<=n-1; j++)
                  { v=i*n+j;
                    a[v]=a[v]-t*a[k*n+j];
                  }
                for (j=0; j<=n-1; j++)
                  { v=j*n+k;
                    a[v]=a[v]+t*a[j*n+i];
                  }
              }
          }
      }
    return;
  }

int qrtt(int n,double *a,double *u,double *v,double eps,int jt)
  { int m,it,i,j,k,l,ii,jj,kk,ll;
    double b,c,w,g,xy,p,q,r,x,s,e,f,z,y;
    it=0; m=n;
    while (m!=0)
      { l=m-1;
        while ((l>0)&&(fabs(a[l*n+l-1])>eps*
	      (fabs(a[(l-1)*n+l-1])+fabs(a[l*n+l])))) l=l-1;
        ii=(m-1)*n+m-1; jj=(m-1)*n+m-2;
        kk=(m-2)*n+m-1; ll=(m-2)*n+m-2;
        if (l==m-1)
          { u[m-1]=a[(m-1)*n+m-1]; v[m-1]=0.0;
            m=m-1; it=0;
          }
        else if (l==m-2)
          { b=-(a[ii]+a[ll]);
            c=a[ii]*a[ll]-a[jj]*a[kk];
            w=b*b-4.0*c;
            y=sqrt(fabs(w));
            if (w>0.0)
              { xy=1.0;
                if (b<0.0) xy=-1.0;
                u[m-1]=(-b-xy*y)/2.0;
                u[m-2]=c/u[m-1];
                v[m-1]=0.0; v[m-2]=0.0;
              }
            else
              { u[m-1]=-b/2.0; u[m-2]=u[m-1];
                v[m-1]=y/2.0; v[m-2]=-v[m-1];
              }
            m=m-2; it=0;
          }
        else
          { if (it>=jt)
              { printf("fail\n");
                return(-1);
              }
            it=it+1;
            for (j=l+2; j<=m-1; j++)
              a[j*n+j-2]=0.0;
            for (j=l+3; j<=m-1; j++)
              a[j*n+j-3]=0.0;
            for (k=l; k<=m-2; k++)
              { if (k!=l)
                  { p=a[k*n+k-1]; q=a[(k+1)*n+k-1];
                    r=0.0;
                    if (k!=m-2) r=a[(k+2)*n+k-1];
                  }
                else
                  { x=a[ii]+a[ll];
                    y=a[ll]*a[ii]-a[kk]*a[jj];
                    ii=l*n+l; jj=l*n+l+1;
                    kk=(l+1)*n+l; ll=(l+1)*n+l+1;
                    p=a[ii]*(a[ii]-x)+a[jj]*a[kk]+y;
                    q=a[kk]*(a[ii]+a[ll]-x);
                    r=a[kk]*a[(l+2)*n+l+1];
                  }
                if ((fabs(p)+fabs(q)+fabs(r))!=0.0)
                  { xy=1.0;
                    if (p<0.0) xy=-1.0;
                    s=xy*sqrt(p*p+q*q+r*r);
                    if (k!=l) a[k*n+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-1; j++)
                      { ii=k*n+j; jj=(k+1)*n+j;
                        p=x*a[ii]+e*a[jj];
                        q=e*a[ii]+y*a[jj];
                        r=f*a[ii]+g*a[jj];
                        if (k!=m-2)
                          { kk=(k+2)*n+j;
                            p=p+f*a[kk];
                            q=q+g*a[kk];
                            r=r+z*a[kk]; a[kk]=r;
                          }
                        a[jj]=q; a[ii]=p;
                      }
                    j=k+3;
                    if (j>=m-1) j=m-1;
                    for (i=l; i<=j; i++)
                      { ii=i*n+k; jj=i*n+k+1;
                        p=x*a[ii]+e*a[jj];
                        q=e*a[ii]+y*a[jj];
                        r=f*a[ii]+g*a[jj];
                        if (k!=m-2)
                          { kk=i*n+k+2;
                            p=p+f*a[kk];
                            q=q+g*a[kk];
                            r=r+z*a[kk]; a[kk]=r;
                          }
                        a[jj]=q; a[ii]=p;
                      }
                  }
              }
          }
      }
    return(1);
  }
void main()
  { int i,j,jt=60;
    double eps=0.000001;
    static double u[4],v[4];
    static double a[4][4]={ {1.0,3.0,2.0,13.0},
                             {7.0,2.0,1.0,-2.0},
                             {9.0,15.0,3.0,-2.0},
                             {-2.0,-2.0,11.0,5.0}};
    hhbg(4,*a);
    printf("\n");
	cout<<"矩阵的全部特征值："<<endl<<endl;
    i=qrtt(4,*a,u,v,eps,jt);
    if (i>0)
      for (i=0; i<=3; i++)
	  {
		  if(v[i]>=0)
		  {
              printf("%13.7f +j %13.7f\n",u[i],v[i]);
		  }
		  else
		  {
			  printf("%13.7f -j %13.7f\n",u[i],abs(v[i]));
		  }
	  }
    printf("\n");
	cout<<"矩阵行列式值："<<endl<<endl;
	cout<<"det(A)=   ";
	double p=1.0,q=0.0,t;
	for(i=0;i<=3;i++)
	{
		t=p;                   //临时变量保存p的数值
		p=u[i]*p-v[i]*q;
		q=v[i]*t+u[i]*q;
	}
	cout<<p<<endl<<endl;
  }