Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
path: root/intern/iksolver/intern/TNT/triang.h
diff options
context:
space:
mode:
Diffstat (limited to 'intern/iksolver/intern/TNT/triang.h')
-rw-r--r--intern/iksolver/intern/TNT/triang.h670
1 files changed, 670 insertions, 0 deletions
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
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-17 06:53:26 +0300