HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | List of all members
UT_MatrixSolverT< T > Class Template Reference

#include <UT_MatrixSolver.h>

Public Member Functions

int LUDecompOld (UT_MatrixT< T > &a, UT_Permutation &index, T &d, T tol=1e-20)
 
template<typename S >
void LUBackSubOld (const UT_MatrixT< T > &a, const UT_Permutation &index, UT_VectorT< S > &b)
 
int LUDecomp (UT_MatrixT< T > &a, UT_Permutation &index, T &d, T tol=1e-20)
 
template<typename S >
void LUBackSub (const UT_MatrixT< T > &a, const UT_Permutation &index, UT_VectorT< S > &b)
 
int LUDecompRect (UT_MatrixT< T > &A, UT_Permutation &pivots, T &detsign, T tol=1e-6)
 
template<typename S >
void LUDecompBackSub (const UT_MatrixT< T > &A, const UT_Permutation &pivots, UT_VectorT< S > &b)
 
int LUDecompCompanion (const UT_MatrixT< T > &P, T a, T b, UT_MatrixT< T > &R, UT_Permutation &index, T tol=1e-20)
 
template<typename S >
void LUBackSubCompanion (const UT_MatrixT< T > &R, T a, T b, const UT_Permutation &index, UT_VectorT< S > &y)
 
template<typename S >
int choleskyDecomp (UT_MatrixT< T > &a, UT_VectorT< S > &d, T tol=1e-5)
 
template<typename S >
int choleskyBackSub (const UT_MatrixT< T > &a, const UT_VectorT< S > &d, const UT_VectorT< S > &b, UT_VectorT< S > &x, T tol=1e-5)
 
template<typename S >
int inversePowerIterate (const UT_MatrixT< T > &C1, const UT_MatrixT< T > &C2, UT_MatrixT< T > &R, UT_MatrixT< T > &tmp, UT_Permutation &index, T s, T &s1, T &s2, UT_VectorT< S > &q, UT_VectorT< S > &Aq, T starttol=1e-2, int startiter=5, T finaltol=1e-5, T doubletol=1e-4, int finaliter=50)
 
void inverse (UT_MatrixT< T > &a, UT_Permutation &index, UT_MatrixT< T > &ia)
 
int findLinIndRows (UT_MatrixT< T > &A, UT_Permutation &index, T tol=1E-6)
 
void fullGaussianElimination (UT_MatrixT< T > &G, UT_MatrixT< T > &Gu, UT_MatrixT< T > &Gv, UT_Permutation &rowperm, UT_Permutation &colperm, T &detsign)
 
void detWithPartials (const UT_MatrixT< T > &G, const UT_MatrixT< T > &Gu, const UT_MatrixT< T > &Gv, const UT_Permutation &rowperm, const UT_Permutation &colperm, T detsign, T &retdet, T &detu, T &detv)
 
T det (UT_MatrixT< T > &a, T d)
 
T det3x3 (const UT_MatrixT< T > &A, const UT_MatrixT< T > &B, const UT_MatrixT< T > &C, int p, int q, int r, int s, int t, int u)
 
template<typename S >
void condEstimate (const UT_MatrixT< T > &A, UT_VectorT< S > &y)
 Estimates condtion number of upper triangular matrix: More...
 
T getConditionLUD (const UT_MatrixT< T > &A, const UT_Permutation &index, T anorm)
 
template<typename S >
void house (const UT_VectorT< S > &x, UT_VectorT< S > &v, T &b)
 Find householder vector of x. More...
 
template<typename S >
void house (const UT_MatrixT< T > &x, UT_VectorT< S > &v, T &b)
 
void findGivens (T a, T b, T &c, T &s, T zerotol=1e-20)
 Finds the givens rotation required to 0 b: More...
 
int findHalfGivens (const UT_MatrixT< T > &A, T &c, T &s, T tol=0.0)
 
void triangularize (UT_MatrixT< T > &A, UT_MatrixT< T > &B, UT_MatrixT< T > *Q=0)
 
void hessentri (UT_MatrixT< T > &A, UT_MatrixT< T > &B)
 Transforms A into upper hessenburg and B into upper triangular. More...
 
