24 #ifndef SVD_EIGEN_HH_HPP
25 #define SVD_EIGEN_HH_HPP
67 template<
typename MatrixA,
typename MatrixUV,
typename VectorS>
69 const Eigen::MatrixBase<MatrixA>& A,
70 Eigen::MatrixBase<MatrixUV>&
U,
71 Eigen::MatrixBase<VectorS>& S,
72 Eigen::MatrixBase<MatrixUV>&
V,
73 Eigen::MatrixBase<VectorS>& tmp,
77 const int rows =
A.rows();
78 const int cols =
A.cols();
82 int i(-1),its(-1),j(-1),jj(-1),k(-1),nm=0;
85 e_scalar maxarg1,maxarg2,anorm(0),
c(0),f(0),h(0),s(0),scale(0),
x(0),
y(0),
z(0),
g(0);
90 for (i=0;i<cols;i++) {
96 for (k=i;k<rows;k++) scale +=
fabs(
U(k,i));
100 for (k=i;k<rows;k++) {
108 for (j=ppi;j<cols;j++) {
110 for (s=0.0,k=i;k<rows;k++) s +=
U(k,i)*
U(k,j);
113 for (k=i;k<rows;k++)
U(k,j) += f*
U(k,i);
115 for (k=i;k<rows;k++)
U(k,i) *= scale;
121 if ((i <rows) && (i+1 != cols)) {
123 for (k=ppi;k<cols;k++) scale +=
fabs(
U(i,k));
125 for (k=ppi;k<cols;k++) {
133 for (k=ppi;k<cols;k++) tmp(k)=
U(i,k)/h;
134 for (j=ppi;j<rows;j++) {
135 for (s=0.0,k=ppi;k<cols;k++) s +=
U(j,k)*
U(i,k);
136 for (k=ppi;k<cols;k++)
U(j,k) += s*tmp(k);
138 for (k=ppi;k<cols;k++)
U(i,k) *= scale;
143 anorm = maxarg1 > maxarg2 ? maxarg1 : maxarg2;
146 for (i=cols-1;i>=0;i--) {
149 for (j=ppi;j<cols;j++)
V(j,i)=(
U(i,j)/
U(i,ppi))/
g;
150 for (j=ppi;j<cols;j++) {
151 for (s=0.0,k=ppi;k<cols;k++) s +=
U(i,k)*
V(k,j);
152 for (k=ppi;k<cols;k++)
V(k,j) += s*
V(k,i);
155 for (j=ppi;j<cols;j++)
V(i,j)=
V(j,i)=0.0;
162 for (i=cols-1<rows-1 ? cols-1:rows-1;i>=0;i--) {
165 for (j=ppi;j<cols;j++)
U(i,j)=0.0;
168 for (j=ppi;j<cols;j++) {
169 for (s=0.0,k=ppi;k<rows;k++) s +=
U(k,i)*
U(k,j);
171 for (k=i;k<rows;k++)
U(k,j) += f*
U(k,i);
173 for (j=i;j<rows;j++)
U(j,i) *=
g;
175 for (j=i;j<rows;j++)
U(j,i)=0.0;
181 for (k=cols-1;k>=0;k--) {
182 for (its=1;its<=maxiter;its++) {
184 for (ppi=k;ppi>=0;ppi--) {
186 if ((
fabs(tmp(ppi))+anorm) == anorm) {
190 if ((
fabs(S(nm)+anorm) == anorm))
break;
195 for (i=ppi;i<=k;i++) {
198 if ((
fabs(f)+anorm) == anorm)
break;
205 for (j=0;j<rows;j++) {
218 for (j=0;j<cols;j++)
V(j,k)=-
V(j,k);
235 for (j=ppi;j<=nm;j++) {
249 for (jj=0;jj<cols;jj++) {
264 for (jj=0;jj<rows;jj++) {
278 for (i=0; i<cols; i++){
282 for (j=i+1; j<cols; j++){
296 U.col(i).swap(
U.col(i_max));
297 V.col(i).swap(
V.col(i_max));
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble z
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint y
ccl_device_inline float2 fabs(const float2 &a)
INLINE Rall1d< T, V, S > sqrt(const Rall1d< T, V, S > &arg)
Scalar PYTHAG(Scalar a, Scalar b)
int svd_eigen_HH(const Eigen::MatrixBase< MatrixA > &A, Eigen::MatrixBase< MatrixUV > &U, Eigen::MatrixBase< VectorS > &S, Eigen::MatrixBase< MatrixUV > &V, Eigen::MatrixBase< VectorS > &tmp, int maxiter=150)
Scalar SIGN(Scalar a, Scalar b)
static const pxr::TfToken b("b", pxr::TfToken::Immortal)
static const pxr::TfToken g("g", pxr::TfToken::Immortal)
CCL_NAMESPACE_BEGIN struct Window V