diff options
Diffstat (limited to 'intern/iksolver/intern/TNT')
24 files changed, 6737 insertions, 0 deletions
diff --git a/intern/iksolver/intern/TNT/cholesky.h b/intern/iksolver/intern/TNT/cholesky.h new file mode 100644 index 00000000000..a3e0118f4e5 --- /dev/null +++ b/intern/iksolver/intern/TNT/cholesky.h @@ -0,0 +1,128 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..80a82417e47 --- /dev/null +++ b/intern/iksolver/intern/TNT/cmat.h @@ -0,0 +1,661 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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> +#include <strstream> +#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); + } + + Matrix(Subscript M, Subscript N, const char *s) + { + initialize(M,N); + std::istrstream ins(s); + + Subscript i, j; + + for (i=0; i<M; i++) + for (j=0; j<N; j++) + ins >> row_[i][j]; + } + + + // 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 new file mode 100644 index 00000000000..877364ad38c --- /dev/null +++ b/intern/iksolver/intern/TNT/fcscmat.h @@ -0,0 +1,197 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..5de9a63813e --- /dev/null +++ b/intern/iksolver/intern/TNT/fmat.h @@ -0,0 +1,609 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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> +#include <strstream> +#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); + } + + + Fortran_Matrix(Subscript M, Subscript N, char *s) + { + initialize(M,N); + std::istrstream ins(s); + + Subscript i, j; + + for (i=1; i<=M; i++) + for (j=1; j<=N; j++) + ins >> (*this)(i,j); + } + + // 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 new file mode 100644 index 00000000000..4f513300588 --- /dev/null +++ b/intern/iksolver/intern/TNT/fortran.h @@ -0,0 +1,99 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..b8dde35b6c6 --- /dev/null +++ b/intern/iksolver/intern/TNT/fspvec.h @@ -0,0 +1,200 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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> +#include <strstream> + +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 new file mode 100644 index 00000000000..885545d11b5 --- /dev/null +++ b/intern/iksolver/intern/TNT/index.h @@ -0,0 +1,115 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..d3556a96071 --- /dev/null +++ b/intern/iksolver/intern/TNT/lapack.h @@ -0,0 +1,225 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..b86027aa386 --- /dev/null +++ b/intern/iksolver/intern/TNT/lu.h @@ -0,0 +1,236 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..074551896b9 --- /dev/null +++ b/intern/iksolver/intern/TNT/qr.h @@ -0,0 +1,261 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..0307ee02c39 --- /dev/null +++ b/intern/iksolver/intern/TNT/region1d.h @@ -0,0 +1,403 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..d612dfdaa79 --- /dev/null +++ b/intern/iksolver/intern/TNT/region2d.h @@ -0,0 +1,499 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..6079ec6f1b1 --- /dev/null +++ b/intern/iksolver/intern/TNT/stopwatch.h @@ -0,0 +1,115 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..4728a79a168 --- /dev/null +++ b/intern/iksolver/intern/TNT/subscript.h @@ -0,0 +1,90 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..9be2d98bd9f --- /dev/null +++ b/intern/iksolver/intern/TNT/svd.h @@ -0,0 +1,507 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +#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). +// work_space is a temporary vector used by this class to compute +// intermediate values during the computation of the SVD. This should +// be of length a.num_cols. This is not checked + +#include "tntmath.h" + +namespace TNT +{ + + +template <class MaTRiX, class VecToR > +void SVD(MaTRiX &a, VecToR &w, MaTRiX &v, VecToR &work_space) { + + int n = a.num_cols(); + int m = a.num_rows(); + + int flag,i,its,j,jj,k,l(0),nm(0); + typename MaTRiX::value_type c,f,h,x,y,z; + typename MaTRiX::value_type anorm(0),g(0),scale(0); + typename MaTRiX::value_type s(0); + + work_space.newsize(n); + + for (i=1;i<=n;i++) { + l=i+1; + work_space(i)=scale*g; + + g = (typename MaTRiX::value_type)0; + + s = (typename MaTRiX::value_type)0; + scale = (typename MaTRiX::value_type)0; + + if (i <= m) { + for (k=i;k<=m;k++) scale += TNT::abs(a(k,i)); + if (scale > (typename MaTRiX::value_type)0) { + for (k=i;k<=m;k++) { + a(k,i) /= scale; + s += a(k,i)*a(k,i); + } + f=a(i,i); + g = -TNT::sign(sqrt(s),f); + h=f*g-s; + a(i,i)=f-g; + if (i != n) { + for (j=l;j<=n;j++) { + s = (typename MaTRiX::value_type)0; + for (k=i;k<=m;k++) s += a(k,i)*a(k,j); + f=s/h; + for (k=i;k<=m;k++) a(k,j) += f*a(k,i); + } + } + for (k=i;k<=m;k++) a(k,i) *= scale; + } + } + w(i)=scale*g; + g = (typename MaTRiX::value_type)0; + s = (typename MaTRiX::value_type)0; + scale = (typename MaTRiX::value_type)0; + if (i <= m && i != n) { + for (k=l;k<=n;k++) scale += TNT::abs(a(i,k)); + if (scale > (typename MaTRiX::value_type)0) { + for (k=l;k<=n;k++) { + a(i,k) /= scale; + s += a(i,k)*a(i,k); + } + f=a(i,l); + g = -TNT::sign(sqrt(s),f); + h=f*g-s; + a(i,l)=f-g; + for (k=l;k<=n;k++) work_space(k)=a(i,k)/h; + if (i != m) { + for (j=l;j<=m;j++) { + s = (typename MaTRiX::value_type)0; + for (k=l;k<=n;k++) s += a(j,k)*a(i,k); + for (k=l;k<=n;k++) a(j,k) += s*work_space(k); + } + } + for (k=l;k<=n;k++) a(i,k) *= scale; + } + } + anorm=TNT::max(anorm,(TNT::abs(w(i))+TNT::abs(work_space(i)))); + } + for (i=n;i>=1;i--) { + if (i < n) { + if (g != (typename MaTRiX::value_type)0) { + for (j=l;j<=n;j++) + v(j,i)=(a(i,j)/a(i,l))/g; + for (j=l;j<=n;j++) { + s = (typename MaTRiX::value_type)0; + for (k=l;k<=n;k++) s += a(i,k)*v(k,j); + for (k=l;k<=n;k++) v(k,j) += s*v(k,i); + } + } + for (j=l;j<=n;j++) v(i,j)=v(j,i)= (typename MaTRiX::value_type)0; + } + v(i,i)= (typename MaTRiX::value_type)1; + g=work_space(i); + l=i; + } + for (i=n;i>=1;i--) { + l=i+1; + g=w(i); + if (i < n) { + for (j=l;j<=n;j++) a(i,j)= (typename MaTRiX::value_type)0; + } + if (g != (typename MaTRiX::value_type)0) { + g= ((typename MaTRiX::value_type)1)/g; + if (i != n) { + for (j=l;j<=n;j++) { + s = (typename MaTRiX::value_type)0; + for (k=l;k<=m;k++) s += a(k,i)*a(k,j); + f=(s/a(i,i))*g; + for (k=i;k<=m;k++) a(k,j) += f*a(k,i); + } + } + for (j=i;j<=m;j++) a(j,i) *= g; + } else { + for (j=i;j<=m;j++) a(j,i)= (typename MaTRiX::value_type)0; + } + ++a(i,i); + } + for (k=n;k>=1;k--) { + for (its=1;its<=30;its++) { + flag=1; + for (l=k;l>=1;l--) { + nm=l-1; + if (TNT::abs(work_space(l))+anorm == anorm) { + flag=0; + break; + } + if (TNT::abs(w(nm))+anorm == anorm) break; + } + if (flag) { + c= (typename MaTRiX::value_type)0; + s= (typename MaTRiX::value_type)1; + for (i=l;i<=k;i++) { + f=s*work_space(i); + if (TNT::abs(f)+anorm != anorm) { + g=w(i); + h= (typename MaTRiX::value_type)TNT::pythag(float(f),float(g)); + w(i)=h; + h= ((typename MaTRiX::value_type)1)/h; + c=g*h; + s=(-f*h); + for (j=1;j<=m;j++) { + y=a(j,nm); + z=a(j,i); + a(j,nm)=y*c+z*s; + a(j,i)=z*c-y*s; + } + } + } + } + z=w(k); + if (l == k) { + if (z < (typename MaTRiX::value_type)0) { + w(k) = -z; + for (j=1;j<=n;j++) v(j,k)=(-v(j,k)); + } + break; + } + + +#if 1 + if (its == 30) + { + TNTException an_exception; + an_exception.i = 0; + throw an_exception; + + return ; + assert(false); + } +#endif + x=w(l); + nm=k-1; + y=w(nm); + g=work_space(nm); + h=work_space(k); + f=((y-z)*(y+z)+(g-h)*(g+h))/(((typename MaTRiX::value_type)2)*h*y); + g=(typename MaTRiX::value_type)TNT::pythag(float(f), float(1)); + f=((x-z)*(x+z)+h*((y/(f+TNT::sign(g,f)))-h))/x; + c = (typename MaTRiX::value_type)1; + s = (typename MaTRiX::value_type)1; + for (j=l;j<=nm;j++) { + i=j+1; + g=work_space(i); + y=w(i); + h=s*g; + g=c*g; + z=(typename MaTRiX::value_type)TNT::pythag(float(f),float(h)); + work_space(j)=z; + c=f/z; + s=h/z; + f=x*c+g*s; + g=g*c-x*s; + h=y*s; + y=y*c; + for (jj=1;jj<=n;jj++) { + x=v(jj,j); + z=v(jj,i); + v(jj,j)=x*c+z*s; + v(jj,i)=z*c-x*s; + } + z=(typename MaTRiX::value_type)TNT::pythag(float(f),float(h)); + w(j)=z; + if (z != (typename MaTRiX::value_type)0) { + z= ((typename MaTRiX::value_type)1)/z; + c=f*z; + s=h*z; + } + f=(c*g)+(s*y); + x=(c*y)-(s*g); + for (jj=1;jj<=m;jj++) { + y=a(jj,j); + z=a(jj,i); + a(jj,j)=y*c+z*s; + a(jj,i)=z*c-y*s; + } + } + work_space(l)= (typename MaTRiX::value_type)0; + work_space(k)=f; + w(k)=x; + } + } +}; + + + +// A is replaced by the column orthogonal matrix U + + +template <class MaTRiX, class VecToR > +void SVD_a( MaTRiX &a, VecToR &w, MaTRiX &v) { + + int n = a.num_cols(); + int m = a.num_rows(); + + int flag,i,its,j,jj,k,l,nm; + typename MaTRiX::value_type anorm,c,f,g,h,s,scale,x,y,z; + + VecToR work_space; + work_space.newsize(n); + + g = scale = anorm = 0.0; + + for (i=1;i <=n;i++) { + l = i+1; + work_space(i) = scale*g; + g = s=scale=0.0; + + if (i <= m) { + for(k=i; k<=m; k++) scale += abs(a(k,i)); + + if (scale) { + for (k = i; k <=m ; k++) { + a(k,i) /= scale; + s += a(k,i)*a(k,i); + } + f = a(i,i); + g = -sign(sqrt(s),f); + h = f*g -s; + a(i,i) = f-g; + + for (j = l; j <=n; j++) { + for (s = 0.0,k =i;k<=m;k++) s += a(k,i)*a(k,j); + f = s/h; + for (k = i; k <= m; k++) a(k,j) += f*a(k,i); + } + for (k = i; k <=m;k++) a(k,i) *= scale; + } + } + + w(i) = scale*g; + g = s = scale = 0.0; + + if (i <=m && i != n) { + for (k = l; k <=n;k++) scale += abs(a(i,k)); + if (scale) { + for(k = l;k <=n;k++) { + a(i,k) /= scale; + s += a(i,k) * a(i,k); + } + + f = a(i,l); + g = -sign(sqrt(s),f); + h= f*g -s; + a(i,l) = f-g; + for(k=l;k<=n;k++) work_space(k) = a(i,k)/h; + for (j=l;j<=m;j++) { + for(s = 0.0,k=l;k<=n;k++) s+= a(j,k)*a(i,k); + for(k=l;k<=n;k++) a(j,k) += s*work_space(k); + } + for(k=l;k<=n;k++) a(i,k)*=scale; + } + } + anorm = max(anorm,(abs(w(i)) + abs(work_space(i)))); + } + for (i=n;i>=1;i--) { + if (i <n) { + if (g) { + for(j=l;j<=n;j++) v(j,i) = (a(i,j)/a(i,l))/g; + for(j=l;j<=n;j++) { + for(s=0.0,k=l;k<=n;k++) s += a(i,k)*v(k,j); + for(k=l; k<=n;k++) v(k,j) +=s*v(k,i); + } + } + for(j=l;j <=n;j++) v(i,j) = v(j,i) = 0.0; + } + v(i,i) = 1.0; + g = work_space(i); + l = i; + } + + for (i = min(m,n);i>=1;i--) { + l = i+1; + g = w(i); + for (j=l;j <=n;j++) a(i,j) = 0.0; + if (g) { + g = 1.0/g; + for (j=l;j<=n;j++) { + for (s = 0.0,k=l;k<=m;k++) s += a(k,i)*a(k,j); + f = (s/a(i,i))*g; + for (k=i;k<=m;k++) a(k,j) += f*a(k,i); + } + for (j=i;j<=m;j++) a(j,i)*=g; + } else { + for (j=i;j<=m;j++) a(j,i) = 0.0; + } + ++a(i,i); + } + + for (k=n;k>=1;k--) { + for (its=1;its<=30;its++) { + flag=1; + for(l=k;l>=1;l--) { + nm = l-1; + if (abs(work_space(l)) + anorm == anorm) { + flag = 0; + break; + } + if (abs(w(nm)) + anorm == anorm) break; + } + if (flag) { + c = 0.0; + s = 1.0; + for (i=l;i<=k;i++) { + f = s*work_space(i); + work_space(i) = c*work_space(i); + if (abs(f) +anorm == anorm) break; + g = w(i); + h = pythag(f,g); + w(i) = h; + h = 1.0/h; + c = g*h; + s = -f*h; + for (j=1;j<=m;j++) { + y= a(j,nm); + z=a(j,i); + a(j,nm) = y*c + z*s; + a(j,i) = z*c - y*s; + } + } + } + z=w(k); + if (l==k) { + if (z <0.0) { + w(k) = -z; + for (j=1;j<=n;j++) v(j,k) = -v(j,k); + } + break; + } + + if (its == 30) assert(false); + + x=w(l); + nm=k-1; + y=w(nm); + g=work_space(nm); + h=work_space(k); + + f= ((y-z)*(y+z) + (g-h)*(g+h))/(2.0*h*y); + g = pythag(f,1.0); + f= ((x-z)*(x+z) + h*((y/(f + sign(g,f)))-h))/x; + c=s=1.0; + + for (j=l;j<=nm;j++) { + i=j+1; + g = work_space(i); + y=w(i); + h=s*g; + g=c*g; + z=pythag(f,h); + work_space(j) = z; + c=f/z; + s=h/z; + f=x*c + g*s; + g= g*c - x*s; + h=y*s; + y*=c; + for(jj=1;jj<=n;jj++) { + x=v(jj,j); + z=v(jj,i); + v(jj,j) = x*c + z*s; + v(jj,i) = z*c- x*s; + } + z=pythag(f,h); + w(j)=z; + if(z) { + z = 1.0/z; + c=f*z; + s=h*z; + } + f=c*g + s*y; + x= c*y-s*g; + + for(jj=1;jj<=m;jj++) { + y=a(jj,j); + z=a(jj,i); + a(jj,j) = y*c+z*s; + a(jj,i) = z*c - y*s; + } + } + + work_space(l) = 0.0; + work_space(k) = f; + w(k) = x; + } + } +} +} +#endif + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/intern/iksolver/intern/TNT/tnt.h b/intern/iksolver/intern/TNT/tnt.h new file mode 100644 index 00000000000..df27079f87a --- /dev/null +++ b/intern/iksolver/intern/TNT/tnt.h @@ -0,0 +1,127 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..3342ff7f5c6 --- /dev/null +++ b/intern/iksolver/intern/TNT/tntmath.h @@ -0,0 +1,177 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 + + + +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 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 new file mode 100644 index 00000000000..c935902554e --- /dev/null +++ b/intern/iksolver/intern/TNT/tntreqs.h @@ -0,0 +1,102 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..47730350acc --- /dev/null +++ b/intern/iksolver/intern/TNT/transv.h @@ -0,0 +1,192 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..883caf138f7 --- /dev/null +++ b/intern/iksolver/intern/TNT/triang.h @@ -0,0 +1,670 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..6a4f6e1a5b3 --- /dev/null +++ b/intern/iksolver/intern/TNT/trisolve.h @@ -0,0 +1,216 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 new file mode 100644 index 00000000000..beb034ce315 --- /dev/null +++ b/intern/iksolver/intern/TNT/vec.h @@ -0,0 +1,535 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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> +#include <strstream> + +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); + } + + Vector(Subscript N, char *s) : v_(0), vm1_(0), n_(0) + { + initialize(N); + std::istrstream ins(s); + + Subscript i; + + for (i=0; i<N; i++) + ins >> v_[i]; + } + + + // 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 << endl; + + for (Subscript i=0; i<N; i++) + s << A[i] << " " << endl; + s << 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 seperate 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 new file mode 100644 index 00000000000..337748c58bd --- /dev/null +++ b/intern/iksolver/intern/TNT/vecadaptor.h @@ -0,0 +1,321 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +/* + +* +* 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 <strstream> +#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*/ char *s) : v_(N) + { + istrstream ins(s); + for (Subscript i=0; i<N; i++) + ins >> v_[i] ; + + 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 new file mode 100644 index 00000000000..c2759e0ce44 --- /dev/null +++ b/intern/iksolver/intern/TNT/version.h @@ -0,0 +1,52 @@ +/** + * $Id$ + * ***** BEGIN GPL/BL DUAL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. The Blender + * Foundation also sells licenses for use in proprietary software under + * the Blender License. See http://www.blender.org/BL/ for information + * about this. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL/BL DUAL LICENSE BLOCK ***** + */ + +// 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 |