UT/UT_MatrixSolver.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * Produced by:
 *      David Wong
 *      Side Effects
 *      477 Richmond Street West
 *      Toronto, Ontario
 *      Canada   M5V 3E7
 *      416-504-9876
 *
 * NAME:        Matrix Solver (C++)
 *
 * COMMENTS:
 *
 */

#ifndef __UT_MatrixSolver_H__
#define __UT_MatrixSolver_H__

#include "UT_API.h"
#include "UT_Assert.h"

#include "UT_Functor.h"
#include "UT_Matrix.h"
#include "UT_Vector.h"

// Set below to 1 in order to use old version of LUDecomp
#define USE_OLD_LUDECOMP 0

class UT_API UT_MatrixSolver {
public:
    /// LU Decomposition of a[1..n][1..n] where a has full rank. 
    /// Output: index[1..n] is the row permutation 
    ///         d = +/-1 indicating row interchanges was even or odd, resp.
    ///         return 0 if the matrix is truly singular 
    ///               1 on success
    ///               2 if the matrix is singular to the precision of the
    ///                 algorithm, where tol is the zero tolerance.
    int         LUDecompOld(UT_Matrix &a, UT_Permutation &index, float &d, 
                            float tol = 1e-20F);

    /// Solve Ax=b
    /// a[1..n][1..n] is a LU Decomposition of matrix A (eg. from LUDecomp())
    /// index describes the row permutations as output from LUDecomp()
    /// b[1..n] is input as rhs, and on output contains the solution x
    void        LUBackSubOld(const UT_Matrix &a, const UT_Permutation &index, 
                             UT_Vector &b);

#if USE_OLD_LUDECOMP
    int         LUDecomp(UT_Matrix &a, UT_Permutation &index, float &d, 
                         float tol = 1e-20F)
    {
        return LUDecompOld(a, index, d, tol);
    }
    void        LUBackSub(const UT_Matrix &a, const UT_Permutation &index, 
                          UT_Vector &b)
    {
        return LUBackSubOld(a, index, b);
    }
#else
    int         LUDecomp(UT_Matrix &a, UT_Permutation &index, float &d, 
                         float tol = 1e-20F)
    {
        return LUDecompRect(a, index, d, tol);
    }
    void        LUBackSub(const UT_Matrix &a, const UT_Permutation &index, 
                          UT_Vector &b)
    {
        return LUDecompBackSub(a, index, b);
    }
#endif

    /// LU decomposition of rectangular A[1..m][1..n] with partial pivoting.
    /// Unlike LUDecomp(), this handles pivot values of 0.
    /// ie. it performs the decomposition of A = P.L.U for unit lower
    ///     triangular L and upper triangular U.
    /// Input:  A[1..m][1..n]
    /// Output: A contains the matrices L and U in compact form, without
    ///           the unit diagonal.
    ///           L is of size min(m,n) by n. U is of size m by min(m,n)
    ///         pivots[1..n] is the row permulation P
    ///         tol is scaled by the largest pivot element
    ///         detsign is the sign of the determinant used by det()
    ///            its +1 if row interchanges was even, -1 if odd
    ///         Returns 0 if the matrix is truly singular 
    ///                 1 on success
    ///                 2 if the matrix is singular to the precision of the
    ///                   algorithm, where tol is the zero tolerance.
    /// This algorithm takes roughly O(m*n^2 - n^3/3) flops
    int         LUDecompRect(UT_Matrix &A, UT_Permutation &pivots,
                             float &detsign,
                             float tol = 1e-6F);

    /// Performs back substitution using the output of LUDecomp() for a square
    /// matrix.
    /// Input:  A[1..n][1..n] and pivots as output from LUDecompRect()
    ///         b[1..n] is the right hand side for solving A.x = b
    /// Output: b is modified to contain the solution x
    void        LUDecompBackSub(const UT_Matrix &A,
                                const UT_Permutation &pivots,
                                UT_Vector &b);