int realSchurFormQZ (UT_MatrixT< T > &oA, UT_MatrixT< T > &oB, T tol=1e-6, int maxIteration=30)
 
void hessenberg (UT_MatrixT< T > &A)
 
void hessenberg (UT_MatrixT< T > &A, UT_MatrixT< T > &P)
 Computes P st A' = P^T A P. More...
 
int bidiagonalize (UT_MatrixT< T > &A)
 Transforms A into a bidiagonal matrix via householder transformations. More...
 
int realSchurForm (UT_MatrixT< T > &A, UT_MatrixT< T > *Q=0, T tol=1E-6, int maxIteration=30)
 
int gkSVD (UT_MatrixT< T > &A, T tol=1e-6, int maxIteration=30)
 
template<typename S >
int SVDDecomp (UT_MatrixT< T > &a, UT_VectorT< S > &w, UT_MatrixT< T > &v, int maxIterstion=30)
 
template<typename S >
int SVDBackSub (UT_MatrixT< T > &u, UT_VectorT< S > &w, UT_MatrixT< T > &v, UT_VectorT< S > &b, UT_VectorT< S > &x, T tol=1e-6)
 
template<typename S >
void findEigenValues (const UT_MatrixT< T > &A, UT_VectorT< S > &vreal, UT_VectorT< S > &vimg)
 Finds eigenvalues from a real schur form matrix: More...
 
template<typename S >
void findEigenValues (const UT_MatrixT< T > &A, const UT_MatrixT< T > &B, UT_VectorT< S > &sreal, UT_VectorT< S > &simg, UT_VectorT< S > &treal, UT_VectorT< S > &timg)
 Finds eigenvalues from generalized system (A - lambdaB) More...
 
template<typename S >
int findUVfromEigenvector (const UT_VectorT< S > &vect, int udeg, int vdeg, int tdeg, T &u, T &v)
 
template<typename S >
void findRightEigenVector (const UT_MatrixT< T > &A, const UT_MatrixT< T > &Q, UT_VectorT< S > &vector, int i)
 

Detailed Description

template<typename T>
class UT_MatrixSolverT< T >

Definition at line 28 of file UT_MatrixSolver.h.

Member Function Documentation

template<typename T >
int UT_MatrixSolverT< T >::bidiagonalize ( UT_MatrixT< T > &  A)

Transforms A into a bidiagonal matrix via householder transformations.

template<typename T >
template<typename S >
int UT_MatrixSolverT< T >::choleskyBackSub ( const UT_MatrixT< T > &  a,
const UT_VectorT< S > &  d,
const UT_VectorT< S > &  b,
UT_VectorT< S > &  x,
T  tol = 1e-5 
)

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.

template<typename T >
template<typename S >
int UT_MatrixSolverT< T >::choleskyDecomp ( UT_MatrixT< T > &  a,
UT_VectorT< S > &  d,
T  tol = 1e-5 
)

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.
template<typename T >
template<typename S >
void UT_MatrixSolverT< T >::condEstimate ( const UT_MatrixT< T > &  A,
UT_VectorT< S > &  y 
)

Estimates condtion number of upper triangular matrix:

template<typename T >
T UT_MatrixSolverT< T >::det ( UT_MatrixT< T > &  a,
T  d 
)

Compute the determinant of a LU decomposed matrix. The d is what one gets from LU decomposition.

template<typename T >
T UT_MatrixSolverT< T >::det3x3 ( const UT_MatrixT< T > &  A,
const UT_MatrixT< T > &  B,
const UT_MatrixT< T > &  C,
int  p,
int  q,
int  r,
int  s,
int  t,
int  u 
)
inline

Find the determiment: (A(p,q), B(r,s), C(t,u)))

Definition at line 255 of file UT_MatrixSolver.h.

template<typename T >
void UT_MatrixSolverT< T >::detWithPartials ( const UT_MatrixT< T > &  G,
const UT_MatrixT< T > &  Gu,
const UT_MatrixT< T > &  Gv,
const UT_Permutation &  rowperm,
const UT_Permutation &  colperm,
T  detsign,
T &  retdet,
T &  detu,
T &  detv 
)

After full Gaussian Elimination, this finds the determinant and partials thereof:

