ludecomposition.txt

基于LU分解的直接求解方法

#define TINY 1.0e-20;

void LU(complex **a, complex *b, int n)

{ int i,imax,j,k;
  complex big,dum,sum,temp;
  complex *vv=new complex[n];
  int *indx=new int[n];

  for (i=0;i<n;i++)
  {
    big=0.0;
    for (j=0;j<n;j++)
    { temp=a[i][j];
      if (abs(temp) >abs( big)) big=temp;
    }

    if (abs(big) == 0.0)      { cout<<"Singular matrix in routine ludcmp"; return; }
    vv[i]=1.0/big;               // Save the scaling.
 }

  for (j=0;j<n;j++)
  {
    for (i=0;i<j;i++)
    {
      sum=a[i][j];
      for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
      a[i][j]=sum;
    }

    big=0.0;      //Initialize for the search for largest pivot element.
    for (i=j;i<n;i++)
    {           //This is i = j of equation (2.3.12) and i = j+1 : ::N of equation (2.3.13).
      sum=a[i][j];
      for (k=0;k<j;k++)  sum -= a[i][k]*a[k][j];
      a[i][j]=sum;
      dum=vv[i]*abs(sum);
      if (abs (dum) >=abs( big))
      { big=dum;
        imax=i;
      }
    }
    if (j != imax)
    {
      for (k=0;k<n;k++)
      {
        dum=a[imax][k];
        a[imax][k]=a[j][k];
        a[j][k]=dum;
      }
     vv[imax]=vv[j];
    }
    indx[j]=imax;
    if (abs(a[j][j])<1.e-8) a[j][j]=TINY;
    if (j != (n-1))
    {
      dum=1.0/(a[j][j]);
      for (i=j+1;i<n;i++) a[i][j] *= dum;
    }
  }
  delete[] vv;
  int ip;
  int ii=-1;

  for (i=0;i<n;i++)
  { ip=indx[i];
    sum=b[ip];
    b[ip]=b[i];

    if(ii>=0) for(j=ii;j<=i-1;j++) sum-=a[i][j]*b[j];
    else if(abs(sum)) ii=i;

    b[i]=sum;
  }
  for (i=n-1;i>=0;i--)
  { sum=b[i];
    for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
    b[i]=sum/a[i][i];
  }
}