diff --git a/Base/quicksort.h b/Base/quicksort.h
index 0650d485d43ba06777922311220462d1c333d710..78f69f558560b2bbd54d4b31e202508b929e9982 100644
--- a/Base/quicksort.h
+++ b/Base/quicksort.h
@@ -14,6 +14,36 @@
// Base
#include "swap.h"
+template <class T>
+void quickSort(T* array, unsigned beg, unsigned end)
+{
+ if (beg < end) {
+ unsigned p = partition_(array, beg, end);
+ quickSort(array, beg, p);
+ quickSort(array, p+1, end);
+ }
+}
+
+template <class T>
+unsigned partition_(T* array, unsigned beg, unsigned end)
+{
+ unsigned i = beg+1;
+ unsigned j = end-1;
+ T m = array[beg];
+
+ for (;;) {
+ while ((i<end) && (array[i] < m)) i++;
+ while ((j>beg) && !(array[j] < m)) j--;
+
+ if (i >= j) break;
+ BASELIB::swap(array[i], array[j]);
+ }
+
+ BASELIB::swap(array[beg], array[j]);
+ return j;
+}
+
+
/**
* version of partition_ that additional updates the permutation vector
* */
diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
index 9dfbbeaeb47d70763f13318d993b46f7e695a25e..cc426573501a7564ecb86a9465baaecbe0c9ac64 100644
--- a/MathLib/CMakeLists.txt
+++ b/MathLib/CMakeLists.txt
@@ -1,49 +1,53 @@
# Source files
-SET( HEADERS
- AnalyticalGeometry.h
- LinearInterpolation.h
- MathTools.h
- Vector3.h
- EarClippingTriangulation.h
- max.h
- sparse.h
- vector_io.h
- LinAlg/MatrixBase.h
- LinAlg/VectorNorms.h
- LinAlg/Dense/Matrix.h
- LinAlg/Preconditioner/generateDiagPrecond.h
- LinAlg/Solvers/LinearSolver.h
- LinAlg/Solvers/DirectLinearSolver.h
- LinAlg/Solvers/DenseDirectLinearSolver.h
- LinAlg/Solvers/GaussAlgorithm.h
- LinAlg/Solvers/TriangularSolve.h
- LinAlg/Solvers/IterativeLinearSolver.h
- LinAlg/Solvers/solver.h
- LinAlg/Solvers/BiCGStab.h
- LinAlg/Solvers/CG.h
- LinAlg/Solvers/GMRes.h
- LinAlg/Sparse/amuxCRS.h
- LinAlg/Sparse/CRSMatrix.h
- LinAlg/Sparse/CRSMatrixPThreads.h
- LinAlg/Sparse/CRSMatrixOpenMP.h
- LinAlg/Sparse/CRSSymMatrix.h
- LinAlg/Sparse/SparseMatrixBase.h
-)
+#SET( HEADERS
+# AnalyticalGeometry.h
+# LinearInterpolation.h
+# MathTools.h
+# Vector3.h
+# EarClippingTriangulation.h
+# max.h
+# sparse.h
+# vector_io.h
+# LinAlg/MatrixBase.h
+# LinAlg/VectorNorms.h
+# LinAlg/Dense/Matrix.h
+# LinAlg/Preconditioner/generateDiagPrecond.h
+# LinAlg/Solvers/LinearSolver.h
+# LinAlg/Solvers/DirectLinearSolver.h
+# LinAlg/Solvers/DenseDirectLinearSolver.h
+# LinAlg/Solvers/GaussAlgorithm.h
+# LinAlg/Solvers/TriangularSolve.h
+# LinAlg/Solvers/IterativeLinearSolver.h
+# LinAlg/Solvers/solver.h
+# LinAlg/Solvers/BiCGStab.h
+# LinAlg/Solvers/CG.h
+# LinAlg/Solvers/GMRes.h
+# LinAlg/Sparse/amuxCRS.h
+# LinAlg/Sparse/CRSMatrix.h
+# LinAlg/Sparse/CRSMatrixPThreads.h
+# LinAlg/Sparse/CRSMatrixOpenMP.h
+# LinAlg/Sparse/CRSSymMatrix.h
+# LinAlg/Sparse/SparseMatrixBase.h
+#)
-SET( SOURCES
- AnalyticalGeometry.cpp
- LinearInterpolation.cpp
- MathTools.cpp
- EarClippingTriangulation.cpp
- LinAlg/Sparse/amuxCRS.cpp
- LinAlg/Solvers/BiCGStab.cpp
- LinAlg/Solvers/CG.cpp
- LinAlg/Solvers/CGParallel.cpp
- LinAlg/Solvers/GMRes.cpp
- LinAlg/Solvers/GaussAlgorithm.cpp
- LinAlg/Solvers/TriangularSolve.cpp
- LinAlg/Preconditioner/generateDiagPrecond.cpp
-)
+FILE(GLOB_RECURSE HEADERS . *.h)
+
+#SET( SOURCES
+# AnalyticalGeometry.cpp
+# LinearInterpolation.cpp
+# MathTools.cpp
+# EarClippingTriangulation.cpp
+# LinAlg/Sparse/amuxCRS.cpp
+# LinAlg/Solvers/BiCGStab.cpp
+# LinAlg/Solvers/CG.cpp
+# LinAlg/Solvers/CGParallel.cpp
+# LinAlg/Solvers/GMRes.cpp
+# LinAlg/Solvers/GaussAlgorithm.cpp
+# LinAlg/Solvers/TriangularSolve.cpp
+# LinAlg/Preconditioner/generateDiagPrecond.cpp
+#)
+
+FILE(GLOB_RECURSE SOURCES . *.cpp)
INCLUDE_DIRECTORIES (
.
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.cpp b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9280eb5bdf118e0ad18946b11ec43b70566b0f5f
--- /dev/null
+++ b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.cpp
@@ -0,0 +1,203 @@
+/*
+ * AdjMat.cpp
+ *
+ * Created on: 02.01.2012
+ * Author: TF
+ */
+
+// Base
+#include "swap.h"
+#include "quicksort.h"
+
+#include "LinAlg/Sparse/NestedDissectionPermutation/AdjMat.h"
+
+namespace MathLib {
+
+AdjMat* AdjMat::getMat(unsigned beg, unsigned end,
+ const unsigned* const op_perm, const unsigned* const po_perm) const
+{
+ const unsigned nsize(end - beg); // size of new matrix
+ unsigned i, c; // row and col idx in permuted matrix
+ unsigned j, idx; // pointer in jA
+ unsigned r; // row idx in original matrix
+
+ unsigned *iAn(new unsigned[nsize + 1]);
+ iAn[0] = 0;
+
+ unsigned *pos(new unsigned[nsize + 1]);
+ for (i = 0; i <= nsize; i++)
+ pos[i] = 0;
+
+ for (i = beg; i < end; i++) {
+ r = op_perm[i];
+ idx = _row_ptr[r + 1];
+ for (j = _row_ptr[r]; j < idx; j++) {
+ c = po_perm[_col_idx[j]];
+ if (beg <= c && c < end)
+ ++pos[i - beg];
+ }
+ }
+ for (i = 0; i < nsize; i++)
+ iAn[i + 1] = iAn[i] + pos[i];
+ for (i = 0; i < nsize; i++)
+ pos[i] = iAn[i];
+
+ unsigned *jAn(new unsigned[iAn[nsize]]);
+ for (i = beg; i < end; i++) {
+ r = op_perm[i];
+ idx = _row_ptr[r + 1];
+ for (j = _row_ptr[r]; j < idx; j++) {
+ c = po_perm[_col_idx[j]];
+ if (beg <= c && c < end)
+ jAn[pos[i - beg]++] = c - beg;
+ }
+ }
+
+ delete[] pos;
+ for (i = 0; i < nsize; ++i)
+ quickSort(jAn, iAn[i], iAn[i + 1]);
+ return new AdjMat(nsize, iAn, jAn, NULL);
+}
+
+/**
+ * generate an adjacency matrix (the upper triangle part)
+ * @param n number of nodes of matrix graph / number of rows/columns of the adjacency matrix
+ * @param iA array of size of the number of rows/columns + 1, array contains pointer into jA array
+ * @param jA array of the length of the number of non-zero entries (edges in the matrix graph)
+ */
+void genAdjMat(unsigned n, unsigned* &iA, unsigned* &jA)
+{
+ unsigned i;
+ // count entries of each row
+ unsigned* iAn = new unsigned[n + 1];
+ for (i = 0; i <= n; ++i)
+ iAn[i] = 0;
+
+ // go through all strictly lower triangular entries (i,j) and check
+ // whether (j,i) exists in the upper triangular part
+
+ // set n pointers to the beginning of each row
+ unsigned* co = new unsigned[n];
+ for (i = 0; i < n; ++i)
+ co[i] = iA[i];
+
+ for (i = 0; i < n; ++i)
+ for (unsigned k = iA[i]; k < iA[i + 1]; ++k) {
+ unsigned j = jA[k];
+ if (i < j)
+ ++iAn[i + 1]; // upper triangular entries count
+ else { // lower triangular only if there is no counter part
+ unsigned k1 = iA[j], k2 = iA[j + 1];
+ if (i < jA[k1] || i > jA[k2 - 1])
+ ++iAn[j + 1]; // i is out of bounds
+ else { // go through all uninspected entries in the jth row
+ while (co[j] < k2 && i > jA[co[j]])
+ ++co[j];
+ if (co[j] == k2 || i < jA[co[j]])
+ ++iAn[j + 1];
+ }
+ }
+ }
+
+ // construct array iAn by summing up the contents of iAn
+ // con is a set of pointer refering to iAn
+ unsigned* con = new unsigned[n];
+ co[0] = con[0] = 0;
+ for (i = 1; i < n; ++i) {
+ co[i] = iA[i];
+ con[i] = iAn[i];
+ iAn[i + 1] += iAn[i];
+ }
+
+ unsigned *jAn = new unsigned[iAn[n]];
+ for (i = 1; i < n; ++i)
+ for (unsigned k = iA[i]; k < iA[i + 1]; ++k) {
+ unsigned j = jA[k];
+ // copy all transposed lower triangular entries and all upper
+ // triangular elements up to that position
+ if (j < i) {
+ while (co[j] < iA[j + 1] && i > jA[co[j]]) {
+ if (jA[co[j]] > j)
+ jAn[con[j]++] = jA[co[j]];
+ ++co[j];
+ }
+
+ if (co[j] == iA[j + 1] || i <= jA[co[j]]) {
+ jAn[con[j]++] = i;
+ ++co[i];
+ if (i == jA[co[j]])
+ ++co[j];
+ }
+ }
+ }
+
+ // finish rows
+ for (i = 0; i < n; ++i)
+ for (unsigned k = co[i]; k < iA[i + 1]; ++k)
+ if (i < jA[k])
+ jAn[con[i]++] = jA[k];
+
+ BASELIB::swap(jA, jAn);
+ BASELIB::swap(iA, iAn);
+
+ delete[] jAn;
+ delete[] con;
+ delete[] co;
+ delete[] iAn;
+}
+
+void genFullAdjMat(unsigned n, unsigned* &iA, unsigned* &jA)
+{
+ unsigned i;
+ // count entries of each column
+ unsigned* cnt = new unsigned[n];
+ for (i = 0; i < n; ++i)
+ cnt[i] = 0;
+
+ for (i = 0; i < n; ++i) {
+ unsigned j = iA[i], idx = iA[i + 1];
+ while (j < idx) {
+ cnt[jA[j]]++;
+ j++;
+ }
+ }
+
+ // summing up entries
+ for (i = 2; i < n; ++i)
+ cnt[i] += cnt[i - 1];
+
+ unsigned* iAn = new unsigned[n + 1]; // VALGRIND meldet hier Fehler
+ iAn[0] = 0;
+ for (i = 1; i <= n; ++i)
+ iAn[i] = iA[i] + cnt[i - 1];
+
+ unsigned *jAn = new unsigned[iAn[n]];
+ for (unsigned k = 0; k < n; k++)
+ cnt[k] = iAn[k];
+
+ for (i = 0; i < n; ++i) {
+ unsigned j = iA[i], idx = iA[i + 1];
+ while (j < idx) {
+ jAn[cnt[i]++] = jA[j];
+ jAn[cnt[jA[j]]++] = i;
+ j++;
+ }
+ }
+
+ BASELIB::swap(jA, jAn);
+ BASELIB::swap(iA, iAn);
+
+ delete[] jAn;
+ delete[] iAn;
+ delete[] cnt;
+}
+
+void AdjMat::makeSymmetric()
+{
+ // store upper triangular mat values
+ genAdjMat(MatrixBase::_n_rows, _row_ptr, _col_idx);
+ // mirror the upper triangular part into lower
+ genFullAdjMat(MatrixBase::_n_rows, _row_ptr, _col_idx);
+}
+
+} // end namespace MathLib
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.h b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.h
new file mode 100644
index 0000000000000000000000000000000000000000..91c8ad7c56bba097e01fe234449341d810bca084
--- /dev/null
+++ b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.h
@@ -0,0 +1,76 @@
+/*
+ * AdjMat.h
+ *
+ * Created on: 02.01.2012
+ * Author: TF
+ */
+
+#ifndef ADJMAT_H_
+#define ADJMAT_H_
+
+#include "LinAlg/Sparse/CRSMatrix.h"
+
+namespace MathLib {
+
+/** AdjMat implements an adjacency matrix used to represent the (possible
+ weighted) edges of a graph. Due to adjacency matrix are mostly sparse
+ the implementation makes use of compressed row storage format.
+*/
+class AdjMat: public CRSMatrix<unsigned, unsigned>
+{
+public:
+ /** constructor with data elements in A
+ * @param s size of the quadratic matrix
+ * @param iA array of length s+1 which holds pointers in jA,
+ * iA[k] points to the first non-zero column-entry of row k,
+ * iA[k]-1 points accordingly to the last non-zero column-entry of row k,
+ * the last entry of iA (iA[s]) takes the number of non zero entries(nnz)
+ * @param jA array of length nnz, each entry is a colum-index
+ * @param A data-array of length nnz of type unsigned (weights)
+ */
+ AdjMat(unsigned s, unsigned *iA, unsigned *jA, unsigned *A = NULL) :
+ CRSMatrix<unsigned, unsigned> (s, iA, jA, A)
+ {}
+
+ /**
+ * destructor
+ */
+ virtual ~AdjMat()
+ {}
+
+ /**
+ * The
+ */
+ void makeSymmetric();
+
+ /** getMat returns the (possibly reducible) block [beg,end-1] x [beg,end-1]
+ * respecting the permutation.
+ * @param beg index of first row/column, it is supposed that 0 <= beg <= n
+ * @param end index one after last row/column, it is supposed that beg <= end <= n
+ * @param op_perm permutation -> original
+ * @param po_perm original -> permutation
+ * @return pointer to an AdjMat object
+ */
+ AdjMat* getMat(unsigned beg, unsigned end, unsigned const* const op_perm,
+ unsigned const* const po_perm) const;
+
+private:
+ /**
+ * copy constructor
+ * @param src another AdjMat with data for initializiation
+ */
+ AdjMat(const AdjMat&);
+
+ /**
+ * assignment operator
+ * @param rhs a instance of AdjacencyCRSMatrix
+ */
+ AdjMat& operator=(AdjMat&)
+ {
+ return *this;
+ }
+};
+
+}
+
+#endif /* ADJMAT_H_ */
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.cpp b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dadbb3920e1a13c44967a81d09896f24920de1ad
--- /dev/null
+++ b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.cpp
@@ -0,0 +1,174 @@
+/*
+ * Cluster.cpp
+ *
+ * Created on: 02.01.2012
+ * Author: TF
+ */
+
+// BaseLib
+#include "swap.h"
+
+#include "LinAlg/Sparse/NestedDissectionPermutation/Cluster.h"
+//#include "blas.h"
+#include "Cluster.h"
+#include "Separator.h"
+#include "AdjMat.h"
+
+
+namespace MathLib {
+
+// METIS function
+extern "C" void METIS_NodeComputeSeparator(unsigned*, unsigned*, unsigned*,
+ unsigned*, unsigned*, unsigned*,
+ unsigned*, unsigned*);
+
+Cluster::Cluster (unsigned n, unsigned* iA, unsigned* jA)
+ : ClusterBase (n, iA, jA)
+{}
+
+
+Cluster::Cluster(ClusterBase* father, unsigned beg, unsigned end,
+ unsigned* op_perm, unsigned* po_perm,
+ AdjMat* global_mat, AdjMat* local_mat)
+ : ClusterBase(father, beg, end, op_perm, po_perm, global_mat, local_mat)
+{}
+
+void Cluster::subdivide(unsigned bmin)
+{
+ const unsigned size(_end - _beg);
+ if (size > bmin) {
+
+ unsigned n_rows(_l_adj_mat->getNRows());
+ unsigned *xadj(const_cast<unsigned*>(_l_adj_mat->getRowPtrArray()));
+ unsigned *adjncy(const_cast<unsigned*>(_l_adj_mat->getColIdxArray()));
+// unsigned nparts = 2;
+ unsigned options = 0; // for METIS
+ unsigned sepsize(0); // for METIS
+ unsigned *vwgt(new unsigned[n_rows + 1]);
+ const unsigned nnz(xadj[n_rows]);
+ unsigned *adjwgt(new unsigned[nnz]);
+ for (unsigned k(0); k < n_rows + 1; k++)
+ vwgt[k] = 1;
+ for (unsigned k(0); k < nnz; k++)
+ adjwgt[k] = 1;
+ unsigned *part(new unsigned[n_rows + 1]);
+
+ // subdivide the index set into three parts employing METIS
+ METIS_NodeComputeSeparator(&n_rows, xadj, adjncy, vwgt, adjwgt, &options,
+ &sepsize, part);
+
+ // create and init local permutations
+ unsigned *l_op_perm(new unsigned[size]);
+ unsigned *l_po_perm(new unsigned[size]);
+ for (unsigned i = 0; i < size; ++i)
+ l_op_perm[i] = l_po_perm[i] = i;
+
+ unsigned isep1, isep2;
+ updatePerm(part, isep1, isep2, l_op_perm, l_po_perm);
+ delete[] part;
+
+ // update global permutation
+ unsigned *t_op_perm = new unsigned[size];
+ for (unsigned k = 0; k < size; ++k)
+ t_op_perm[k] = _g_op_perm[_beg + l_op_perm[k]];
+
+ for (unsigned k = _beg; k < _end; ++k) {
+ _g_op_perm[k] = t_op_perm[k - _beg];
+ _g_po_perm[_g_op_perm[k]] = k;
+ }
+ delete[] t_op_perm;
+
+ // next recursion step
+ if ((isep1 >= bmin) && (isep2 - isep1 >= bmin)) {
+ // construct adj matrices for [0, isep1), [isep1,isep2), [isep2, _end)
+ AdjMat *l_adj0(_l_adj_mat->getMat(0, isep1, l_op_perm, l_po_perm));
+ AdjMat *l_adj1(_l_adj_mat->getMat(isep1, isep2, l_op_perm, l_po_perm));
+ AdjMat *l_adj2(_l_adj_mat->getMat(isep2, size, l_op_perm, l_po_perm));
+
+ delete[] l_op_perm;
+ delete[] l_po_perm;
+ delete _l_adj_mat;
+ _l_adj_mat = NULL;
+
+ _n_sons = 3;
+ _sons = new ClusterBase*[_n_sons];
+
+ isep1 += _beg;
+ isep2 += _beg;
+
+ // constructing child nodes for index cluster tree
+ _sons[0] = new Cluster(this, _beg, isep1, _g_op_perm, _g_po_perm, _g_adj_mat, l_adj0);
+ _sons[1] = new Cluster(this, isep1, isep2, _g_op_perm, _g_po_perm, _g_adj_mat, l_adj1);
+ _sons[2] = new Separator(this, isep2, _end, _g_op_perm, _g_po_perm, _g_adj_mat, l_adj2);
+
+ dynamic_cast<Cluster*>(_sons[0])->subdivide(bmin);
+ dynamic_cast<Cluster*>(_sons[1])->subdivide(bmin);
+
+ } else {
+ delete _l_adj_mat;
+ _l_adj_mat = NULL;
+ } // end if next recursion step
+ } // end if ( connected && size () > bmin )
+}
+
+
+void Cluster::updatePerm(unsigned* reordering, unsigned &isep0,
+ unsigned &isep1, unsigned* l_op_perm, unsigned* l_po_perm)
+{
+ unsigned beg = 0, end = _end - _beg;
+ while (beg < end) {
+ if (reordering[beg] >= 1) {
+ --end;
+ while (beg < end && reordering[end] >= 1)
+ --end;
+ // local permutation
+ BASELIB::swap(l_op_perm[beg], l_op_perm[end]);
+ BASELIB::swap(l_po_perm[l_op_perm[beg]], l_po_perm[l_op_perm[end]]);
+ BASELIB::swap(reordering[beg], reordering[end]);
+ }
+ ++beg;
+ }
+ if (beg > end)
+ isep0 = beg - 1;
+ else
+ isep0 = end;
+
+ beg = isep0, end = _end - _beg;
+ while (beg < end) {
+ if (reordering[beg] == 2) {
+ --end;
+ while (beg < end && reordering[end] == 2)
+ --end;
+ // local permutation
+ BASELIB::swap(l_op_perm[beg], l_op_perm[end]);
+ BASELIB::swap(l_po_perm[l_op_perm[beg]], l_po_perm[l_op_perm[end]]);
+ BASELIB::swap(reordering[beg], reordering[end]);
+ }
+ ++beg;
+ }
+ if (beg > end)
+ isep1 = beg - 1;
+ else
+ isep1 = end;
+}
+
+
+void Cluster::createClusterTree(unsigned* op_perm, unsigned* po_perm,
+ unsigned bmin)
+{
+ _g_op_perm = op_perm;
+ _g_po_perm = po_perm;
+
+ // *** 1 create local problem
+ unsigned n = _g_adj_mat->getNRows();
+ unsigned *l_op_perm = new unsigned[n];
+ unsigned *l_po_perm = new unsigned[n];
+ for (unsigned k = 0; k < n; ++k)
+ l_op_perm[k] = l_po_perm[k] = k;
+ _l_adj_mat = _l_adj_mat->getMat(0, n, l_op_perm, l_po_perm);
+
+ // *** 2 create cluster tree
+ subdivide(bmin);
+}
+
+} // end namespace MathLib
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.h b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.h
new file mode 100644
index 0000000000000000000000000000000000000000..1ee767c2b454e4e8dba4afad2e43c589cc068eb5
--- /dev/null
+++ b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.h
@@ -0,0 +1,79 @@
+/*
+ * Cluster.h
+ *
+ * Created on: 02.01.2012
+ * Author: TF
+ */
+
+#ifndef CLUSTER_H_
+#define CLUSTER_H_
+
+#include "ClusterBase.h"
+
+namespace MathLib {
+
+/** \brief class for storing clusters of degrees of freedom without
+ geometric information. This class stores clusters (not separators).
+*/
+
+class Cluster: public ClusterBase
+{
+public:
+ /**
+ * Constructor creates the root of the cluster tree
+ * @param n
+ * @param jA
+ * @param iA
+ * @return
+ */
+ Cluster(unsigned n, unsigned* iA, unsigned* jA);
+
+ virtual void subdivide(unsigned bmin);
+
+ /** Method returns the status of this ClusterBase object. In this case
+ * instances of this class are "normal" Clusters.
+ * @return false
+ */
+ virtual bool isSeparator() const
+ {
+ return false;
+ }
+
+ /** Destructor. */
+ virtual ~Cluster() {}
+
+ /**
+ * Method creates recursively the cluster tree, i.e. changes the permutation
+ * op_perm and po_perm and create child cluster trees. For this task only the
+ * adjacency matrix is used.
+ * @param op_perm permutation: original_idx = op_perm[permutated_idx]
+ * @param po_perm reverse permutation: permutated_idx = po_perm[original_idx]
+ * @param bmin threshold value for stopping further refinement
+ * @return a cluster tree
+ */
+ virtual void createClusterTree(unsigned* op_perm, unsigned* po_perm,
+ unsigned bmin = 50);
+
+protected:
+ /** \brief Constructor
+ \param father parent node in cluster tree
+ \param beg beginning index of the cluster
+ \param end beginning index of the next cluster
+ \param op_perm permutation
+ \param po_perm permutation
+ \param global_mat reference to adjacency matrix of the matrix graph in
+ crs format
+ \param local_mat pointer to the local adjacency matrix of the matrix
+ graph in crs format
+ */
+ Cluster(ClusterBase* father, unsigned beg, unsigned end, unsigned* op_perm,
+ unsigned* po_perm, AdjMat* global_mat, AdjMat* local_mat);
+
+private:
+ /** update perm */
+ void updatePerm(unsigned* reordering, unsigned &isep0, unsigned &isep1, unsigned* l_op_perm, unsigned* l_po_perm);
+};
+
+} // end namespace MathLib
+
+#endif /* CLUSTER_H_ */
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/ClusterBase.cpp b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/ClusterBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..215672da6d6f443b4f3717ea0065fadaef71ff64
--- /dev/null
+++ b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/ClusterBase.cpp
@@ -0,0 +1,48 @@
+/*
+ * ClusterBase.cpp
+ *
+ * Created on: 02.01.2012
+ * Author: TF
+ */
+
+//#include "blas.h"
+#include "AdjMat.h"
+#include "LinAlg/Sparse/NestedDissectionPermutation/ClusterBase.h"
+
+namespace MathLib {
+
+ClusterBase::ClusterBase(unsigned n, unsigned const* const iA,
+ unsigned const*const jA) :
+ _beg(0), _end(n), _n_sons(0), _sons(NULL), _parent(NULL), _g_op_perm(NULL),
+ _g_po_perm(NULL), _l_adj_mat(NULL)
+{
+ const unsigned nnz = iA[n];
+
+ // create adjacency matrix
+ unsigned *row_ptr = new unsigned[n + 1];
+ for (unsigned k = 0; k <= n; ++k)
+ row_ptr[k] = iA[k];
+ unsigned *col_idx = new unsigned[nnz];
+ for (unsigned k = 0; k < nnz; ++k)
+ col_idx[k] = jA[k];
+
+ _l_adj_mat = new AdjMat(n, row_ptr, col_idx);
+ _l_adj_mat->makeSymmetric();
+}
+
+ClusterBase::ClusterBase(ClusterBase *father, unsigned beg, unsigned end,
+ unsigned* op_perm, unsigned* po_perm, AdjMat* global_mat, AdjMat* local_mat) :
+ _beg(beg), _end(end), _n_sons(0), _sons(NULL), _parent(father),
+ _g_op_perm(op_perm), _g_po_perm(po_perm), _g_adj_mat(global_mat),
+ _l_adj_mat(local_mat)
+{
+}
+
+ClusterBase::~ClusterBase()
+{
+ if (_parent == NULL)
+ delete _g_adj_mat;
+ delete _l_adj_mat;
+}
+
+} // end namespace MathLib
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/ClusterBase.h b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/ClusterBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..291df29ab15faa4e8cc9c843403d0589d2aa656b
--- /dev/null
+++ b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/ClusterBase.h
@@ -0,0 +1,102 @@
+/*
+ * ClusterBase.h
+ *
+ * Created on: 02.01.2012
+ * Author: TF
+ */
+
+#ifndef CLUSTERBASE_H_
+#define CLUSTERBASE_H_
+
+namespace MathLib {
+
+class AdjMat;
+
+/** \brief Base class for storing clusters of degrees of freedom without
+ geometric information.
+ */
+class ClusterBase
+{
+public:
+ /**
+ * Constructor creates the root of the cluster tree
+ * @param n number of rows/columns
+ * @param iA row pointer array
+ * @param jA column index array
+ * @return
+ */
+ ClusterBase(unsigned n, unsigned const*const iA, unsigned const*const jA);
+ /*!
+ \brief Constructor
+ \param father parent node in cluster tree
+ \param beg beginning index of the cluster
+ \param end beginning index of the next cluster
+ \param op_perm global permutation array (original_idx = op_perm[permuted_idx])
+ \param po_perm global permutation array (permuted_idx = po_perm[original_idx])
+ \param global_mat reference to the global adjacency matrix of the matrix graph in crs format
+ \param local_mat pointer to the local adjacency matrix of the matrix graph in crs format
+ */
+ ClusterBase(ClusterBase* father, unsigned beg, unsigned end,
+ unsigned* op_perm, unsigned* po_perm, AdjMat* global_mat, AdjMat* local_mat);
+
+ /** \brief Destructor.
+ * Destructor frees all form the objects allocated memory.
+ * */
+ virtual ~ClusterBase();
+
+ virtual bool isSeparator() const = 0;
+
+protected:
+ /** \brief Method returns the pointer to the parent cluster.
+ \returns parent cluster */
+ ClusterBase* getParent() const
+ {
+ return _parent;
+ }
+ /**
+ * beginning index in the global permutation arrays
+ */
+ const unsigned _beg;
+
+ /**
+ * beginning index of next next cluster in the global permutation arrays
+ */
+ const unsigned _end;
+
+ /**
+ * number of sons, _n_sons==0 iff this is a leaf
+ */
+ unsigned _n_sons;
+
+ /**
+ * the array of sons of this cluster in the cluster tree
+ */
+ ClusterBase** _sons;
+
+ /**
+ * pointer to parent
+ */
+ ClusterBase *_parent;
+ /**
+ * pointer global permutation array (original_idx = op_perm[permuted_idx])
+ */
+ unsigned* _g_op_perm;
+ /**
+ * global permutation: permutation <- po_perm <- original
+ */
+ unsigned* _g_po_perm;
+ /**
+ * pointer to global adjacency matrix
+ * The attribute _g_adj_mat stores the set of edges of the matrix graph $G = (V,E)$.
+ * (see class AdjMat)
+ */
+ AdjMat* _g_adj_mat;
+ /**
+ * local adjacency matrix
+ */
+ AdjMat* _l_adj_mat;
+};
+
+} // end namespace MathLib
+
+#endif /* CLUSTERBASE_H_ */
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Separator.cpp b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Separator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7d7d211ab0b67f51cc7cb7d6763d0e6e8db18998
--- /dev/null
+++ b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Separator.cpp
@@ -0,0 +1,45 @@
+/*
+ * Separator.cpp
+ *
+ * Created on: 02.01.2012
+ * Author: TF
+ */
+
+// BaseLib
+#include "swap.h"
+
+#include "LinAlg/Sparse/NestedDissectionPermutation/Separator.h"
+
+namespace MathLib {
+
+extern "C" void METIS_PartGraphRecursive(unsigned*, unsigned*, unsigned*,
+ unsigned*, unsigned*, unsigned*,
+ unsigned*, unsigned*, unsigned*,
+ unsigned*, unsigned*);
+
+Separator::Separator(ClusterBase* father, unsigned beg, unsigned end,
+ unsigned* op_perm, unsigned* po_perm,
+ AdjMat* global_mat, AdjMat* local_mat)
+ : ClusterBase (father, beg, end, op_perm, po_perm, global_mat, local_mat)
+{}
+
+Separator::~Separator()
+{}
+
+unsigned Separator::updatePerm(unsigned* reordering, unsigned* l_op_perm, unsigned* l_po_perm)
+{
+ unsigned beg = 0, end = _end-_beg;
+ while (beg < end) {
+ if (reordering[beg] == 1) {
+ --end;
+ while (beg < end && reordering[end] == 1) --end;
+ // local permutation
+ BASELIB::swap(l_op_perm [beg], l_op_perm[end]);
+ BASELIB::swap(l_po_perm[l_op_perm [beg]], l_po_perm[l_op_perm[end]]);
+ }
+ ++beg;
+ }
+ return ((beg>end) ? beg-1 : end);
+}
+
+} // end namespace MathLib
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Separator.h b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Separator.h
new file mode 100644
index 0000000000000000000000000000000000000000..a90ba000139b037cb21646eb461d06b685ce963e
--- /dev/null
+++ b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Separator.h
@@ -0,0 +1,57 @@
+/*
+ * Separator.h
+ *
+ * Created on: 02.01.2012
+ * Author: TF
+ */
+
+#ifndef SEPARATOR_H_
+#define SEPARATOR_H_
+
+#include "LinAlg/Sparse/NestedDissectionPermutation/ClusterBase.h"
+
+namespace MathLib {
+
+class Cluster;
+class AdjMat;
+
+/** \brief class for storing clusters of degrees of freedom without
+ geometric information.
+*/
+class Separator: public ClusterBase
+{
+public:
+ /** brief Constructor builds a initial object for clustering
+ \param father pointer to the father node in cluster tree
+ \param beg index in op_perm and po_perm which describes the begin of the index set of the Separator
+ \param end index in op_perm and po_perm which describes the begin of the index set of the next
+ ClusterBase
+ \param op_perm permutation
+ \param po_perm inverse permutation
+ \param global_mat reference to adjacency matrix of the matrix graph in crs format
+ \param local_mat reference to the local adjacency matrix of the matrix graph in crs format
+ */
+ Separator(ClusterBase* father, unsigned beg, unsigned end,
+ unsigned* op_perm, unsigned* po_perm, AdjMat* global_mat,
+ AdjMat* local_mat);
+
+ /** Destructor. */
+ virtual ~Separator();
+
+ /** Method returns the status of this ClusterAlg object. Instances
+ of this class are Separators.
+ \returns true
+ */
+ virtual bool isSeparator() const
+ {
+ return true;
+ }
+
+private:
+ /** update perm */
+ unsigned updatePerm(unsigned *reordering, unsigned* l_op_perm, unsigned* l_po_perm);
+};
+
+} // end namespace MathLib
+
+#endif /* SEPARATOR_H_ */
diff --git a/MathLib/LinAlg/Sparse/SparseMatrixBase.h b/MathLib/LinAlg/Sparse/SparseMatrixBase.h
index 57f16b088b65e7b79f9086aa6e9c6b901489eee6..6607f9d8ec9589a36f2dfb0da07b96021a7a7856 100644
--- a/MathLib/LinAlg/Sparse/SparseMatrixBase.h
+++ b/MathLib/LinAlg/Sparse/SparseMatrixBase.h
@@ -10,7 +10,13 @@ template<typename FP_TYPE, typename IDX_TYPE> class SparseMatrixBase : public Ma
public:
SparseMatrixBase(IDX_TYPE n1, IDX_TYPE n2) : MatrixBase (n1,n2) {}
SparseMatrixBase() : MatrixBase () {}
- virtual void amux(FP_TYPE d, FP_TYPE const * const __restrict__ x, FP_TYPE * __restrict__ y) const = 0; // y +=d*Ax
+ /**
+ * y = d * A * x
+ * @param d scalar factor
+ * @param x vector to multiply with
+ * @param y result vector
+ */
+ virtual void amux(FP_TYPE d, FP_TYPE const * const __restrict__ x, FP_TYPE * __restrict__ y) const = 0;
virtual ~SparseMatrixBase() {};
};