HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
UT_MatrixSolver.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: Matrix Solver (C++)
7  *
8  * COMMENTS:
9  *
10  */
11 
12 #ifndef __UT_MatrixSolver_H__
13 #define __UT_MatrixSolver_H__
14 
15 #include "UT_API.h"
16 #include "UT_Assert.h"
17 #include "UT_Functor.h"
18 #include "UT_Matrix.h"
19 #include "UT_Vector.h"
20 #include "UT_SparseMatrix.h"
21 #include "UT_Interrupt.h"
22 #include <SYS/SYS_Types.h>
23 
24 // Set below to 1 in order to use old version of LUDecomp
25 #define USE_OLD_LUDECOMP 0
26 
27 template <typename T>
29 {
30 public:
31  /// LU Decomposition of a[1..n][1..n] where a has full rank.
32  /// Output: index[1..n] is the row permutation
33  /// d = +/-1 indicating row interchanges was even or odd, resp.
34  /// return 0 if the matrix is truly singular
35  /// 1 on success
36  /// 2 if the matrix is singular to the precision of the
37  /// algorithm, where tol is the zero tolerance.
38  int LUDecompOld(UT_MatrixT<T> &a, UT_Permutation &index, T &d,
39  T tol = 1e-20);
40 
41  /// Solve Ax=b
42  /// a[1..n][1..n] is a LU Decomposition of matrix A (eg. from LUDecomp())
43  /// index describes the row permutations as output from LUDecomp()
44  /// b[1..n] is input as rhs, and on output contains the solution x
45  template <typename S>
46  void LUBackSubOld(const UT_MatrixT<T> &a,
47  const UT_Permutation &index,
48  UT_VectorT<S> &b);
49 
50 #if USE_OLD_LUDECOMP
51  int LUDecomp(UT_MatrixT<T> &a, UT_Permutation &index, T &d,
52  T tol = 1e-20)
53  {
54  return LUDecompOld(a, index, d, tol);
55  }
56  template <typename S>
57  void LUBackSub(const UT_MatrixT<T> &a,
58  const UT_Permutation &index,
59  UT_VectorT<S> &b)
60  {
61  return LUBackSubOld(a, index, b);
62  }
63 #else
64  int LUDecomp(UT_MatrixT<T> &a, UT_Permutation &index, T &d,
65  T tol = 1e-20)
66  {
67  return LUDecompRect(a, index, d, tol);
68  }
69  template <typename S>
70  void LUBackSub(const UT_MatrixT<T> &a, const UT_Permutation &index,
71  UT_VectorT<S> &b)
72  {
73  return LUDecompBackSub(a, index, b);
74  }
75 #endif
76 
77  /// LU decomposition of rectangular A[1..m][1..n] with partial pivoting.
78  /// Unlike LUDecomp(), this handles pivot values of 0.
79  /// ie. it performs the decomposition of A = P.L.U for unit lower
80  /// triangular L and upper triangular U.
81  /// Input: A[1..m][1..n]
82  /// Output: A contains the matrices L and U in compact form, without
83  /// the unit diagonal.
84  /// L is of size min(m,n) by n. U is of size m by min(m,n)
85  /// pivots[1..n] is the row permulation P
86  /// tol is scaled by the largest pivot element
87  /// detsign is the sign of the determinant used by det()
88  /// its +1 if row interchanges was even, -1 if odd
89  /// Returns 0 if the matrix is truly singular
90  /// 1 on success
91  /// 2 if the matrix is singular to the precision of the
92  /// algorithm, where tol is the zero tolerance.
93  /// This algorithm takes roughly O(m*n^2 - n^3/3) flops
94  int LUDecompRect(UT_MatrixT<T> &A, UT_Permutation &pivots,
95  T &detsign,
96  T tol = 1e-6);
97 
98  /// Performs back substitution using the output of LUDecomp() for a square
99  /// matrix.
100  /// Input: A[1..n][1..n] and pivots as output from LUDecompRect()
101  /// b[1..n] is the right hand side for solving A.x = b
102  /// Output: b is modified to contain the solution x
103  template <typename S>
104  void LUDecompBackSub(const UT_MatrixT<T> &A,
105  const UT_Permutation &pivots,
106  UT_VectorT<S> &b);
107 
108  /// LU decompose a companion matrix
109  /// The matrix given is the base of the companion matrix:
110  /// [aI bI 0 0 ]
111  /// [0 aI bI 0 ]
112  /// [0 0 aI bI]
113  /// [P1 P2 P3 P4]
114  /// ie: [P1 P2 P3 P4]
115  /// The resulting LU decompostion is written into the R matrix and
116  /// is:
117  /// [aI 0 0 0] [I b/aI 0 0 ]
118  /// [0 aI 0 0] [0 I b/aI 0 ]
119  /// [0 0 aI 0] [0 0 I b/aI]
120  /// [R1 R2 R3 L] [0 0 0 U ]
121  /// LU is stored in the last portion of the matrix & index is the
122  /// permuation that resulted in it. Thus, index should be the height
123  /// of P, NOT the entire companion matrix.
124  ///
125  /// Return 0 on success, -1 for singular, and -2 for singular to
126  /// machine precision (As judged by tol), and -3 for a being near zero.
127  int LUDecompCompanion(const UT_MatrixT<T> &P, T a, T b,
128  UT_MatrixT<T> &R, UT_Permutation &index,
129  T tol=1e-20);
130 
131  /// This takes the result of the previous operation and does a
132  /// back substitution.
133  template <typename S>
134  void LUBackSubCompanion(const UT_MatrixT<T> &R, T a, T b,
135  const UT_Permutation &index,
136  UT_VectorT<S> &y);
137 
138  ///
139  /// Cholesky factorization
140  ///
141  /// To solve a full rank system A . x = b where A can be undertermined
142  /// (ie. A has less rows than columns), we can make use of the Cholesky
143  /// factorization functions below by way of the normal equations like so:
144  /// 1. Left multiply both sides by the tranpose of A, A^T to get
145  /// (A^T . A) . x = (A^T . b)
146  /// 2. Now solve for x using P . x = c where P = A^T . A and c = A^T . b
147  /// This is done by first calling choleskyDecomp() using P since it
148  /// will be symmetric non-negative definite. Then use choleskyBackSub()
149  /// to obtain x for particular values of c.
150 
151  /// Performs a Cholesky decomposition of the symmetric non-negative definite
152  /// matrix A such that A = L . L^T where L^T is tranpose of the lower
153  /// triangular matrix L.
154  /// Input: Symmetric non-negative definite matrix a[1..n][1..n].
155  /// Only the upper triangle of a is read from.
156  /// tol is the tolerance check for when a sum is deemed 0
157  /// Output: Cholesky factor L is returned in the lower triangle of A, except
158  /// for its diagonal elements which is returned in d[1..n]
159  /// Returns the number of articial zero's placed into d. If we
160  /// encounter negative diagonal values, they are zeroed out and we
161  /// return -N where N is the number of negative values found.
162  template <typename S>
163  int choleskyDecomp(UT_MatrixT<T> &a, UT_VectorT<S> &d, T tol=1e-5);
164 
165  /// Solves the set of n linear equations, A.x = b, where a is a nonnegative
166  /// definite symmetric matrix.
167  /// Input: a[1..n][1..n] and d[1..n] as output from choleskyDecomp() above.
168  /// Only the lower triangle of a is read from.
169  /// b[1..n] is the right-hand side of A.x = b to solve for
170  /// tol is the tolerance check for when an element of d is deemed 0
171  /// Output: x[1..n], the solution of A.x = b
172  /// Returns the number of articial zero's placed into x.
173  template <typename S>
174  int choleskyBackSub(const UT_MatrixT<T> &a, const UT_VectorT<S> &d,
175  const UT_VectorT<S> &b, UT_VectorT<S> &x,
176  T tol=1e-5);
177 
178  /// Performs inverse iteration to find the closest eigenvalue(s)
179  /// to s.
180  /// Does a two pass method, first using single iterations until
181  /// starttol is reached, in which case singles are used to
182  /// achieve final tolerance. If starttol is not reached in
183  /// startiter iterations, uses double method to find using
184  /// final tolerance/iterations.
185  ///
186  /// The vector q is updated with the newest eigenvector.
187  /// Source matrices are:
188  /// [1 0 0] [0 1 0]
189  /// C1' = [0 1 0] C2' = [0 0 1]
190  /// [0 0 C1] [ C2 ]
191  /// Thus, the initial sizes of the matrices/vectors are (row*col):
192  /// C1: n*n, C2: n*(n*deg), R: n*(n*deg), tmp: n*(n*deg),
193  /// index: n, q: n*deg, Aq: n*deg
194  /// The result code determines how s1 & s2 are to be interpreted:
195  /// Result:
196  /// 0 - One eigenvalue found at s1, eigenvector q.
197  /// 1 - Two real eigenvalues found: s1 & s2.
198  /// 2 - Complex conjugate eigenvalues found: s1 +- s2i.
199  /// -1 - Failed to converge to final tolerance using single method.
200  /// -2 - Failed to converge to final tolerance using double method.
201  template <typename S>
202  int inversePowerIterate(const UT_MatrixT<T> &C1,
203  const UT_MatrixT<T> &C2,
204  UT_MatrixT<T> &R, UT_MatrixT<T> &tmp,
205  UT_Permutation &index,
206  T s,
207  T &s1, T &s2,
209  T starttol = 1e-2, int startiter = 5,
210  T finaltol = 1e-5,
211  T doubletol = 1e-4, int finaliter = 50);
212 
213  /// Compute the inverse of a LU decomposed matrix.
214  /// If you want to solve systems of equations AX=B,
215  /// where B is a matrix. It is better (faster) to do a
216  /// LU decomposition and use back substituation for columns of B.
217  void inverse(UT_MatrixT<T> &a, UT_Permutation &index,
218  UT_MatrixT<T> &ia);
219 
220  /// Finds all linearly independent rows in A, does NOT require
221  /// A to be any particular shape. Returns the number of independent
222  /// rows & first entries of index are their indices in the original
223  /// matrix.
224  /// Destroys A.
225  int findLinIndRows(UT_MatrixT<T> &A, UT_Permutation &index,
226  T tol = 1E-6);
227 
228  /// Similar to findLinIndRows, performs gaussian elimination
229  /// with full pivoting on G. Updates as well Gu, and Gv, the
230  /// partial derivitive matrix of G(u, v). After this is done,
231  /// the determinant of G & partials of the determinant can
232  /// be computed.
233  void fullGaussianElimination(UT_MatrixT<T> &G, UT_MatrixT<T> &Gu,
234  UT_MatrixT<T> &Gv,
235  UT_Permutation &rowperm,
236  UT_Permutation &colperm,
237  T &detsign);
238 
239  /// After full Gaussian Elimination, this finds the determinant
240  /// and partials thereof:
241  void detWithPartials(const UT_MatrixT<T> &G, const UT_MatrixT<T> &Gu,
242  const UT_MatrixT<T> &Gv,
243  const UT_Permutation &rowperm,
244  const UT_Permutation &colperm,
245  T detsign,
246  T &retdet, T &detu, T &detv);
247 
248 
249  /// Compute the determinant of a LU decomposed matrix.
250  /// The d is what one gets from LU decomposition.
251  T det(UT_MatrixT<T> &a, T d);
252 
253  /// Find the determiment:
254  /// (A(p,q), B(r,s), C(t,u)))
255  inline T det3x3(const UT_MatrixT<T> &A, const UT_MatrixT<T> &B,
256  const UT_MatrixT<T> &C, int p, int q, int r, int s,
257  int t, int u)
258  {
259  return (A(p,q) * (B(r,s) * C(t,u) - B(t,u) * C(r,s)) -
260  B(p,q) * (A(r,s) * C(t,u) - A(t,u) * C(r,s)) +
261  C(p,q) * (A(r,s) * B(t,u) - A(t,u) * B(r,s)) );
262  }
263 
264  /// Estimates condtion number of upper triangular matrix:
265  template <typename S>
266  void condEstimate(const UT_MatrixT<T> &A, UT_VectorT<S> &y);
267 
268  /// Finds the condition number of a LUDecomposed matrix.
269  /// This is just an approximation
270  /// You must provide the infinite norm of A before LUDecomposition.
271  T getConditionLUD(const UT_MatrixT<T> &A,
272  const UT_Permutation &index, T anorm);
273 
274  /// Find householder vector of x.
275  template <typename S>
276  void house(const UT_VectorT<S> &x, UT_VectorT<S> &v, T &b);
277  template <typename S>
278  void house(const UT_MatrixT<T> &x, UT_VectorT<S> &v, T &b);
279 
280  /// Finds the givens rotation required to 0 b:
281  void findGivens(T a, T b, T &c, T &s,
282  T zerotol = 1e-20);
283 
284  /// Find half givens, return 0 if success, otherwise the required
285  /// rotation is impossible (ie:complex eigenvalue)
286  int findHalfGivens(const UT_MatrixT<T> &A, T &c, T &s,
287  T tol = 0.0);
288 
289  /// Upper triangulizes A with householder transforms, applying
290  /// simultaneously to A and optionally Q.
291  void triangularize(UT_MatrixT<T> &A, UT_MatrixT<T> &B,
292  UT_MatrixT<T> *Q = 0);
293 
294  /// Transforms A into upper hessenburg and B into upper triangular.
295  void hessentri(UT_MatrixT<T> &A, UT_MatrixT<T> &B);
296 
297  /// Converts A into quasi upper triangular & B into upper triangular
298  /// form using the QZ algorithm.
299  int realSchurFormQZ(UT_MatrixT<T> &oA, UT_MatrixT<T> &oB,
300  T tol=1e-6,
301  int maxIteration = 30);
302 
303  /// Transform A into a hessenberg matrix through a series of
304  /// householder reflections.
305  void hessenberg(UT_MatrixT<T> &A);
306  /// Computes P st A' = P^T A P
307  void hessenberg(UT_MatrixT<T> &A, UT_MatrixT<T> &P);
308 
309  /// Transforms A into a bidiagonal matrix via householder transformations.
310  int bidiagonalize(UT_MatrixT<T> &A);
311 
312  /// Reduces the matrix to real Schur form:
313  /// Orthogonal transform responible for getting to real schur form
314  /// is accumulated in Q, if it is non-zero.
315  int realSchurForm(UT_MatrixT<T> &A, UT_MatrixT<T> *Q = 0,
316  T tol = 1E-6,
317  int maxIteration = 30);
318 
319  /// Golub-Kahan SVD algorithm. Resulting SVD matrix found in
320  /// diaganol of A.
321  int gkSVD(UT_MatrixT<T> &A,
322  T tol = 1e-6,
323  int maxIteration = 30);
324 
325  /// Use Numerical Recipes method for SVD decomposition, with results
326  /// that are suitable for doing back substitution.
327  /// NOTE: The UT_MatrixF should start with 1,1!
328  ///
329  /// Decomposition a = u * w * v^T
330  /// a[1..m][1..n] also stores the result u[1..m][1..n]
331  /// w[1..n][1..n] is a diagonal matrix, stored as a vector
332  /// v[1..n][1..n]
333  template <typename S>
334  int SVDDecomp(UT_MatrixT<T> &a, UT_VectorT<S> &w, UT_MatrixT<T> &v,
335  int maxIterstion = 30);
336  /// Perform back substitution on an SVD decomposed matrix. Solves
337  /// equation Ax = b.
338  template <typename S>
339  int SVDBackSub(UT_MatrixT<T> &u, UT_VectorT<S> &w, UT_MatrixT<T> &v,
340  UT_VectorT<S> &b, UT_VectorT<S> &x, T tol = 1e-6);
341 
342  /// Finds eigenvalues from a real schur form matrix:
343  template <typename S>
344  void findEigenValues(const UT_MatrixT<T> &A,
345  UT_VectorT<S> &vreal, UT_VectorT<S> &vimg);
346  /// Finds eigenvalues from generalized system (A - lambdaB)
347  template <typename S>
348  void findEigenValues(const UT_MatrixT<T> &A, const UT_MatrixT<T> &B,
349  UT_VectorT<S> &sreal, UT_VectorT<S> &simg,
350  UT_VectorT<S> &treal, UT_VectorT<S> &timg);
351 
352  /// Finds the u & v eigenvalues from an eigenvector, whose udeg,
353  /// vdeg, and tdeg are as specified. Returns 0 if fails.
354  /// NB: Does the s/(1+s) transformation!
355  template <typename S>
356  int findUVfromEigenvector(const UT_VectorT<S> &vect,
357  int udeg, int vdeg, int tdeg, T &u, T &v);
358 
359  /// Finds a specified eigenvector from a real schur form matrix.
360  /// The eigenvector must be real!
361  template <typename S>
362  void findRightEigenVector(const UT_MatrixT<T> &A,
363  const UT_MatrixT<T> &Q, UT_VectorT<S> &vector, int i);
364 
365 
366 private:
367  /// Pushes down zeros, causing A to be a reduced hessenberg.
368  void chaseZerosQZ(UT_MatrixT<T> &A, UT_MatrixT<T> &B, int zeropos);
369 
370  /// Pushes across zeros generated in Golub-Kahan SVD method
371  void chaseZerosGK(UT_MatrixT<T> &B, int zeropos);
372 
373  /// Applies one iteration of the QZ algorithm to unreduced
374  /// upper hessenberg A & non-singular upper triangular B.
375  void QZStep(UT_MatrixT<T> &A, UT_MatrixT<T> &B);
376 
377  /// Performs the francisQRStep on H.
378  void francisQRStep(UT_MatrixT<T> &H);
379 
380  // Performs francisQRStep, updating H12 = H12*Z,
381  // Q = Q*Z, H23 = ZT * H23 simultaneously with finding H = ZT H Z
382  void francisQRStep(
383  UT_MatrixT<T> &H,
384  UT_MatrixT<T> &H12,
385  UT_MatrixT<T> &H23,
386  UT_MatrixT<T> &Q);
387 
388  /// Performs one step in the Golub-Kahan SVD algorithm:
389  /// B is a bidiagonal matrix with no zeros on diagonal nor supradiagonal,
390  void golubKahanSVDStep(UT_MatrixT<T> &B);
391 
392  /// Solves (A - lambda I) x = gamma B for x, requires A is quasi triangular.
393  /// B is column.
394  template <typename S>
395  void solveQuasiUpperTriangular(const UT_MatrixT<T> &A, T lambda,
396  T gamma,
397  const UT_MatrixT<T> &B, UT_VectorT<S> &x,
398  T tol = 1e-6);
399 
400  /// Performs inverse iteration until convergence or maxiter is reached.
401  /// The vector is updated with the newest eigenvector.
402  /// Call with findDecomp = 0 to use R directly as opposed to decomposing
403  /// it. Source matrices are:
404  /// [1 0 0] [0 1 0]
405  /// C1' = [0 1 0] C2' = [0 0 1]
406  /// [0 0 C1] [ C2 ]
407  /// Returns:
408  /// 1: Decomp found the pencil to be degenerate, perfect solution.
409  /// 0: Successfully converged
410  /// <0: Failed to converge
411  template <typename S>
412  int inversePowerSingle( const UT_MatrixT<T> &C1,
413  const UT_MatrixT<T> &C2,
414  UT_MatrixT<T> &R, UT_MatrixT<T> &tmp,
415  UT_Permutation &index,
416  T &s, T lasts,
418  int findDecomp = 1,
419  T tol = 1e-5, int maxiter = 50);
420 
421  /// Performs inverse iteration until convergence or maxiter is reached.
422  /// The vector is updated with the newest eigenvector.
423  /// Call with findDecomp = 0 to use R directly as opposed to decomposing
424  /// it. Source matrices are:
425  /// [1 0 0] [0 1 0]
426  /// C1' = [0 1 0] C2' = [0 0 1]
427  /// [0 0 C1] [ C2 ]
428  /// Computes two solutions, namely P & Q such that
429  /// A = (C2 - sC1), and
430  /// lim i->inf (A^i+2 - pA^i+1 - qA^i) u = 0
431  /// Returns:
432  /// 1: Decomp found the pencil to be degenerate, perfect solution.
433  /// 0: Successfully converged
434  /// <0: Failed to converge
435  template <typename S>
436  int inversePowerDouble( const UT_MatrixT<T> &C1,
437  const UT_MatrixT<T> &C2,
438  UT_MatrixT<T> &R, UT_MatrixT<T> &tmp,
439  UT_Permutation &index,
440  T s,
441  T &p, T &q,
443  int findDecomp = 1,
444  T tol = 1e-5, int maxiter = 50);
445 
446 };
447 
451 typedef UT_MatrixSolverT<fpreal64> UT_MatrixSolver; // Matches UT_MatrixF
452 
453 // UT_MatrixIterSolver is mainly for sparse matrices.
454 // The sparse matrices may be defined explicitly with UT_SparseMatrix,
455 // or implicitly with functions.
456 template <typename T>
458 {
459 public:
465 
466  // Creates a new iterative matrix solver.
467  // Inputs:
468  // 1. maxIterations: The upper bound of iterations at which the iterative
469  // solver will automatically stop.
470  // If this is < 0 then the default will be 3*N where N
471  // is the length of the solving vector
472  // 2. tolerance: The maximum error at which the iteration will stop
473  // 3. useAbsoluteError: If true, the tolerance will be used as an absolute
474  // tolerance. Otherwise, the tolerance is divided by
475  // the vector norm to produce a relative tolerance.
476  UT_MatrixIterSolverT(exint maxIterations = -1,
477  UT_PreciseT<T> tolerance = 1e-5, bool useAbsoluteError = false);
478 
479  // Sets the operation for A*x using a custom functor
480  void setMultVec(const MatrixOp &op);
481 
482  // Sets the operation for A*x using UT_SparseMatrix::multVec.
483  // Also sets myMultVecAndDotOp to UT_SparseMatrix::multVecAndDot to speed
484  // up solvers that can use it.
485  void setMultVec(const SparseMatrixType &matrix);
486 
487  // Sets the operation for At*x using a custom functor
488  void setTransposeMultVec(const MatrixOp &op);
489 
490  // Sets the operation for At*x using UT_SparseMatrix::multVec
491  void setTransposeMultVec(const SparseMatrixType &transposematrix);
492 
493  // Sets the operation for inv(M)*x where M is the preconditioner using a
494  // custom functor
495  void setPreconditioner(const MatrixOp &op);
496 
497  // Sets the operation for inv(M)*x using an inverted diagonal vector and
498  // UT_Vector::multAndSet
499  void setPreconditioner(const UT_VectorT<T> &invertedDiagonal);
500 
501  // Sets the operation for inv(M)*x using lower and upper triangular matrices
502  // For r = inv(LU)*x, this solves Lz = x, Ur = z.
503  // L and U can refer to the same matrix if the matrix is storing both.
504  // Set LhasUnitDiagonal/UhasUnitDiagonal to true if either L/U should have
505  // its diagonal assumed to be all ones.
506  void setPreconditioner(const SparseMatrixType &L, const SparseMatrixType &U,
507  bool LhasUnitDiagonal = false, bool UhasUnitDiagonal = false);
508 
509  // Sets the operation for inv(Mt)*x. The usage is identical to
510  // setPreconditioner.
511  void setTransposePreconditioner(const MatrixOp &op);
512  void setTransposePreconditioner(const UT_VectorT<T> &invertedDiagonal);
513 
514  // Sets the operation to determine if iteration should prematurely stop.
515  // The operation is called every iteration and will break the iteration if
516  // it returns true.
517  void setStopCondition(const StopConditionOp &op);
518 
519  // Solves Ax = b using the Preconditioned Conjugate Gradient method,
520  // where A is a symmetric positive definite matrix.
521  // Requires MultVec and PreconditionerOp to be set.
522  // myMultVecAndDotOp will be used if set for a performance increase.
523  void PCG(UT_VectorT<T> &x, const UT_VectorT<T> &b,
524  exint *resultIters = NULL, UT_PreciseT<T> *resultError = NULL,
525  UT_Interrupt *boss = NULL) const;
526 
527  // Solves Ax = b using the Preconditioned Conjugate Gradient method
528  // for least squares. i.e. Computes vector x that minimizes the Euclidean
529  // 2-norm ||b - Ax||^2. A does not need to be square.
530  // Given A is an MxN matrix, this requires length(x) == N, length(b) == M.
531  // Requires MultVec, TransposeMultVec and PreconditionerOp to be set.
532  void PCGLS(UT_VectorT<T> &x, const UT_VectorT<T> &b,
533  exint *resultIters = NULL, UT_PreciseT<T> *resultError = NULL,
534  UT_Interrupt *boss = NULL) const;
535 
536  // Solve Ax = b using the Preconditioned Biconjugate Gradient method.
537  // A must be square but does not need to be symmetric or self-adjoint.
538  // Requires MultVec, TransposeMultVec, PreconditionerOp, and
539  // TransposePreconditionerOp to be set.
540  void PBiCG(UT_VectorT<T> &x, const UT_VectorT<T> &b,
541  exint *resultIters = NULL, UT_PreciseT<T> *resultError = NULL,
542  UT_Interrupt *boss = NULL) const;
543 
544  // Solves Ax = b using the Preconditioned Biconjugate Gradient Stabilized
545  // method. A must be square but does not need to be symmetric/self-adjoint.
546  // Requires MultVec and PreconditionerOp to be set.
547  void PBiCGStab(UT_VectorT<T> &x, const UT_VectorT<T> &b,
548  exint *resultIters = NULL, UT_PreciseT<T> *resultError = NULL,
549  UT_Interrupt *boss = NULL) const;
550 
551 private:
552 
553  exint myMaxIterations;
554  UT_PreciseT<T> myTolerance;
555  bool myUseAbsError;
556 
557  MatrixOp myMultVecOp;
558  MatrixOp myTransposeMultVecOp;
559  MatrixOp myPreconditionerOp;
560  MatrixOp myTransposePreconditionerOp;
561  StopConditionOp myStopConditionOp;
562  MatrixDotOp myMultVecAndDotOp;
563 
564 public:
565  // TODO: Deprecate these
566 
567  /// Preconditioned Conjugate Gradient method for solving Ax=b,
568  /// given the symmetric positive definite matrix A[0..n-1][0..n-1]
569  /// indirectly through AMult (ie. Ax),
570  /// a preconditioner for A through ASolve,
571  /// the vector b[0..n-1] and an initial guess x[0..n-1].
572  /// maxIter defaults to 2n.
573  /// error |b-Ax|/|b| < tol.
574  /// normType indicates the type of norm used for error:
575  /// 0 L-infinity norm (ie. max)
576  /// 1 L1-norm (ie. sum of abs)
577  /// 2 L2-norm (ie. Euclidean distance)
578  /// Output: x[0..n-1] and error.
579  T PCG(void (*AMult)(const UT_VectorT<T> &x, UT_VectorT<T> &result),
580  void (*ASolve)(const UT_VectorT<T> &b, UT_VectorT<T> &x),
581  int n, UT_VectorT<T> &x, const UT_VectorT<T> &b,
582  T tol=1e-3, int normType=2, int maxIter=-1);
583 
584  /// Preconditioned Conjugate Gradient method for solving Ax=b with
585  /// the symmetric positive definite matrix A, preconditioner, and iteration
586  /// conditions implicitly defined by functors.
587  static void PCG(UT_VectorT<T> &x, const UT_VectorT<T> &b,
588  const UT_Functor2<void, const UT_VectorT<T> &, UT_VectorT<T> &> &AMult,
589  const UT_Functor2<void, const UT_VectorT<T> &, UT_VectorT<T> &> &ASolve,
590  const UT_Functor2<bool, int, const UT_VectorT<T> &> &iterateTest);
591 
592  /// Preconditioned Conjugate Gradient method for solving least squares Ax=b,
593  /// given the matrix A[0..m-1][0..n-1] indirectly through AMult
594  /// (ie. Ax if int = 0 or A^tx if int = 1),
595  /// a preconditioner for A^tA through AtASolve,
596  /// the vector b[0..m-1] and an initial guess x[0..n-1].
597  /// maxIter defaults to 2n.
598  /// error |A^t(b-Ax)|/|A^tb| < tol.
599  /// normType indicates the type of norm used for error:
600  /// 0 L-infinity norm (ie. max)
601  /// 1 L1-norm (ie. sum of abs)
602  /// 2 L2-norm (ie. Euclidean distance)
603  /// Output: x[0..n-1] and error.
604  T PCGLS(void (*AMult)(const UT_VectorT<T> &x, UT_VectorT<T> &result,
605  int transpose),
606  void (*AtASolve)(const UT_VectorT<T> &b, UT_VectorT<T> &x),
607  int m, int n,
608  UT_VectorT<T> &x, const UT_VectorT<T> &b,
609  T tol=1e-3, int normType=2, int maxIter=-1);
610 
611  ///
612  /// Modified interface for the above Preconditioned Conjugate Gradient
613  /// for solving Ax=b. The template parameter class C must implement
614  /// the following two methods:
615  ///
616  ///
617  /// 1) A method to compute either Av or A^tv depending on transpose
618  /// being 0 or 1.
619  ///
620  /// void AMult(const UT_VectorT<T> &v, UT_VectorT<T> &result,
621  /// bool transpose);
622  ///
623  ///
624  /// 2) The preconditioner as above.
625  ///
626  /// void AtASolve(const UT_VectorT<fpreal> &b, UT_VectorT<fpreal> &x)
627  ///
628  /// The rel_tol and abs_tol respectively represent relative and absolute
629  /// tolerance. Set either to a positive number to stop the solver as soon
630  /// as the corresponding error drops below the given threshold. If both
631  /// positive, the solver stops as soon as one of the errors meets its
632  /// corresponding threshold.
633  ///
634  template <typename AMult, typename AtASolve>
635  T PCGLS(const AMult &amult, const AtASolve &atasolve,
636  int rows, int cols, UT_VectorT<T> &x, const UT_VectorT<T> &b,
637  fpreal tol = -1.0, int norm_type = 2, int max_iter = -1);
638 
639  /// Preconditioned Conjugate Gradient method for solving Ax=b,
640  /// given the matrix A[0..m-1][0..n-1] indirectly through AMult
641  /// (ie. Ax if int = 0 or A^tx if int = 1),
642  /// a preconditioner for A through ASolve,
643  /// the vector b[0..m-1] and an initial guess x[0..n-1].
644  /// maxIter defaults to 2n.
645  /// error |b-Ax|/|b| < tol.
646  /// normType indicates the type of norm used for error:
647  /// 0 L-infinity norm (ie. max)
648  /// 1 L1-norm (ie. sum of abs)
649  /// 2 L2-norm (ie. Euclidean distance)
650  /// Output: x[0..n-1] and error.
651  T biPCG(void (*AMult)(const UT_VectorT<T> &x, UT_VectorT<T> &result,
652  int transpose),
653  void (*ASolve)(const UT_VectorT<T> &b, UT_VectorT<T> &x,
654  int transpose),
655  int n,
656  UT_VectorT<T> &x, const UT_VectorT<T> &b,
657  T tol=1e-3, int normType=2, int maxIter=-1);
658 };
659 
664 
665 #include "UT_MatrixSolverImpl.h"
666 
667 #endif
UT_MatrixSolverT< fpreal > UT_MatrixSolverR
void LUBackSub(const UT_MatrixT< T > &a, const UT_Permutation &index, UT_VectorT< S > &b)
png_voidp s1
Definition: png.h:2193
const GLdouble * v
Definition: glcorearb.h:836
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
UT_MatrixIterSolverT< fpreal64 > UT_MatrixIterSolverD
#define UT_API
Definition: UT_API.h:12
GLint y
Definition: glcorearb.h:102
UT_Functor2< void, const UT_VectorT< T > &, UT_VectorT< T > & > MatrixOp
UT_SparseMatrixCSRT< T > SparseMatrixType
UT_MatrixSolverT< fpreal32 > UT_MatrixSolverF
png_uint_32 i
Definition: png.h:2877
GLsizei GLboolean transpose
Definition: glcorearb.h:831
GLdouble n
Definition: glcorearb.h:2007
int64 exint
Definition: SYS_Types.h:109
UT_Functor2< bool, exint, T > StopConditionOp
UT_MatrixSolverT< fpreal64 > UT_MatrixSolverD
typename UT_TypePromoteT< T >::PreciseType UT_PreciseT
UT_MatrixIterSolverT< fpreal > UT_MatrixIterSolverR
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
UT_MatrixIterSolverT< fpreal64 > UT_MatrixIterSolver
int LUDecomp(UT_MatrixT< T > &a, UT_Permutation &index, T &d, T tol=1e-20)
double fpreal
Definition: SYS_Types.h:263
png_voidp png_voidp s2
Definition: png.h:2193
UT_MatrixIterSolverT< fpreal32 > UT_MatrixIterSolverF
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)
GLuint index
Definition: glcorearb.h:785
UT_Functor2< UT_PreciseT< T >, const UT_VectorT< T > &, UT_VectorT< T > & > MatrixDotOp
GLint GLenum GLint x
Definition: glcorearb.h:408
GLubyte GLubyte GLubyte GLubyte w
Definition: glcorearb.h:856
GLboolean r
Definition: glcorearb.h:1221
UT_MatrixSolverT< fpreal64 > UT_MatrixSolver