template<typename T >
template<typename S >
void UT_MatrixSolverT< T >::findEigenValues ( const UT_MatrixT< T > &  A,
UT_VectorT< S > &  vreal,
UT_VectorT< S > &  vimg 
)

Finds eigenvalues from a real schur form matrix:

template<typename T >
template<typename S >
void UT_MatrixSolverT< T >::findEigenValues ( const UT_MatrixT< T > &  A,
const UT_MatrixT< T > &  B,
UT_VectorT< S > &  sreal,
UT_VectorT< S > &  simg,
UT_VectorT< S > &  treal,
UT_VectorT< S > &  timg 
)

Finds eigenvalues from generalized system (A - lambdaB)

template<typename T >
void UT_MatrixSolverT< T >::findGivens ( T  a,
T  b,
T &  c,
T &  s,
T  zerotol = 1e-20 
)

Finds the givens rotation required to 0 b:

template<typename T >
int UT_MatrixSolverT< T >::findHalfGivens ( const UT_MatrixT< T > &  A,
T &  c,
T &  s,
T  tol = 0.0 
)

Find half givens, return 0 if success, otherwise the required rotation is impossible (ie:complex eigenvalue)

template<typename T >
int UT_MatrixSolverT< T >::findLinIndRows ( UT_MatrixT< T > &  A,
UT_Permutation &  index,
T  tol = 1E-6 
)

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.

template<typename T >
template<typename S >
void UT_MatrixSolverT< T >::findRightEigenVector ( const UT_MatrixT< T > &  A,
const UT_MatrixT< T > &  Q,
UT_VectorT< S > &  vector,
int  i 
)

Finds a specified eigenvector from a real schur form matrix. The eigenvector must be real!

template<typename T >
template<typename S >
int UT_MatrixSolverT< T >::findUVfromEigenvector ( const UT_VectorT< S > &  vect,
int  udeg,
int  vdeg,
int  tdeg,
T &  u,
T &  v 
)

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!

template<typename T >
void UT_MatrixSolverT< T >::fullGaussianElimination ( UT_MatrixT< T > &  G,
UT_MatrixT< T > &  Gu,
UT_MatrixT< T > &  Gv,
UT_Permutation &  rowperm,
UT_Permutation &  colperm,
T &  detsign 
)

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.

template<typename T >
T UT_MatrixSolverT< T >::getConditionLUD ( const UT_MatrixT< T > &  A,
const UT_Permutation &  index,
T  anorm 
)

Finds the condition number of a LUDecomposed matrix. This is just an approximation You must provide the infinite norm of A before LUDecomposition.

template<typename T >
int UT_MatrixSolverT< T >::gkSVD ( UT_MatrixT< T > &  A,
T  tol = 1e-6,
int  maxIteration = 30 
)

Golub-Kahan SVD algorithm. Resulting SVD matrix found in diaganol of A.

template<typename T >
void UT_MatrixSolverT< T >::hessenberg ( UT_MatrixT< T > &  A)

Transform A into a hessenberg matrix through a series of householder reflections.

template<typename T >
void UT_MatrixSolverT< T >::hessenberg ( UT_MatrixT< T > &  A,
UT_MatrixT< T > &  P 
)

Computes P st A' = P^T A P.

template<typename T >
void UT_MatrixSolverT< T >::hessentri ( UT_MatrixT< T > &  A,
UT_MatrixT< T > &  B 
)

Transforms A into upper hessenburg and B into upper triangular.

template<typename T >
template<typename S >
void UT_MatrixSolverT< T >::house ( const UT_VectorT< S > &  x,
UT_VectorT< S > &  v,
T &  b 
)

Find householder vector of x.

template<typename T >
template<typename S >
void UT_MatrixSolverT< T >::house ( const UT_MatrixT< T > &  x,
UT_VectorT< S > &  v,
T &  b 
)
template<typename T >
void UT_MatrixSolverT< T >::inverse ( UT_MatrixT< T > &  a,
UT_Permutation &  index,
UT_MatrixT< T > &  ia 
)

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.

