inv.h

电力系统计算中的矩阵求逆程序

void MatInv(double *pa, int n)
    {
      int *is,*js,i,j,k,l,u,v;
      double d,p;
      is=new int[n];
      js=new int[n];

      for (k=0; k<=n-1; k++)
	{ 
	 d=0.0;
	 for (i=k; i<=n-1; i++)
	    for (j=k; j<=n-1; j++)
	       { l=i*n+j; p=fabs(*(pa+l));
	         if (p>d) { d=p; is[k]=i; js[k]=j;}
	       }
	    if (d+1.0==1.0)
	       { delete is; delete js; printf("err**not inv\n");}
	    if (is[k]!=k)
		for (j=0; j<=n-1; j++)
		   { u=k*n+j; v=is[k]*n+j;
		      p=*(pa+u);*(pa+u)=*(pa+v);*(pa+v)=p;
		   }
	    if (js[k]!=k)
	        for (i=0; i<=n-1; i++)
		    { u=i*n+k; v=i*n+js[k];
		       p=*(pa+u);*(pa+u)=*(pa+v);*(pa+v)=p;
		    }
	    l=k*n+k;
	    *(pa+l)=1.0/(*(pa+l));
            for (j=0; j<=n-1; j++)
		if (j!=k)
		    { u=k*n+j; *(pa+u)=*(pa+u)**(pa+l);}
	    for (i=0; i<=n-1; i++)
	        if (i!=k)
		   for (j=0; j<=n-1; j++)
			if (j!=k)
			   { u=i*n+j;
			     *(pa+u)=*(pa+u)-*(pa+i*n+k)**(pa+k*n+j);
			   }
	    for (i=0; i<=n-1; i++)
		if (i!=k)
		   { u=i*n+k; *(pa+u)=-*(pa+u)**(pa+l);}
	 }
	for (k=n-1; k>=0; k--)
	   { if (js[k]!=k)
	      for (j=0; j<=n-1; j++)
	          { u=k*n+j; v=js[k]*n+j;
	            p=*(pa+u);*(pa+u)=*(pa+v);*(pa+v)=p;
	          }
	     if (is[k]!=k)
	       for (i=0; i<=n-1; i++)
		  { u=i*n+k; v=i*n+is[k];
	            p=*(pa+u);*(pa+u)=*(pa+v);*(pa+v)=p;
		  }
	   }
	delete is; delete js;
}