    /// LU decompose a companion matrix
    /// The matrix given is the base of the companion matrix:
    /// [aI bI 0  0 ]
    /// [0  aI bI 0 ]
    /// [0  0  aI bI]
    /// [P1 P2 P3 P4]
    /// ie: [P1 P2 P3 P4]
    /// The resulting LU decompostion is written into the R matrix and
    /// is:
    /// [aI 0   0 0] [I b/aI  0    0  ]
    /// [0  aI  0 0] [0  I   b/aI  0  ]
    /// [0  0  aI 0] [0  0    I   b/aI]
    /// [R1 R2 R3 L] [0  0    0    U  ]
    /// LU is  stored in the last portion of the matrix & index is the
    /// permuation that resulted in it.  Thus, index should be the height
    /// of P, NOT the entire companion matrix.
    ///
    /// Return 0 on success, -1 for singular, and -2 for singular to
    /// machine precision (As judged by tol), and -3 for a being near zero.
    int         LUDecompCompanion(const UT_Matrix &P, float a, float b,
                                    UT_Matrix &R, UT_Permutation &index,
                                    float tol=1e-20F);

    /// This takes the result of the previous operation and does a 
    /// back substitution.
    void        LUBackSubCompanion(const UT_Matrix &R, float a, float b, 
                                   const UT_Permutation &index,
                                   UT_Vector &y);

    ///
    /// Cholesky factorization
    ///
    /// To solve a full rank system A . x = b where A can be undertermined
    /// (ie. A has less rows than columns), we can make use of the Cholesky
    /// factorization functions below by way of the normal equations like so:
    /// 1. Left multiply both sides by the tranpose of A, A^T to get
    ///         (A^T . A) . x = (A^T . b)
    /// 2. Now solve for x using P . x = c where P = A^T . A and c = A^T . b
    ///    This is done by first calling choleskyDecomp() using P since it
    ///    will be symmetric non-negative definite. Then use choleskyBackSub()
    ///    to obtain x for particular values of c.

    /// Performs a Cholesky decomposition of the symmetric non-negative definite
    /// matrix A such that A = L . L^T where L^T is tranpose of the lower
    /// triangular matrix L.
    /// Input:  Symmetric non-negative definite matrix a[1..n][1..n].
    ///         Only the upper triangle of a is read from.
    ///         tol is the tolerance check for when a sum is deemed 0
    /// Output: Cholesky factor L is returned in the lower triangle of A, except
    ///         for its diagonal elements which is returned in d[1..n]
    ///         Returns the number of articial zero's placed into d. If we
    ///         encounter negative diagonal values, they are zeroed out and we
    ///         return -N where N is the number of negative values found.
    int         choleskyDecomp(UT_Matrix &a, UT_Vector &d, float tol=1e-5F);

    /// Solves the set of n linear equations, A.x = b, where a is a nonnegative
    /// definite symmetric matrix.
    /// Input:  a[1..n][1..n] and d[1..n] as output from choleskyDecomp() above.
    ///         Only the lower triangle of a is read from.
    ///        b[1..n] is the right-hand side of A.x = b to solve for
    ///         tol is the tolerance check for when an element of d is deemed 0
    /// Output: x[1..n], the solution of A.x = b
    ///         Returns the number of articial zero's placed into x.
    int         choleskyBackSub(const UT_Matrix &a, const UT_Vector &d,
                                const UT_Vector &b, UT_Vector &x,
                                float tol=1e-5F);
    
