C, Numerical Recipes in C, solution of linear algebraic equations, LU decomposition source program

 

《Numerical Recipes in C++》 Summary of the original ：

 

 

 

Make up the number of words in the nonsense translation ：

In linear algebra , LU decompose (LU Factorization) It's a kind of matrix factorization , A matrix can be decomposed into the product of a lower triangular matrix and an upper triangular matrix （ Sometimes it's the product of them and a permutation matrix ）.LU Decomposition is mainly used in numerical analysis , To solve linear equations 、 Find the inverse matrix or calculate the determinant .
The coefficient matrix A Into two equivalent matrices L and U The product of the , among L and U They are unit lower triangular matrix and upper triangular matrix . When A All the sequential principal and sub expressions of are not 0 when , matrix A Can be uniquely decomposed into A=LU. among L It's a lower triangular matrix ,U It's an upper triangular matrix .
LU Decomposition is essentially an expression of Gauss elimination method . In essence, I will A It is transformed into an upper triangular matrix by elementary row transformation , Its transformation matrix is a unit lower triangular matrix . This is the so-called durritt algorithm （Doolittle algorithm）： From bottom to top, the matrix A Do elementary line transformation , Change the element at the bottom left of the diagonal to zero , Then it is proved that the effect of these row transformations is equivalent to left multiplying a series of unit lower triangular matrices , The inverse of the product of this series of unit lower triangular matrices is L matrix , It is also a unit lower triangular matrix . The complexity of this kind of algorithm is generally ( Two thirds of n Third power ) about .

Algorithm improvements

（i）Doolittle decompose
For nonsingular matrices （ ren n The order principal and subordinate formulas are not all 0） Matrix of A, All can be done Doolittle decompose , obtain A=LU, among L Is the unit lower triangular matrix ,U Is the upper triangular matrix ; there Doolittle Decomposition is actually Gauss Transformation ;
（ii）Crout decompose
For nonsingular matrices （ ren n The order principal and subordinate formulas are not all 0） Matrix of A, All can be done Crout decompose , obtain A=LU, among L Is the lower triangular matrix ,U Is the unit upper triangular matrix ;
（iii） Column principal component triangular decomposition
For the square matrix of nonsingular matrix A, Use column principal element triangular decomposition , obtain PA=LU, among P Is a permutation matrix ,L,U And Doolittle The decomposition rules are the same ;
（iv） Triangular decomposition of all principal components
For the square matrix of nonsingular matrix A, Use full principal component triangular decomposition , obtain PAQ=LU, among P,Q Is a permutation matrix ,L,U And Doolittle The decomposition rules are the same ;
（v） Direct triangular decomposition
For the square matrix of nonsingular matrix A, The formula derived from direct triangular decomposition （Doolittle Decomposition formula or Crout Decomposition formula ）, You can perform recursive operations , To facilitate computer programming ;
（vi）“ Catch up method ”
The catch-up method is for banded matrices （ Especially tridiagonal matrix ） The special structure of this large sparse matrix , Derivation of a band preserving decomposition formula obtained , The substantive result is also LU decompose ; Because large sparse matrix is widely used in Engineering , So this part of the content needs to be mastered .
（vii）Cholesky Decomposition （ Square root method ） And the improved square root method
Cholesky The decomposition method is aimed at the decomposition of positive definite matrix , As a result, A=LDLT=LD(1/2)D(1/2)LT=L1L1T. How to get L1, In fact, the recursive formula is also given .
The improved square root method is Cholesky An improvement of decomposition . To avoid square in the formula , And what you get is A=LDLT=TLT, It also gives the solution T,L Formula .
Summary ：
（1） from （i）~（iv） It is the basic method of manual calculation ,（v）~（vi） It is guided by the algorithm formula of computer-aided calculation ;
（2） The purpose of these methods is to get the solution of linear equations , The essence is Gauss elimination .
 

 

