Compact LU Factorization / Decomposition schemes of Crout, Doolittle and Cholesky

Tested C++ code for the compact LU factorization / decomposition schemes of Crout, Doolittle and Cholesky

LU Factorization or Decomposition is an efficient and common method for directly solving linear systems like Ax = b. The principles of LU decomposition are described in most numerical methods texts. The code on this page implements C / C++ versions of compact LU decomposition schemes that are useful for practical solutions of linear systems where high speed and low storage requirements are important.

These compact schemes are meant to perform all operations in-place so that the matrix A is overwritten by the desired LU decomposition. This is desirable for some problems, but not others, so the code on this page (below) uses separate source S and destination D matrices even though this is not required. By separating S and D the source matrix may be saved for future use without having to copy it manually. But if in-place computation is preferred then a naive re-naming of all D matrices to S (or vice versa) will work just fine.

The d by d matrices S and D are stored as single 1D arrays. C-style array indices apply here so the first element of an array is S[0] and the last element is S[N-1] where N is the total number of elements in the array.

Each compact scheme must be used with its own forward substitution function. This is required because the non-unit diagonal elements could belong to either the upper or lower triangular matrix and the forward substitution function must know which is which.

The best description that I've seen of these algorithms is: Dahlquist, Bjorck and Anderson "Numerical Methods" (1974). I've been dissatisfied with the treatment of LU decomposition by several popular undergraduate numerical methods texts because they:

often present pseudo-code with arrays going from 1 to N
often present pseudo-code that is simply incorrect
rarely mention that different forward substitution functions must be used for different schemes
rarely discuss the various LU decomposition options such as Crout, Doolittle and Cholesky

A single-file code that implements and tests these methods can be downloaded here. In gcc just compile with something like "g++ -O3 -ffastmath LU.cpp" and then run and results for decompositions will print out.

Some helpful discussion is online at: http://math.fullerton.edu/mathews/n2003/CholeskyMod.html. The results printed by the single-file code can be compared to the correct answers on this site.

I'd be happy to hear any comments or suggestions you might have; here is how to contact me.

// Crout uses unit diagonals for the upper triangle
void Crout(int d,double*S,double*D){
   for(int k=0;k<d;++k){
      for(int i=k;i<d;++i){
         double sum=0.;
         for(int p=0;p<k;++p)sum+=D[i*d+p]*D[p*d+k];
         D[i*d+k]=S[i*d+k]-sum; // not dividing by diagonals
      }
      for(int j=k+1;j<d;++j){
         double sum=0.;
         for(int p=0;p<k;++p)sum+=D[k*d+p]*D[p*d+j];
         D[k*d+j]=(S[k*d+j]-sum)/D[k*d+k];
      }
   }
}
void solveCrout(int d,double*LU,double*b,double*x){
   double y[d];
   for(int i=0;i<d;++i){
      double sum=0.;
      for(int k=0;k<i;++k)sum+=LU[i*d+k]*y[k];
      y[i]=(b[i]-sum)/LU[i*d+i];
   }
   for(int i=d-1;i>=0;--i){
      double sum=0.;
      for(int k=i+1;k<d;++k)sum+=LU[i*d+k]*x[k];
      x[i]=(y[i]-sum); // not dividing by diagonals
   }
}



// Doolittle uses unit diagonals for the lower triangle
void Doolittle(int d,double*S,double*D){
   for(int k=0;k<d;++k){
      for(int j=k;j<d;++j){
         double sum=0.;
         for(int p=0;p<k;++p)sum+=D[k*d+p]*D[p*d+j];
         D[k*d+j]=(S[k*d+j]-sum); // not dividing by diagonals
      }
      for(int i=k+1;i<d;++i){
         double sum=0.;
         for(int p=0;p<k;++p)sum+=D[i*d+p]*D[p*d+k];
         D[i*d+k]=(S[i*d+k]-sum)/D[k*d+k];
      }
   }
}
void solveDoolittle(int d,double*LU,double*b,double*x){
   double y[d];
   for(int i=0;i<d;++i){
      double sum=0.;
      for(int k=0;k<i;++k)sum+=LU[i*d+k]*y[k];
      y[i]=(b[i]-sum); // not dividing by diagonals
   }
   for(int i=d-1;i>=0;--i){
      double sum=0.;
      for(int k=i+1;k<d;++k)sum+=LU[i*d+k]*x[k];
      x[i]=(y[i]-sum)/LU[i*d+i];
   }
}



// Cholesky requires the matrix to be symmetric positive-definite
void Cholesky(int d,double*S,double*D){
   for(int k=0;k<d;++k){
      double sum=0.;
      for(int p=0;p<k;++p)sum+=D[k*d+p]*D[k*d+p];
      D[k*d+k]=sqrt(S[k*d+k]-sum);
      for(int i=k+1;i<d;++i){
         double sum=0.;
         for(int p=0;p<k;++p)sum+=D[i*d+p]*D[k*d+p];
         D[i*d+k]=(S[i*d+k]-sum)/D[k*d+k];
      }
   }
}
// This version could be more efficient on some architectures
// Use solveCholesky for both Cholesky decompositions
void CholeskyRow(int d,double*S,double*D){
   for(int k=0;k<d;++k){
      for(int j=0;j<d;++j){
         double sum=0.;
         for(int p=0;p<j;++p)sum+=D[k*d+p]*D[j*d+p];
         D[k*d+j]=(S[k*d+j]-sum)/D[j*d+j];
      }
      double sum=0.;
      for(int p=0;p<k;++p)sum+=D[k*d+p]*D[k*d+p];
      D[k*d+k]=sqrt(S[k*d+k]-sum);
   }
}
void solveCholesky(int d,double*LU,double*b,double*x){
   double y[d];
   for(int i=0;i<d;++i){
      double sum=0.;
      for(int k=0;k<i;++k)sum+=LU[i*d+k]*y[k];
      y[i]=(b[i]-sum)/LU[i*d+i];
   }
   for(int i=d-1;i>=0;--i){
      double sum=0.;
      for(int k=i+1;k<d;++k)sum+=LU[k*d+i]*x[k];
      x[i]=(y[i]-sum)/LU[i*d+i];
   }
}

Home