    /// Performs inverse iteration to find the closest eigenvalue(s)
    /// to s.
    /// Does a two pass method, first using single iterations until
    /// starttol is reached, in which case singles are used to
    /// achieve final tolerance.  If starttol is not reached in
    /// startiter iterations, uses double method to find using
    /// final tolerance/iterations.
    ///
    /// The vector q is updated with the newest eigenvector.
    /// Source matrices are:
    ///        [1  0  0]        [0 1 0]
    ///  C1' = [0  1  0]  C2' = [0 0 1]
    ///        [0  0 C1]        [ C2  ]
    /// Thus, the initial sizes of the matrices/vectors are (row*col):
    /// C1: n*n, C2: n*(n*deg), R: n*(n*deg), tmp: n*(n*deg),
    /// index: n, q: n*deg, Aq: n*deg
    /// The result code determines how s1 & s2 are to be interpreted:
    /// Result:
    ///  0 - One eigenvalue found at s1, eigenvector q.
    ///  1 - Two real eigenvalues found: s1 & s2.
    ///  2 - Complex conjugate eigenvalues found: s1 +- s2i.
    /// -1 - Failed to converge to final tolerance using single method.
    /// -2 - Failed to converge to final tolerance using double method.
    int         inversePowerIterate(const UT_Matrix &C1, 
                                const UT_Matrix &C2, 
                                UT_Matrix &R, UT_Matrix &tmp, 
                                UT_Permutation &index, 
                                float s, 
                                float &s1, float &s2,
                                UT_Vector &q, UT_Vector &Aq,
                                float starttol = 1e-2F, int startiter = 5,
                                float finaltol = 1e-5F, 
                                float doubletol = 1e-4F, int finaliter = 50);

    /// Compute the inverse of a LU decomposed matrix.
    /// If you want to solve systems of equations AX=B,
    /// where B is a matrix. It is better (faster) to do a 
    /// LU decomposition and use back substituation for columns of B.
    void        inverse(UT_Matrix &a, UT_Permutation &index, UT_Matrix &ia);

    /// Finds all linearly independent rows in A, does NOT require
    /// A to be any particular shape.  Returns the number of independent
    /// rows & first entries of index are their indices in the original
    /// matrix.
    /// Destroys A.
    int         findLinIndRows(UT_Matrix &A, UT_Permutation &index, 
                                    float tol = 1E-6F);

    /// Similar to findLinIndRows, performs gaussian elimination
    /// with full pivoting on G.  Updates as well Gu, and Gv, the
    /// partial derivitive matrix of G(u, v).  After this is done,
    /// the determinant of G & partials of the determinant can
    /// be computed.
    void        fullGaussianElimination(UT_Matrix &G, UT_Matrix &Gu,
                                        UT_Matrix &Gv,
                                        UT_Permutation &rowperm,
                                        UT_Permutation &colperm,
                                        float &detsign);

    /// After full Gaussian Elimination, this finds the determinant
    /// and partials thereof:
    void        detWithPartials(const UT_Matrix &G, const UT_Matrix &Gu,
                                const UT_Matrix &Gv, 
                                const UT_Permutation &rowperm,
                                const UT_Permutation &colperm,
                                float detsign,
                                float &retdet, float &detu, float &detv);


    /// Compute the determinant of a LU decomposed matrix.
    /// The d is what one gets from LU decomposition.
    float       det(UT_Matrix &a, float d);

    /// Find the determiment:
    /// (A(p,q), B(r,s), C(t,u)))
    inline float        det3x3(const UT_Matrix &A, const UT_Matrix &B, 
                    const UT_Matrix &C, int p, int q, int r, int s,
                    int t, int u)
            {
                return (A(p,q) * (B(r,s) * C(t,u) - B(t,u) * C(r,s)) -
                        B(p,q) * (A(r,s) * C(t,u) - A(t,u) * C(r,s)) +
                        C(p,q) * (A(r,s) * B(t,u) - A(t,u) * B(r,s)) );
            }

    /// Estimates condtion number of upper triangular matrix:
    void        condEstimate(const UT_Matrix &A, UT_Vector &y);

    /// Finds the condition number of a LUDecomposed matrix.
    /// This is just an approximation
    /// You must provide the infinite norm of A before LUDecomposition.
    float       getConditionLUD(const UT_Matrix &A, 
                                const UT_Permutation &index, float anorm);

