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 
80  for (exint iter = 0; iter < max_iter; ++iter)
81  {
82  if (!(iter % 256) && boss.wasInterrupted())
83  break;
84 
85  amult(p, tmp, 0);
86  alpha = abs_new / dot(tmp, tmp);
87  x.addScaledVec(alpha, p);
88  r.addScaledVec(-alpha, tmp);
89  amult(r, s, 1); // s = A^t * r
90 
91  if (norm_type == 2)
92  residual_norm2 = dot(s, s);
93  else
94  {
95  residual_norm2 = s.norm(norm_type);
96  residual_norm2 *= residual_norm2;
97  }
98  if (residual_norm2 < threshold)
99  break;
100 
101  atasolve(s, z);
102  fpreal64 absOld = abs_new;
103  abs_new = dot(s, z);
104  p.scaleAddVec(abs_new / absOld, z);
105  }
106 
107  return residual_norm2 / rhs_norm2;
108 }
109 
110 
111 #endif // __UT_MATRIXSOLVERIMPL_H_INCLUDED__
void zero()
Initialize to zeros.
Definition: UT_Vector.h:57
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
int64 exint
Definition: SYS_Types.h:125
GLint GLenum GLint x
Definition: glcorearb.h:409
double fpreal64
Definition: SYS_Types.h:201
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
Definition: CE_Vector.h:130
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:848
GLfloat GLfloat p
Definition: glew.h:16656
void PCGLS(UT_VectorT< T > &x, const UT_VectorT< T > &b, exint *resultIters=NULL, UT_PreciseT< T > *resultError=NULL, UT_Interrupt *boss=NULL) const
fpreal64 fpreal
Definition: SYS_Types.h:277
GLboolean r
Definition: glcorearb.h:1222
GLfloat GLfloat GLfloat alpha
Definition: glcorearb.h:112
GLdouble s
Definition: glew.h:1395