diff options
Diffstat (limited to 'intern/iksolver')
28 files changed, 84 insertions, 6068 deletions
diff --git a/intern/iksolver/CMakeLists.txt b/intern/iksolver/CMakeLists.txt index 4e5699f1abd..67e7aa6abd3 100644 --- a/intern/iksolver/CMakeLists.txt +++ b/intern/iksolver/CMakeLists.txt @@ -44,30 +44,6 @@ set(SRC intern/IK_QJacobianSolver.h intern/IK_QSegment.h intern/IK_QTask.h - intern/TNT/cholesky.h - intern/TNT/cmat.h - intern/TNT/fcscmat.h - intern/TNT/fmat.h - intern/TNT/fortran.h - intern/TNT/fspvec.h - intern/TNT/index.h - intern/TNT/lapack.h - intern/TNT/lu.h - intern/TNT/qr.h - intern/TNT/region1d.h - intern/TNT/region2d.h - intern/TNT/stopwatch.h - intern/TNT/subscript.h - intern/TNT/svd.h - intern/TNT/tnt.h - intern/TNT/tntmath.h - intern/TNT/tntreqs.h - intern/TNT/transv.h - intern/TNT/triang.h - intern/TNT/trisolve.h - intern/TNT/vec.h - intern/TNT/vecadaptor.h - intern/TNT/version.h ) blender_add_lib(bf_intern_iksolver "${SRC}" "${INC}" "${INC_SYS}") diff --git a/intern/iksolver/intern/IK_Math.h b/intern/iksolver/intern/IK_Math.h index ba0eb047c8a..78d99837c43 100644 --- a/intern/iksolver/intern/IK_Math.h +++ b/intern/iksolver/intern/IK_Math.h @@ -33,9 +33,11 @@ #include <cmath> +using Eigen::Affine3d; using Eigen::Matrix3d; +using Eigen::MatrixXd; using Eigen::Vector3d; -using Eigen::Affine3d; +using Eigen::VectorXd; static const double IK_EPSILON = 1e-20; diff --git a/intern/iksolver/intern/IK_QJacobian.cpp b/intern/iksolver/intern/IK_QJacobian.cpp index cd77d5d662d..2925cadf435 100644 --- a/intern/iksolver/intern/IK_QJacobian.cpp +++ b/intern/iksolver/intern/IK_QJacobian.cpp @@ -32,7 +32,6 @@ #include "IK_QJacobian.h" -#include "TNT/svd.h" IK_QJacobian::IK_QJacobian() : m_sdls(true), m_min_damp(1.0) @@ -48,59 +47,51 @@ void IK_QJacobian::ArmMatrices(int dof, int task_size) m_dof = dof; m_task_size = task_size; - m_jacobian.newsize(task_size, dof); - m_jacobian = 0; + m_jacobian.resize(task_size, dof); + m_jacobian.setZero(); - m_alpha.newsize(dof); - m_alpha = 0; + m_alpha.resize(dof); + m_alpha.setZero(); - m_null.newsize(dof, dof); + m_nullspace.resize(dof, dof); - m_d_theta.newsize(dof); - m_d_theta_tmp.newsize(dof); - m_d_norm_weight.newsize(dof); + m_d_theta.resize(dof); + m_d_theta_tmp.resize(dof); + m_d_norm_weight.resize(dof); - m_norm.newsize(dof); - m_norm = 0.0; + m_norm.resize(dof); + m_norm.setZero(); - m_beta.newsize(task_size); + m_beta.resize(task_size); - m_weight.newsize(dof); - m_weight_sqrt.newsize(dof); - m_weight = 1.0; - m_weight_sqrt = 1.0; + m_weight.resize(dof); + m_weight_sqrt.resize(dof); + m_weight.setOnes(); + m_weight_sqrt.setOnes(); if (task_size >= dof) { m_transpose = false; - m_jacobian_tmp.newsize(task_size, dof); + m_jacobian_tmp.resize(task_size, dof); - m_svd_u.newsize(task_size, dof); - m_svd_v.newsize(dof, dof); - m_svd_w.newsize(dof); + m_svd_u.resize(task_size, dof); + m_svd_v.resize(dof, dof); + m_svd_w.resize(dof); - m_work1.newsize(task_size); - m_work2.newsize(dof); - - m_svd_u_t.newsize(dof, task_size); - m_svd_u_beta.newsize(dof); + m_svd_u_beta.resize(dof); } else { // use the SVD of the transpose jacobian, it works just as well // as the original, and often allows using smaller matrices. m_transpose = true; - m_jacobian_tmp.newsize(dof, task_size); - - m_svd_u.newsize(task_size, task_size); - m_svd_v.newsize(dof, task_size); - m_svd_w.newsize(task_size); + m_jacobian_tmp.resize(dof, task_size); - m_work1.newsize(dof); - m_work2.newsize(task_size); + m_svd_u.resize(task_size, task_size); + m_svd_v.resize(dof, task_size); + m_svd_w.resize(task_size); - m_svd_u_t.newsize(task_size, task_size); - m_svd_u_beta.newsize(task_size); + m_svd_u_beta.resize(task_size); } } @@ -113,9 +104,9 @@ void IK_QJacobian::SetBetas(int id, int, const Vector3d& v) void IK_QJacobian::SetDerivatives(int id, int dof_id, const Vector3d& v, double norm_weight) { - m_jacobian[id + 0][dof_id] = v.x() * m_weight_sqrt[dof_id]; - m_jacobian[id + 1][dof_id] = v.y() * m_weight_sqrt[dof_id]; - m_jacobian[id + 2][dof_id] = v.z() * m_weight_sqrt[dof_id]; + m_jacobian(id + 0, dof_id) = v.x() * m_weight_sqrt[dof_id]; + m_jacobian(id + 1, dof_id) = v.y() * m_weight_sqrt[dof_id]; + m_jacobian(id + 2, dof_id) = v.z() * m_weight_sqrt[dof_id]; m_d_norm_weight[dof_id] = norm_weight; } @@ -125,14 +116,18 @@ void IK_QJacobian::Invert() if (m_transpose) { // SVD will decompose Jt into V*W*Ut with U,V orthogonal and W diagonal, // so J = U*W*Vt and Jinv = V*Winv*Ut - TNT::transpose(m_jacobian, m_jacobian_tmp); - TNT::SVD(m_jacobian_tmp, m_svd_v, m_svd_w, m_svd_u, m_work1, m_work2); + Eigen::JacobiSVD<MatrixXd> svd(m_jacobian.transpose(), Eigen::ComputeThinU | Eigen::ComputeThinV); + m_svd_u = svd.matrixV(); + m_svd_w = svd.singularValues(); + m_svd_v = svd.matrixU(); } else { // SVD will decompose J into U*W*Vt with U,V orthogonal and W diagonal, // so Jinv = V*Winv*Ut - m_jacobian_tmp = m_jacobian; - TNT::SVD(m_jacobian_tmp, m_svd_u, m_svd_w, m_svd_v, m_work1, m_work2); + Eigen::JacobiSVD<MatrixXd> svd(m_jacobian, Eigen::ComputeThinU | Eigen::ComputeThinV); + m_svd_u = svd.matrixU(); + m_svd_w = svd.singularValues(); + m_svd_v = svd.matrixV(); } if (m_sdls) @@ -154,26 +149,24 @@ bool IK_QJacobian::ComputeNullProjection() if (rank < m_task_size) return false; - TMatrix basis(m_svd_v.num_rows(), rank); - TMatrix basis_t(rank, m_svd_v.num_rows()); + MatrixXd basis(m_svd_v.rows(), rank); int b = 0; for (i = 0; i < m_svd_w.size(); i++) if (m_svd_w[i] > epsilon) { - for (j = 0; j < m_svd_v.num_rows(); j++) - basis[j][b] = m_svd_v[j][i]; + for (j = 0; j < m_svd_v.rows(); j++) + basis(j, b) = m_svd_v(j, i); b++; } - TNT::transpose(basis, basis_t); - TNT::matmult(m_null, basis, basis_t); + m_nullspace = basis * basis.transpose(); - for (i = 0; i < m_null.num_rows(); i++) - for (j = 0; j < m_null.num_cols(); j++) + for (i = 0; i < m_nullspace.rows(); i++) + for (j = 0; j < m_nullspace.cols(); j++) if (i == j) - m_null[i][j] = 1.0 - m_null[i][j]; + m_nullspace(i, j) = 1.0 - m_nullspace(i, j); else - m_null[i][j] = -m_null[i][j]; + m_nullspace(i, j) = -m_nullspace(i, j); return true; } @@ -184,7 +177,7 @@ void IK_QJacobian::SubTask(IK_QJacobian& jacobian) return; // restrict lower priority jacobian - jacobian.Restrict(m_d_theta, m_null); + jacobian.Restrict(m_d_theta, m_nullspace); // add angle update from lower priority jacobian.Invert(); @@ -197,19 +190,15 @@ void IK_QJacobian::SubTask(IK_QJacobian& jacobian) m_d_theta[i] = m_d_theta[i] + /*m_min_damp * */ jacobian.AngleUpdate(i); } -void IK_QJacobian::Restrict(TVector& d_theta, TMatrix& null) +void IK_QJacobian::Restrict(VectorXd& d_theta, MatrixXd& nullspace) { // subtract part already moved by higher task from beta - TVector beta_sub(m_beta.size()); - - TNT::matmult(beta_sub, m_jacobian, d_theta); - m_beta = m_beta - beta_sub; + m_beta = m_beta - m_jacobian * d_theta; // note: should we be using the norm of the unrestricted jacobian for SDLS? // project jacobian on to null space of higher priority task - TMatrix jacobian_copy(m_jacobian); - TNT::matmult(m_jacobian, jacobian_copy, null); + m_jacobian = m_jacobian * nullspace; } void IK_QJacobian::InvertSDLS() @@ -234,16 +223,16 @@ void IK_QJacobian::InvertSDLS() double epsilon = 1e-10; int i, j; - m_d_theta = 0; + m_d_theta.setZero(); m_min_damp = 1.0; for (i = 0; i < m_dof; i++) { m_norm[i] = 0.0; for (j = 0; j < m_task_size; j += 3) { double n = 0.0; - n += m_jacobian[j][i] * m_jacobian[j][i]; - n += m_jacobian[j + 1][i] * m_jacobian[j + 1][i]; - n += m_jacobian[j + 2][i] * m_jacobian[j + 2][i]; + n += m_jacobian(j, i) * m_jacobian(j, i); + n += m_jacobian(j + 1, i) * m_jacobian(j + 1, i); + n += m_jacobian(j + 2, i) * m_jacobian(j + 2, i); m_norm[i] += sqrt(n); } } @@ -257,17 +246,17 @@ void IK_QJacobian::InvertSDLS() double N = 0.0; // compute alpha and N - for (j = 0; j < m_svd_u.num_rows(); j += 3) { - alpha += m_svd_u[j][i] * m_beta[j]; - alpha += m_svd_u[j + 1][i] * m_beta[j + 1]; - alpha += m_svd_u[j + 2][i] * m_beta[j + 2]; + for (j = 0; j < m_svd_u.rows(); j += 3) { + alpha += m_svd_u(j, i) * m_beta[j]; + alpha += m_svd_u(j + 1, i) * m_beta[j + 1]; + alpha += m_svd_u(j + 2, i) * m_beta[j + 2]; // note: for 1 end effector, N will always be 1, since U is // orthogonal, .. so could be optimized double tmp; - tmp = m_svd_u[j][i] * m_svd_u[j][i]; - tmp += m_svd_u[j + 1][i] * m_svd_u[j + 1][i]; - tmp += m_svd_u[j + 2][i] * m_svd_u[j + 2][i]; + tmp = m_svd_u(j, i) * m_svd_u(j, i); + tmp += m_svd_u(j + 1, i) * m_svd_u(j + 1, i); + tmp += m_svd_u(j + 2, i) * m_svd_u(j + 2, i); N += sqrt(tmp); } alpha *= wInv; @@ -277,7 +266,7 @@ void IK_QJacobian::InvertSDLS() double max_dtheta = 0.0, abs_dtheta; for (j = 0; j < m_d_theta.size(); j++) { - double v = m_svd_v[j][i]; + double v = m_svd_v(j, i); M += fabs(v) * m_norm[j]; // compute tmporary dTheta's @@ -355,7 +344,7 @@ void IK_QJacobian::InvertDLS() double epsilon = 1e-10; double max_angle_change = 0.1; - double x_length = sqrt(TNT::dot_prod(m_beta, m_beta)); + double x_length = sqrt(m_beta.dot(m_beta)); int i, j; double w_min = std::numeric_limits<double>::max(); @@ -386,10 +375,9 @@ void IK_QJacobian::InvertDLS() // rather than matrix*matrix products // compute Ut*Beta - TNT::transpose(m_svd_u, m_svd_u_t); - TNT::matmult(m_svd_u_beta, m_svd_u_t, m_beta); + m_svd_u_beta = m_svd_u.transpose() * m_beta; - m_d_theta = 0.0; + m_d_theta.setZero(); for (i = 0; i < m_svd_w.size(); i++) { if (m_svd_w[i] > epsilon) { @@ -399,7 +387,7 @@ void IK_QJacobian::InvertDLS() m_svd_u_beta[i] *= wInv; for (j = 0; j < m_d_theta.size(); j++) - m_d_theta[j] += m_svd_v[j][i] * m_svd_u_beta[i]; + m_d_theta[j] += m_svd_v(j, i) * m_svd_u_beta[i]; } } @@ -412,8 +400,8 @@ void IK_QJacobian::Lock(int dof_id, double delta) int i; for (i = 0; i < m_task_size; i++) { - m_beta[i] -= m_jacobian[i][dof_id] * delta; - m_jacobian[i][dof_id] = 0.0; + m_beta[i] -= m_jacobian(i, dof_id) * delta; + m_jacobian(i, dof_id) = 0.0; } m_norm[dof_id] = 0.0; // unneeded diff --git a/intern/iksolver/intern/IK_QJacobian.h b/intern/iksolver/intern/IK_QJacobian.h index 2a5bac49373..f541866c6a7 100644 --- a/intern/iksolver/intern/IK_QJacobian.h +++ b/intern/iksolver/intern/IK_QJacobian.h @@ -33,16 +33,11 @@ #pragma once -#include "TNT/cmat.h" -#include <vector> #include "IK_Math.h" class IK_QJacobian { public: - typedef TNT::Matrix<double> TMatrix; - typedef TNT::Vector<double> TVector; - IK_QJacobian(); ~IK_QJacobian(); @@ -65,7 +60,7 @@ public: // Secondary task bool ComputeNullProjection(); - void Restrict(TVector& d_theta, TMatrix& null); + void Restrict(VectorXd& d_theta, MatrixXd& nullspace); void SubTask(IK_QJacobian& jacobian); private: @@ -77,39 +72,35 @@ private: bool m_transpose; // the jacobian matrix and it's null space projector - TMatrix m_jacobian, m_jacobian_tmp; - TMatrix m_null; + MatrixXd m_jacobian, m_jacobian_tmp; + MatrixXd m_nullspace; /// the vector of intermediate betas - TVector m_beta; + VectorXd m_beta; /// the vector of computed angle changes - TVector m_d_theta; - TVector m_d_norm_weight; + VectorXd m_d_theta; + VectorXd m_d_norm_weight; /// space required for SVD computation + VectorXd m_svd_w; + MatrixXd m_svd_v; + MatrixXd m_svd_u; - TVector m_svd_w; - TMatrix m_svd_v; - TMatrix m_svd_u; - TVector m_work1; - TVector m_work2; - - TMatrix m_svd_u_t; - TVector m_svd_u_beta; + VectorXd m_svd_u_beta; // space required for SDLS bool m_sdls; - TVector m_norm; - TVector m_d_theta_tmp; + VectorXd m_norm; + VectorXd m_d_theta_tmp; double m_min_damp; // null space task vector - TVector m_alpha; + VectorXd m_alpha; // dof weighting - TVector m_weight; - TVector m_weight_sqrt; + VectorXd m_weight; + VectorXd m_weight_sqrt; }; diff --git a/intern/iksolver/intern/TNT/cholesky.h b/intern/iksolver/intern/TNT/cholesky.h deleted file mode 100644 index 89c7dc64d5b..00000000000 --- a/intern/iksolver/intern/TNT/cholesky.h +++ /dev/null @@ -1,98 +0,0 @@ -/* - */ - -/* -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - -#ifndef CHOLESKY_H -#define CHOLESKY_H - -#include <cmath> - -// index method - -namespace TNT -{ - - -// -// Only upper part of A is used. Cholesky factor is returned in -// lower part of L. Returns 0 if successful, 1 otherwise. -// -template <class SPDMatrix, class SymmMatrix> -int Cholesky_upper_factorization(SPDMatrix &A, SymmMatrix &L) -{ - Subscript M = A.dim(1); - Subscript N = A.dim(2); - - assert(M == N); // make sure A is square - - // readjust size of L, if necessary - - if (M != L.dim(1) || N != L.dim(2)) - L = SymmMatrix(N,N); - - Subscript i,j,k; - - - typename SPDMatrix::element_type dot=0; - - - for (j=1; j<=N; j++) // form column j of L - { - dot= 0; - - for (i=1; i<j; i++) // for k= 1 TO j-1 - dot = dot + L(j,i)*L(j,i); - - L(j,j) = A(j,j) - dot; - - for (i=j+1; i<=N; i++) - { - dot = 0; - for (k=1; k<j; k++) - dot = dot + L(i,k)*L(j,k); - L(i,j) = A(j,i) - dot; - } - - if (L(j,j) <= 0.0) return 1; - - L(j,j) = sqrt( L(j,j) ); - - for (i=j+1; i<=N; i++) - L(i,j) = L(i,j) / L(j,j); - - } - - return 0; -} - - - - -} -// namespace TNT - -#endif -// CHOLESKY_H - diff --git a/intern/iksolver/intern/TNT/cmat.h b/intern/iksolver/intern/TNT/cmat.h deleted file mode 100644 index fd3a1851262..00000000000 --- a/intern/iksolver/intern/TNT/cmat.h +++ /dev/null @@ -1,614 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -// C compatible matrix: row-oriented, 0-based [i][j] and 1-based (i,j) indexing -// - -#ifndef CMAT_H -#define CMAT_H - -#include "subscript.h" -#include "vec.h" -#include <stdlib.h> -#include <assert.h> -#include <iostream> -#ifdef TNT_USE_REGIONS -#include "region2d.h" -#endif - -namespace TNT -{ - -template <class T> -class Matrix -{ - - - public: - - typedef Subscript size_type; - typedef T value_type; - typedef T element_type; - typedef T* pointer; - typedef T* iterator; - typedef T& reference; - typedef const T* const_iterator; - typedef const T& const_reference; - - Subscript lbound() const { return 1;} - - protected: - Subscript m_; - Subscript n_; - Subscript mn_; // total size - T* v_; - T** row_; - T* vm1_ ; // these point to the same data, but are 1-based - T** rowm1_; - - // internal helper function to create the array - // of row pointers - - void initialize(Subscript M, Subscript N) - { - mn_ = M*N; - m_ = M; - n_ = N; - - v_ = new T[mn_]; - row_ = new T*[M]; - rowm1_ = new T*[M]; - - assert(v_ != NULL); - assert(row_ != NULL); - assert(rowm1_ != NULL); - - T* p = v_; - vm1_ = v_ - 1; - for (Subscript i=0; i<M; i++) - { - row_[i] = p; - rowm1_[i] = p-1; - p += N ; - - } - - rowm1_ -- ; // compensate for 1-based offset - } - - void copy(const T* v) - { - Subscript N = m_ * n_; - Subscript i; - -#ifdef TNT_UNROLL_LOOPS - Subscript Nmod4 = N & 3; - Subscript N4 = N - Nmod4; - - for (i=0; i<N4; i+=4) - { - v_[i] = v[i]; - v_[i+1] = v[i+1]; - v_[i+2] = v[i+2]; - v_[i+3] = v[i+3]; - } - - for (i=N4; i< N; i++) - v_[i] = v[i]; -#else - - for (i=0; i< N; i++) - v_[i] = v[i]; -#endif - } - - void set(const T& val) - { - Subscript N = m_ * n_; - Subscript i; - -#ifdef TNT_UNROLL_LOOPS - Subscript Nmod4 = N & 3; - Subscript N4 = N - Nmod4; - - for (i=0; i<N4; i+=4) - { - v_[i] = val; - v_[i+1] = val; - v_[i+2] = val; - v_[i+3] = val; - } - - for (i=N4; i< N; i++) - v_[i] = val; -#else - - for (i=0; i< N; i++) - v_[i] = val; - -#endif - } - - - - void destroy() - { - /* do nothing, if no memory has been previously allocated */ - if (v_ == NULL) return ; - - /* if we are here, then matrix was previously allocated */ - if (v_ != NULL) delete [] (v_); - if (row_ != NULL) delete [] (row_); - - /* return rowm1_ back to original value */ - rowm1_ ++; - if (rowm1_ != NULL ) delete [] (rowm1_); - } - - - public: - - operator T**(){ return row_; } - operator T**() const { return row_; } - - - Subscript size() const { return mn_; } - - // constructors - - Matrix() : m_(0), n_(0), mn_(0), v_(0), row_(0), vm1_(0), rowm1_(0) {} - - Matrix(const Matrix<T> &A) - { - initialize(A.m_, A.n_); - copy(A.v_); - } - - Matrix(Subscript M, Subscript N, const T& value = T()) - { - initialize(M,N); - set(value); - } - - Matrix(Subscript M, Subscript N, const T* v) - { - initialize(M,N); - copy(v); - } - - - // destructor - // - ~Matrix() - { - destroy(); - } - - - // reallocating - // - Matrix<T>& newsize(Subscript M, Subscript N) - { - if (num_rows() == M && num_cols() == N) - return *this; - - destroy(); - initialize(M,N); - - return *this; - } - - void - diagonal(Vector<T> &diag) - { - int sz = diag.dim(); - newsize(sz,sz); - set(0); - - Subscript i; - for (i = 0; i < sz; i++) { - row_[i][i] = diag[i]; - } - } - - - - // assignments - // - Matrix<T>& operator=(const Matrix<T> &A) - { - if (v_ == A.v_) - return *this; - - if (m_ == A.m_ && n_ == A.n_) // no need to re-alloc - copy(A.v_); - - else - { - destroy(); - initialize(A.m_, A.n_); - copy(A.v_); - } - - return *this; - } - - Matrix<T>& operator=(const T& scalar) - { - set(scalar); - return *this; - } - - - Subscript dim(Subscript d) const - { -#ifdef TNT_BOUNDS_CHECK - assert( d >= 1); - assert( d <= 2); -#endif - return (d==1) ? m_ : ((d==2) ? n_ : 0); - } - - Subscript num_rows() const { return m_; } - Subscript num_cols() const { return n_; } - - - - - inline T* operator[](Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(0<=i); - assert(i < m_) ; -#endif - return row_[i]; - } - - inline const T* operator[](Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(0<=i); - assert(i < m_) ; -#endif - return row_[i]; - } - - inline reference operator()(Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= mn_) ; -#endif - return vm1_[i]; - } - - inline const_reference operator()(Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= mn_) ; -#endif - return vm1_[i]; - } - - - - inline reference operator()(Subscript i, Subscript j) - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= m_) ; - assert(1<=j); - assert(j <= n_); -#endif - return rowm1_[i][j]; - } - - - - inline const_reference operator() (Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= m_) ; - assert(1<=j); - assert(j <= n_); -#endif - return rowm1_[i][j]; - } - - - -#ifdef TNT_USE_REGIONS - - typedef Region2D<Matrix<T> > Region; - - - Region operator()(const Index1D &I, const Index1D &J) - { - return Region(*this, I,J); - } - - - typedef const_Region2D< Matrix<T> > const_Region; - const_Region operator()(const Index1D &I, const Index1D &J) const - { - return const_Region(*this, I,J); - } - -#endif - - -}; - - -/* *************************** I/O ********************************/ - -template <class T> -std::ostream& operator<<(std::ostream &s, const Matrix<T> &A) -{ - Subscript M=A.num_rows(); - Subscript N=A.num_cols(); - - s << M << " " << N << "\n"; - - for (Subscript i=0; i<M; i++) - { - for (Subscript j=0; j<N; j++) - { - s << A[i][j] << " "; - } - s << "\n"; - } - - - return s; -} - -template <class T> -std::istream& operator>>(std::istream &s, Matrix<T> &A) -{ - - Subscript M, N; - - s >> M >> N; - - if ( !(M == A.num_rows() && N == A.num_cols() )) - { - A.newsize(M,N); - } - - - for (Subscript i=0; i<M; i++) - for (Subscript j=0; j<N; j++) - { - s >> A[i][j]; - } - - - return s; -} - -// *******************[ basic matrix algorithms ]*************************** - -template <class T> -Matrix<T> operator+(const Matrix<T> &A, - const Matrix<T> &B) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(M==B.num_rows()); - assert(N==B.num_cols()); - - Matrix<T> tmp(M,N); - Subscript i,j; - - for (i=0; i<M; i++) - for (j=0; j<N; j++) - tmp[i][j] = A[i][j] + B[i][j]; - - return tmp; -} - -template <class T> -Matrix<T> operator-(const Matrix<T> &A, - const Matrix<T> &B) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(M==B.num_rows()); - assert(N==B.num_cols()); - - Matrix<T> tmp(M,N); - Subscript i,j; - - for (i=0; i<M; i++) - for (j=0; j<N; j++) - tmp[i][j] = A[i][j] - B[i][j]; - - return tmp; -} - -template <class T> -Matrix<T> mult_element(const Matrix<T> &A, - const Matrix<T> &B) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(M==B.num_rows()); - assert(N==B.num_cols()); - - Matrix<T> tmp(M,N); - Subscript i,j; - - for (i=0; i<M; i++) - for (j=0; j<N; j++) - tmp[i][j] = A[i][j] * B[i][j]; - - return tmp; -} - -template <class T> -void transpose(const Matrix<T> &A, Matrix<T> &S) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(M==S.num_cols()); - assert(N==S.num_rows()); - - Subscript i, j; - - for (i=0; i<M; i++) - for (j=0; j<N; j++) - S[j][i] = A[i][j]; - -} - - -template <class T> -inline void matmult(Matrix<T>& C, const Matrix<T> &A, - const Matrix<T> &B) -{ - - assert(A.num_cols() == B.num_rows()); - - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - Subscript K = B.num_cols(); - - C.newsize(M,K); - - T sum; - - const T* row_i; - const T* col_k; - - for (Subscript i=0; i<M; i++) - for (Subscript k=0; k<K; k++) - { - row_i = &(A[i][0]); - col_k = &(B[0][k]); - sum = 0; - for (Subscript j=0; j<N; j++) - { - sum += *row_i * *col_k; - row_i++; - col_k += K; - } - C[i][k] = sum; - } - -} - -template <class T> -void matmult(Vector<T> &y, const Matrix<T> &A, const Vector<T> &x) -{ - -#ifdef TNT_BOUNDS_CHECK - assert(A.num_cols() == x.dim()); - assert(A.num_rows() == y.dim()); -#endif - - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - T sum; - - for (Subscript i=0; i<M; i++) - { - sum = 0; - const T* rowi = A[i]; - for (Subscript j=0; j<N; j++) - sum = sum + rowi[j] * x[j]; - - y[i] = sum; - } -} - -template <class T> -inline void matmultdiag( - Matrix<T>& C, - const Matrix<T> &A, - const Vector<T> &diag -){ -#ifdef TNT_BOUNDS_CHECK - assert(A.num_cols() ==A.num_rows()== diag.dim()); -#endif - - Subscript M = A.num_rows(); - Subscript K = diag.dim(); - - C.newsize(M,K); - - const T* row_i; - const T* col_k; - - for (Subscript i=0; i<M; i++) { - for (Subscript k=0; k<K; k++) - { - C[i][k] = A[i][k] * diag[k]; - } - } -} - - -template <class T> -inline void matmultdiag( - Matrix<T>& C, - const Vector<T> &diag, - const Matrix<T> &A -){ -#ifdef TNT_BOUNDS_CHECK - assert(A.num_cols() ==A.num_rows()== diag.dim()); -#endif - - Subscript M = A.num_rows(); - Subscript K = diag.dim(); - - C.newsize(M,K); - - for (Subscript i=0; i<M; i++) { - - const T diag_element = diag[i]; - - for (Subscript k=0; k<K; k++) - { - C[i][k] = A[i][k] * diag_element; - } - } -} - -} // namespace TNT - -#endif // CMAT_H - diff --git a/intern/iksolver/intern/TNT/fcscmat.h b/intern/iksolver/intern/TNT/fcscmat.h deleted file mode 100644 index 8865dc7a039..00000000000 --- a/intern/iksolver/intern/TNT/fcscmat.h +++ /dev/null @@ -1,167 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -// Templated compressed sparse column matrix (Fortran conventions). -// uses 1-based offsets in storing row indices. -// Used primarily to interface with Fortran sparse matrix libaries. -// (CANNOT BE USED AS AN STL CONTAINER.) - - -#ifndef FCSCMAT_H -#define FCSCMAT_H - -#include <iostream> -#include <cassert> -#include "tnt.h" -#include "vec.h" - -using namespace std; - -namespace TNT -{ - -template <class T> -class Fortran_Sparse_Col_Matrix -{ - - protected: - - Vector<T> val_; // data values (nz_ elements) - Vector<Subscript> rowind_; // row_ind (nz_ elements) - Vector<Subscript> colptr_; // col_ptr (n_+1 elements) - - int nz_; // number of nonzeros - Subscript m_; // global dimensions - Subscript n_; - - public: - - - Fortran_Sparse_Col_Matrix(void); - Fortran_Sparse_Col_Matrix(const Fortran_Sparse_Col_Matrix<T> &S) - : val_(S.val_), rowind_(S.rowind_), colptr_(S.colptr_), nz_(S.nz_), - m_(S.m_), n_(S.n_) {}; - Fortran_Sparse_Col_Matrix(Subscript M, Subscript N, - Subscript nz, const T *val, const Subscript *r, - const Subscript *c) : val_(nz, val), rowind_(nz, r), - colptr_(N+1, c), nz_(nz), m_(M), n_(N) {}; - - Fortran_Sparse_Col_Matrix(Subscript M, Subscript N, - Subscript nz, char *val, char *r, - char *c) : val_(nz, val), rowind_(nz, r), - colptr_(N+1, c), nz_(nz), m_(M), n_(N) {}; - - Fortran_Sparse_Col_Matrix(Subscript M, Subscript N, - Subscript nz, const T *val, Subscript *r, Subscript *c) - : val_(nz, val), rowind_(nz, r), colptr_(N+1, c), nz_(nz), - m_(M), n_(N) {}; - - ~Fortran_Sparse_Col_Matrix() {}; - - - T & val(Subscript i) { return val_(i); } - const T & val(Subscript i) const { return val_(i); } - - Subscript & row_ind(Subscript i) { return rowind_(i); } - const Subscript & row_ind(Subscript i) const { return rowind_(i); } - - Subscript col_ptr(Subscript i) { return colptr_(i);} - const Subscript col_ptr(Subscript i) const { return colptr_(i);} - - - Subscript num_cols() const { return m_;} - Subscript num_rows() const { return n_; } - - Subscript dim(Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert( 1 <= i ); - assert( i <= 2 ); -#endif - if (i==1) return m_; - else if (i==2) return m_; - else return 0; - } - - Subscript num_nonzeros() const {return nz_;}; - Subscript lbound() const {return 1;} - - - - Fortran_Sparse_Col_Matrix& operator=(const - Fortran_Sparse_Col_Matrix &C) - { - val_ = C.val_; - rowind_ = C.rowind_; - colptr_ = C.colptr_; - nz_ = C.nz_; - m_ = C.m_; - n_ = C.n_; - - return *this; - } - - Fortran_Sparse_Col_Matrix& newsize(Subscript M, Subscript N, - Subscript nz) - { - val_.newsize(nz); - rowind_.newsize(nz); - colptr_.newsize(N+1); - return *this; - } -}; - -template <class T> -ostream& operator<<(ostream &s, const Fortran_Sparse_Col_Matrix<T> &A) -{ - Subscript M=A.num_rows(); - Subscript N=A.num_cols(); - - s << M << " " << N << " " << A.num_nonzeros() << endl; - - - for (Subscript k=1; k<=N; k++) - { - Subscript start = A.col_ptr(k); - Subscript end = A.col_ptr(k+1); - - for (Subscript i= start; i<end; i++) - { - s << A.row_ind(i) << " " << k << " " << A.val(i) << endl; - } - } - - return s; -} - - -} // namespace TNT - -#endif /* FCSCMAT_H */ - diff --git a/intern/iksolver/intern/TNT/fmat.h b/intern/iksolver/intern/TNT/fmat.h deleted file mode 100644 index edb64003143..00000000000 --- a/intern/iksolver/intern/TNT/fmat.h +++ /dev/null @@ -1,569 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -// Fortran-compatible matrix: column oriented, 1-based (i,j) indexing - -#ifndef FMAT_H -#define FMAT_H - -#include "subscript.h" -#include "vec.h" -#include <cstdlib> -#include <cassert> -#include <iostream> -#ifdef TNT_USE_REGIONS -#include "region2d.h" -#endif - -// simple 1-based, column oriented Matrix class - -namespace TNT -{ - -template <class T> -class Fortran_Matrix -{ - - - public: - - typedef T value_type; - typedef T element_type; - typedef T* pointer; - typedef T* iterator; - typedef T& reference; - typedef const T* const_iterator; - typedef const T& const_reference; - - Subscript lbound() const { return 1;} - - protected: - T* v_; // these are adjusted to simulate 1-offset - Subscript m_; - Subscript n_; - T** col_; // these are adjusted to simulate 1-offset - - // internal helper function to create the array - // of row pointers - - void initialize(Subscript M, Subscript N) - { - // adjust col_[] pointers so that they are 1-offset: - // col_[j][i] is really col_[j-1][i-1]; - // - // v_[] is the internal contiguous array, it is still 0-offset - // - v_ = new T[M*N]; - col_ = new T*[N]; - - assert(v_ != NULL); - assert(col_ != NULL); - - - m_ = M; - n_ = N; - T* p = v_ - 1; - for (Subscript i=0; i<N; i++) - { - col_[i] = p; - p += M ; - - } - col_ --; - } - - void copy(const T* v) - { - Subscript N = m_ * n_; - Subscript i; - -#ifdef TNT_UNROLL_LOOPS - Subscript Nmod4 = N & 3; - Subscript N4 = N - Nmod4; - - for (i=0; i<N4; i+=4) - { - v_[i] = v[i]; - v_[i+1] = v[i+1]; - v_[i+2] = v[i+2]; - v_[i+3] = v[i+3]; - } - - for (i=N4; i< N; i++) - v_[i] = v[i]; -#else - - for (i=0; i< N; i++) - v_[i] = v[i]; -#endif - } - - void set(const T& val) - { - Subscript N = m_ * n_; - Subscript i; - -#ifdef TNT_UNROLL_LOOPS - Subscript Nmod4 = N & 3; - Subscript N4 = N - Nmod4; - - for (i=0; i<N4; i+=4) - { - v_[i] = val; - v_[i+1] = val; - v_[i+2] = val; - v_[i+3] = val; - } - - for (i=N4; i< N; i++) - v_[i] = val; -#else - - for (i=0; i< N; i++) - v_[i] = val; - -#endif - } - - - - void destroy() - { - /* do nothing, if no memory has been previously allocated */ - if (v_ == NULL) return ; - - /* if we are here, then matrix was previously allocated */ - delete [] (v_); - col_ ++; // changed back to 0-offset - delete [] (col_); - } - - - public: - - T* begin() { return v_; } - const T* begin() const { return v_;} - - T* end() { return v_ + m_*n_; } - const T* end() const { return v_ + m_*n_; } - - - // constructors - - Fortran_Matrix() : v_(0), m_(0), n_(0), col_(0) {}; - Fortran_Matrix(const Fortran_Matrix<T> &A) - { - initialize(A.m_, A.n_); - copy(A.v_); - } - - Fortran_Matrix(Subscript M, Subscript N, const T& value = T()) - { - initialize(M,N); - set(value); - } - - Fortran_Matrix(Subscript M, Subscript N, const T* v) - { - initialize(M,N); - copy(v); - } - - - // destructor - ~Fortran_Matrix() - { - destroy(); - } - - - // assignments - // - Fortran_Matrix<T>& operator=(const Fortran_Matrix<T> &A) - { - if (v_ == A.v_) - return *this; - - if (m_ == A.m_ && n_ == A.n_) // no need to re-alloc - copy(A.v_); - - else - { - destroy(); - initialize(A.m_, A.n_); - copy(A.v_); - } - - return *this; - } - - Fortran_Matrix<T>& operator=(const T& scalar) - { - set(scalar); - return *this; - } - - - Subscript dim(Subscript d) const - { -#ifdef TNT_BOUNDS_CHECK - assert( d >= 1); - assert( d <= 2); -#endif - return (d==1) ? m_ : ((d==2) ? n_ : 0); - } - - Subscript num_rows() const { return m_; } - Subscript num_cols() const { return n_; } - - Fortran_Matrix<T>& newsize(Subscript M, Subscript N) - { - if (num_rows() == M && num_cols() == N) - return *this; - - destroy(); - initialize(M,N); - - return *this; - } - - - - // 1-based element access - // - inline reference operator()(Subscript i, Subscript j) - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= m_) ; - assert(1<=j); - assert(j <= n_); -#endif - return col_[j][i]; - } - - inline const_reference operator() (Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= m_) ; - assert(1<=j); - assert(j <= n_); -#endif - return col_[j][i]; - } - - -#ifdef TNT_USE_REGIONS - - typedef Region2D<Fortran_Matrix<T> > Region; - typedef const_Region2D< Fortran_Matrix<T> > const_Region; - - Region operator()(const Index1D &I, const Index1D &J) - { - return Region(*this, I,J); - } - - const_Region operator()(const Index1D &I, const Index1D &J) const - { - return const_Region(*this, I,J); - } - -#endif - - -}; - - -/* *************************** I/O ********************************/ - -template <class T> -std::ostream& operator<<(std::ostream &s, const Fortran_Matrix<T> &A) -{ - Subscript M=A.num_rows(); - Subscript N=A.num_cols(); - - s << M << " " << N << "\n"; - - for (Subscript i=1; i<=M; i++) - { - for (Subscript j=1; j<=N; j++) - { - s << A(i,j) << " "; - } - s << "\n"; - } - - - return s; -} - -template <class T> -std::istream& operator>>(std::istream &s, Fortran_Matrix<T> &A) -{ - - Subscript M, N; - - s >> M >> N; - - if ( !(M == A.num_rows() && N == A.num_cols())) - { - A.newsize(M,N); - } - - - for (Subscript i=1; i<=M; i++) - for (Subscript j=1; j<=N; j++) - { - s >> A(i,j); - } - - - return s; -} - -// *******************[ basic matrix algorithms ]*************************** - - -template <class T> -Fortran_Matrix<T> operator+(const Fortran_Matrix<T> &A, - const Fortran_Matrix<T> &B) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(M==B.num_rows()); - assert(N==B.num_cols()); - - Fortran_Matrix<T> tmp(M,N); - Subscript i,j; - - for (i=1; i<=M; i++) - for (j=1; j<=N; j++) - tmp(i,j) = A(i,j) + B(i,j); - - return tmp; -} - -template <class T> -Fortran_Matrix<T> operator-(const Fortran_Matrix<T> &A, - const Fortran_Matrix<T> &B) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(M==B.num_rows()); - assert(N==B.num_cols()); - - Fortran_Matrix<T> tmp(M,N); - Subscript i,j; - - for (i=1; i<=M; i++) - for (j=1; j<=N; j++) - tmp(i,j) = A(i,j) - B(i,j); - - return tmp; -} - -// element-wise multiplication (use matmult() below for matrix -// multiplication in the linear algebra sense.) -// -// -template <class T> -Fortran_Matrix<T> mult_element(const Fortran_Matrix<T> &A, - const Fortran_Matrix<T> &B) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(M==B.num_rows()); - assert(N==B.num_cols()); - - Fortran_Matrix<T> tmp(M,N); - Subscript i,j; - - for (i=1; i<=M; i++) - for (j=1; j<=N; j++) - tmp(i,j) = A(i,j) * B(i,j); - - return tmp; -} - - -template <class T> -Fortran_Matrix<T> transpose(const Fortran_Matrix<T> &A) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - Fortran_Matrix<T> S(N,M); - Subscript i, j; - - for (i=1; i<=M; i++) - for (j=1; j<=N; j++) - S(j,i) = A(i,j); - - return S; -} - - - -template <class T> -inline Fortran_Matrix<T> matmult(const Fortran_Matrix<T> &A, - const Fortran_Matrix<T> &B) -{ - -#ifdef TNT_BOUNDS_CHECK - assert(A.num_cols() == B.num_rows()); -#endif - - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - Subscript K = B.num_cols(); - - Fortran_Matrix<T> tmp(M,K); - T sum; - - for (Subscript i=1; i<=M; i++) - for (Subscript k=1; k<=K; k++) - { - sum = 0; - for (Subscript j=1; j<=N; j++) - sum = sum + A(i,j) * B(j,k); - - tmp(i,k) = sum; - } - - return tmp; -} - -template <class T> -inline Fortran_Matrix<T> operator*(const Fortran_Matrix<T> &A, - const Fortran_Matrix<T> &B) -{ - return matmult(A,B); -} - -template <class T> -inline int matmult(Fortran_Matrix<T>& C, const Fortran_Matrix<T> &A, - const Fortran_Matrix<T> &B) -{ - - assert(A.num_cols() == B.num_rows()); - - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - Subscript K = B.num_cols(); - - C.newsize(M,K); // adjust shape of C, if necessary - - - T sum; - - const T* row_i; - const T* col_k; - - for (Subscript i=1; i<=M; i++) - { - for (Subscript k=1; k<=K; k++) - { - row_i = &A(i,1); - col_k = &B(1,k); - sum = 0; - for (Subscript j=1; j<=N; j++) - { - sum += *row_i * *col_k; - row_i += M; - col_k ++; - } - - C(i,k) = sum; - } - - } - - return 0; -} - - -template <class T> -Vector<T> matmult(const Fortran_Matrix<T> &A, const Vector<T> &x) -{ - -#ifdef TNT_BOUNDS_CHECK - assert(A.num_cols() == x.dim()); -#endif - - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - Vector<T> tmp(M); - T sum; - - for (Subscript i=1; i<=M; i++) - { - sum = 0; - for (Subscript j=1; j<=N; j++) - sum = sum + A(i,j) * x(j); - - tmp(i) = sum; - } - - return tmp; -} - -template <class T> -inline Vector<T> operator*(const Fortran_Matrix<T> &A, const Vector<T> &x) -{ - return matmult(A,x); -} - -template <class T> -inline Fortran_Matrix<T> operator*(const Fortran_Matrix<T> &A, const T &x) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - Subscript MN = M*N; - - Fortran_Matrix<T> res(M,N); - const T* a = A.begin(); - T* t = res.begin(); - T* tend = res.end(); - - for (t=res.begin(); t < tend; t++, a++) - *t = *a * x; - - return res; -} - -} // namespace TNT - -#endif // FMAT_H - diff --git a/intern/iksolver/intern/TNT/fortran.h b/intern/iksolver/intern/TNT/fortran.h deleted file mode 100644 index c58a859741b..00000000000 --- a/intern/iksolver/intern/TNT/fortran.h +++ /dev/null @@ -1,69 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -// Header file to define C/Fortran conventions (Platform specific) - -#ifndef FORTRAN_H -#define FORTRAN_H - -// help map between C/C++ data types and Fortran types - -typedef int Fortran_integer; -typedef float Fortran_float; -typedef double Fortran_double; - - -typedef Fortran_double *fda_; // (in/out) double precision array -typedef const Fortran_double *cfda_; // (in) double precsion array - -typedef Fortran_double *fd_; // (in/out) single double precision -typedef const Fortran_double *cfd_; // (in) single double precision - -typedef Fortran_float *ffa_; // (in/out) float precision array -typedef const Fortran_float *cffa_; // (in) float precsion array - -typedef Fortran_float *ff_; // (in/out) single float precision -typedef const Fortran_float *cff_; // (in) single float precision - -typedef Fortran_integer *fia_; // (in/out) single integer array -typedef const Fortran_integer *cfia_; // (in) single integer array - -typedef Fortran_integer *fi_; // (in/out) single integer -typedef const Fortran_integer *cfi_; // (in) single integer - -typedef char *fch_; // (in/out) single character -typedef char *cfch_; // (in) single character - - -#ifndef TNT_SUBSCRIPT_TYPE -#define TNT_SUBSCRIPT_TYPE TNT::Fortran_integer -#endif - -#endif // FORTRAN_H - diff --git a/intern/iksolver/intern/TNT/fspvec.h b/intern/iksolver/intern/TNT/fspvec.h deleted file mode 100644 index 8b14fa89216..00000000000 --- a/intern/iksolver/intern/TNT/fspvec.h +++ /dev/null @@ -1,171 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - -// Templated sparse vector (Fortran conventions). -// Used primarily to interface with Fortran sparse matrix libaries. -// (CANNOT BE USED AS AN STL CONTAINER.) - -#ifndef FSPVEC_H -#define FSPVEC_H - -#include "tnt.h" -#include "vec.h" -#include <cstdlib> -#include <cassert> -#include <iostream> - -using namespace std; - -namespace TNT -{ - -template <class T> -class Fortran_Sparse_Vector -{ - - - public: - - typedef Subscript size_type; - typedef T value_type; - typedef T element_type; - typedef T* pointer; - typedef T* iterator; - typedef T& reference; - typedef const T* const_iterator; - typedef const T& const_reference; - - Subscript lbound() const { return 1;} - - protected: - Vector<T> val_; - Vector<Subscript> index_; - Subscript dim_; // prescribed dimension - - - public: - - // size and shape information - - Subscript dim() const { return dim_; } - Subscript num_nonzeros() const { return val_.dim(); } - - // access - - T& val(Subscript i) { return val_(i); } - const T& val(Subscript i) const { return val_(i); } - - Subscript &index(Subscript i) { return index_(i); } - const Subscript &index(Subscript i) const { return index_(i); } - - // constructors - - Fortran_Sparse_Vector() : val_(), index_(), dim_(0) {}; - Fortran_Sparse_Vector(Subscript N, Subscript nz) : val_(nz), - index_(nz), dim_(N) {}; - Fortran_Sparse_Vector(Subscript N, Subscript nz, const T *values, - const Subscript *indices): val_(nz, values), index_(nz, indices), - dim_(N) {} - - Fortran_Sparse_Vector(const Fortran_Sparse_Vector<T> &S): - val_(S.val_), index_(S.index_), dim_(S.dim_) {} - - // initialize from string, e.g. - // - // Fortran_Sparse_Vector<T> A(N, 2, "1.0 2.1", "1 3"); - // - Fortran_Sparse_Vector(Subscript N, Subscript nz, char *v, - char *ind) : val_(nz, v), index_(nz, ind), dim_(N) {} - - // assignments - - Fortran_Sparse_Vector<T> & newsize(Subscript N, Subscript nz) - { - val_.newsize(nz); - index_.newsize(nz); - dim_ = N; - return *this; - } - - Fortran_Sparse_Vector<T> & operator=( const Fortran_Sparse_Vector<T> &A) - { - val_ = A.val_; - index_ = A.index_; - dim_ = A.dim_; - - return *this; - } - - // methods - - - -}; - - -/* *************************** I/O ********************************/ - -template <class T> -ostream& operator<<(ostream &s, const Fortran_Sparse_Vector<T> &A) -{ - // output format is : N nz val1 ind1 val2 ind2 ... - Subscript nz=A.num_nonzeros(); - - s << A.dim() << " " << nz << endl; - - for (Subscript i=1; i<=nz; i++) - s << A.val(i) << " " << A.index(i) << endl; - s << endl; - - return s; -} - - -template <class T> -istream& operator>>(istream &s, Fortran_Sparse_Vector<T> &A) -{ - // output format is : N nz val1 ind1 val2 ind2 ... - - Subscript N; - Subscript nz; - - s >> N >> nz; - - A.newsize(N, nz); - - for (Subscript i=1; i<=nz; i++) - s >> A.val(i) >> A.index(i); - - - return s; -} - -} // namespace TNT - -#endif // FSPVEC_H - diff --git a/intern/iksolver/intern/TNT/index.h b/intern/iksolver/intern/TNT/index.h deleted file mode 100644 index 1abe20ba729..00000000000 --- a/intern/iksolver/intern/TNT/index.h +++ /dev/null @@ -1,87 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -// Vector/Matrix/Array Index Module - -#ifndef INDEX_H -#define INDEX_H - -#include "subscript.h" - -namespace TNT -{ - -class Index1D -{ - Subscript lbound_; - Subscript ubound_; - - public: - - Subscript lbound() const { return lbound_; } - Subscript ubound() const { return ubound_; } - - Index1D(const Index1D &D) : lbound_(D.lbound_), ubound_(D.ubound_) {} - Index1D(Subscript i1, Subscript i2) : lbound_(i1), ubound_(i2) {} - - Index1D & operator=(const Index1D &D) - { - lbound_ = D.lbound_; - ubound_ = D.ubound_; - return *this; - } - -}; - -inline Index1D operator+(const Index1D &D, Subscript i) -{ - return Index1D(i+D.lbound(), i+D.ubound()); -} - -inline Index1D operator+(Subscript i, const Index1D &D) -{ - return Index1D(i+D.lbound(), i+D.ubound()); -} - - - -inline Index1D operator-(Index1D &D, Subscript i) -{ - return Index1D(D.lbound()-i, D.ubound()-i); -} - -inline Index1D operator-(Subscript i, Index1D &D) -{ - return Index1D(i-D.lbound(), i-D.ubound()); -} - -} // namespace TNT - -#endif - diff --git a/intern/iksolver/intern/TNT/lapack.h b/intern/iksolver/intern/TNT/lapack.h deleted file mode 100644 index e015fa9b617..00000000000 --- a/intern/iksolver/intern/TNT/lapack.h +++ /dev/null @@ -1,189 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -// Header file for Fortran Lapack - -#ifndef LAPACK_H -#define LAPACK_H - -// This file incomplete and included here to only demonstrate the -// basic framework for linking with the Fortran Lapack routines. - -#include "fortran.h" -#include "vec.h" -#include "fmat.h" - - -#define F77_DGESV dgesv_ -#define F77_DGELS dgels_ -#define F77_DSYEV dsyev_ -#define F77_DGEEV dgeev_ - -extern "C" -{ - - // linear equations (general) using LU factorizaiton - // - void F77_DGESV(cfi_ N, cfi_ nrhs, fda_ A, cfi_ lda, - fia_ ipiv, fda_ b, cfi_ ldb, fi_ info); - - // solve linear least squares using QR or LU factorization - // - void F77_DGELS(cfch_ trans, cfi_ M, - cfi_ N, cfi_ nrhs, fda_ A, cfi_ lda, fda_ B, cfi_ ldb, fda_ work, - cfi_ lwork, fi_ info); - - // solve symmetric eigenvalues - // - void F77_DSYEV( cfch_ jobz, cfch_ uplo, cfi_ N, fda_ A, cfi_ lda, - fda_ W, fda_ work, cfi_ lwork, fi_ info); - - // solve unsymmetric eigenvalues - // - void F77_DGEEV(cfch_ jobvl, cfch_ jobvr, cfi_ N, fda_ A, cfi_ lda, - fda_ wr, fda_ wi, fda_ vl, cfi_ ldvl, fda_ vr, - cfi_ ldvr, fda_ work, cfi_ lwork, fi_ info); - -} - -// solve linear equations using LU factorization - -using namespace TNT; - -Vector<double> Lapack_LU_linear_solve(const Fortran_Matrix<double> &A, - const Vector<double> &b) -{ - const Fortran_integer one=1; - Subscript M=A.num_rows(); - Subscript N=A.num_cols(); - - Fortran_Matrix<double> Tmp(A); - Vector<double> x(b); - Vector<Fortran_integer> index(M); - Fortran_integer info = 0; - - F77_DGESV(&N, &one, &Tmp(1,1), &M, &index(1), &x(1), &M, &info); - - if (info != 0) return Vector<double>(0); - else - return x; -} - -// solve linear least squares problem using QR factorization -// -Vector<double> Lapack_LLS_QR_linear_solve(const Fortran_Matrix<double> &A, - const Vector<double> &b) -{ - const Fortran_integer one=1; - Subscript M=A.num_rows(); - Subscript N=A.num_cols(); - - Fortran_Matrix<double> Tmp(A); - Vector<double> x(b); - Fortran_integer info = 0; - - char transp = 'N'; - Fortran_integer lwork = 5 * (M+N); // temporary work space - Vector<double> work(lwork); - - F77_DGELS(&transp, &M, &N, &one, &Tmp(1,1), &M, &x(1), &M, &work(1), - &lwork, &info); - - if (info != 0) return Vector<double>(0); - else - return x; -} - -// *********************** Eigenvalue problems ******************* - -// solve symmetric eigenvalue problem (eigenvalues only) -// -Vector<double> Upper_symmetric_eigenvalue_solve(const Fortran_Matrix<double> &A) -{ - char jobz = 'N'; - char uplo = 'U'; - Subscript N = A.num_rows(); - - assert(N == A.num_cols()); - - Vector<double> eigvals(N); - Fortran_integer worksize = 3*N; - Fortran_integer info = 0; - Vector<double> work(worksize); - Fortran_Matrix<double> Tmp = A; - - F77_DSYEV(&jobz, &uplo, &N, &Tmp(1,1), &N, eigvals.begin(), work.begin(), - &worksize, &info); - - if (info != 0) return Vector<double>(); - else - return eigvals; -} - - -// solve unsymmetric eigenvalue problems -// -int eigenvalue_solve(const Fortran_Matrix<double> &A, - Vector<double> &wr, Vector<double> &wi) -{ - char jobvl = 'N'; - char jobvr = 'N'; - - Fortran_integer N = A.num_rows(); - - - assert(N == A.num_cols()); - - if (N<1) return 1; - - Fortran_Matrix<double> vl(1,N); /* should be NxN ? **** */ - Fortran_Matrix<double> vr(1,N); - Fortran_integer one = 1; - - Fortran_integer worksize = 5*N; - Fortran_integer info = 0; - Vector<double> work(worksize, 0.0); - Fortran_Matrix<double> Tmp = A; - - wr.newsize(N); - wi.newsize(N); - -// void F77_DGEEV(cfch_ jobvl, cfch_ jobvr, cfi_ N, fda_ A, cfi_ lda, -// fda_ wr, fda_ wi, fda_ vl, cfi_ ldvl, fda_ vr, -// cfi_ ldvr, fda_ work, cfi_ lwork, fi_ info); - - F77_DGEEV(&jobvl, &jobvr, &N, &Tmp(1,1), &N, &(wr(1)), - &(wi(1)), &(vl(1,1)), &one, &(vr(1,1)), &one, - &(work(1)), &worksize, &info); - - return (info==0 ? 0: 1); -} - -#endif // LAPACK_H - diff --git a/intern/iksolver/intern/TNT/lu.h b/intern/iksolver/intern/TNT/lu.h deleted file mode 100644 index f64172bd31f..00000000000 --- a/intern/iksolver/intern/TNT/lu.h +++ /dev/null @@ -1,208 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -#ifndef LU_H -#define LU_H - -// Solve system of linear equations Ax = b. -// -// Typical usage: -// -// Matrix(double) A; -// Vector(Subscript) ipiv; -// Vector(double) b; -// -// 1) LU_Factor(A,ipiv); -// 2) LU_Solve(A,ipiv,b); -// -// Now b has the solution x. Note that both A and b -// are overwritten. If these values need to be preserved, -// one can make temporary copies, as in -// -// O) Matrix(double) T = A; -// 1) LU_Factor(T,ipiv); -// 1a) Vector(double) x=b; -// 2) LU_Solve(T,ipiv,x); -// -// See details below. -// - - -// for fabs() -// -#include <cmath> - -// right-looking LU factorization algorithm (unblocked) -// -// Factors matrix A into lower and upper triangular matrices -// (L and U respectively) in solving the linear equation Ax=b. -// -// -// Args: -// -// A (input/output) Matrix(1:n, 1:n) In input, matrix to be -// factored. On output, overwritten with lower and -// upper triangular factors. -// -// indx (output) Vector(1:n) Pivot vector. Describes how -// the rows of A were reordered to increase -// numerical stability. -// -// Return value: -// -// int (0 if successful, 1 otherwise) -// -// - - -namespace TNT -{ - -template <class MaTRiX, class VecToRSubscript> -int LU_factor( MaTRiX &A, VecToRSubscript &indx) -{ - assert(A.lbound() == 1); // currently for 1-offset - assert(indx.lbound() == 1); // vectors and matrices - - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - if (M == 0 || N==0) return 0; - if (indx.dim() != M) - indx.newsize(M); - - Subscript i=0,j=0,k=0; - Subscript jp=0; - - typename MaTRiX::element_type t; - - Subscript minMN = (M < N ? M : N) ; // min(M,N); - - for (j=1; j<= minMN; j++) - { - - // find pivot in column j and test for singularity. - - jp = j; - t = fabs(A(j,j)); - for (i=j+1; i<=M; i++) - if ( fabs(A(i,j)) > t) - { - jp = i; - t = fabs(A(i,j)); - } - - indx(j) = jp; - - // jp now has the index of maximum element - // of column j, below the diagonal - - if ( A(jp,j) == 0 ) - return 1; // factorization failed because of zero pivot - - - if (jp != j) // swap rows j and jp - for (k=1; k<=N; k++) - { - t = A(j,k); - A(j,k) = A(jp,k); - A(jp,k) =t; - } - - if (j<M) // compute elements j+1:M of jth column - { - // note A(j,j), was A(jp,p) previously which was - // guarranteed not to be zero (Label #1) - // - typename MaTRiX::element_type recp = 1.0 / A(j,j); - - for (k=j+1; k<=M; k++) - A(k,j) *= recp; - } - - - if (j < minMN) - { - // rank-1 update to trailing submatrix: E = E - x*y; - // - // E is the region A(j+1:M, j+1:N) - // x is the column vector A(j+1:M,j) - // y is row vector A(j,j+1:N) - - Subscript ii,jj; - - for (ii=j+1; ii<=M; ii++) - for (jj=j+1; jj<=N; jj++) - A(ii,jj) -= A(ii,j)*A(j,jj); - } - } - - return 0; -} - - - - -template <class MaTRiX, class VecToR, class VecToRSubscripts> -int LU_solve(const MaTRiX &A, const VecToRSubscripts &indx, VecToR &b) -{ - assert(A.lbound() == 1); // currently for 1-offset - assert(indx.lbound() == 1); // vectors and matrices - assert(b.lbound() == 1); - - Subscript i,ii=0,ip,j; - Subscript n = b.dim(); - typename MaTRiX::element_type sum = 0.0; - - for (i=1;i<=n;i++) - { - ip=indx(i); - sum=b(ip); - b(ip)=b(i); - if (ii) - for (j=ii;j<=i-1;j++) - sum -= A(i,j)*b(j); - else if (sum) ii=i; - b(i)=sum; - } - for (i=n;i>=1;i--) - { - sum=b(i); - for (j=i+1;j<=n;j++) - sum -= A(i,j)*b(j); - b(i)=sum/A(i,i); - } - - return 0; -} - -} // namespace TNT - -#endif // LU_H - diff --git a/intern/iksolver/intern/TNT/qr.h b/intern/iksolver/intern/TNT/qr.h deleted file mode 100644 index 9df34e2d7cc..00000000000 --- a/intern/iksolver/intern/TNT/qr.h +++ /dev/null @@ -1,233 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - -#ifndef QR_H -#define QR_H - -// Classical QR factorization example, based on Stewart[1973]. -// -// -// This algorithm computes the factorization of a matrix A -// into a product of an orthognal matrix (Q) and an upper triangular -// matrix (R), such that QR = A. -// -// Parameters: -// -// A (in): Matrix(1:N, 1:N) -// -// Q (output): Matrix(1:N, 1:N), collection of Householder -// column vectors Q1, Q2, ... QN -// -// R (output): upper triangular Matrix(1:N, 1:N) -// -// Returns: -// -// 0 if successful, 1 if A is detected to be singular -// - - -#include <cmath> //for sqrt() & fabs() -#include "tntmath.h" // for sign() - -// Classical QR factorization, based on Stewart[1973]. -// -// -// This algorithm computes the factorization of a matrix A -// into a product of an orthognal matrix (Q) and an upper triangular -// matrix (R), such that QR = A. -// -// Parameters: -// -// A (in/out): On input, A is square, Matrix(1:N, 1:N), that represents -// the matrix to be factored. -// -// On output, Q and R is encoded in the same Matrix(1:N,1:N) -// in the following manner: -// -// R is contained in the upper triangular section of A, -// except that R's main diagonal is in D. The lower -// triangular section of A represents Q, where each -// column j is the vector Qj = I - uj*uj'/pi_j. -// -// C (output): vector of Pi[j] -// D (output): main diagonal of R, i.e. D(i) is R(i,i) -// -// Returns: -// -// 0 if successful, 1 if A is detected to be singular -// - -namespace TNT -{ - -template <class MaTRiX, class Vector> -int QR_factor(MaTRiX &A, Vector& C, Vector &D) -{ - assert(A.lbound() == 1); // ensure these are all - assert(C.lbound() == 1); // 1-based arrays and vectors - assert(D.lbound() == 1); - - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(M == N); // make sure A is square - - Subscript i,j,k; - typename MaTRiX::element_type eta, sigma, sum; - - // adjust the shape of C and D, if needed... - - if (N != C.size()) C.newsize(N); - if (N != D.size()) D.newsize(N); - - for (k=1; k<N; k++) - { - // eta = max |M(i,k)|, for k <= i <= n - // - eta = 0; - for (i=k; i<=N; i++) - { - double absA = fabs(A(i,k)); - eta = ( absA > eta ? absA : eta ); - } - - if (eta == 0) // matrix is singular - { - cerr << "QR: k=" << k << "\n"; - return 1; - } - - // form Qk and premiltiply M by it - // - for(i=k; i<=N; i++) - A(i,k) = A(i,k) / eta; - - sum = 0; - for (i=k; i<=N; i++) - sum = sum + A(i,k)*A(i,k); - sigma = sign(A(k,k)) * sqrt(sum); - - - A(k,k) = A(k,k) + sigma; - C(k) = sigma * A(k,k); - D(k) = -eta * sigma; - - for (j=k+1; j<=N; j++) - { - sum = 0; - for (i=k; i<=N; i++) - sum = sum + A(i,k)*A(i,j); - sum = sum / C(k); - - for (i=k; i<=N; i++) - A(i,j) = A(i,j) - sum * A(i,k); - } - - D(N) = A(N,N); - } - - return 0; -} - -// modified form of upper triangular solve, except that the main diagonal -// of R (upper portion of A) is in D. -// -template <class MaTRiX, class Vector> -int R_solve(const MaTRiX &A, /*const*/ Vector &D, Vector &b) -{ - assert(A.lbound() == 1); // ensure these are all - assert(D.lbound() == 1); // 1-based arrays and vectors - assert(b.lbound() == 1); - - Subscript i,j; - Subscript N = A.num_rows(); - - assert(N == A.num_cols()); - assert(N == D.dim()); - assert(N == b.dim()); - - typename MaTRiX::element_type sum; - - if (D(N) == 0) - return 1; - - b(N) = b(N) / - D(N); - - for (i=N-1; i>=1; i--) - { - if (D(i) == 0) - return 1; - sum = 0; - for (j=i+1; j<=N; j++) - sum = sum + A(i,j)*b(j); - b(i) = ( b(i) - sum ) / - D(i); - } - - return 0; -} - - -template <class MaTRiX, class Vector> -int QR_solve(const MaTRiX &A, const Vector &c, /*const*/ Vector &d, - Vector &b) -{ - assert(A.lbound() == 1); // ensure these are all - assert(c.lbound() == 1); // 1-based arrays and vectors - assert(d.lbound() == 1); - - Subscript N=A.num_rows(); - - assert(N == A.num_cols()); - assert(N == c.dim()); - assert(N == d.dim()); - assert(N == b.dim()); - - Subscript i,j; - typename MaTRiX::element_type sum, tau; - - for (j=1; j<N; j++) - { - // form Q'*b - sum = 0; - for (i=j; i<=N; i++) - sum = sum + A(i,j)*b(i); - if (c(j) == 0) - return 1; - tau = sum / c(j); - for (i=j; i<=N; i++) - b(i) = b(i) - tau * A(i,j); - } - return R_solve(A, d, b); // solve Rx = Q'b -} - -} // namespace TNT - -#endif // QR_H - diff --git a/intern/iksolver/intern/TNT/region1d.h b/intern/iksolver/intern/TNT/region1d.h deleted file mode 100644 index 8acaac3ae23..00000000000 --- a/intern/iksolver/intern/TNT/region1d.h +++ /dev/null @@ -1,375 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - - -#ifndef REGION1D_H -#define REGION1D_H - - -#include "subscript.h" -#include "index.h" -#include <iostream> -#include <cassert> - -namespace TNT -{ - -template <class Array1D> -class const_Region1D; - -template <class Array1D> -class Region1D -{ - protected: - - Array1D & A_; - Subscript offset_; // 0-based - Subscript dim_; - - typedef typename Array1D::element_type T; - - public: - const Array1D & array() const { return A_; } - - Subscript offset() const { return offset_;} - Subscript dim() const { return dim_; } - - Subscript offset(Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(i==TNT_BASE_OFFSET); -#endif - return offset_; - } - - Subscript dim(Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(i== TNT_BASE_OFFSET); -#endif - return offset_; - } - - - Region1D(Array1D &A, Subscript i1, Subscript i2) : A_(A) - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET <= i1 ); - assert(i2 <= A.dim() + (TNT_BASE_OFFSET-1)); - assert(i1 <= i2); -#endif - offset_ = i1 - TNT_BASE_OFFSET; - dim_ = i2-i1 + 1; - } - - Region1D(Array1D &A, const Index1D &I) : A_(A) - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET <=I.lbound()); - assert(I.ubound() <= A.dim() + (TNT_BASE_OFFSET-1)); - assert(I.lbound() <= I.ubound()); -#endif - offset_ = I.lbound() - TNT_BASE_OFFSET; - dim_ = I.ubound() - I.lbound() + 1; - } - - Region1D(Region1D<Array1D> &A, Subscript i1, Subscript i2) : - A_(A.A_) - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET <= i1 ); - assert(i2 <= A.dim() + (TNT_BASE_OFFSET - 1)); - assert(i1 <= i2); -#endif - // (old-offset) (new-offset) - // - offset_ = (i1 - TNT_BASE_OFFSET) + A.offset_; - dim_ = i2-i1 + 1; - } - - Region1D<Array1D> operator()(Subscript i1, Subscript i2) - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET <= i1); - assert(i2 <= dim() + (TNT_BASE_OFFSET -1)); - assert(i1 <= i2); -#endif - // offset_ is 0-based, so no need for - // ( - TNT_BASE_OFFSET) - // - return Region1D<Array1D>(A_, i1+offset_, - offset_ + i2); - } - - - Region1D<Array1D> operator()(const Index1D &I) - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET<=I.lbound()); - assert(I.ubound() <= dim() + (TNT_BASE_OFFSET-1)); - assert(I.lbound() <= I.ubound()); -#endif - return Region1D<Array1D>(A_, I.lbound()+offset_, - offset_ + I.ubound()); - } - - - - - T & operator()(Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET <= i); - assert(i <= dim() + (TNT_BASE_OFFSET-1)); -#endif - return A_(i+offset_); - } - - const T & operator() (Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET <= i); - assert(i <= dim() + (TNT_BASE_OFFSET-1)); -#endif - return A_(i+offset_); - } - - - Region1D<Array1D> & operator=(const Region1D<Array1D> &R) - { - // make sure both sides conform - assert(dim() == R.dim()); - - Subscript N = dim(); - Subscript i; - Subscript istart = TNT_BASE_OFFSET; - Subscript iend = istart + N-1; - - for (i=istart; i<=iend; i++) - (*this)(i) = R(i); - - return *this; - } - - - - Region1D<Array1D> & operator=(const const_Region1D<Array1D> &R) - { - // make sure both sides conform - assert(dim() == R.dim()); - - Subscript N = dim(); - Subscript i; - Subscript istart = TNT_BASE_OFFSET; - Subscript iend = istart + N-1; - - for (i=istart; i<=iend; i++) - (*this)(i) = R(i); - - return *this; - - } - - - Region1D<Array1D> & operator=(const T& t) - { - Subscript N=dim(); - Subscript i; - Subscript istart = TNT_BASE_OFFSET; - Subscript iend = istart + N-1; - - for (i=istart; i<= iend; i++) - (*this)(i) = t; - - return *this; - - } - - - Region1D<Array1D> & operator=(const Array1D &R) - { - // make sure both sides conform - Subscript N = dim(); - assert(dim() == R.dim()); - - Subscript i; - Subscript istart = TNT_BASE_OFFSET; - Subscript iend = istart + N-1; - - for (i=istart; i<=iend; i++) - (*this)(i) = R(i); - - return *this; - - } - -}; - -template <class Array1D> -std::ostream& operator<<(std::ostream &s, Region1D<Array1D> &A) -{ - Subscript N=A.dim(); - Subscript istart = TNT_BASE_OFFSET; - Subscript iend = N - 1 + TNT_BASE_OFFSET; - - for (Subscript i=istart; i<=iend; i++) - s << A(i) << endl; - - return s; -} - - -/* --------- class const_Region1D ------------ */ - -template <class Array1D> -class const_Region1D -{ - protected: - - const Array1D & A_; - Subscript offset_; // 0-based - Subscript dim_; - typedef typename Array1D::element_type T; - - public: - const Array1D & array() const { return A_; } - - Subscript offset() const { return offset_;} - Subscript dim() const { return dim_; } - - Subscript offset(Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(i==TNT_BASE_OFFSET); -#endif - return offset_; - } - - Subscript dim(Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(i== TNT_BASE_OFFSET); -#endif - return offset_; - } - - - const_Region1D(const Array1D &A, Subscript i1, Subscript i2) : A_(A) - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET <= i1 ); - assert(i2 <= A.dim() + (TNT_BASE_OFFSET-1)); - assert(i1 <= i2); -#endif - offset_ = i1 - TNT_BASE_OFFSET; - dim_ = i2-i1 + 1; - } - - const_Region1D(const Array1D &A, const Index1D &I) : A_(A) - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET <=I.lbound()); - assert(I.ubound() <= A.dim() + (TNT_BASE_OFFSET-1)); - assert(I.lbound() <= I.ubound()); -#endif - offset_ = I.lbound() - TNT_BASE_OFFSET; - dim_ = I.ubound() - I.lbound() + 1; - } - - const_Region1D(const_Region1D<Array1D> &A, Subscript i1, Subscript i2) : - A_(A.A_) - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET <= i1 ); - assert(i2 <= A.dim() + (TNT_BASE_OFFSET - 1)); - assert(i1 <= i2); -#endif - // (old-offset) (new-offset) - // - offset_ = (i1 - TNT_BASE_OFFSET) + A.offset_; - dim_ = i2-i1 + 1; - } - - const_Region1D<Array1D> operator()(Subscript i1, Subscript i2) - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET <= i1); - assert(i2 <= dim() + (TNT_BASE_OFFSET -1)); - assert(i1 <= i2); -#endif - // offset_ is 0-based, so no need for - // ( - TNT_BASE_OFFSET) - // - return const_Region1D<Array1D>(A_, i1+offset_, - offset_ + i2); - } - - - const_Region1D<Array1D> operator()(const Index1D &I) - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET<=I.lbound()); - assert(I.ubound() <= dim() + (TNT_BASE_OFFSET-1)); - assert(I.lbound() <= I.ubound()); -#endif - return const_Region1D<Array1D>(A_, I.lbound()+offset_, - offset_ + I.ubound()); - } - - - const T & operator() (Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(TNT_BASE_OFFSET <= i); - assert(i <= dim() + (TNT_BASE_OFFSET-1)); -#endif - return A_(i+offset_); - } - - - - -}; - -template <class Array1D> -std::ostream& operator<<(std::ostream &s, const_Region1D<Array1D> &A) -{ - Subscript N=A.dim(); - - for (Subscript i=1; i<=N; i++) - s << A(i) << endl; - - return s; -} - - -} // namespace TNT - -#endif // const_Region1D_H - diff --git a/intern/iksolver/intern/TNT/region2d.h b/intern/iksolver/intern/TNT/region2d.h deleted file mode 100644 index 6af50262382..00000000000 --- a/intern/iksolver/intern/TNT/region2d.h +++ /dev/null @@ -1,471 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - -// 2D Regions for arrays and matrices - -#ifndef REGION2D_H -#define REGION2D_H - -#include "index.h" -#include <iostream> -#include <cassert> - -namespace TNT -{ - -template <class Array2D> -class const_Region2D; - - -template <class Array2D> -class Region2D -{ - protected: - - Array2D & A_; - Subscript offset_[2]; // 0-offset internally - Subscript dim_[2]; - - public: - typedef typename Array2D::value_type T; - typedef Subscript size_type; - typedef T value_type; - typedef T element_type; - typedef T* pointer; - typedef T* iterator; - typedef T& reference; - typedef const T* const_iterator; - typedef const T& const_reference; - - Array2D & array() { return A_; } - const Array2D & array() const { return A_; } - Subscript lbound() const { return A_.lbound(); } - Subscript num_rows() const { return dim_[0]; } - Subscript num_cols() const { return dim_[1]; } - Subscript offset(Subscript i) const // 1-offset - { -#ifdef TNT_BOUNDS_CHECK - assert( A_.lbound() <= i); - assert( i<= dim_[0] + A_.lbound()-1); -#endif - return offset_[i-A_.lbound()]; - } - - Subscript dim(Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert( A_.lbound() <= i); - assert( i<= dim_[0] + A_.lbound()-1); -#endif - return dim_[i-A_.lbound()]; - } - - - - Region2D(Array2D &A, Subscript i1, Subscript i2, Subscript j1, - Subscript j2) : A_(A) - { -#ifdef TNT_BOUNDS_CHECK - assert( i1 <= i2 ); - assert( j1 <= j2); - assert( A.lbound() <= i1); - assert( i2<= A.dim(A.lbound()) + A.lbound()-1); - assert( A.lbound() <= j1); - assert( j2<= A.dim(A.lbound()+1) + A.lbound()-1 ); -#endif - - - offset_[0] = i1-A.lbound(); - offset_[1] = j1-A.lbound(); - dim_[0] = i2-i1+1; - dim_[1] = j2-j1+1; - } - - Region2D(Array2D &A, const Index1D &I, const Index1D &J) : A_(A) - { -#ifdef TNT_BOUNDS_CHECK - assert( I.lbound() <= I.ubound() ); - assert( J.lbound() <= J.ubound() ); - assert( A.lbound() <= I.lbound()); - assert( I.ubound()<= A.dim(A.lbound()) + A.lbound()-1); - assert( A.lbound() <= J.lbound()); - assert( J.ubound() <= A.dim(A.lbound()+1) + A.lbound()-1 ); -#endif - - offset_[0] = I.lbound()-A.lbound(); - offset_[1] = J.lbound()-A.lbound(); - dim_[0] = I.ubound() - I.lbound() + 1; - dim_[1] = J.ubound() - J.lbound() + 1; - } - - Region2D(Region2D<Array2D> &A, Subscript i1, Subscript i2, - Subscript j1, Subscript j2) : A_(A.A_) - { -#ifdef TNT_BOUNDS_CHECK - assert( i1 <= i2 ); - assert( j1 <= j2); - assert( A.lbound() <= i1); - assert( i2<= A.dim(A.lbound()) + A.lbound()-1); - assert( A.lbound() <= j1); - assert( j2<= A.dim(A.lbound()+1) + A.lbound()-1 ); -#endif - offset_[0] = (i1 - A.lbound()) + A.offset_[0]; - offset_[1] = (j1 - A.lbound()) + A.offset_[1]; - dim_[0] = i2-i1 + 1; - dim_[1] = j2-j1+1; - } - - Region2D<Array2D> operator()(Subscript i1, Subscript i2, - Subscript j1, Subscript j2) - { -#ifdef TNT_BOUNDS_CHECK - assert( i1 <= i2 ); - assert( j1 <= j2); - assert( A_.lbound() <= i1); - assert( i2<= dim_[0] + A_.lbound()-1); - assert( A_.lbound() <= j1); - assert( j2<= dim_[1] + A_.lbound()-1 ); -#endif - return Region2D<Array2D>(A_, - i1+offset_[0], offset_[0] + i2, - j1+offset_[1], offset_[1] + j2); - } - - - Region2D<Array2D> operator()(const Index1D &I, - const Index1D &J) - { -#ifdef TNT_BOUNDS_CHECK - assert( I.lbound() <= I.ubound() ); - assert( J.lbound() <= J.ubound() ); - assert( A_.lbound() <= I.lbound()); - assert( I.ubound()<= dim_[0] + A_.lbound()-1); - assert( A_.lbound() <= J.lbound()); - assert( J.ubound() <= dim_[1] + A_.lbound()-1 ); -#endif - - return Region2D<Array2D>(A_, I.lbound()+offset_[0], - offset_[0] + I.ubound(), offset_[1]+J.lbound(), - offset_[1] + J.ubound()); - } - - inline T & operator()(Subscript i, Subscript j) - { -#ifdef TNT_BOUNDS_CHECK - assert( A_.lbound() <= i); - assert( i<= dim_[0] + A_.lbound()-1); - assert( A_.lbound() <= j); - assert( j<= dim_[1] + A_.lbound()-1 ); -#endif - return A_(i+offset_[0], j+offset_[1]); - } - - inline const T & operator() (Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert( A_.lbound() <= i); - assert( i<= dim_[0] + A_.lbound()-1); - assert( A_.lbound() <= j); - assert( j<= dim_[1] + A_.lbound()-1 ); -#endif - return A_(i+offset_[0], j+offset_[1]); - } - - - Region2D<Array2D> & operator=(const Region2D<Array2D> &R) - { - Subscript M = num_rows(); - Subscript N = num_cols(); - - // make sure both sides conform - assert(M == R.num_rows()); - assert(N == R.num_cols()); - - - Subscript start = R.lbound(); - Subscript Mend = start + M - 1; - Subscript Nend = start + N - 1; - for (Subscript i=start; i<=Mend; i++) - for (Subscript j=start; j<=Nend; j++) - (*this)(i,j) = R(i,j); - - return *this; - } - - Region2D<Array2D> & operator=(const const_Region2D<Array2D> &R) - { - Subscript M = num_rows(); - Subscript N = num_cols(); - - // make sure both sides conform - assert(M == R.num_rows()); - assert(N == R.num_cols()); - - - Subscript start = R.lbound(); - Subscript Mend = start + M - 1; - Subscript Nend = start + N - 1; - for (Subscript i=start; i<=Mend; i++) - for (Subscript j=start; j<=Nend; j++) - (*this)(i,j) = R(i,j); - - return *this; - } - - Region2D<Array2D> & operator=(const Array2D &R) - { - Subscript M = num_rows(); - Subscript N = num_cols(); - - // make sure both sides conform - assert(M == R.num_rows()); - assert(N == R.num_cols()); - - - Subscript start = R.lbound(); - Subscript Mend = start + M - 1; - Subscript Nend = start + N - 1; - for (Subscript i=start; i<=Mend; i++) - for (Subscript j=start; j<=Nend; j++) - (*this)(i,j) = R(i,j); - - return *this; - } - - Region2D<Array2D> & operator=(const T &scalar) - { - Subscript start = lbound(); - Subscript Mend = lbound() + num_rows() - 1; - Subscript Nend = lbound() + num_cols() - 1; - - - for (Subscript i=start; i<=Mend; i++) - for (Subscript j=start; j<=Nend; j++) - (*this)(i,j) = scalar; - - return *this; - } - -}; - -//************************ - -template <class Array2D> -class const_Region2D -{ - protected: - - const Array2D & A_; - Subscript offset_[2]; // 0-offset internally - Subscript dim_[2]; - - public: - typedef typename Array2D::value_type T; - typedef T value_type; - typedef T element_type; - typedef const T* const_iterator; - typedef const T& const_reference; - - const Array2D & array() const { return A_; } - Subscript lbound() const { return A_.lbound(); } - Subscript num_rows() const { return dim_[0]; } - Subscript num_cols() const { return dim_[1]; } - Subscript offset(Subscript i) const // 1-offset - { -#ifdef TNT_BOUNDS_CHECK - assert( TNT_BASE_OFFSET <= i); - assert( i<= dim_[0] + TNT_BASE_OFFSET-1); -#endif - return offset_[i-TNT_BASE_OFFSET]; - } - - Subscript dim(Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert( TNT_BASE_OFFSET <= i); - assert( i<= dim_[0] + TNT_BASE_OFFSET-1); -#endif - return dim_[i-TNT_BASE_OFFSET]; - } - - - const_Region2D(const Array2D &A, Subscript i1, Subscript i2, - Subscript j1, Subscript j2) : A_(A) - { -#ifdef TNT_BOUNDS_CHECK - assert( i1 <= i2 ); - assert( j1 <= j2); - assert( TNT_BASE_OFFSET <= i1); - assert( i2<= A.dim(TNT_BASE_OFFSET) + TNT_BASE_OFFSET-1); - assert( TNT_BASE_OFFSET <= j1); - assert( j2<= A.dim(TNT_BASE_OFFSET+1) + TNT_BASE_OFFSET-1 ); -#endif - - offset_[0] = i1-TNT_BASE_OFFSET; - offset_[1] = j1-TNT_BASE_OFFSET; - dim_[0] = i2-i1+1; - dim_[1] = j2-j1+1; - } - - const_Region2D(const Array2D &A, const Index1D &I, const Index1D &J) - : A_(A) - { -#ifdef TNT_BOUNDS_CHECK - assert( I.lbound() <= I.ubound() ); - assert( J.lbound() <= J.ubound() ); - assert( TNT_BASE_OFFSET <= I.lbound()); - assert( I.ubound()<= A.dim(TNT_BASE_OFFSET) + TNT_BASE_OFFSET-1); - assert( TNT_BASE_OFFSET <= J.lbound()); - assert( J.ubound() <= A.dim(TNT_BASE_OFFSET+1) + TNT_BASE_OFFSET-1 ); -#endif - - offset_[0] = I.lbound()-TNT_BASE_OFFSET; - offset_[1] = J.lbound()-TNT_BASE_OFFSET; - dim_[0] = I.ubound() - I.lbound() + 1; - dim_[1] = J.ubound() - J.lbound() + 1; - } - - - const_Region2D(const_Region2D<Array2D> &A, Subscript i1, - Subscript i2, - Subscript j1, Subscript j2) : A_(A.A_) - { -#ifdef TNT_BOUNDS_CHECK - assert( i1 <= i2 ); - assert( j1 <= j2); - assert( TNT_BASE_OFFSET <= i1); - assert( i2<= A.dim(TNT_BASE_OFFSET) + TNT_BASE_OFFSET-1); - assert( TNT_BASE_OFFSET <= j1); - assert( j2<= A.dim(TNT_BASE_OFFSET+1) + TNT_BASE_OFFSET-1 ); -#endif - offset_[0] = (i1 - TNT_BASE_OFFSET) + A.offset_[0]; - offset_[1] = (j1 - TNT_BASE_OFFSET) + A.offset_[1]; - dim_[0] = i2-i1 + 1; - dim_[1] = j2-j1+1; - } - - const_Region2D<Array2D> operator()(Subscript i1, Subscript i2, - Subscript j1, Subscript j2) - { -#ifdef TNT_BOUNDS_CHECK - assert( i1 <= i2 ); - assert( j1 <= j2); - assert( TNT_BASE_OFFSET <= i1); - assert( i2<= dim_[0] + TNT_BASE_OFFSET-1); - assert( TNT_BASE_OFFSET <= j1); - assert( j2<= dim_[0] + TNT_BASE_OFFSET-1 ); -#endif - return const_Region2D<Array2D>(A_, - i1+offset_[0], offset_[0] + i2, - j1+offset_[1], offset_[1] + j2); - } - - - const_Region2D<Array2D> operator()(const Index1D &I, - const Index1D &J) - { -#ifdef TNT_BOUNDS_CHECK - assert( I.lbound() <= I.ubound() ); - assert( J.lbound() <= J.ubound() ); - assert( TNT_BASE_OFFSET <= I.lbound()); - assert( I.ubound()<= dim_[0] + TNT_BASE_OFFSET-1); - assert( TNT_BASE_OFFSET <= J.lbound()); - assert( J.ubound() <= dim_[1] + TNT_BASE_OFFSET-1 ); -#endif - - return const_Region2D<Array2D>(A_, I.lbound()+offset_[0], - offset_[0] + I.ubound(), offset_[1]+J.lbound(), - offset_[1] + J.ubound()); - } - - - inline const T & operator() (Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert( TNT_BASE_OFFSET <= i); - assert( i<= dim_[0] + TNT_BASE_OFFSET-1); - assert( TNT_BASE_OFFSET <= j); - assert( j<= dim_[1] + TNT_BASE_OFFSET-1 ); -#endif - return A_(i+offset_[0], j+offset_[1]); - } - -}; - - -// ************** std::ostream algorithms ******************************* - -template <class Array2D> -std::ostream& operator<<(std::ostream &s, const const_Region2D<Array2D> &A) -{ - Subscript start = A.lbound(); - Subscript Mend=A.lbound()+ A.num_rows() - 1; - Subscript Nend=A.lbound() + A.num_cols() - 1; - - - s << A.num_rows() << " " << A.num_cols() << "\n"; - for (Subscript i=start; i<=Mend; i++) - { - for (Subscript j=start; j<=Nend; j++) - { - s << A(i,j) << " "; - } - s << "\n"; - } - - - return s; -} - -template <class Array2D> -std::ostream& operator<<(std::ostream &s, const Region2D<Array2D> &A) -{ - Subscript start = A.lbound(); - Subscript Mend=A.lbound()+ A.num_rows() - 1; - Subscript Nend=A.lbound() + A.num_cols() - 1; - - - s << A.num_rows() << " " << A.num_cols() << "\n"; - for (Subscript i=start; i<=Mend; i++) - { - for (Subscript j=start; j<=Nend; j++) - { - s << A(i,j) << " "; - } - s << "\n"; - } - - - return s; - -} - -} // namespace TNT - -#endif // REGION2D_H - diff --git a/intern/iksolver/intern/TNT/stopwatch.h b/intern/iksolver/intern/TNT/stopwatch.h deleted file mode 100644 index 55cd532d6a4..00000000000 --- a/intern/iksolver/intern/TNT/stopwatch.h +++ /dev/null @@ -1,83 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - -#ifndef STPWATCH_H -#define STPWATCH_H - -// for clock() and CLOCKS_PER_SEC -#include <ctime> - -namespace TNT -{ - -/* Simple stopwatch object: - - void start() : start timing - double stop() : stop timing - void reset() : set elapsed time to 0.0 - double read() : read elapsed time (in seconds) - -*/ - -inline double seconds(void) -{ - static const double secs_per_tick = 1.0 / CLOCKS_PER_SEC; - return ( (double) clock() ) * secs_per_tick; -} - - -class stopwatch { - private: - int running; - double last_time; - double total; - - public: - stopwatch() : running(0), last_time(0.0), total(0.0) {} - void reset() { running = 0; last_time = 0.0; total=0.0; } - void start() { if (!running) { last_time = seconds(); running = 1;}} - double stop() { if (running) - { - total += seconds() - last_time; - running = 0; - } - return total; - } - double read() { if (running) - { - total+= seconds() - last_time; - last_time = seconds(); - } - return total; - } - -}; - -} // namespace TNT - -#endif - diff --git a/intern/iksolver/intern/TNT/subscript.h b/intern/iksolver/intern/TNT/subscript.h deleted file mode 100644 index 75555e1c0b6..00000000000 --- a/intern/iksolver/intern/TNT/subscript.h +++ /dev/null @@ -1,63 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - -#ifndef SUBSCRPT_H -#define SUBSCRPT_H - - -//--------------------------------------------------------------------- -// This definition describes the default TNT data type used for -// indexing into TNT matrices and vectors. The data type should -// be wide enough to index into large arrays. It defaults to an -// "int", but can be overriden at compile time redefining TNT_SUBSCRIPT_TYPE, -// e.g. -// -// g++ -DTNT_SUBSCRIPT_TYPE='unsigned int' ... -// -//--------------------------------------------------------------------- -// - -#ifndef TNT_SUBSCRIPT_TYPE -#define TNT_SUBSCRIPT_TYPE int -#endif - -namespace TNT -{ - typedef TNT_SUBSCRIPT_TYPE Subscript; -} - - -// () indexing in TNT means 1-offset, i.e. x(1) and A(1,1) are the -// first elements. This offset is left as a macro for future -// purposes, but should not be changed in the current release. -// -// -#define TNT_BASE_OFFSET (1) - -#endif - diff --git a/intern/iksolver/intern/TNT/svd.h b/intern/iksolver/intern/TNT/svd.h deleted file mode 100644 index cfff73caa3f..00000000000 --- a/intern/iksolver/intern/TNT/svd.h +++ /dev/null @@ -1,435 +0,0 @@ -/** - */ - -#ifndef SVD_H - -#define SVD_H - -// Compute the Single Value Decomposition of an arbitrary matrix A -// That is compute the 3 matrices U,W,V with U column orthogonal (m,n) -// ,W a diagonal matrix and V an orthogonal square matrix s.t. -// A = U.W.Vt. From this decomposition it is trivial to compute the -// inverse of A as Ainv = V.Winv.tranpose(U). -// -// s = diagonal elements of W -// work1 = workspace, length must be A.num_rows -// work2 = workspace, length must be A.num_cols - -#include "tntmath.h" - -#define SVD_MAX_ITER 200 - -namespace TNT -{ - -template <class MaTRiX, class VecToR > -void SVD(MaTRiX &A, MaTRiX &U, VecToR &s, MaTRiX &V, VecToR &work1, VecToR &work2, int maxiter=SVD_MAX_ITER) { - - int m = A.num_rows(); - int n = A.num_cols(); - int nu = min(m,n); - - VecToR& work = work1; - VecToR& e = work2; - - U = 0; - s = 0; - - int i=0, j=0, k=0; - - // Reduce A to bidiagonal form, storing the diagonal elements - // in s and the super-diagonal elements in e. - - int nct = min(m-1,n); - int nrt = max(0,min(n-2,m)); - for (k = 0; k < max(nct,nrt); k++) { - if (k < nct) { - - // Compute the transformation for the k-th column and - // place the k-th diagonal in s[k]. - // Compute 2-norm of k-th column without under/overflow. - s[k] = 0; - for (i = k; i < m; i++) { - s[k] = hypot(s[k],A[i][k]); - } - if (s[k] != 0.0) { - if (A[k][k] < 0.0) { - s[k] = -s[k]; - } - for (i = k; i < m; i++) { - A[i][k] /= s[k]; - } - A[k][k] += 1.0; - } - s[k] = -s[k]; - } - for (j = k+1; j < n; j++) { - if ((k < nct) && (s[k] != 0.0)) { - - // Apply the transformation. - - typename MaTRiX::value_type t = 0; - for (i = k; i < m; i++) { - t += A[i][k]*A[i][j]; - } - t = -t/A[k][k]; - for (i = k; i < m; i++) { - A[i][j] += t*A[i][k]; - } - } - - // Place the k-th row of A into e for the - // subsequent calculation of the row transformation. - - e[j] = A[k][j]; - } - if (k < nct) { - - // Place the transformation in U for subsequent back - // multiplication. - - for (i = k; i < m; i++) - U[i][k] = A[i][k]; - } - if (k < nrt) { - - // Compute the k-th row transformation and place the - // k-th super-diagonal in e[k]. - // Compute 2-norm without under/overflow. - e[k] = 0; - for (i = k+1; i < n; i++) { - e[k] = hypot(e[k],e[i]); - } - if (e[k] != 0.0) { - if (e[k+1] < 0.0) { - e[k] = -e[k]; - } - for (i = k+1; i < n; i++) { - e[i] /= e[k]; - } - e[k+1] += 1.0; - } - e[k] = -e[k]; - if ((k+1 < m) & (e[k] != 0.0)) { - - // Apply the transformation. - - for (i = k+1; i < m; i++) { - work[i] = 0.0; - } - for (j = k+1; j < n; j++) { - for (i = k+1; i < m; i++) { - work[i] += e[j]*A[i][j]; - } - } - for (j = k+1; j < n; j++) { - typename MaTRiX::value_type t = -e[j]/e[k+1]; - for (i = k+1; i < m; i++) { - A[i][j] += t*work[i]; - } - } - } - - // Place the transformation in V for subsequent - // back multiplication. - - for (i = k+1; i < n; i++) - V[i][k] = e[i]; - } - } - - // Set up the final bidiagonal matrix or order p. - - int p = min(n,m+1); - if (nct < n) { - s[nct] = A[nct][nct]; - } - if (m < p) { - s[p-1] = 0.0; - } - if (nrt+1 < p) { - e[nrt] = A[nrt][p-1]; - } - e[p-1] = 0.0; - - // If required, generate U. - - for (j = nct; j < nu; j++) { - for (i = 0; i < m; i++) { - U[i][j] = 0.0; - } - U[j][j] = 1.0; - } - for (k = nct-1; k >= 0; k--) { - if (s[k] != 0.0) { - for (j = k+1; j < nu; j++) { - typename MaTRiX::value_type t = 0; - for (i = k; i < m; i++) { - t += U[i][k]*U[i][j]; - } - t = -t/U[k][k]; - for (i = k; i < m; i++) { - U[i][j] += t*U[i][k]; - } - } - for (i = k; i < m; i++ ) { - U[i][k] = -U[i][k]; - } - U[k][k] = 1.0 + U[k][k]; - for (i = 0; i < k-1; i++) { - U[i][k] = 0.0; - } - } else { - for (i = 0; i < m; i++) { - U[i][k] = 0.0; - } - U[k][k] = 1.0; - } - } - - // If required, generate V. - - for (k = n-1; k >= 0; k--) { - if ((k < nrt) & (e[k] != 0.0)) { - for (j = k+1; j < nu; j++) { - typename MaTRiX::value_type t = 0; - for (i = k+1; i < n; i++) { - t += V[i][k]*V[i][j]; - } - t = -t/V[k+1][k]; - for (i = k+1; i < n; i++) { - V[i][j] += t*V[i][k]; - } - } - } - for (i = 0; i < n; i++) { - V[i][k] = 0.0; - } - V[k][k] = 1.0; - } - - // Main iteration loop for the singular values. - - int pp = p-1; - int iter = 0; - typename MaTRiX::value_type eps = pow(2.0,-52.0); - while (p > 0) { - int kase=0; - k=0; - - // Test for maximum iterations to avoid infinite loop - if(maxiter == 0) - break; - maxiter--; - - // This section of the program inspects for - // negligible elements in the s and e arrays. On - // completion the variables kase and k are set as follows. - - // kase = 1 if s(p) and e[k-1] are negligible and k<p - // kase = 2 if s(k) is negligible and k<p - // kase = 3 if e[k-1] is negligible, k<p, and - // s(k), ..., s(p) are not negligible (qr step). - // kase = 4 if e(p-1) is negligible (convergence). - - for (k = p-2; k >= -1; k--) { - if (k == -1) { - break; - } - if (TNT::abs(e[k]) <= eps*(TNT::abs(s[k]) + TNT::abs(s[k+1]))) { - e[k] = 0.0; - break; - } - } - if (k == p-2) { - kase = 4; - } else { - int ks; - for (ks = p-1; ks >= k; ks--) { - if (ks == k) { - break; - } - typename MaTRiX::value_type t = (ks != p ? TNT::abs(e[ks]) : 0.) + - (ks != k+1 ? TNT::abs(e[ks-1]) : 0.); - if (TNT::abs(s[ks]) <= eps*t) { - s[ks] = 0.0; - break; - } - } - if (ks == k) { - kase = 3; - } else if (ks == p-1) { - kase = 1; - } else { - kase = 2; - k = ks; - } - } - k++; - - // Perform the task indicated by kase. - - switch (kase) { - - // Deflate negligible s(p). - - case 1: { - typename MaTRiX::value_type f = e[p-2]; - e[p-2] = 0.0; - for (j = p-2; j >= k; j--) { - typename MaTRiX::value_type t = hypot(s[j],f); - typename MaTRiX::value_type cs = s[j]/t; - typename MaTRiX::value_type sn = f/t; - s[j] = t; - if (j != k) { - f = -sn*e[j-1]; - e[j-1] = cs*e[j-1]; - } - - for (i = 0; i < n; i++) { - t = cs*V[i][j] + sn*V[i][p-1]; - V[i][p-1] = -sn*V[i][j] + cs*V[i][p-1]; - V[i][j] = t; - } - } - } - break; - - // Split at negligible s(k). - - case 2: { - typename MaTRiX::value_type f = e[k-1]; - e[k-1] = 0.0; - for (j = k; j < p; j++) { - typename MaTRiX::value_type t = hypot(s[j],f); - typename MaTRiX::value_type cs = s[j]/t; - typename MaTRiX::value_type sn = f/t; - s[j] = t; - f = -sn*e[j]; - e[j] = cs*e[j]; - - for (i = 0; i < m; i++) { - t = cs*U[i][j] + sn*U[i][k-1]; - U[i][k-1] = -sn*U[i][j] + cs*U[i][k-1]; - U[i][j] = t; - } - } - } - break; - - // Perform one qr step. - - case 3: { - - // Calculate the shift. - - typename MaTRiX::value_type scale = max(max(max(max( - TNT::abs(s[p-1]),TNT::abs(s[p-2])),TNT::abs(e[p-2])), - TNT::abs(s[k])),TNT::abs(e[k])); - typename MaTRiX::value_type sp = s[p-1]/scale; - typename MaTRiX::value_type spm1 = s[p-2]/scale; - typename MaTRiX::value_type epm1 = e[p-2]/scale; - typename MaTRiX::value_type sk = s[k]/scale; - typename MaTRiX::value_type ek = e[k]/scale; - typename MaTRiX::value_type b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0; - typename MaTRiX::value_type c = (sp*epm1)*(sp*epm1); - typename MaTRiX::value_type shift = 0.0; - if ((b != 0.0) || (c != 0.0)) { - shift = sqrt(b*b + c); - if (b < 0.0) { - shift = -shift; - } - shift = c/(b + shift); - } - typename MaTRiX::value_type f = (sk + sp)*(sk - sp) + shift; - typename MaTRiX::value_type g = sk*ek; - - // Chase zeros. - - for (j = k; j < p-1; j++) { - typename MaTRiX::value_type t = hypot(f,g); - /* division by zero checks added to avoid NaN (brecht) */ - typename MaTRiX::value_type cs = (t == 0.0f)? 0.0f: f/t; - typename MaTRiX::value_type sn = (t == 0.0f)? 0.0f: g/t; - if (j != k) { - e[j-1] = t; - } - f = cs*s[j] + sn*e[j]; - e[j] = cs*e[j] - sn*s[j]; - g = sn*s[j+1]; - s[j+1] = cs*s[j+1]; - - for (i = 0; i < n; i++) { - t = cs*V[i][j] + sn*V[i][j+1]; - V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1]; - V[i][j] = t; - } - - t = hypot(f,g); - /* division by zero checks added to avoid NaN (brecht) */ - cs = (t == 0.0f)? 0.0f: f/t; - sn = (t == 0.0f)? 0.0f: g/t; - s[j] = t; - f = cs*e[j] + sn*s[j+1]; - s[j+1] = -sn*e[j] + cs*s[j+1]; - g = sn*e[j+1]; - e[j+1] = cs*e[j+1]; - if (j < m-1) { - for (i = 0; i < m; i++) { - t = cs*U[i][j] + sn*U[i][j+1]; - U[i][j+1] = -sn*U[i][j] + cs*U[i][j+1]; - U[i][j] = t; - } - } - } - e[p-2] = f; - iter = iter + 1; - } - break; - - // Convergence. - - case 4: { - - // Make the singular values positive. - - if (s[k] <= 0.0) { - s[k] = (s[k] < 0.0 ? -s[k] : 0.0); - - for (i = 0; i <= pp; i++) - V[i][k] = -V[i][k]; - } - - // Order the singular values. - - while (k < pp) { - if (s[k] >= s[k+1]) { - break; - } - typename MaTRiX::value_type t = s[k]; - s[k] = s[k+1]; - s[k+1] = t; - if (k < n-1) { - for (i = 0; i < n; i++) { - t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t; - } - } - if (k < m-1) { - for (i = 0; i < m; i++) { - t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t; - } - } - k++; - } - iter = 0; - p--; - } - break; - } - } -} - -} - -#endif - diff --git a/intern/iksolver/intern/TNT/tnt.h b/intern/iksolver/intern/TNT/tnt.h deleted file mode 100644 index c15ae27c547..00000000000 --- a/intern/iksolver/intern/TNT/tnt.h +++ /dev/null @@ -1,93 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - -#ifndef TNT_H -#define TNT_H - -//--------------------------------------------------------------------- -// tnt.h TNT general header file. Defines default types -// and conventions. -//--------------------------------------------------------------------- - -//--------------------------------------------------------------------- -// Include current version -//--------------------------------------------------------------------- -#include "version.h" - -//--------------------------------------------------------------------- -// Define the data type used for matrix and vector Subscripts. -// This will default to "int", but it can be overriden at compile time, -// e.g. -// -// g++ -DTNT_SUBSCRIPT_TYPE='unsinged long' ... -// -// See subscript.h for details. -//--------------------------------------------------------------------- - -#include "subscript.h" - - - -//--------------------------------------------------------------------- -// Define this macro if you want TNT to ensure all refernces -// are within the bounds of the array. This encurs a run-time -// overhead, of course, but is recommended while developing -// code. It can be turned off for production runs. -// -// #define TNT_BOUNDS_CHECK -//--------------------------------------------------------------------- -// -#define TNT_BOUNDS_CHECK -#ifdef TNT_NO_BOUNDS_CHECK -#undef TNT_BOUNDS_CHECK -#endif - -//--------------------------------------------------------------------- -// Define this macro if you want to utilize matrix and vector -// regions. This is typically on, but you can save some -// compilation time by turning it off. If you do this and -// attempt to use regions you will get an error message. -// -// #define TNT_USE_REGIONS -//--------------------------------------------------------------------- -// -#define TNT_USE_REGIONS - -//--------------------------------------------------------------------- -// -//--------------------------------------------------------------------- -// if your system doesn't have abs() min(), and max() uncoment the following -//--------------------------------------------------------------------- -// -// -//#define __NEED_ABS_MIN_MAX_ - -#include "tntmath.h" - -#endif // TNT_H - diff --git a/intern/iksolver/intern/TNT/tntmath.h b/intern/iksolver/intern/TNT/tntmath.h deleted file mode 100644 index be72796da59..00000000000 --- a/intern/iksolver/intern/TNT/tntmath.h +++ /dev/null @@ -1,154 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -// Header file for scalar math functions - -#ifndef TNTMATH_H -#define TNTMATH_H - -// conventional functions required by several matrix algorithms - -#if defined(_MSC_VER) && (_MSC_VER < 1800) -#define hypot _hypot -#endif - -namespace TNT -{ - -struct TNTException { - int i; -}; - - -inline double abs(double t) -{ - return ( t > 0 ? t : -t); -} - -inline double min(double a, double b) -{ - return (a < b ? a : b); -} - -inline double max(double a, double b) -{ - return (a > b ? a : b); -} - -inline float abs(float t) -{ - return ( t > 0 ? t : -t); -} - -inline float min(float a, float b) -{ - return (a < b ? a : b); -} - -inline int min(int a, int b) -{ - return (a < b ? a : b); -} - -inline int max(int a, int b) -{ - return (a > b ? a : b); -} - -inline float max(float a, float b) -{ - return (a > b ? a : b); -} - -inline double sign(double a) -{ - return (a > 0 ? 1.0 : -1.0); -} - -inline double sign(double a,double b) { - return (b >= 0.0 ? TNT::abs(a) : -TNT::abs(a)); -} - -inline float sign(float a,float b) { - return (b >= 0.0f ? TNT::abs(a) : -TNT::abs(a)); -} - -inline float sign(float a) -{ - return (a > 0.0 ? 1.0f : -1.0f); -} - -inline float pythag(float a, float b) -{ - float absa,absb; - absa = abs(a); - absb = abs(b); - - if (absa > absb) { - float sqr = absb/absa; - sqr *= sqr; - return absa * float(sqrt(1 + sqr)); - } else { - if (absb > float(0)) { - float sqr = absa/absb; - sqr *= sqr; - return absb * float(sqrt(1 + sqr)); - } else { - return float(0); - } - } -} - -inline double pythag(double a, double b) -{ - double absa,absb; - absa = abs(a); - absb = abs(b); - - if (absa > absb) { - double sqr = absb/absa; - sqr *= sqr; - return absa * double(sqrt(1 + sqr)); - } else { - - if (absb > double(0)) { - double sqr = absa/absb; - sqr *= sqr; - return absb * double(sqrt(1 + sqr)); - } else { - return double(0); - } - } -} - - -} /* namespace TNT */ - -#endif /* TNTMATH_H */ - diff --git a/intern/iksolver/intern/TNT/tntreqs.h b/intern/iksolver/intern/TNT/tntreqs.h deleted file mode 100644 index 2ee124d5940..00000000000 --- a/intern/iksolver/intern/TNT/tntreqs.h +++ /dev/null @@ -1,73 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - -// The requirements for a bare-bones vector class: -// -// -// o) must have 0-based [] indexing for const and -// non-const objects (i.e. operator[] defined) -// -// o) must have size() method to denote the number of -// elements -// o) must clean up after itself when destructed -// (i.e. no memory leaks) -// -// -) must have begin() and end() methods (The begin() -// method is necessary, because relying on -// &v_[0] may not work on a empty vector (i.e. v_ is NULL.) -// -// o) must be templated -// o) must have X::value_type defined to be the types of elements -// o) must have X::X(const &x) copy constructor (by *value*) -// o) must have X::X(int N) constructor to N-length vector -// (NOTE: this constructor need *NOT* initalize elements) -// -// -) must have X::X(int N, T scalar) constructor to initalize -// elements to value of "scalar". -// -// ( removed, because valarray<> class uses (scalar, N) rather -// than (N, scalar) ) -// -) must have X::X(int N, const T* scalars) constructor to copy from -// any C linear array -// -// ( removed, because of same reverse order of valarray<> ) -// -// o) must have assignment A=B, by value -// -// NOTE: this class is *NOT* meant to be derived from, -// so its methods (particularly indexing) need not be -// declared virtual. -// -// -// Some things it *DOES NOT* need to do are -// -// o) bounds checking -// o) array referencing (e.g. reference counting) -// o) support () indexing -// o) I/O -// - diff --git a/intern/iksolver/intern/TNT/transv.h b/intern/iksolver/intern/TNT/transv.h deleted file mode 100644 index 5500cecff40..00000000000 --- a/intern/iksolver/intern/TNT/transv.h +++ /dev/null @@ -1,164 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -// Matrix Transpose Views - -#ifndef TRANSV_H -#define TRANSV_H - -#include <iostream> -#include <cassert> -#include "vec.h" - -namespace TNT -{ - -template <class Array2D> -class Transpose_View -{ - protected: - - const Array2D & A_; - - public: - - typedef typename Array2D::element_type T; - typedef T value_type; - typedef T element_type; - typedef T* pointer; - typedef T* iterator; - typedef T& reference; - typedef const T* const_iterator; - typedef const T& const_reference; - - - const Array2D & array() const { return A_; } - Subscript num_rows() const { return A_.num_cols();} - Subscript num_cols() const { return A_.num_rows();} - Subscript lbound() const { return A_.lbound(); } - Subscript dim(Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert( A_.lbound() <= i); - assert( i<= A_.lbound()+1); -#endif - if (i== A_.lbound()) - return num_rows(); - else - return num_cols(); - } - - - Transpose_View(const Transpose_View<Array2D> &A) : A_(A.A_) {}; - Transpose_View(const Array2D &A) : A_(A) {}; - - - inline const typename Array2D::element_type & operator()( - Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(lbound()<=i); - assert(i<=A_.num_cols() + lbound() - 1); - assert(lbound()<=j); - assert(j<=A_.num_rows() + lbound() - 1); -#endif - - return A_(j,i); - } - - -}; - -template <class Matrix> -Transpose_View<Matrix> Transpose_view(const Matrix &A) -{ - return Transpose_View<Matrix>(A); -} - -template <class Matrix, class T> -Vector<T> matmult( - const Transpose_View<Matrix> & A, - const Vector<T> &B) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(B.dim() == N); - - Vector<T> x(N); - - Subscript i, j; - T tmp = 0; - - for (i=1; i<=M; i++) - { - tmp = 0; - for (j=1; j<=N; j++) - tmp += A(i,j) * B(j); - x(i) = tmp; - } - - return x; -} - -template <class Matrix, class T> -inline Vector<T> operator*(const Transpose_View<Matrix> & A, const Vector<T> &B) -{ - return matmult(A,B); -} - - -template <class Matrix> -std::ostream& operator<<(std::ostream &s, const Transpose_View<Matrix> &A) -{ - Subscript M=A.num_rows(); - Subscript N=A.num_cols(); - - Subscript start = A.lbound(); - Subscript Mend = M + A.lbound() - 1; - Subscript Nend = N + A.lbound() - 1; - - s << M << " " << N << endl; - for (Subscript i=start; i<=Mend; i++) - { - for (Subscript j=start; j<=Nend; j++) - { - s << A(i,j) << " "; - } - s << endl; - } - - - return s; -} - -} // namespace TNT - -#endif // TRANSV_H - diff --git a/intern/iksolver/intern/TNT/triang.h b/intern/iksolver/intern/TNT/triang.h deleted file mode 100644 index 10ec2cb18c3..00000000000 --- a/intern/iksolver/intern/TNT/triang.h +++ /dev/null @@ -1,637 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -// Triangular Matrices (Views and Adpators) - -#ifndef TRIANG_H -#define TRIANG_H - -// default to use lower-triangular portions of arrays -// for symmetric matrices. - -namespace TNT -{ - -template <class MaTRiX> -class LowerTriangularView -{ - protected: - - - const MaTRiX &A_; - const typename MaTRiX::element_type zero_; - - public: - - - typedef typename MaTRiX::const_reference const_reference; - typedef const typename MaTRiX::element_type element_type; - typedef const typename MaTRiX::element_type value_type; - typedef element_type T; - - Subscript dim(Subscript d) const { return A_.dim(d); } - Subscript lbound() const { return A_.lbound(); } - Subscript num_rows() const { return A_.num_rows(); } - Subscript num_cols() const { return A_.num_cols(); } - - - // constructors - - LowerTriangularView(/*const*/ MaTRiX &A) : A_(A), zero_(0) {} - - - inline const_reference get(Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(lbound()<=i); - assert(i<=A_.num_rows() + lbound() - 1); - assert(lbound()<=j); - assert(j<=A_.num_cols() + lbound() - 1); -#endif - if (i<j) - return zero_; - else - return A_(i,j); - } - - - inline const_reference operator() (Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(lbound()<=i); - assert(i<=A_.num_rows() + lbound() - 1); - assert(lbound()<=j); - assert(j<=A_.num_cols() + lbound() - 1); -#endif - if (i<j) - return zero_; - else - return A_(i,j); - } - -#ifdef TNT_USE_REGIONS - - typedef const_Region2D< LowerTriangularView<MaTRiX> > - const_Region; - - const_Region operator()(/*const*/ Index1D &I, - /*const*/ Index1D &J) const - { - return const_Region(*this, I, J); - } - - const_Region operator()(Subscript i1, Subscript i2, - Subscript j1, Subscript j2) const - { - return const_Region(*this, i1, i2, j1, j2); - } - - - -#endif -// TNT_USE_REGIONS - -}; - - -/* *********** Lower_triangular_view() algorithms ****************** */ - -template <class MaTRiX, class VecToR> -VecToR matmult(/*const*/ LowerTriangularView<MaTRiX> &A, VecToR &x) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(N == x.dim()); - - Subscript i, j; - typename MaTRiX::element_type sum=0.0; - VecToR result(M); - - Subscript start = A.lbound(); - Subscript Mend = M + A.lbound() -1 ; - - for (i=start; i<=Mend; i++) - { - sum = 0.0; - for (j=start; j<=i; j++) - sum = sum + A(i,j)*x(j); - result(i) = sum; - } - - return result; -} - -template <class MaTRiX, class VecToR> -inline VecToR operator*(/*const*/ LowerTriangularView<MaTRiX> &A, VecToR &x) -{ - return matmult(A,x); -} - -template <class MaTRiX> -class UnitLowerTriangularView -{ - protected: - - const MaTRiX &A_; - const typename MaTRiX::element_type zero; - const typename MaTRiX::element_type one; - - public: - - typedef typename MaTRiX::const_reference const_reference; - typedef typename MaTRiX::element_type element_type; - typedef typename MaTRiX::element_type value_type; - typedef element_type T; - - Subscript lbound() const { return 1; } - Subscript dim(Subscript d) const { return A_.dim(d); } - Subscript num_rows() const { return A_.num_rows(); } - Subscript num_cols() const { return A_.num_cols(); } - - - // constructors - - UnitLowerTriangularView(/*const*/ MaTRiX &A) : A_(A), zero(0), one(1) {} - - - inline const_reference get(Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i<=A_.dim(1)); - assert(1<=j); - assert(j<=A_.dim(2)); - assert(0<=i && i<A_.dim(0) && 0<=j && j<A_.dim(1)); -#endif - if (i>j) - return A_(i,j); - else if (i==j) - return one; - else - return zero; - } - - - inline const_reference operator() (Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i<=A_.dim(1)); - assert(1<=j); - assert(j<=A_.dim(2)); -#endif - if (i>j) - return A_(i,j); - else if (i==j) - return one; - else - return zero; - } - - -#ifdef TNT_USE_REGIONS - // These are the "index-aware" features - - typedef const_Region2D< UnitLowerTriangularView<MaTRiX> > - const_Region; - - const_Region operator()(/*const*/ Index1D &I, - /*const*/ Index1D &J) const - { - return const_Region(*this, I, J); - } - - const_Region operator()(Subscript i1, Subscript i2, - Subscript j1, Subscript j2) const - { - return const_Region(*this, i1, i2, j1, j2); - } -#endif -// TNT_USE_REGIONS -}; - -template <class MaTRiX> -LowerTriangularView<MaTRiX> Lower_triangular_view( - /*const*/ MaTRiX &A) -{ - return LowerTriangularView<MaTRiX>(A); -} - - -template <class MaTRiX> -UnitLowerTriangularView<MaTRiX> Unit_lower_triangular_view( - /*const*/ MaTRiX &A) -{ - return UnitLowerTriangularView<MaTRiX>(A); -} - -template <class MaTRiX, class VecToR> -VecToR matmult(/*const*/ UnitLowerTriangularView<MaTRiX> &A, VecToR &x) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(N == x.dim()); - - Subscript i, j; - typename MaTRiX::element_type sum=0.0; - VecToR result(M); - - Subscript start = A.lbound(); - Subscript Mend = M + A.lbound() -1 ; - - for (i=start; i<=Mend; i++) - { - sum = 0.0; - for (j=start; j<i; j++) - sum = sum + A(i,j)*x(j); - result(i) = sum + x(i); - } - - return result; -} - -template <class MaTRiX, class VecToR> -inline VecToR operator*(/*const*/ UnitLowerTriangularView<MaTRiX> &A, VecToR &x) -{ - return matmult(A,x); -} - - -//********************** Algorithms ************************************* - - - -template <class MaTRiX> -std::ostream& operator<<(std::ostream &s, const LowerTriangularView<MaTRiX>&A) -{ - Subscript M=A.num_rows(); - Subscript N=A.num_cols(); - - s << M << " " << N << endl; - - for (Subscript i=1; i<=M; i++) - { - for (Subscript j=1; j<=N; j++) - { - s << A(i,j) << " "; - } - s << endl; - } - - - return s; -} - -template <class MaTRiX> -std::ostream& operator<<(std::ostream &s, - const UnitLowerTriangularView<MaTRiX>&A) -{ - Subscript M=A.num_rows(); - Subscript N=A.num_cols(); - - s << M << " " << N << endl; - - for (Subscript i=1; i<=M; i++) - { - for (Subscript j=1; j<=N; j++) - { - s << A(i,j) << " "; - } - s << endl; - } - - - return s; -} - - - -// ******************* Upper Triangular Section ************************** - -template <class MaTRiX> -class UpperTriangularView -{ - protected: - - - /*const*/ MaTRiX &A_; - /*const*/ typename MaTRiX::element_type zero_; - - public: - - - typedef typename MaTRiX::const_reference const_reference; - typedef /*const*/ typename MaTRiX::element_type element_type; - typedef /*const*/ typename MaTRiX::element_type value_type; - typedef element_type T; - - Subscript dim(Subscript d) const { return A_.dim(d); } - Subscript lbound() const { return A_.lbound(); } - Subscript num_rows() const { return A_.num_rows(); } - Subscript num_cols() const { return A_.num_cols(); } - - - // constructors - - UpperTriangularView(/*const*/ MaTRiX &A) : A_(A), zero_(0) {} - - - inline const_reference get(Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(lbound()<=i); - assert(i<=A_.num_rows() + lbound() - 1); - assert(lbound()<=j); - assert(j<=A_.num_cols() + lbound() - 1); -#endif - if (i>j) - return zero_; - else - return A_(i,j); - } - - - inline const_reference operator() (Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(lbound()<=i); - assert(i<=A_.num_rows() + lbound() - 1); - assert(lbound()<=j); - assert(j<=A_.num_cols() + lbound() - 1); -#endif - if (i>j) - return zero_; - else - return A_(i,j); - } - -#ifdef TNT_USE_REGIONS - - typedef const_Region2D< UpperTriangularView<MaTRiX> > - const_Region; - - const_Region operator()(const Index1D &I, - const Index1D &J) const - { - return const_Region(*this, I, J); - } - - const_Region operator()(Subscript i1, Subscript i2, - Subscript j1, Subscript j2) const - { - return const_Region(*this, i1, i2, j1, j2); - } - - - -#endif -// TNT_USE_REGIONS - -}; - - -/* *********** Upper_triangular_view() algorithms ****************** */ - -template <class MaTRiX, class VecToR> -VecToR matmult(/*const*/ UpperTriangularView<MaTRiX> &A, VecToR &x) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(N == x.dim()); - - Subscript i, j; - typename VecToR::element_type sum=0.0; - VecToR result(M); - - Subscript start = A.lbound(); - Subscript Mend = M + A.lbound() -1 ; - - for (i=start; i<=Mend; i++) - { - sum = 0.0; - for (j=i; j<=N; j++) - sum = sum + A(i,j)*x(j); - result(i) = sum; - } - - return result; -} - -template <class MaTRiX, class VecToR> -inline VecToR operator*(/*const*/ UpperTriangularView<MaTRiX> &A, VecToR &x) -{ - return matmult(A,x); -} - -template <class MaTRiX> -class UnitUpperTriangularView -{ - protected: - - const MaTRiX &A_; - const typename MaTRiX::element_type zero; - const typename MaTRiX::element_type one; - - public: - - typedef typename MaTRiX::const_reference const_reference; - typedef typename MaTRiX::element_type element_type; - typedef typename MaTRiX::element_type value_type; - typedef element_type T; - - Subscript lbound() const { return 1; } - Subscript dim(Subscript d) const { return A_.dim(d); } - Subscript num_rows() const { return A_.num_rows(); } - Subscript num_cols() const { return A_.num_cols(); } - - - // constructors - - UnitUpperTriangularView(/*const*/ MaTRiX &A) : A_(A), zero(0), one(1) {} - - - inline const_reference get(Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i<=A_.dim(1)); - assert(1<=j); - assert(j<=A_.dim(2)); - assert(0<=i && i<A_.dim(0) && 0<=j && j<A_.dim(1)); -#endif - if (i<j) - return A_(i,j); - else if (i==j) - return one; - else - return zero; - } - - - inline const_reference operator() (Subscript i, Subscript j) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i<=A_.dim(1)); - assert(1<=j); - assert(j<=A_.dim(2)); -#endif - if (i<j) - return A_(i,j); - else if (i==j) - return one; - else - return zero; - } - - -#ifdef TNT_USE_REGIONS - // These are the "index-aware" features - - typedef const_Region2D< UnitUpperTriangularView<MaTRiX> > - const_Region; - - const_Region operator()(const Index1D &I, - const Index1D &J) const - { - return const_Region(*this, I, J); - } - - const_Region operator()(Subscript i1, Subscript i2, - Subscript j1, Subscript j2) const - { - return const_Region(*this, i1, i2, j1, j2); - } -#endif -// TNT_USE_REGIONS -}; - -template <class MaTRiX> -UpperTriangularView<MaTRiX> Upper_triangular_view( - /*const*/ MaTRiX &A) -{ - return UpperTriangularView<MaTRiX>(A); -} - - -template <class MaTRiX> -UnitUpperTriangularView<MaTRiX> Unit_upper_triangular_view( - /*const*/ MaTRiX &A) -{ - return UnitUpperTriangularView<MaTRiX>(A); -} - -template <class MaTRiX, class VecToR> -VecToR matmult(/*const*/ UnitUpperTriangularView<MaTRiX> &A, VecToR &x) -{ - Subscript M = A.num_rows(); - Subscript N = A.num_cols(); - - assert(N == x.dim()); - - Subscript i, j; - typename VecToR::element_type sum=0.0; - VecToR result(M); - - Subscript start = A.lbound(); - Subscript Mend = M + A.lbound() -1 ; - - for (i=start; i<=Mend; i++) - { - sum = x(i); - for (j=i+1; j<=N; j++) - sum = sum + A(i,j)*x(j); - result(i) = sum + x(i); - } - - return result; -} - -template <class MaTRiX, class VecToR> -inline VecToR operator*(/*const*/ UnitUpperTriangularView<MaTRiX> &A, VecToR &x) -{ - return matmult(A,x); -} - - -//********************** Algorithms ************************************* - - - -template <class MaTRiX> -std::ostream& operator<<(std::ostream &s, - /*const*/ UpperTriangularView<MaTRiX>&A) -{ - Subscript M=A.num_rows(); - Subscript N=A.num_cols(); - - s << M << " " << N << endl; - - for (Subscript i=1; i<=M; i++) - { - for (Subscript j=1; j<=N; j++) - { - s << A(i,j) << " "; - } - s << endl; - } - - - return s; -} - -template <class MaTRiX> -std::ostream& operator<<(std::ostream &s, - /*const*/ UnitUpperTriangularView<MaTRiX>&A) -{ - Subscript M=A.num_rows(); - Subscript N=A.num_cols(); - - s << M << " " << N << endl; - - for (Subscript i=1; i<=M; i++) - { - for (Subscript j=1; j<=N; j++) - { - s << A(i,j) << " "; - } - s << endl; - } - - - return s; -} - -} // namespace TNT - -#endif //TRIANG_H - diff --git a/intern/iksolver/intern/TNT/trisolve.h b/intern/iksolver/intern/TNT/trisolve.h deleted file mode 100644 index e6e2a0afe3a..00000000000 --- a/intern/iksolver/intern/TNT/trisolve.h +++ /dev/null @@ -1,188 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -// Triangular Solves - -#ifndef TRISLV_H -#define TRISLV_H - - -#include "triang.h" - -namespace TNT -{ - -template <class MaTriX, class VecToR> -VecToR Lower_triangular_solve(/*const*/ MaTriX &A, /*const*/ VecToR &b) -{ - Subscript N = A.num_rows(); - - // make sure matrix sizes agree; A must be square - - assert(A.num_cols() == N); - assert(b.dim() == N); - - VecToR x(N); - - Subscript i; - for (i=1; i<=N; i++) - { - typename MaTriX::element_type tmp=0; - - for (Subscript j=1; j<i; j++) - tmp = tmp + A(i,j)*x(j); - - x(i) = (b(i) - tmp)/ A(i,i); - } - - return x; -} - - -template <class MaTriX, class VecToR> -VecToR Unit_lower_triangular_solve(/*const*/ MaTriX &A, /*const*/ VecToR &b) -{ - Subscript N = A.num_rows(); - - // make sure matrix sizes agree; A must be square - - assert(A.num_cols() == N); - assert(b.dim() == N); - - VecToR x(N); - - Subscript i; - for (i=1; i<=N; i++) - { - - typename MaTriX::element_type tmp=0; - - for (Subscript j=1; j<i; j++) - tmp = tmp + A(i,j)*x(j); - - x(i) = b(i) - tmp; - } - - return x; -} - - -template <class MaTriX, class VecToR> -VecToR linear_solve(/*const*/ LowerTriangularView<MaTriX> &A, - /*const*/ VecToR &b) -{ - return Lower_triangular_solve(A, b); -} - -template <class MaTriX, class VecToR> -VecToR linear_solve(/*const*/ UnitLowerTriangularView<MaTriX> &A, - /*const*/ VecToR &b) -{ - return Unit_lower_triangular_solve(A, b); -} - - - -//********************** Upper triangular section **************** - -template <class MaTriX, class VecToR> -VecToR Upper_triangular_solve(/*const*/ MaTriX &A, /*const*/ VecToR &b) -{ - Subscript N = A.num_rows(); - - // make sure matrix sizes agree; A must be square - - assert(A.num_cols() == N); - assert(b.dim() == N); - - VecToR x(N); - - Subscript i; - for (i=N; i>=1; i--) - { - - typename MaTriX::element_type tmp=0; - - for (Subscript j=i+1; j<=N; j++) - tmp = tmp + A(i,j)*x(j); - - x(i) = (b(i) - tmp)/ A(i,i); - } - - return x; -} - - -template <class MaTriX, class VecToR> -VecToR Unit_upper_triangular_solve(/*const*/ MaTriX &A, /*const*/ VecToR &b) -{ - Subscript N = A.num_rows(); - - // make sure matrix sizes agree; A must be square - - assert(A.num_cols() == N); - assert(b.dim() == N); - - VecToR x(N); - - Subscript i; - for (i=N; i>=1; i--) - { - - typename MaTriX::element_type tmp=0; - - for (Subscript j=i+1; j<i; j++) - tmp = tmp + A(i,j)*x(j); - - x(i) = b(i) - tmp; - } - - return x; -} - - -template <class MaTriX, class VecToR> -VecToR linear_solve(/*const*/ UpperTriangularView<MaTriX> &A, - /*const*/ VecToR &b) -{ - return Upper_triangular_solve(A, b); -} - -template <class MaTriX, class VecToR> -VecToR linear_solve(/*const*/ UnitUpperTriangularView<MaTriX> &A, - /*const*/ VecToR &b) -{ - return Unit_upper_triangular_solve(A, b); -} - - -} // namespace TNT - -#endif // TRISLV_H - diff --git a/intern/iksolver/intern/TNT/vec.h b/intern/iksolver/intern/TNT/vec.h deleted file mode 100644 index 5d4ef14bf73..00000000000 --- a/intern/iksolver/intern/TNT/vec.h +++ /dev/null @@ -1,491 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - -// Basic TNT numerical vector (0-based [i] AND 1-based (i) indexing ) -// - -#ifndef VEC_H -#define VEC_H - -#include "subscript.h" -#include <stdlib.h> -#include <assert.h> -#include <iostream> - -namespace TNT -{ - -template <class T> -class Vector -{ - - - public: - - typedef Subscript size_type; - typedef T value_type; - typedef T element_type; - typedef T* pointer; - typedef T* iterator; - typedef T& reference; - typedef const T* const_iterator; - typedef const T& const_reference; - - Subscript lbound() const { return 1;} - - protected: - T* v_; - T* vm1_; // pointer adjustment for optimzied 1-offset indexing - Subscript n_; - - // internal helper function to create the array - // of row pointers - - void initialize(Subscript N) - { - // adjust pointers so that they are 1-offset: - // v_[] is the internal contiguous array, it is still 0-offset - // - assert(v_ == NULL); - v_ = new T[N]; - assert(v_ != NULL); - vm1_ = v_-1; - n_ = N; - } - - void copy(const T* v) - { - Subscript N = n_; - Subscript i; - -#ifdef TNT_UNROLL_LOOPS - Subscript Nmod4 = N & 3; - Subscript N4 = N - Nmod4; - - for (i=0; i<N4; i+=4) - { - v_[i] = v[i]; - v_[i+1] = v[i+1]; - v_[i+2] = v[i+2]; - v_[i+3] = v[i+3]; - } - - for (i=N4; i< N; i++) - v_[i] = v[i]; -#else - - for (i=0; i< N; i++) - v_[i] = v[i]; -#endif - } - - void set(const T& val) - { - Subscript N = n_; - Subscript i; - -#ifdef TNT_UNROLL_LOOPS - Subscript Nmod4 = N & 3; - Subscript N4 = N - Nmod4; - - for (i=0; i<N4; i+=4) - { - v_[i] = val; - v_[i+1] = val; - v_[i+2] = val; - v_[i+3] = val; - } - - for (i=N4; i< N; i++) - v_[i] = val; -#else - - for (i=0; i< N; i++) - v_[i] = val; - -#endif - } - - - - void destroy() - { - /* do nothing, if no memory has been previously allocated */ - if (v_ == NULL) return ; - - /* if we are here, then matrix was previously allocated */ - delete [] (v_); - - v_ = NULL; - vm1_ = NULL; - } - - - public: - - // access - - iterator begin() { return v_;} - iterator end() { return v_ + n_; } - const iterator begin() const { return v_;} - const iterator end() const { return v_ + n_; } - - // destructor - - ~Vector() - { - destroy(); - } - - // constructors - - Vector() : v_(0), vm1_(0), n_(0) {} - - Vector(const Vector<T> &A) : v_(0), vm1_(0), n_(0) - { - initialize(A.n_); - copy(A.v_); - } - - Vector(Subscript N, const T& value = T()) : v_(0), vm1_(0), n_(0) - { - initialize(N); - set(value); - } - - Vector(Subscript N, const T* v) : v_(0), vm1_(0), n_(0) - { - initialize(N); - copy(v); - } - - - // methods - // - Vector<T>& newsize(Subscript N) - { - if (n_ == N) return *this; - - destroy(); - initialize(N); - - return *this; - } - - - // assignments - // - Vector<T>& operator=(const Vector<T> &A) - { - if (v_ == A.v_) - return *this; - - if (n_ == A.n_) // no need to re-alloc - copy(A.v_); - - else - { - destroy(); - initialize(A.n_); - copy(A.v_); - } - - return *this; - } - - Vector<T>& operator=(const T& scalar) - { - set(scalar); - return *this; - } - - inline Subscript dim() const - { - return n_; - } - - inline Subscript size() const - { - return n_; - } - - - inline reference operator()(Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= n_) ; -#endif - return vm1_[i]; - } - - inline const_reference operator() (Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i <= n_) ; -#endif - return vm1_[i]; - } - - inline reference operator[](Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(0<=i); - assert(i < n_) ; -#endif - return v_[i]; - } - - inline const_reference operator[](Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(0<=i) ; - assert(i < n_) ; -#endif - return v_[i]; - } - - - -}; - - -/* *************************** I/O ********************************/ - -template <class T> -std::ostream& operator<<(std::ostream &s, const Vector<T> &A) -{ - Subscript N=A.dim(); - - s << N << std::endl; - - for (Subscript i=0; i<N; i++) - s << A[i] << " " << std::endl; - s << std::endl; - - return s; -} - -template <class T> -std::istream & operator>>(std::istream &s, Vector<T> &A) -{ - - Subscript N; - - s >> N; - - if ( !(N == A.size() )) - { - A.newsize(N); - } - - - for (Subscript i=0; i<N; i++) - s >> A[i]; - - - return s; -} - -// *******************[ basic matrix algorithms ]*************************** - - -template <class T> -Vector<T> operator+(const Vector<T> &A, - const Vector<T> &B) -{ - Subscript N = A.dim(); - - assert(N==B.dim()); - - Vector<T> tmp(N); - Subscript i; - - for (i=0; i<N; i++) - tmp[i] = A[i] + B[i]; - - return tmp; -} - -template <class T> -Vector<T> operator-(const Vector<T> &A, - const Vector<T> &B) -{ - Subscript N = A.dim(); - - assert(N==B.dim()); - - Vector<T> tmp(N); - Subscript i; - - for (i=0; i<N; i++) - tmp[i] = A[i] - B[i]; - - return tmp; -} - -template <class T> -Vector<T> operator*(const Vector<T> &A, - const Vector<T> &B) -{ - Subscript N = A.dim(); - - assert(N==B.dim()); - - Vector<T> tmp(N); - Subscript i; - - for (i=0; i<N; i++) - tmp[i] = A[i] * B[i]; - - return tmp; -} - -template <class T> -Vector<T> operator*(const Vector<T> &A, - const T &B) -{ - Subscript N = A.dim(); - - Vector<T> tmp(N); - Subscript i; - - for (i=0; i<N; i++) - tmp[i] = A[i] * B; - - return tmp; -} - - -template <class T> -T dot_prod(const Vector<T> &A, const Vector<T> &B) -{ - Subscript N = A.dim(); - assert(N == B.dim()); - - Subscript i; - T sum = 0; - - for (i=0; i<N; i++) - sum += A[i] * B[i]; - - return sum; -} - -// inplace versions of the above template functions - -// A = A + B - -template <class T> -void vectoradd( - Vector<T> &A, - const Vector<T> &B) -{ - const Subscript N = A.dim(); - assert(N==B.dim()); - Subscript i; - - for (i=0; i<N; i++) - A[i] += B[i]; -} - -// same with separate output vector - -template <class T> -void vectoradd( - Vector<T> &C, - const Vector<T> &A, - const Vector<T> &B) -{ - const Subscript N = A.dim(); - assert(N==B.dim()); - Subscript i; - - for (i=0; i<N; i++) - C[i] = A[i] + B[i]; -} - -// A = A - B - -template <class T> -void vectorsub( - Vector<T> &A, - const Vector<T> &B) -{ - const Subscript N = A.dim(); - assert(N==B.dim()); - Subscript i; - - for (i=0; i<N; i++) - A[i] -= B[i]; -} - -template <class T> -void vectorsub( - Vector<T> &C, - const Vector<T> &A, - const Vector<T> &B) -{ - const Subscript N = A.dim(); - assert(N==B.dim()); - Subscript i; - - for (i=0; i<N; i++) - C[i] = A[i] - B[i]; -} - -template <class T> -void vectorscale( - Vector<T> &C, - const Vector<T> &A, - const T &B) -{ - const Subscript N = A.dim(); - Subscript i; - - for (i=0; i<N; i++) - C[i] = A[i] * B; -} - -template <class T> -void vectorscale( - Vector<T> &A, - const T &B) -{ - const Subscript N = A.dim(); - Subscript i; - - for (i=0; i<N; i++) - A[i] *= B; -} - -} /* namespace TNT */ - -#endif // VEC_H - diff --git a/intern/iksolver/intern/TNT/vecadaptor.h b/intern/iksolver/intern/TNT/vecadaptor.h deleted file mode 100644 index 77d4e2c05a4..00000000000 --- a/intern/iksolver/intern/TNT/vecadaptor.h +++ /dev/null @@ -1,284 +0,0 @@ -/** - */ - -/* - -* -* Template Numerical Toolkit (TNT): Linear Algebra Module -* -* Mathematical and Computational Sciences Division -* National Institute of Technology, -* Gaithersburg, MD USA -* -* -* This software was developed at the National Institute of Standards and -* Technology (NIST) by employees of the Federal Government in the course -* of their official duties. Pursuant to title 17 Section 105 of the -* United States Code, this software is not subject to copyright protection -* and is in the public domain. The Template Numerical Toolkit (TNT) is -* an experimental system. NIST assumes no responsibility whatsoever for -* its use by other parties, and makes no guarantees, expressed or implied, -* about its quality, reliability, or any other characteristic. -* -* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -* see http://math.nist.gov/tnt for latest updates. -* -*/ - - - -#ifndef VECADAPTOR_H -#define VECADAPTOR_H - -#include <cstdlib> -#include <iostream> -#include <cassert> - -#include "subscript.h" - -#ifdef TNT_USE_REGIONS -#include "region1d.h" -#endif - -namespace TNT -{ - -// see "tntreq.h" for TNT requirements for underlying vector -// class. This need NOT be the STL vector<> class, but a subset -// that provides minimal services. -// -// This is a container adaptor that provides the following services. -// -// o) adds 1-offset operator() access ([] is always 0 offset) -// o) adds TNT_BOUNDS_CHECK to () and [] -// o) adds initialization from strings, e.g. "1.0 2.0 3.0"; -// o) adds newsize(N) function (does not preserve previous values) -// o) adds dim() and dim(1) -// o) adds free() function to release memory used by vector -// o) adds regions, e.g. A(Index(1,10)) = ... -// o) add getVector() method to return adapted container -// o) adds simple I/O for ostreams - -template <class BBVec> -class Vector_Adaptor -{ - - public: - typedef typename BBVec::value_type T; - typedef T value_type; - typedef T element_type; - typedef T* pointer; - typedef T* iterator; - typedef T& reference; - typedef const T* const_iterator; - typedef const T& const_reference; - - Subscript lbound() const { return 1; } - - protected: - BBVec v_; - T* vm1_; - - public: - - Subscript size() const { return v_.size(); } - - // These were removed so that the ANSI C++ valarray class - // would work as a possible storage container. - // - // - //iterator begin() { return v_.begin();} - //iterator begin() { return &v_[0];} - // - //iterator end() { return v_.end(); } - //iterator end() { return &v_[0] + v_.size(); } - // - //const_iterator begin() const { return v_.begin();} - //const_iterator begin() const { return &v_[0];} - // - //const_iterator end() const { return v_.end(); } - //const_iterator end() const { return &v_[0] + v_.size(); } - - BBVec& getVector() { return v_; } - Subscript dim() const { return v_.size(); } - Subscript dim(Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(i==TNT_BASE_OFFSET); -#endif - return (i==TNT_BASE_OFFSET ? v_.size() : 0 ); - } - Vector_Adaptor() : v_() {}; - Vector_Adaptor(const Vector_Adaptor<BBVec> &A) : v_(A.v_) - { - vm1_ = ( v_.size() > 0 ? &(v_[0]) -1 : NULL); - - } - - Vector_Adaptor(Subscript N, const T& value = T()) : v_(N) - { - for (Subscript i=0; i<N; i++) - v_[i] = value; - - vm1_ = ( v_.size() > 0 ? &(v_[0]) -1 : NULL); - } - - Vector_Adaptor(Subscript N, const T* values) : v_(N) - { - for (Subscript i=0; i<N; i++) - v_[i] = values[i]; - vm1_ = ( v_.size() > 0 ? &(v_[0]) -1 : NULL); - } - Vector_Adaptor(const BBVec & A) : v_(A) - { - vm1_ = ( v_.size() > 0 ? &(v_[0]) -1 : NULL); - } - - // NOTE: this assumes that BBVec(0) constructor creates an - // null vector that does not take up space... It would be - // great to require that BBVec have a corresponding free() - // function, but in particular STL vectors do not. - // - Vector_Adaptor<BBVec>& free() - { - return *this = Vector_Adaptor<BBVec>(0); - } - - Vector_Adaptor<BBVec>& operator=(const Vector_Adaptor<BBVec> &A) - { - v_ = A.v_ ; - vm1_ = ( v_.size() > 0 ? &(v_[0]) -1 : NULL); - return *this; - } - - Vector_Adaptor<BBVec>& newsize(Subscript N) - { - // NOTE: this is not as efficient as it could be - // but to retain compatiblity with STL interface - // we cannot assume underlying implementation - // has a newsize() function. - - return *this = Vector_Adaptor<BBVec>(N); - - } - - Vector_Adaptor<BBVec>& operator=(const T &a) - { - Subscript i; - Subscript N = v_.size(); - for (i=0; i<N; i++) - v_[i] = a; - - return *this; - } - - Vector_Adaptor<BBVec>& resize(Subscript N) - { - if (N == size()) return *this; - - Vector_Adaptor<BBVec> tmp(N); - Subscript n = (N < size() ? N : size()); // min(N, size()); - Subscript i; - - for (i=0; i<n; i++) - tmp[i] = v_[i]; - - - return (*this = tmp); - - } - - - reference operator()(Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i<=dim()); -#endif - return vm1_[i]; - } - - const_reference operator()(Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(1<=i); - assert(i<=dim()); -#endif - return vm1_[i]; - } - - reference operator[](Subscript i) - { -#ifdef TNT_BOUNDS_CHECK - assert(0<=i); - assert(i<dim()); -#endif - return v_[i]; - } - - const_reference operator[](Subscript i) const - { -#ifdef TNT_BOUNDS_CHECK - assert(0<=i); - assert(i<dim()); -#endif - return v_[i]; - } - - -#ifdef TNT_USE_REGIONS - // "index-aware" features, all of these are 1-based offsets - - typedef Region1D<Vector_Adaptor<BBVec> > Region; - - typedef const_Region1D< Vector_Adaptor<BBVec> > const_Region; - - Region operator()(const Index1D &I) - { return Region(*this, I); } - - Region operator()(const Subscript i1, Subscript i2) - { return Region(*this, i1, i2); } - - const_Region operator()(const Index1D &I) const - { return const_Region(*this, I); } - - const_Region operator()(const Subscript i1, Subscript i2) const - { return const_Region(*this, i1, i2); } -#endif -// TNT_USE_REGIONS - - -}; - -#include <iostream> - -template <class BBVec> -std::ostream& operator<<(std::ostream &s, const Vector_Adaptor<BBVec> &A) -{ - Subscript M=A.size(); - - s << M << endl; - for (Subscript i=1; i<=M; i++) - s << A(i) << endl; - return s; -} - -template <class BBVec> -std::istream& operator>>(std::istream &s, Vector_Adaptor<BBVec> &A) -{ - Subscript N; - - s >> N; - - A.resize(N); - - for (Subscript i=1; i<=N; i++) - s >> A(i); - - return s; -} - -} // namespace TNT - -#endif - diff --git a/intern/iksolver/intern/TNT/version.h b/intern/iksolver/intern/TNT/version.h deleted file mode 100644 index 58017dd80f9..00000000000 --- a/intern/iksolver/intern/TNT/version.h +++ /dev/null @@ -1,25 +0,0 @@ -/** - */ - -// Template Numerical Toolkit (TNT) for Linear Algebra - -// -// BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE -// Please see http://math.nist.gov/tnt for updates -// -// R. Pozo -// Mathematical and Computational Sciences Division -// National Institute of Standards and Technology - - -#ifndef TNT_VERSION_H -#define TNT_VERSION_H - - -#define TNT_MAJOR_VERSION '0' -#define TNT_MINOR_VERSION '9' -#define TNT_SUBMINOR_VERSION '4' -#define TNT_VERSION_STRING "0.9.4" - -#endif - |