    /// Find householder vector of x.
    void        house(const UT_Vector &x, UT_Vector &v, float &b);
    void        house(const UT_Matrix &x, UT_Vector &v, float &b);

    /// Finds the givens rotation required to 0 b:
    void        findGivens(float a, float b, float &c, float &s,
                            float zerotol = 1e-20F);

    /// Find half givens, return 0 if success, otherwise the required
    /// rotation is impossible (ie:complex eigenvalue)
    int         findHalfGivens(const UT_Matrix &A, float &c, float &s,
                                float tol = 0.0F);

    /// Upper triangulizes A with householder transforms, applying
    /// simultaneously to A and optionally Q.
    void        triangularize(UT_Matrix &A, UT_Matrix &B, UT_Matrix *Q = 0);

    /// Transforms A into upper hessenburg and B into upper triangular.
    void        hessentri(UT_Matrix &A, UT_Matrix &B);

    /// Converts A into quasi upper triangular & B into upper triangular
    /// form using the QZ algorithm.
    int         realSchurFormQZ(UT_Matrix &oA, UT_Matrix &oB, float tol=1e-6F,
                                int maxIteration = 30);

    /// Transform A into a hessenberg matrix through a series of 
    /// householder reflections.
    void        hessenberg(UT_Matrix &A);
    /// Computes P st A' = P^T A P
    void        hessenberg(UT_Matrix &A, UT_Matrix &P);

    /// Transforms A into a bidiagonal matrix via householder transformations.
    int         bidiagonalize(UT_Matrix &A);

    /// Reduces the matrix to real Schur form:
    /// Orthogonal transform responible for getting to real schur form
    /// is accumulated in Q, if it is non-zero.
    int         realSchurForm(UT_Matrix &A, UT_Matrix *Q = 0, 
                                float tol = 1E-6F,
                                int maxIteration = 30);

    /// Golub-Kahan SVD algorithm.  Resulting SVD matrix found in
    /// diaganol of A.
    int         gkSVD(UT_Matrix &A, float tol = 1e-6F,
                                int maxIteration = 30);

    /// Use Numerical Recipes method for SVD decomposition, with results
    /// that are suitable for doing back substitution.
    /// NOTE: The UT_Matrix should start with 1,1!
    ///
    /// Decomposition a = u * w * v^T
    ///   a[1..m][1..n] also stores the result u[1..m][1..n]
    ///   w[1..n][1..n] is a diagonal matrix, stored as a vector
    ///   v[1..n][1..n]
    int         SVDDecomp(UT_Matrix &a, UT_Vector &w, UT_Matrix &v,
                                int maxIterstion = 30);
    /// Perform back substitution on an SVD decomposed matrix. Solves
    /// equation Ax = b.
    int         SVDBackSub(UT_Matrix &u, UT_Vector &w, UT_Matrix &v,
                                UT_Vector &b, UT_Vector &x, float tol = 1e-6F);

    /// Finds eigenvalues from a real schur form matrix:
    void        findEigenValues(const UT_Matrix &A, 
                                UT_Vector &vreal, UT_Vector &vimg);
    /// Finds eigenvalues from generalized system (A - lambdaB)
    void        findEigenValues(const UT_Matrix &A, const UT_Matrix &B,
                                UT_Vector &sreal, UT_Vector &simg,
                                UT_Vector &treal, UT_Vector &timg);

    /// Finds the u & v eigenvalues from an eigenvector, whose udeg,
    /// vdeg, and tdeg are as specified.  Returns 0 if fails.
    /// NB: Does the s/(1+s) transformation!
    int         findUVfromEigenvector(const UT_Vector &vect,
                        int udeg, int vdeg, int tdeg, float &u, float &v);

