www.pudn.com > cghost.rar > AGJDN.C

 
  #include "stdlib.h" 
  #include "math.h" 
  #include "stdio.h" 
  int agjdn(a,b,n,m) 
  int n,m; 
  double a[],b[]; 
  { int *js,l,k,i,j,is,p,q; 
    double d,t; 
    js=malloc(n*sizeof(int)); 
    l=1; 
    for (k=0;k<=n-1;k++) 
      { d=0.0; 
        for (i=k;i<=n-1;i++) 
          for (j=k;j<=n-1;j++) 
            { t=fabs(a[i*n+j]); 
              if (t>d) { d=t; js[k]=j; is=i;} 
            } 
        if (d+1.0==1.0) l=0; 
        else 
          { if (js[k]!=k) 
              for (i=0;i<=n-1;i++) 
                { p=i*n+k; q=i*n+js[k]; 
                  t=a[p]; a[p]=a[q]; a[q]=t; 
                } 
            if (is!=k) 
              { for (j=k;j<=n-1;j++) 
                  { p=k*n+j; q=is*n+j; 
                    t=a[p]; a[p]=a[q]; a[q]=t; 
                  } 
                for (j=0;j<=m-1;j++) 
                  { p=k*m+j; q=is*m+j; 
                    t=b[p]; b[p]=b[q]; b[q]=t; 
                  } 
              } 
          } 
        if (l==0) 
          { free(js); printf("fail\n"); 
            return(0); 
          } 
        d=a[k*n+k]; 
        for (j=k+1;j<=n-1;j++) 
          { p=k*n+j; a[p]=a[p]/d;} 
        for (j=0;j<=m-1;j++) 
          { p=k*m+j; b[p]=b[p]/d;} 
        for (j=k+1;j<=n-1;j++) 
          for (i=0;i<=n-1;i++) 
            { p=i*n+j; 
              if (i!=k) 
                a[p]=a[p]-a[i*n+k]*a[k*n+j]; 
            } 
        for (j=0;j<=m-1;j++) 
        for (i=0;i<=n-1;i++) 
          { p=i*m+j; 
            if (i!=k) 
              b[p]=b[p]-a[i*n+k]*b[k*m+j]; 
          } 
      } 
    for (k=n-1;k>=0;k--) 
      if (js[k]!=k) 
        for (j=0;j<=m-1;j++) 
          { p=k*m+j; q=js[k]*m+j; 
            t=b[p]; b[p]=b[q]; b[q]=t; 
          } 
    free(js); 
    return(1); 
  }