编译无误的全选主元高斯消去法

算法来源于<<数值分析与算法>>   徐士良 著
但,我发现他的算法在传递二维数组时,怎么也不能把数据传入.可能是我的水平有限,没有搞明白为什么.但问题总得解决,所以我做了一点点小小的改动.能通过编译,我用matlalab难计算结果.结果正确,
我用的编译器是devc++.现贴出原码,希望他人少走弯路.

///正确的高斯全选项主元消去法

#include <cstdlib>

#include <iostream>

#include <math.h>

#include <fstream>

using namespace std;

int mygauss(double *allm,int n,double * b)

{

    ofstream resultfile("result.out");

    int * js,l,k,i,j,is,p,q;

    double d,t;

    js=(int *)malloc(n*sizeof(int));

    l=1;//置非奇异标志

    //cout<<"hui";

    for(k=0;k<=n-2;k++)

    {

        d=0.0;

        for(i=k;i<=n-1;i++)

            for(j=k;j<=n-1;j++)

                 {

                 

                    t=fabs(allm[i*n+j]);

                    if(t>d){d=t;js[k]=j;is=i;};//记忆行列交换信息,,,js[k]保存列      

                  }  

               //cout<<d<<'#';

              //cout<<js[k]<<'@';  

              // cout<<"i"<<is;

         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=allm[p];allm[p]=allm[q];allm[q]=t;

                        }   

                   }

                   

                  if(is!=k)//行交换 

                    {

                       for(j=k;j<=n-1;j++)

                          { 

                              p=k*n+j;q=is*n+j;

                              t=allm[p];allm[p]=allm[q];allm[q]=t;

                          }   

                        t=b[k];b[k]=b[is];b[is]=t;   

                    }           

            } 

            

            

            if(l==0)//奇异返回

             {

                free(js);printf("fail/n"); return 0;              

             } 

             

          d=allm[k*n+k];

          

          for(j=k+1;j<=n-1;j++)//归一化

             {p=k*n+j;allm[p]=allm[p]/d;}

             b[k]=b[k]/d;

          for(i=k+1;i<=n-1;i++)//消元 

             {

                for(j=k+1;j<=n-1;j++)

                   {

                      p=i*n+j;

                      allm[p]=allm[p]-allm[i*n+k]*allm[k*n+j];            

                   }  

                b[i]=b[i] -allm[i*n+k]*b[k];                      

             }  

       

  

            //for(int w=0;w<n;w++)

               //{

               //for(int y=0;y<n;y++)

                //{

                //     cout<<allm[w*n+y]<<'#'; 

                 // }

                 

                 // cout<<endl;

               //}         

    }

    

    

    

      d=allm[(n-1)*n+n-1];  

             

       

       if(fabs(d)+1.0==1.0)//奇异返回

         {

              free(js);printf("fail/n");

              return 0;                       

         }  

         

         b[n-1]=b[n-1]/d; 

         for(i=n-2;i>=0;i--)

         {

              t=0.0;

              for(j=i+1;j<=n-1;j++)

                 t=t+allm[i*n+j]*b[j];

                 b[i]=b[i]-t;

         }

         

         js[n-1]=n-1;

         

         for(k=n-1;k>=0;k--)

            {

               if(js[k]!=k)

                  {

                       t=b[k];b[k]=b[js[k]];b[js[k]]=t;    

                  }             

            }

            for(int pos=0;pos<n;pos++)

            {

                //cout<<b[pos]<<"flag";  

                resultfile<<b[pos];

                resultfile<<'/n';  

            }

            free(js);

            resultfile.close();

    return 1;

    

    

}

int main(int argc, char *argv[])

{

    double * dall;

    double a[16]={0.2368,0.2471,0.2568,1.2671,  

                  0.1968,0.2071,1.2168,0.2271,  

                  0.1581,1.1675,0.1768,0.1871,  

                  1.1161,0.1254,0.1397,0.1490};

    dall=a; 

    

    double * djg;

    double jga[4]={1.8471,1.7471,1.6471,1.5471};

    

    djg=jga;

    mygauss(dall,4,djg);

             

                        

                  

    system("PAUSE");

    return EXIT_SUCCESS;

}

//result.out

1.04058

0.987051

0.93504

0.881282