    /// Finds a specified eigenvector from a real schur form matrix.
    /// The eigenvector must be real!
    void        findRightEigenVector(const UT_Matrix &A, const UT_Matrix &Q,
                                UT_Vector &vector, int i);
                                

private:
    /// Pushes down zeros, causing A to be a reduced hessenberg.
    void        chaseZerosQZ(UT_Matrix &A, UT_Matrix &B, int zeropos);

    /// Pushes across zeros generated in Golub-Kahan SVD method
    void        chaseZerosGK(UT_Matrix &B, int zeropos);

    /// Applies one iteration of the QZ algorithm to unreduced
    /// upper hessenberg A & non-singular upper triangular B.
    void        QZStep(UT_Matrix &A, UT_Matrix &B);

    /// Performs the francisQRStep on H.
    void        francisQRStep(UT_Matrix &H);

    // Performs francisQRStep, updating H12 = H12*Z,
    // Q = Q*Z, H23 = ZT * H23 simultaneously with finding H = ZT H Z
    void        francisQRStep(UT_Matrix &H, 
                            UT_Matrix &H12, UT_Matrix &H23, UT_Matrix &Q);

    /// Performs one step in the Golub-Kahan SVD algorithm:
    /// B is a bidiagonal matrix with no zeros on diagonal nor supradiagonal,
    void        golubKahanSVDStep(UT_Matrix &B);

    /// Solves (A - lambda I) x = gamma B for x, requires A is quasi triangular.
    /// B is column.
    void        solveQuasiUpperTriangular(const UT_Matrix &A, float lambda,
                                    float gamma, 
                                    const UT_Matrix &B, UT_Vector &x,
                                    float tol = 1e-6F);

    /// Performs inverse iteration until convergence or maxiter is reached.
    /// The vector is updated with the newest eigenvector.
    /// Call with findDecomp = 0 to use R directly as opposed to decomposing
    /// it.  Source matrices are:
    ///        [1  0  0]        [0 1 0]
    ///  C1' = [0  1  0]  C2' = [0 0 1]
    ///        [0  0 C1]        [ C2  ]
    /// Returns:
    ///  1: Decomp found the pencil to be degenerate, perfect solution.
    ///  0: Successfully converged
    /// <0: Failed to converge
    int         inversePowerSingle( const UT_Matrix &C1, 
                                    const UT_Matrix &C2, 
                                    UT_Matrix &R, UT_Matrix &tmp, 
                                    UT_Permutation &index, 
                                    float &s, float lasts,
                                    UT_Vector &q, UT_Vector &Aq,
                                    int findDecomp = 1, 
                                    float tol = 1e-5F, int maxiter = 50);

    /// Performs inverse iteration until convergence or maxiter is reached.
    /// The vector is updated with the newest eigenvector.
    /// Call with findDecomp = 0 to use R directly as opposed to decomposing
    /// it.  Source matrices are:
    ///        [1  0  0]        [0 1 0]
    ///  C1' = [0  1  0]  C2' = [0 0 1]
    ///        [0  0 C1]        [ C2  ]
    /// Computes two solutions, namely P & Q such that
    /// A = (C2 - sC1), and
    /// lim i->inf (A^i+2 - pA^i+1 - qA^i) u = 0
    /// Returns:
    ///  1: Decomp found the pencil to be degenerate, perfect solution.
    ///  0: Successfully converged
    /// <0: Failed to converge
    int         inversePowerDouble( const UT_Matrix &C1, 
                                    const UT_Matrix &C2, 
                                    UT_Matrix &R, UT_Matrix &tmp, 
                                    UT_Permutation &index, 
                                    float s, 
                                    float &p, float &q,
                                    UT_Vector &u, UT_Vector &Au,
                                    int findDecomp = 1, 
                                    float tol = 1e-5F, int maxiter = 50);
};


// UT_MatrixIterSolver class is mainly for sparse matrices.
// The sparse matrices may be defined explicitly with UT_SparseMatrix,
// or implicitly with functions.