Source program （POWER BY 315SOFT.COM）

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// LU decomposition - PTC
    /// Iterative Improvement of a Solution to Linear Equations
    /// </summary>
    public class LUdcmp
    {
        private int n { get; set; }
        /// <summary>
        /// Stores the decomposition
        /// </summary>
        private double[,] lu { get; set; }
        /// <summary>
        /// Stores the permutation.
        /// </summary>
        private int[] indx { get; set; }
        /// <summary>
        /// Used by det.
        /// </summary>
        private double d { get; set; }
        private double[,] aref { get; set; }

        /// <summary>
        /// Given a matrix a[0..n - 1][0..n-1], this routine replaces it by the LU decomposition of a
        /// rowwise permutation of itself.a is input.On output, it is arranged as in equation(2.3.14)
        /// above; indx[0..n - 1] is an output vector that records the row permutation effected by the
        /// partial pivoting; d is output as +/-1 depending on whether the number of row interchanges
        /// was even or odd, respectively. This routine is used in combination with solve to solve linear
        /// equations or invert a matrix.
        /// </summary>
        /// <param name="a"></param>
        /// <exception cref="Exception"></exception>
        public LUdcmp(double[,] a)
        {
            this.n = a.GetLength(0);
            this.lu = a;
            this.aref = a;
            this.indx = new int[n];
            const double TINY = 1.0e-40;

            double[] vv = new double[n];
            d = 1.0;
            for (int i = 0; i < n; i++)
            {
                double big = 0.0;
                for (int j = 0; j < n; j++)
                {
                    double temp = Math.Abs(lu[i, j]);
                    if ((temp) > big)
                    {
                        big = temp;
                    }
                }
                //if (big == 0.0)
                if (Math.Abs(big) <= float.Epsilon)
                {
                    throw new Exception("Singular matrix in LUdcmp");
                }
                vv[i] = 1.0 / big;
            }
            for (int k = 0; k < n; k++)
            {
                double big = 0.0;
                int imax = k;
                for (int i = k; i < n; i++)
                {
                    double temp = vv[i] * Math.Abs(lu[i, k]);
                    if (temp > big)
                    {
                        big = temp;
                        imax = i;
                    }
                }
                if (k != imax)
                {
                    for (int j = 0; j < n; j++)
                    {
                        double temp = lu[imax, j];
                        lu[imax, j] = lu[k, j];
                        lu[k, j] = temp;
                    }
                    d = -d;
                    vv[imax] = vv[k];
                }
                indx[k] = imax;
                //if (lu[k, k] == 0.0)
                if (Math.Abs(lu[k, k]) <= float.Epsilon)
                {
                    lu[k, k] = TINY;
                }
                for (int i = k + 1; i < n; i++)
                {
                    double temp = lu[i, k] /= lu[k, k];
                    for (int j = k + 1; j < n; j++)
                    {
                        lu[i, j] -= temp * lu[k, j];
                    }
                }
            }
        }

        /// <summary>
        /// Solves the set of n linear equations A*x = b using the stored LU decomposition of A.
        /// b[0..n - 1] is input as the right-hand side vector b, while x returns the solution vector x; b and
        /// x may reference the same vector, in which case the solution overwrites the input.This routine
        /// takes into account the possibility that b will begin with many zero elements, so it is efficient for
        /// use in matrix inversion.
        /// </summary>
        /// <param name="b"></param>
        /// <param name="x"></param>
        /// <exception cref="Exception"></exception>
        public void solve(double[] b, double[] x)
        {
            int ii = 0;
            if (b.Length != n || x.Length != n)
            {
                throw new Exception("LUdcmp::solve bad sizes");
            }
            for (int i = 0; i < n; i++)
            {
                x[i] = b[i];
            }
            for (int i = 0; i < n; i++)
            {
                int ip = indx[i];
                double sum = x[ip];
                x[ip] = x[i];
                if (ii != 0)
                {
                    for (int j = ii - 1; j < i; j++)
                    {
                        sum -= lu[i, j] * x[j];
                    }
                }
                else if (sum != 0.0)
                {
                    ii = i + 1;
                }
                x[i] = sum;
            }
            for (int i = n - 1; i >= 0; i--)
            {
                double sum = x[i];
                for (int j = i + 1; j < n; j++)
                {
                    sum -= lu[i, j] * x[j];
                }
                x[i] = sum / lu[i, i];
            }
        }

        /// <summary>
        /// Solves m sets of n linear equations A*X = B using the stored LU decomposition of A.The
        /// matrix b[0..n - 1][0..m - 1] inputs the right-hand sides, while x[0..n - 1][0..m - 1] returns the
        /// solution A^-1* B.b and x may reference the same matrix, in which case the solution overwrites
        /// the input.
        /// </summary>
        /// <param name="b"></param>
        /// <param name="x"></param>
        /// <exception cref="Exception"></exception>
        public void solve(double[,] b, double[,] x)
        {
            int m = b.GetLength(1);
            if (b.GetLength(0) != n || x.GetLength(0) != n || b.GetLength(1) != x.GetLength(1))
            {
                throw new Exception("LUdcmp::solve bad sizes");
            }
            double[] xx = new double[n];
            for (int j = 0; j < m; j++)
            {
                for (int i = 0; i < n; i++)
                {
                    xx[i] = b[i, j];
                }
                solve( xx,  xx);
                for (int i = 0; i < n; i++)
                {
                    x[i, j] = xx[i];
                }
            }
        }

        /// <summary>
        /// Using the stored LU decomposition, 
        /// return in ainv the matrix inverse
        /// </summary>
        /// <param name="ainv"></param>
        public void inverse(double[,] ainv)
        {
            //ainv.resize(n, n);
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    ainv[i, j] = 0.0;
                }
                ainv[i, i] = 1.0;
            }
            solve( ainv,  ainv);
        }

        /// <summary>
        /// Using the stored LU decomposition, 
        /// return the determinant of the matrix A.
        /// </summary>
        /// <returns></returns>
        public double det()
        {
            double dd = d;
            for (int i = 0; i < n; i++)
            {
                dd *= lu[i, i];
            }
            return dd;
        }

        /// <summary>
        /// Improves a solution vector x[0..n - 1] of the linear set of equations A*x= b.
        /// The vectors b[0..n - 1] and x[0..n - 1] are input.On output, x[0..n - 1] is
        /// modified, to an improved set of values.
        /// </summary>
        /// <param name="b"></param>
        /// <param name="x"></param>
        public void mprove(double[] b, double[] x)
        {
            double[] r = new double[n];
            for (int i = 0; i < n; i++)
            {
                double sdp = -b[i];
                for (int j = 0; j < n; j++)
                {
                    sdp += (double)aref[i, j] * (double)x[j];
                }
                r[i] = sdp;
            }
            solve( r,  r);
            for (int i = 0; i < n; i++)
            {
                x[i] -= r[i];
            }
        }
    }
}