IK Solver: replace TNT math library with Eigen.

Performance is about the same or slightly better for typical IK chains.
In extreme cases with many bones and multiple targets, of which some are
unreachable, I've seen 2x speedups.

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
-