class UT_API UT_MatrixIterSolver {
public:
    /// Preconditioned Conjugate Gradient method for solving Ax=b,
    /// given the symmetric positive definite matrix A[0..n-1][0..n-1] 
    /// indirectly through AMult (ie. Ax),
    /// a preconditioner for A through ASolve, 
    /// the vector b[0..n-1] and an initial guess x[0..n-1].
    /// maxIter defaults to 2n.
    /// error |b-Ax|/|b| < tol.
    /// normType indicates the type of norm used for error:
    ///     0    L-infinity norm (ie. max)
    ///     1    L1-norm         (ie. sum of abs)
    ///     2    L2-norm         (ie. Euclidean distance)
    /// Output: x[0..n-1] and error.
    float       PCG(void (*AMult)(const UT_VectorD &x, UT_VectorD &result),
                    void (*ASolve)(const UT_VectorD &b, UT_VectorD &x),
                    int n, UT_VectorD &x, const UT_VectorD &b,
                    float tol=1e-3F, int normType=2, int maxIter=-1);
    float       PCG(void (*AMult)(const UT_VectorF &x, UT_VectorF &result),
                    void (*ASolve)(const UT_VectorF &b, UT_VectorF &x),
                    int n, UT_VectorF &x, const UT_VectorF &b,
                    float tol=1e-3F, int normType=2, int maxIter=-1);

    /// Preconditioned Conjugate Gradient method for solving Ax=b with
    /// the symmetric positive definite matrix A, preconditioner, and iteration
    /// conditions implicitly defined by functors.
    static void PCG(UT_Vector &x, const UT_Vector &b,
                const UT_Functor2<void, const UT_Vector &, UT_Vector &> &AMult,
                const UT_Functor2<void, const UT_Vector &, UT_Vector &> &ASolve,
                const UT_Functor2<bool, int, const UT_Vector &> &iterateTest);

    /// Preconditioned Conjugate Gradient method for solving least squares Ax=b,
    /// given the matrix A[0..m-1][0..n-1] indirectly through AMult 
    /// (ie. Ax if int = 0 or A^tx if int = 1),
    /// a preconditioner for A^tA through AtASolve, 
    /// the vector b[0..m-1] and an initial guess x[0..n-1].
    /// maxIter defaults to 2n.
    /// error |A^t(b-Ax)|/|A^tb| < tol.
    /// normType indicates the type of norm used for error:
    ///     0    L-infinity norm (ie. max)
    ///     1    L1-norm         (ie. sum of abs)
    ///     2    L2-norm         (ie. Euclidean distance)
    /// Output: x[0..n-1] and error.
    float       PCGLS(void (*AMult)(const UT_Vector &x, UT_Vector &result, 
                                    int transpose),
                      void (*AtASolve)(const UT_Vector &b, UT_Vector &x),
                      int m, int n,
                      UT_Vector &x, const UT_Vector &b,
                      float tol=1e-3F, int normType=2, int maxIter=-1);

    /// Preconditioned Conjugate Gradient method for solving Ax=b,
    /// given the matrix A[0..m-1][0..n-1] indirectly through AMult 
    /// (ie. Ax if int = 0 or A^tx if int = 1),
    /// a preconditioner for A through ASolve, 
    /// the vector b[0..m-1] and an initial guess x[0..n-1].
    /// maxIter defaults to 2n.
    /// error |b-Ax|/|b| < tol.
    /// normType indicates the type of norm used for error:
    ///     0    L-infinity norm (ie. max)
    ///     1    L1-norm         (ie. sum of abs)
    ///     2    L2-norm         (ie. Euclidean distance)
    /// Output: x[0..n-1] and error.
    float       biPCG(void (*AMult)(const UT_Vector &x, UT_Vector &result, 
                                    int transpose),
                      void (*ASolve)(const UT_Vector &b, UT_Vector &x, 
                                    int transpose),
                      int n,
                      UT_Vector &x, const UT_Vector &b,
                      float tol=1e-3F, int normType=2, int maxIter=-1);
};


#endif

---