template<typename T >
template<typename S >
int UT_MatrixSolverT< T >::inversePowerIterate ( const UT_MatrixT< T > &  C1,
const UT_MatrixT< T > &  C2,
UT_MatrixT< T > &  R,
UT_MatrixT< T > &  tmp,
UT_Permutation &  index,
T  s,
T &  s1,
T &  s2,
UT_VectorT< S > &  q,
UT_VectorT< S > &  Aq,
T  starttol = 1e-2,
int  startiter = 5,
T  finaltol = 1e-5,
T  doubletol = 1e-4,
int  finaliter = 50 
)

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.

template<typename T >
template<typename S >
void UT_MatrixSolverT< T >::LUBackSub ( const UT_MatrixT< T > &  a,
const UT_Permutation &  index,
UT_VectorT< S > &  b 
)
inline

Definition at line 70 of file UT_MatrixSolver.h.

template<typename T >
template<typename S >
void UT_MatrixSolverT< T >::LUBackSubCompanion ( const UT_MatrixT< T > &  R,
T  a,
T  b,
const UT_Permutation &  index,
UT_VectorT< S > &  y 
)

This takes the result of the previous operation and does a back substitution.

template<typename T >
template<typename S >
void UT_MatrixSolverT< T >::LUBackSubOld ( const UT_MatrixT< T > &  a,
const UT_Permutation &  index,
UT_VectorT< S > &  b 
)

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

template<typename T >
int UT_MatrixSolverT< T >::LUDecomp ( UT_MatrixT< T > &  a,
UT_Permutation &  index,
T &  d,
T  tol = 1e-20 
)
inline

Definition at line 64 of file UT_MatrixSolver.h.

template<typename T >
template<typename S >
void UT_MatrixSolverT< T >::LUDecompBackSub ( const UT_MatrixT< T > &  A,
const UT_Permutation &  pivots,
UT_VectorT< S > &  b 
)

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

template<typename T >
int UT_MatrixSolverT< T >::LUDecompCompanion ( const UT_MatrixT< T > &  P,
T  a,
T  b,
UT_MatrixT< T > &  R,
UT_Permutation &  index,
T  tol = 1e-20 
)

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.

template<typename T >
int UT_MatrixSolverT< T >::LUDecompOld ( UT_MatrixT< T > &  a,
UT_Permutation &  index,
T &  d,
T  tol = 1e-20 
)

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.

template<typename T >
int UT_MatrixSolverT< T >::LUDecompRect ( UT_MatrixT< T > &  A,
UT_Permutation &  pivots,
T &  detsign,
T  tol = 1e-6 
)

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

template<typename T >
int UT_MatrixSolverT< T >::realSchurForm ( UT_MatrixT< T > &  A,
UT_MatrixT< T > *  Q = 0,
T  tol = 1E-6,
int  maxIteration = 30 
)

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.

template<typename T >
int UT_MatrixSolverT< T >::realSchurFormQZ ( UT_MatrixT< T > &  oA,
UT_MatrixT< T > &  oB,
T  tol = 1e-6,
int  maxIteration = 30 
)

Converts A into quasi upper triangular & B into upper triangular form using the QZ algorithm.

template<typename T >
template<typename S >
int UT_MatrixSolverT< T >::SVDBackSub ( UT_MatrixT< T > &  u,
UT_VectorT< S > &  w,
UT_MatrixT< T > &  v,
UT_VectorT< S > &  b,
UT_VectorT< S > &  x,
T  tol = 1e-6 
)

Perform back substitution on an SVD decomposed matrix. Solves equation Ax = b.

template<typename T >
template<typename S >
int UT_MatrixSolverT< T >::SVDDecomp ( UT_MatrixT< T > &  a,
UT_VectorT< S > &  w,
UT_MatrixT< T > &  v,
int  maxIterstion = 30 
)

Use Numerical Recipes method for SVD decomposition, with results that are suitable for doing back substitution. NOTE: The UT_MatrixF 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]

template<typename T >
void UT_MatrixSolverT< T >::triangularize ( UT_MatrixT< T > &  A,
UT_MatrixT< T > &  B,
UT_MatrixT< T > *  Q = 0 
)

Upper triangulizes A with householder transforms, applying simultaneously to A and optionally Q.

The documentation for this class was generated from the following file: