LUDecomposition.hpp

Go to the documentation of this file.

/**--------------------------------------------------------------------------------------------
Copyright (c) 2019, NDEVR LLC
tyler.parke@ndevr.org
  __    __   ____     _____ __     __  _______
 |  \  |  | |  __ \  |  ___|\ \   / / |   __  \
 |   \ |  | | |  \ \ | |___  \ \ / /  |  |__)  |
 |  . \|  | | |__/ / | |___   \ V /   |   _   /
 |  |\    |_|_____/__|_____|___\_/____|  | \  \
 |__| \__________________________________|  \__\
 
Subject to the terms of the Enterprise+ Agreement, NDEVR hereby grants 
Licensee a limited, non-exclusive, non-transferable, royalty-free license 
(without the right to sublicense) to use the API solely for the purpose of 
Licensee's internal development efforts to develop applications for which 
the API was provided.
 
The above copyright notice and this permission notice shall be included in all 
copies or substantial portions of the Software.
 
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE
FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
 
Library: Base
File: LUDecomposition
Included in API: True
Author(s): Tyler Parke
 *-----------------------------------------------------------------------------------------**/
#pragma once
#include <NDEVR/Matrix.h>
namespace NDEVR
{
    template<class t_type, uint01 t_m, uint01 t_n>
    class LUDecomposition
    {
    public:
        LUDecomposition(const Matrix<t_type, t_m, t_n>& A)
            : LU(A)
        {
 
            // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
 
            for (uint01 i = 0; i < t_m; i++)
            {
                piv[i] = i;
            }
            pivsign = 1;
            Vector<t_m, t_type> LUrowi(0.0);
            Vector<t_m, t_type> LUcolj(0.0);
 
            // Outer loop.
 
            for (uint01 j = 0; j < t_n; j++) 
            {
 
                // Make a copy of the j-th column to localize references.
                for (uint01 i = 0; i < t_m; i++) 
                    LUcolj[i] = LU[i][j];
 
                // Apply previous transformations.
 
                for (uint01 i = 0; i < t_m; i++) 
                {
                    LUrowi = LU[i];
 
                    // Most of the time is spent in the following dot product.
 
                    uint01 kmax = getMin(i, j);
                    t_type s(0);
                    for (uint01 k = 0; k < kmax; k++)
                        s += LUrowi[k] * LUcolj[k];
 
                    LUrowi[j] = LUcolj[i] -= s;
                }
 
                // Find pivot and exchange if necessary.
 
                uint01 p = j;
                for (uint01 i = j + 1; i < t_m; i++) 
                {
                    if (abs(LUcolj[i]) > abs(LUcolj[p]))
                        p = i;
                }
                if (p != j) 
                {
                    for (uint01 k = 0; k < t_n; k++) 
                    {
                        std::swap(LU[p][k], LU[j][k]);
                    }
                    std::swap(piv[p], piv[j]);
                    pivsign = -pivsign;
                }
 
                // Compute multipliers.
 
                if (j < t_m && LU[j][j] != 0.0) 
                {
                    for (uint01 i = j + 1; i < t_m; i++) {
                        LU[i][j] /= LU[j][j];
                    }
                }
            }
        }
 
         /* ------------------------
            Public Methods
          * ------------------------ */
 
          /** Is the matrix nonsingular?
          @return     true if U, and hence A, is nonsingular.
          */
 
        bool isNonsingular() const 
        {
            for (uint01 j = 0; j < t_n; j++) 
            {
                if (LU[j][j] == 0)
                    return false;
            }
            return true;
        }
 
        /** Return lower triangular factor
        @return     L
        */
 
        Matrix<t_type, t_m, t_n> getL() const
        {
            Matrix<t_type, t_m, t_n> mat(0.0);
            for (uint01 i = 0; i < t_m; i++) 
            {
                for (uint01 j = 0; j < t_n; j++) 
                {
                    if (i > j)
                        mat[i][j] = LU[i][j];
                    else if (i == j)
                        mat[i][j] = t_type(1);
                    else
                        mat[i][j] = t_type(0);
                }
            }
            return mat;
        }
 
        /** Return upper triangular factor
        @return     U
        */
 
        Matrix<t_type, t_m, t_n> getU() const
        {
            Matrix<t_type, t_n, t_n> mat(0.0);
            for (uint01 i = 0; i < t_n; i++)
            {
                for (uint01 j = 0; j < t_n; j++)
                {
                    if (i <= j)
                        mat[i][j] = LU[i][j];
                    else
                        mat[i][j] = t_type(0);
                }
            }
            return mat;
        }
 
        /** Return pivot permutation vector
        @return     piv
        */
        template<class t_pivot_type = sint04>
        Vector<t_m, t_pivot_type> getPivot() const
        {
            return piv.template as<t_m, t_pivot_type>();
        }
 
        /** Determinant
        @return     det(A)
        @exception  IllegalArgumentException  Matrix must be square
        */
 
        t_type det() const
        {
            static_assert(t_m == t_n, "Matrix must be square.");
            t_type d = cast<t_type>(pivsign);
            for (uint01 j = 0; j < t_n; j++)
                d *= LU[j][j];
            return d;
        }
 
        /** Solve A*X = B
        @param  B   A Matrix with as many rows as A and any number of columns.
        @return     X so that L*U*X = B(piv,:)
        @exception  IllegalArgumentException Matrix row dimensions must agree.
        @exception  RuntimeException  Matrix is singular.
        */
        template<uint01 t_nx>
        Matrix<t_type, t_m, t_nx> solve(const Matrix<t_type, t_m, t_nx>& B)
        {
            lib_assert(isNonsingular(), "Matrix cannot be singular.");
 
            // Copy right hand side with pivoting
            Matrix<t_type, t_m, t_nx> mat = B.template subMatrix<0, t_nx - 1>(piv);
 
            // Solve L*Y = B(piv,:)
            for (uint01 k = 0; k < t_n; k++) 
            {
                for (uint01 i = k + 1; i < t_n; i++) 
                {
                    for (uint01 j = 0; j < t_nx; j++) 
                    {
                        mat[i][j] -= mat[k][j] * LU[i][k];
                    }
                }
            }
            // Solve U*X = Y;
            for (uint01 k = t_n - 1; !isNaN(k); k--)
            {
                for (uint01 j = 0; j < t_nx; j++)
                    mat[k][j] /= LU[k][k];
                for (uint01 i = 0; i < k; i++) 
                {
                    for (uint01 j = 0; j < t_nx; j++)
                        mat[i][j] -= mat[k][j] * LU[i][k];
                }
            }
            return mat;
        }
 
        private:
            Matrix<t_type, t_m, t_n> LU;
 
            /** Row and column dimensions, and pivot sign.
            @serial column dimension.
            @serial row dimension.
            @serial pivot sign.
            */
            sint04 pivsign;
 
            /** Internal storage of pivot vector.
            @serial pivot vector.
            */
            Vector<t_m, uint01> piv;
    };
}