HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UT_MatrixSolverImpl.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: Utility Library (C++)
7  *
8  * COMMENTS:
9  * This is meant to be included by UT_MatrixSolver.h and includes
10  * the template implementations needed by external code.
11  */
12 
13 #pragma once
14 
15 #ifndef __UT_MATRIXSOLVERIMPL_H_INCLUDED__
16 #define __UT_MATRIXSOLVERIMPL_H_INCLUDED__
17 
18 #include "UT_Interrupt.h"
19 
20 
21 template <typename T>
22 template <typename AMult, typename AtASolve>
23 T
24 UT_MatrixIterSolverT<T>::PCGLS(const AMult &amult, const AtASolve &atasolve,
25  int rows, int cols,
26  UT_VectorT<T> &x, const UT_VectorT<T> &b,
27  fpreal tol, int norm_type, int max_iter)
28 {
29  UT_AutoInterrupt boss("Solving PCGLS");
30 
31  fpreal64 alpha, abs_new, residual_norm2, rhs_norm2;
32 
33  UT_VectorT<T> r(0, rows - 1), tmp(0, rows - 1);
34  UT_VectorT<T> p(0, cols - 1), s(0, cols - 1), z(0, cols - 1);
35 
36 
37  // s is the normal residual.
38  amult(b, s, 1); // s = A^t * b
39 
40  if (norm_type == 2)
41  rhs_norm2 = dot(s, s);
42  else
43  {
44  rhs_norm2 = s.norm(norm_type);
45  rhs_norm2 *= rhs_norm2;
46  }
47 
48  if (rhs_norm2 == 0.0)
49  {
50  x.zero();
51  return 0;
52  }
53 
54  fpreal64 threshold = tol * tol * rhs_norm2;
55 
56  if (max_iter < 0)
57  max_iter = 2 * (rows + cols);
58 
59  amult(x, r, 0); // Residual r = b - Ax
60  r.negPlus(b);
61 
62  amult(r, s, 1); // s = A^t * r
63 
64  if (norm_type == 2)
65  residual_norm2 = dot(s, s);
66  else
67  {
68  residual_norm2 = s.norm(norm_type);
69  residual_norm2 *= residual_norm2;
70  }
71 
72  if (residual_norm2 < threshold)
73  {
74  return SYSsqrt(residual_norm2 / rhs_norm2);
75  }
76 
77  atasolve(s, p);
78  abs_new = dot(s, p);
79  p = p;
80 
81  for (exint iter = 0; iter < max_iter; ++iter)
82  {
83  if (!(iter % 256) && boss.wasInterrupted())
84  break;
85 
86  amult(p, tmp, 0);
87  alpha = abs_new / dot(tmp, tmp);
88  x.addScaledVec(alpha, p);
89  r.addScaledVec(-alpha, tmp);
90  amult(r, s, 1); // s = A^t * r
91 
92  if (norm_type == 2)
93  residual_norm2 = dot(s, s);
94  else
95  {
96  residual_norm2 = s.norm(norm_type);
97  residual_norm2 *= residual_norm2;
98  }
99  if (residual_norm2 < threshold)
100  break;
101 
102  atasolve(s, z);
103  fpreal64 absOld = abs_new;
104  abs_new = dot(s, z);
105  p.scaleAddVec(abs_new / absOld, z);
106  }
107 
108  return residual_norm2 / rhs_norm2;
109 }
110 
111 
112 #endif // __UT_MATRIXSOLVERIMPL_H_INCLUDED__
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:847
int64 exint
Definition: SYS_Types.h:116
double fpreal64
Definition: SYS_Types.h:192
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
Definition: CE_Vector.h:218
GLfloat GLfloat GLfloat alpha
Definition: glcorearb.h:111
void PCGLS(UT_VectorT< T > &x, const UT_VectorT< T > &b, exint *resultIters=NULL, UT_PreciseT< T > *resultError=NULL, UT_Interrupt *boss=NULL) const
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
double fpreal
Definition: SYS_Types.h:270
void zero(exint nl, exint nh)
Initialize to zeros.
Definition: UT_Vector.C:158
GLint GLenum GLint x
Definition: glcorearb.h:408
GLboolean r
Definition: glcorearb.h:1221