diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6fc4e08d77a279e6e3e991b46a5ad96a3ffa0af3..e8cfeabcf761c4a894ece6e23f6ba98152d922be 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -212,9 +212,6 @@ if( OGS_BUILD_TESTS AND NOT IS_SUBPROJECT )
else()
add_subdirectory( SimpleTests/MatrixTests )
add_subdirectory( SimpleTests/MeshTests )
- if(NOT MSVC AND BLAS_FOUND AND LAPACK_FOUND)
- add_subdirectory( SimpleTests/SolverTests )
- endif()
endif()
endif() # OGS_BUILD_TESTS
diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
index cc22a818570e27649669d52c397066f689ad2b64..dd6028af94d6df408ae4e2d1dd041ca4c1252a3d 100644
--- a/MathLib/CMakeLists.txt
+++ b/MathLib/CMakeLists.txt
@@ -14,15 +14,9 @@ set(SOURCES ${SOURCES} ${SOURCES_LINALG})
GET_SOURCE_FILES(SOURCES_LINALG_DENSE LinAlg/Dense)
set(SOURCES ${SOURCES} ${SOURCES_LINALG_DENSE})
-GET_SOURCE_FILES(SOURCES_LINALG_SPARSE LinAlg/Sparse)
-set(SOURCES ${SOURCES} ${SOURCES_LINALG_SPARSE})
-
GET_SOURCE_FILES(SOURCES_LINALG_SOLVERS LinAlg/Solvers)
set(SOURCES ${SOURCES} ${SOURCES_LINALG_SOLVERS})
-GET_SOURCE_FILES(SOURCES_LINALG_PRECOND LinAlg/Preconditioner)
-set(SOURCES ${SOURCES} ${SOURCES_LINALG_PRECOND})
-
GET_SOURCE_FILES(SOURCES_ODE ODE)
set(SOURCES ${SOURCES} ${SOURCES_ODE})
@@ -46,17 +40,9 @@ if(OGS_USE_PETSC)
set(SOURCES ${SOURCES} ${SOURCES_LINALG_PETSC})
endif()
-if(METIS_FOUND)
- GET_SOURCE_FILES(SOURCES_LINALG_SPARSE_NESTEDDISSECTION LinAlg/Sparse/NestedDissectionPermutation)
- set(SOURCES ${SOURCES} ${SOURCES_LINALG_SPARSE_NESTEDDISSECTION})
-endif ()
-
GET_SOURCE_FILES(SOURCES_NONLINEAR Nonlinear)
set(SOURCES ${SOURCES} ${SOURCES_NONLINEAR})
-if(METIS_FOUND)
- include_directories(${METIS_INCLUDE_DIR})
-endif()
# Create the library
add_library(MathLib STATIC ${SOURCES})
@@ -71,10 +57,6 @@ if (CVODE_FOUND)
target_link_libraries(MathLib ${CVODE_LIBRARIES})
endif()
-if(METIS_FOUND)
- target_link_libraries(MathLib ${METIS_LIBRARIES})
-endif()
-
if(LAPACK_FOUND)
target_link_libraries(MathLib ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES})
endif()
diff --git a/MathLib/LinAlg/ApplyKnownSolution.h b/MathLib/LinAlg/ApplyKnownSolution.h
index 09c71d7311a5488d9f9e57adeadfba31f364ce21..94b7e9255266822db920f752d4e3e381619d933d 100644
--- a/MathLib/LinAlg/ApplyKnownSolution.h
+++ b/MathLib/LinAlg/ApplyKnownSolution.h
@@ -15,10 +15,6 @@
#include "MathLib/LinAlg/Eigen/EigenTools.h"
#endif // OGS_USE_EIGEN
-#ifdef USE_LIS
-#include "MathLib/LinAlg/Lis/LisTools.h"
-#endif // USE_LIS
-
#ifdef USE_PETSC
#include "MathLib/LinAlg/PETSc/PETScTools.h"
#endif // USE_PETSC
diff --git a/MathLib/LinAlg/Lis/LisTools.cpp b/MathLib/LinAlg/Lis/LisTools.cpp
deleted file mode 100644
index 8a3e024ccc4af2b46d52e19278ac7bdcc6f76a4c..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Lis/LisTools.cpp
+++ /dev/null
@@ -1,112 +0,0 @@
-/**
- * \file
- * \author Norihiro Watanabe
- * \date 2013-04-16
- * \brief Implementation of Lis utility functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "LisTools.h"
-
-#include <cassert>
-
-#include "logog/include/logog.hpp"
-
-#include "LisMatrix.h"
-#include "LisVector.h"
-
-#include "BaseLib/quicksort.h"
-#include "MathLib/LinAlg/Sparse/CRSMatrix.h"
-#include "MathLib/LinAlg/Sparse/CRSTools.h"
-
-namespace MathLib
-{
-
-namespace detail
-{
-/// Converts the internal format of a lis matrix into the compressed row storage
-/// format. Since the column entries of a row in the internal lis format aren't
-/// sorted lis2crs not only transfers the indices and entries, it also
-/// sorts the columns and values, accordingly.
-MathLib::CRSMatrix<double, typename LisMatrix::IndexType>* lis2crs(LisMatrix &a)
-{
- using IndexType = LisMatrix::IndexType;
-
- LIS_MATRIX &A = a.getRawMatrix();
-
- IndexType const n_rows(A->n); // number of rows
- IndexType *iA(new IndexType[n_rows+1]); // row ptr array
- iA[0] = 0;
- for (LIS_INT k=1; k<n_rows+1; ++k) {
- iA[k] = iA[k-1] + A->w_row[k-1 - A->is];
- }
-
- IndexType *jA(new IndexType[iA[n_rows]]); // column indices array
- double *entries(new double[iA[n_rows]]);
- for (IndexType r(0); r<n_rows; ++r) {
- IndexType const beg_idx(iA[r]);
- IndexType const end_idx(iA[r+1]);
- for (IndexType j(beg_idx); j<end_idx; ++j) {
- jA[j] = A->w_index[r-A->is][j-beg_idx];
- entries[j] = A->w_value[r-A->is][j-beg_idx];
- }
- }
-
- for (IndexType r(0); r<n_rows; ++r) {
- IndexType const beg_idx(iA[r]);
- IndexType const end_idx(iA[r+1]);
- // sort the column entries of the row
- BaseLib::quicksort(jA, beg_idx, end_idx, entries);
- }
-
- return new MathLib::CRSMatrix<double,IndexType>(A->n, iA, jA, entries);
-}
-
-// This function resets the the column indices and the entries, respectively.
-// The LIS_MATRIX must have reserved enough memory for each row already!
-void crs2lis(
- MathLib::CRSMatrix<double, typename LisMatrix::IndexType> const& mat,
- LIS_MATRIX &A)
-{
- LisMatrix::IndexType const*const jA(mat.getColIdxArray());
- double * entries(const_cast<double*>(mat.getEntryArray()));
-
- // reset the entries in the lis matrix
- LisMatrix::IndexType cnt(0);
- for (LIS_INT row_i = 0; row_i < A->n; ++row_i) {
- for (LIS_INT j = 0; j < A->w_row[row_i - A->is]; ++j) {
- A->w_index[row_i-A->is][j] = jA[cnt];
- A->w_value[row_i-A->is][j] = entries[cnt];
- cnt++;
- }
- }
-}
-} // end namespace detail
-
-void applyKnownSolution(LisMatrix &eqsA, LisVector &eqsRHS, LisVector &/*eqsX*/,
- const std::vector<LisMatrix::IndexType> &input_rows,
- const std::vector<double> &input_vals)
-{
- // unfortunatly the input is not sorted => copy and sort
- std::vector<LisMatrix::IndexType> rows(input_rows);
- std::vector<double> vals(input_vals);
- BaseLib::quicksort(rows, 0, rows.size(), vals);
-
- MathLib::CRSMatrix<double, typename LisMatrix::IndexType> *crs_mat(
- MathLib::detail::lis2crs(eqsA));
-
- // The following function is defined in CRSTools-impl.h
- applyKnownSolution(crs_mat, eqsRHS, input_rows, input_vals);
-
- detail::crs2lis(*crs_mat, eqsA.getRawMatrix());
-
- delete crs_mat;
-}
-
-} // MathLib
diff --git a/MathLib/LinAlg/Lis/LisTools.h b/MathLib/LinAlg/Lis/LisTools.h
deleted file mode 100644
index c20b6f85def0dd67e6d7bab150bd6e31f583f90d..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Lis/LisTools.h
+++ /dev/null
@@ -1,45 +0,0 @@
-/**
- * \file
- * \author Norihiro Watanabe
- * \date 2013-04-16
- * \brief Definition of Lis utility functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef LISTOOLS_H_
-#define LISTOOLS_H_
-
-#include <vector>
-
-#include "LisMatrix.h" // for LisMatrix::IndexType
-
-namespace MathLib
-{
-class LisVector;
-
-/**
- * Integrate Dirichlet boundary conditions into a system of linear equations.
- *
- * This function introduces the constants into the system by setting
- * appropriated row and column entries of the matrix to zero (except the
- * diagonal entries) and modifying values within the right hand side vector.
- *
- * @param eqsA Coefficient matrix
- * @param eqsRHS RHS vector
- * @param rows a vector of known solution entry IDs
- * @param vals a vector of known solutions
- */
-void applyKnownSolution(LisMatrix &eqsA, LisVector &eqsRHS, LisVector &/*eqsX*/,
- const std::vector<LisMatrix::IndexType> &rows,
- const std::vector<double> &vals);
-
-} // MathLib
-
-#endif //LISTOOLS_H_
-
diff --git a/MathLib/LinAlg/Preconditioner/generateDiagPrecond.cpp b/MathLib/LinAlg/Preconditioner/generateDiagPrecond.cpp
deleted file mode 100644
index eda88c34bb4a7af4632ffc90de5c89dbe2cfc35e..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Preconditioner/generateDiagPrecond.cpp
+++ /dev/null
@@ -1,88 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date no date
- * \brief Implementation of the generateDiagPrecond functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <limits>
-#include <cmath>
-#include <iostream>
-
-namespace MathLib {
-
-bool generateDiagPrecond (unsigned n, unsigned const*const iA, unsigned const*const jA,
- double const*const A, double* diag)
-{
- unsigned idx; // first idx of next row
- unsigned c; // column
- unsigned j;
- bool has_no_diag;
-
- for (unsigned r(0); r<n; ++r) {
- idx=iA[r+1];
- has_no_diag=true;
- for (j=iA[r]; j<idx && has_no_diag; ++j) {
- c=jA[j];
- if (c==r) {
- has_no_diag=false;
- diag[r] = 1.0/A[j];
- }
- }
- if (j==idx && has_no_diag) {
- std::cout << "row " << r << " has no diagonal element " << std::endl;
- return false;
- }
- }
- return true;
-}
-
-bool generateDiagPrecondRowSum(unsigned n, unsigned const*const iA, double const*const A, double* diag)
-{
- unsigned idx; // first idx of next row
- unsigned j;
-
- for (unsigned r(0); r<n; ++r) {
- diag[r] = 0.0;
- idx=iA[r+1];
- for (j=iA[r]; j<idx; ++j) {
- diag[r] += fabs(A[j]);
- }
- if (fabs(diag[r]) < std::numeric_limits<double>::epsilon()) {
- std::cout << "row " << r << " has only very small entries" << std::endl;
- return false;
- }
- diag[r] = 1.0/diag[r];
- }
- return true;
-}
-
-bool generateDiagPrecondRowMax(unsigned n, unsigned const*const iA, double const*const A, double* diag)
-{
- unsigned idx; // first idx of next row
- unsigned j;
-
- for (unsigned r(0); r<n; ++r) {
- idx=iA[r+1];
- diag[r] = A[idx];
- for (j=iA[r]; j<idx; ++j) {
- if (A[j] > diag[r])
- diag[r] = A[j];
- }
- if (fabs(diag[r]) < std::numeric_limits<double>::epsilon()) {
- std::cout << "the maximum entry of row " << r << " has only very small value" << std::endl;
- return false;
- }
- diag[r] = 1.0/diag[r];
- }
- return true;
-}
-
-} // end namespace MathLib
diff --git a/MathLib/LinAlg/Preconditioner/generateDiagPrecond.h b/MathLib/LinAlg/Preconditioner/generateDiagPrecond.h
deleted file mode 100644
index 99747a8a2c81cf744f0601a7315836d5e4745243..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Preconditioner/generateDiagPrecond.h
+++ /dev/null
@@ -1,56 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-09-28
- * \brief Definition of the generateDiagPrecond functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef GENERATEDIAGPRECOND_H_
-#define GENERATEDIAGPRECOND_H_
-
-namespace MathLib {
-
-/**
- * diagonal preconditioner \f$P_{ii} = a_{ii}^{-1}\f$ associated with \f$n \times n\f$ matrix \f$A\f$
- * @param n number of rows / columns
- * @param iA row pointer of compressed row storage format
- * @param jA column index of compressed row storage format
- * @param A data entries of compressed row storage format
- * @param diag inverse entries of the diagonal
- * @return true, if all diagonal entries are distinct from zero, else false
- */
-bool generateDiagPrecond(unsigned n, unsigned const*const iA, unsigned const*const jA,
- double const*const A, double* diag);
-
-/**
- * diagonal preconditioner \f$P_{ii} = \left(\sum_{j} |a_{ij}|\right)^{-1}\f$ associated with \f$n \times n\f$ matrix \f$A\f$
- * @param n number of rows / columns
- * @param iA row pointer of compressed row storage format
- * @param A data entries of compressed row storage format
- * @param diag inverse entries of the diagonal
- * @return true, if all row sums are distinct from zero, else false
- */
-bool generateDiagPrecondRowSum(unsigned n, unsigned const*const iA, double const*const A,
- double* diag);
-
-/**
- * diagonal preconditioner \f$P_{ii} = \left(\max_{j} a_{ij}\right)^{-1}\f$ associated with \f$n \times n\f$ matrix \f$A\f$
- * @param n number of rows / columns
- * @param iA row pointer of compressed row storage format
- * @param A data entries of compressed row storage format
- * @param diag inverse entries of the diagonal
- * @return true, if all row sums are distinct from zero, else false
- */
-bool generateDiagPrecondRowMax(unsigned n, unsigned const*const iA, double const*const A,
- double* diag);
-
-} // end namespace MathLib
-
-#endif /* PRECONDITIONER_H_ */
diff --git a/MathLib/LinAlg/Solvers/BiCGStab.cpp b/MathLib/LinAlg/Solvers/BiCGStab.cpp
deleted file mode 100644
index a7013bc1d98be7f235f0695fb1b3c34d51f3781b..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/BiCGStab.cpp
+++ /dev/null
@@ -1,151 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-10-04
- * \brief Implementation of the BiCGStab function.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "BiCGStab.h"
-#include "MathLib/LinAlg/Sparse/CRSMatrix.h"
-
-
-#include "MathTools.h"
-#include "blas.h"
-
-namespace MathLib {
-
-unsigned BiCGStab(CRSMatrix<double, unsigned> const& A, double* const b, double* const x,
- double& eps, unsigned& nsteps)
-{
- const unsigned N(A.getNRows());
- double *v (new double[8* N]);
- double *p (v + N);
- double *phat (p + N);
- double *s (phat + N);
- double *shat (s + N);
- double *t (shat + N);
- double *r (t + N);
- double *r0 (r + N);
- double resid;
-
- // normb = |b|
- double nrmb = blas::nrm2(N, b);
- if (nrmb < D_PREC) nrmb = D_ONE;
-
- // r = r0 = b - A x0
- blas::copy(N, b, r0);
- A.amux(D_MONE, x, r0);
- blas::copy(N, r0, r);
-
- resid = blas::nrm2(N, r) / nrmb;
-
- if (resid < eps) {
- eps = resid;
- nsteps = 0;
- delete[] v;
- return 0;
- }
-
- double alpha = D_ZERO, omega = D_ZERO, rho2 = D_ZERO;
-
- for (unsigned l = 1; l <= nsteps; ++l) {
- // rho1 = r0 * r
- const double rho1 = blas::scpr(N, r0, r);
- if (fabs(rho1) < D_PREC) {
- eps = blas::nrm2(N, r) / nrmb;
- delete[] v;
- return 2;
- }
-
- if (l == 1)
- blas::copy(N, r, p); // p = r
- else {
-// blas::axpy(N, -omega, v, p); // p = (p-omega v)*beta+r
- const double beta = rho1 * alpha / (rho2 * omega);
-// blas::scal(N, beta, p);
-// blas::axpy(N, D_ONE, r, p);
- // p = (p-omega v)*beta+r
- for (unsigned k(0); k<N; k++) {
- p[k] = (p[k] - omega * v[k]) * beta + r[k];
- }
- }
-
- // p^ = C p
- blas::copy(N, p, phat);
- A.precondApply(phat);
- // v = A p^
- blas::setzero(N, v);
- A.amux(D_ONE, phat, v);
-
- alpha = rho1 / blas::scpr(N, r0, v);
-
- // s = r - alpha v
-// blas::copy(N, r, s);
-// blas::axpy(N, -alpha, v, s);
- for (unsigned k(0); k<N; k++) {
- s[k] = r[k] - alpha * v[k];
- }
-
- resid = blas::nrm2(N, s) / nrmb;
-#ifndef NDEBUG
- std::cout << "Step " << l << ", resid=" << resid << std::endl;
-#endif
- if (resid < eps) {
- // x += alpha p^
- blas::axpy(N, alpha, phat, x);
- eps = resid;
- nsteps = l;
- delete[] v;
- return 0;
- }
-
- // s^ = C s
- blas::copy(N, s, shat);
- A.precondApply(shat);
-
- // t = A s^
- blas::setzero(N, t);
- A.amux(D_ONE, shat, t);
-
- // omega = t*s / t*t
- omega = blas::scpr(N, t, s) / blas::scpr(N, t, t);
-
- // x += alpha p^ + omega s^
- blas::axpy(N, alpha, phat, x);
- blas::axpy(N, omega, shat, x);
-
- // r = s - omega t
- blas::copy(N, s, r);
- blas::axpy(N, -omega, t, r);
-
- rho2 = rho1;
-
- resid = blas::nrm2(N, r) / nrmb;
-
- if (resid < eps) {
- eps = resid;
- nsteps = l;
- delete[] v;
- return 0;
- }
-
- if (fabs(omega) < D_PREC) {
- eps = resid;
- delete[] v;
- return 3;
- }
- }
-
- eps = resid;
- delete[] v;
- return 1;
-}
-
-} // end namespace MathLib
diff --git a/MathLib/LinAlg/Solvers/BiCGStab.h b/MathLib/LinAlg/Solvers/BiCGStab.h
deleted file mode 100644
index 290c888cf18999821c78354e37ad1f9b70f5fd69..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/BiCGStab.h
+++ /dev/null
@@ -1,31 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-10-04
- * \brief Definition of the BiCGStab function.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef BICGSTAB_H_
-#define BICGSTAB_H_
-
-namespace MathLib
-{
-template <typename FP_TYPE, typename IDX_TYPE>
-class CRSMatrix;
-}
-
-namespace MathLib {
-
-unsigned BiCGStab(CRSMatrix<double, unsigned> const& A, double* const b, double* const x,
- double& eps, unsigned& nsteps);
-
-} // end namespace MathLib
-
-#endif /* BICGSTAB_H_ */
diff --git a/MathLib/LinAlg/Solvers/CG.cpp b/MathLib/LinAlg/Solvers/CG.cpp
deleted file mode 100644
index e93f6db443c892e59ae1a49596ceb0c12778e6a0..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/CG.cpp
+++ /dev/null
@@ -1,121 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-09-27
- * \brief Implementation of the CG functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <limits>
-
-#include "MathTools.h"
-#include "blas.h"
-#include "../Sparse/CRSMatrix.h"
-#include "../Sparse/CRSMatrixDiagPrecond.h"
-
-// CG solves the symmetric positive definite linear
-// system Ax=b using the Conjugate Gradient method.
-//
-// The return value indicates convergence within max_iter (input)
-// iterations (0), or no convergence within max_iter iterations (1).
-//
-// Upon successful return, output arguments have the following values:
-//
-// x -- approximate solution to Ax = b
-// nsteps -- the number of iterations performed before the
-// tolerance was reached
-// eps -- the residual after the final iteration
-
-namespace MathLib {
-
-extern
-unsigned CG(CRSMatrix<double,unsigned> const * mat, double const * const b,
- double* const x, double& eps, unsigned& nsteps)
-{
- unsigned N = mat->getNRows();
- double *p, *q, *r, *rhat, rho, rho1 = 0.0;
-
- p = new double[4* N];
- q = p + N;
- r = q + N;
- rhat = r + N;
-
- double nrmb = sqrt(scalarProduct(b, b, N));
- if (nrmb < std::numeric_limits<double>::epsilon()) {
- blas::setzero(N, x);
- eps = 0.0;
- nsteps = 0;
- delete[] p;
- return 0;
- }
-
- // r0 = b - Ax0
- mat->amux(D_MONE, x, r);
- for (unsigned k(0); k < N; k++) {
- r[k] = b[k] - r[k];
- }
-
- double resid = blas::nrm2(N, r);
- if (resid <= eps * nrmb) {
- eps = resid / nrmb;
- nsteps = 0;
- delete[] p;
- return 0;
- }
-
- for (unsigned l = 1; l <= nsteps; ++l) {
-#ifndef NDEBUG
- std::cout << "Step " << l << ", resid=" << resid / nrmb << std::endl;
-#endif
- // r^ = C r
- blas::copy(N, r, rhat);
- mat->precondApply(rhat);
-
- // rho = r * r^;
- rho = scalarProduct(r, rhat, N); // num_threads);
-
- if (l > 1) {
- double beta = rho / rho1;
- // p = r^ + beta * p
- unsigned k;
- for (k = 0; k < N; k++) {
- p[k] = rhat[k] + beta * p[k];
- }
- } else blas::copy(N, rhat, p);
-
- // q = Ap
- blas::setzero(N, q);
- mat->amux(D_ONE, p, q);
-
- // alpha = rho / p*q
- double alpha = rho / scalarProduct(p, q, N);
-
- // x += alpha * p
- blas::axpy(N, alpha, p, x);
-
- // r -= alpha * q
- blas::axpy(N, -alpha, q, r);
-
- resid = sqrt(scalarProduct(r, r, N));
-
- if (resid <= eps * nrmb) {
- eps = resid / nrmb;
- nsteps = l;
- delete[] p;
- return 0;
- }
-
- rho1 = rho;
- }
- eps = resid / nrmb;
- delete[] p;
- return 1;
-}
-
-} // end namespace MathLib
diff --git a/MathLib/LinAlg/Solvers/CG.h b/MathLib/LinAlg/Solvers/CG.h
deleted file mode 100644
index aed3d6eebc30764911add7e6ead3f4ba3d278eae..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/CG.h
+++ /dev/null
@@ -1,33 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-09-27
- * \brief Definition of the CG functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef CG_H_
-#define CG_H_
-
-namespace MathLib {
-
-// forward declaration
-template <typename PF_TYPE, typename IDX_TYPE> class CRSMatrix;
-
-unsigned CG(CRSMatrix<double,unsigned> const * mat, double const * const b,
- double* const x, double& eps, unsigned& nsteps);
-
-#ifdef _OPENMP
-unsigned CGParallel(CRSMatrix<double,unsigned> const * mat, double const * const b,
- double* const x, double& eps, unsigned& nsteps);
-#endif
-
-} // end namespace MathLib
-
-#endif /* SOLVER_H_ */
diff --git a/MathLib/LinAlg/Solvers/CGParallel.cpp b/MathLib/LinAlg/Solvers/CGParallel.cpp
deleted file mode 100644
index 1ba30c35f3e7e35a4b8eca7fcaf9f9b18c128a49..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/CGParallel.cpp
+++ /dev/null
@@ -1,165 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-12-02
- * \brief Implementation of the CGParallel functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <limits>
-
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
-#include "MathTools.h"
-#include "blas.h"
-#include "../Sparse/CRSMatrix.h"
-#include "../Sparse/CRSMatrixDiagPrecond.h"
-
-// CG solves the symmetric positive definite linear
-// system Ax=b using the Conjugate Gradient method.
-//
-// The return value indicates convergence within max_iter (input)
-// iterations (0), or no convergence within max_iter iterations (1).
-//
-// Upon successful return, output arguments have the following values:
-//
-// x -- approximate solution to Ax = b
-// nsteps -- the number of iterations performed before the
-// tolerance was reached
-// eps -- the residual after the final iteration
-
-namespace MathLib {
-
-#ifdef _OPENMP
-unsigned CGParallel(CRSMatrix<double,unsigned> const * mat, double const * const b,
- double* const x, double& eps, unsigned& nsteps)
-{
-#ifdef WIN32
-#pragma warning ( push )
-#pragma warning ( disable: 4018 )
-#endif
- const unsigned N(mat->getNRows());
- double * __restrict__ p(new double[N]);
- double * __restrict__ q(new double[N]);
- double * __restrict__ r(new double[N]);
- double * __restrict__ rhat(new double[N]);
- double rho, rho1 = 0.0;
-
- double nrmb = sqrt(scalarProduct(b, b, N));
-
- if (nrmb < std::numeric_limits<double>::epsilon()) {
- blas::setzero(N, x);
- eps = 0.0;
- nsteps = 0;
- delete[] p;
- return 0;
- }
-
- // r0 = b - Ax0
- mat->amux(D_MONE, x, r);
- for (unsigned k(0); k < N; k++) {
- r[k] = b[k] - r[k];
- }
-
- double resid = blas::nrm2(N, r);
- if (resid <= eps * nrmb) {
- eps = resid / nrmb;
- nsteps = 0;
- delete[] p;
- delete[] q;
- delete[] r;
- delete[] rhat;
- return 0;
- }
-
- OPENMP_LOOP_TYPE k;
- for (unsigned l = 1; l <= nsteps; ++l) {
-#ifndef NDEBUG
- std::cout << "Step " << l << ", resid=" << resid / nrmb << std::endl;
-#endif
-
- // r^ = C r
- // rhat = r
-// blas::copy(N, r, rhat);
-#pragma omp parallel for
- for (k = 0; k < N; k++) {
- rhat[k] = r[k];
- }
- mat->precondApply(rhat);
-
- // rho = r * r^;
- rho = scalarProduct(r, rhat, N);
-
- if (l > 1) {
- double beta = rho / rho1;
- // p = r^ + beta * p
-#pragma omp parallel for
- for (k = 0; k < N; k++) {
- p[k] = rhat[k] + beta * p[k];
- }
- } else {
-// blas::copy(N, rhat, p);
- #pragma omp parallel for
- for (k = 0; k < N; k++) {
- p[k] = rhat[k];
- }
- }
-
- // q = Ap
- mat->amux(D_ONE, p, q);
-
- // alpha = rho / p*q
- double alpha = rho / scalarProduct(p, q, N);
-
- #pragma omp parallel
- {
- // x += alpha * p
- #pragma omp for nowait
- for (k = 0; k < N; k++) {
- x[k] += alpha * p[k];
- }
-
- // r -= alpha * q
- #pragma omp for nowait
- for (k = 0; k < N; k++) {
- r[k] -= alpha * q[k];
- }
-
- #pragma omp barrier
- } // end #pragma omp parallel
-
- resid = sqrt(scalarProduct(r, r, N));
-
- if (resid <= eps * nrmb) {
- eps = resid / nrmb;
- nsteps = l;
- delete[] p;
- delete[] q;
- delete[] r;
- delete[] rhat;
- return 0;
- }
-
- rho1 = rho;
- }
- eps = resid / nrmb;
- delete[] p;
- delete[] q;
- delete[] r;
- delete[] rhat;
- return 1;
-#ifdef WIN32
-#pragma warning ( pop )
-#endif
-}
-#endif
-
-} // end of namespace MathLib
diff --git a/MathLib/LinAlg/Solvers/GMRes.cpp b/MathLib/LinAlg/Solvers/GMRes.cpp
deleted file mode 100644
index 21208584704a1e852881cbfaae5ca832fa65ce19..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/GMRes.cpp
+++ /dev/null
@@ -1,174 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-10-04
- * \brief Implementation of the GMRes functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "GMRes.h"
-
-#include <cmath>
-#include <limits>
-#include "blas.h"
-
-namespace MathLib {
-
-static void genPlRot(double dx, double dy, double& cs, double& sn)
-{
- if (dy <= std::numeric_limits<double>::epsilon()) {
- cs = 1.0;
- sn = 0.0;
- } else if (fabs(dy) > fabs(dx)) {
- const double tmp = dx / dy;
- sn = 1.0 / sqrt(1.0 + tmp * tmp);
- cs = tmp * sn;
- } else {
- const double tmp = dy / dx;
- cs = 1.0 / sqrt(1.0 + tmp * tmp);
- sn = tmp * cs;
- }
-}
-
-inline void applPlRot(double& dx, double& dy, double cs, double sn)
-{
- const double tmp = cs * dx + sn * dy;
- dy = cs * dy - sn * dx;
- dx = tmp;
-}
-
-// solve H y = s and update x += MVy
-static void update(const CRSMatrix<double,unsigned>& A, unsigned k, double* H,
- unsigned ldH, double* s, double* V, double* x)
-{
- const unsigned n(static_cast<unsigned>(A.getNRows()));
- double *y = new double[k];
- double *xh = new double[n];
- blas::copy(k, s, y);
- int inf;
-
- dtrtrs_(JOB_STR + 5, JOB_STR, JOB_STR, &k, &N_ONE, H, &ldH, y, &k, &inf);
- assert(inf == 0);
-
- // x += M V y
- blas::setzero(n, xh);
- blas::gemva(n, k, D_ONE, V, y, xh);
- A.precondApply(xh);
- blas::add(n, xh, x);
-
- delete[] xh;
- delete[] y;
-}
-
-unsigned GMRes(const CRSMatrix<double,unsigned>& A, double* const b, double* const x,
- double& eps, unsigned m, unsigned& nsteps)
-{
- double resid;
- unsigned j = 1;
-
- const unsigned n (static_cast<unsigned>(A.getNRows()));
-
- double *r = new double[2*n + (n + m + 4) * (m + 1)]; // n
- double *V = r + n; // n x (m+1)
- double *H = V + n * (m + 1); // m+1 x m
- double *cs = H + (m + 1) * m; // m+1
- double *sn = cs + m + 1; // m+1
- double *s = sn + m + 1; // m+1
- double *xh = s + m + 1; // m+1
-
- // normb = norm(b)
- double normb = blas::nrm2(n, b);
- if (normb == 0.0) {
- blas::setzero(n, x);
- eps = 0.0;
- nsteps = 0;
- delete[] r;
- return 0;
- }
-
- // r = b - Ax
- blas::copy(n, b, r);
- A.amux(D_MONE, x, r);
-
- double beta = blas::nrm2(n, r);
-
- if ((resid = beta / normb) <= eps) {
- eps = resid;
- nsteps = 0;
- delete[] r;
- return 0;
- }
-
- while (j <= nsteps) {
- blas::copy(n, r, V); // v0 first orthonormal vector
- blas::scal(n, 1.0 / beta, V);
-
- s[0] = beta;
- blas::setzero(m, s + 1);
-
- for (unsigned i = 0; i < m && j <= nsteps; i++, j++) {
-
- // w = A M * v[i];
- blas::copy(n, V + i * n, xh);
- A.precondApply(xh);
- blas::setzero(n, V + (i + 1) * n);
- A.amux(D_ONE, xh, V + (i + 1) * n);
-
- for (unsigned k = 0; k <= i; k++) {
- H[k + i * (m + 1)] = blas::scpr(n, V + (i + 1) * n, V + k * n);
- blas::axpy(n, -H[k + i * (m + 1)], V + k * n, V + (i + 1) * n);
- }
-
- H[i * (m + 2) + 1] = blas::nrm2(n, V + (i + 1) * n);
- blas::scal(n, 1.0 / H[i * (m + 2) + 1], V + (i + 1) * n);
-
- // apply old Givens rotations to the last column in H
- for (unsigned k = 0; k < i; k++)
- applPlRot(H[k + i * (m + 1)], H[k + 1 + i * (m + 1)], cs[k],
- sn[k]);
-
- // generate new Givens rotation which eleminates H[i*(m+2)+1]
- genPlRot(H[i * (m + 2)], H[i * (m + 2) + 1], cs[i], sn[i]);
- // apply it to H and s
- applPlRot(H[i * (m + 2)], H[i * (m + 2) + 1], cs[i], sn[i]);
- applPlRot(s[i], s[i + 1], cs[i], sn[i]);
-
- if ((resid = fabs(s[i + 1] / normb)) < eps) {
- update(A, i + 1, H, m + 1, s, V, x);
- eps = resid;
- nsteps = j;
- delete[] r;
- return 0;
- }
-#ifndef NDEBUG
- std::cout << "Step " << j << ", resid=" << resid << std::endl;
-#endif
- }
-
- update(A, m, H, m + 1, s, V, x);
-
- // r = b - A x;
- blas::copy(n, b, r);
- A.amux(D_MONE, x, r);
- beta = blas::nrm2(n, r);
-
- if ((resid = beta / normb) < eps) {
- eps = resid;
- nsteps = j;
- delete[] r;
- return 0;
- }
- }
-
- eps = resid;
- delete[] r;
- return 1;
-}
-
-} // end namespace MathLib
diff --git a/MathLib/LinAlg/Solvers/GMRes.h b/MathLib/LinAlg/Solvers/GMRes.h
deleted file mode 100644
index 68e74e3a0e65504042253c2a8dc7e5703c36032a..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/GMRes.h
+++ /dev/null
@@ -1,27 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-10-04
- * \brief Definition of the GMRes function.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef GMRES_H_
-#define GMRES_H_
-
-#include "../Sparse/CRSMatrix.h"
-
-namespace MathLib {
-
-unsigned GMRes(const CRSMatrix<double,unsigned>& mat, double* const b, double* const x,
- double& eps, unsigned m, unsigned& steps);
-
-} // end namespace MathLib
-
-#endif /* GMRES_H_ */
diff --git a/MathLib/LinAlg/Solvers/IterativeLinearSolver.h b/MathLib/LinAlg/Solvers/IterativeLinearSolver.h
deleted file mode 100644
index e3de3943c6ca9c23394e1b34748cba47a1217a94..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/IterativeLinearSolver.h
+++ /dev/null
@@ -1,30 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-01-07
- * \brief Definition of the IterativeLinearSolver class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef ITERATIVELINEARSOLVER_H_
-#define ITERATIVELINEARSOLVER_H_
-
-#include <LinearSolver.h>
-
-namespace MathLib {
-
-class IterativeLinearSolver: public MathLib::LinearSolver {
-public:
- IterativeLinearSolver() {};
- virtual ~IterativeLinearSolver() {};
-};
-
-}
-
-#endif /* ITERATIVELINEARSOLVER_H_ */
diff --git a/MathLib/LinAlg/Solvers/blas.h b/MathLib/LinAlg/Solvers/blas.h
deleted file mode 100644
index 77298c12192703e9bf740a78eda38d4c6133afac..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/blas.h
+++ /dev/null
@@ -1,1332 +0,0 @@
-#ifndef BLAS_H
-#define BLAS_H
-
-#include <cassert>
-#include <cstdlib>
-#include <cmath>
-
-
-inline unsigned GE(unsigned i, unsigned j, unsigned n)
-{
- return i+j*n;
-}
-
-const double D_ZERO = 0.0;
-const double D_ONE = 1.0;
-const double D_MONE = -1.0;
-const float S_ZERO = 0.0;
-const float S_ONE = 1.0;
-const float S_MONE = -1.0;
-
-const double D_PREC = 1e-16;
-
-const unsigned N_ONE = 1;
-const int N_MONE = -1;
-const char JOB_STR[] = "NTOSVULCRA";
-
-extern "C"
-{
- /******************************************************************/
- //double precision real
- /******************************************************************/
- unsigned idamax_(const unsigned*, const double*, const unsigned*);
- void dcopy_(const unsigned*, const double*, const unsigned*,
- double*, const unsigned*);
- void daxpy_(const unsigned*, const double*, const double*,
- const unsigned*, double*, const unsigned*);
- void dscal_(const unsigned*, const double*, const double*,
- const unsigned*);
- double ddot_(const unsigned*, const double*, const unsigned*,
- const double*, const unsigned*);
- double dnrm2_(const unsigned*, const double* const, const unsigned*);
-
- void dgtsv_(const unsigned*, const unsigned*, const double*,
- const double*, const double*, const double*, const unsigned*,
- const int*);
- void dgemm_(const char*, const char*, const unsigned*, const unsigned*,
- const unsigned*, const double*, const double*, const unsigned*,
- const double*, const unsigned*, const double*, double*,
- const unsigned*);
- void dger_(const unsigned*, const unsigned*, const double*, const double*,
- const unsigned*, double*, const unsigned*, const double*,
- const unsigned*);
- void dgemv_(const char*, const unsigned*, const unsigned*, const double*,
- const double*, const unsigned*, const double*, const unsigned*,
- const double*, double*, const unsigned*);
- void dorgqr_(const unsigned*, const unsigned*, const unsigned*,
- double*, const unsigned*, double*, double*,
- const unsigned*, int*);
- void dormqr_(const char*, const char*, const unsigned*, const unsigned*,
- const unsigned*, double*, const unsigned*, double*,
- double*, const unsigned*, double*, const unsigned*,
- int*);
- void dsyev_(const char*, const char*, const unsigned*, double*,
- const unsigned*, double*, double*, const unsigned*, int*);
- void dgeqrf_(const unsigned*, const unsigned*, double*, const unsigned*,
- double*, double*, const unsigned*, int*);
- void dgeqp3_(const unsigned*, const unsigned*, const double*,
- const unsigned*, const unsigned*, const double*, const double*,
- const unsigned*, int*);
- void dgeqpf_(const unsigned*, const unsigned*, const double*,
- const unsigned*, const unsigned*, const double*, const double*,
- const unsigned*, int*);
- void dgesvd_(const char*, const char*, const unsigned*, const unsigned*,
- double*, const unsigned*, double*, double*, const unsigned*,
- double*, const unsigned*, double*, const unsigned*, int*);
- void dgetrf_(const unsigned*, const unsigned*, double*, const unsigned*,
- unsigned*, int*);
- void dgetrs_(const char*, const unsigned*, const unsigned*, double*,
- const unsigned*, const unsigned*, double*, const unsigned*,
- int*);
- void dgetri_(const unsigned*, double*, const unsigned*, unsigned*, double*,
- const unsigned*, int*);
- void dspmv_(const char*, const unsigned*, const double*,
- const double*, const double*, const unsigned*,
- const double*, double*, const unsigned*);
- void dsptrf_(const char*, const unsigned*, const double*, int*, int*);
- void dsptri_(const char*, const unsigned*, const double*, int*, double*, int*);
- void dpotrf_(const char*, const unsigned*, const double*, const unsigned*,
- int*);
- void dpotri_(const char*, const unsigned*, const double*, const unsigned*,
- int*);
- void dpptrf_(const char*, const unsigned*, const double*, int*);
- void dpptri_(const char*, const unsigned*, const double*, int*);
- void dtptrs_(const char*, const char*, const char*, const unsigned*,
-
- const unsigned*, double*, const unsigned*, double*,
- double*, const unsigned*, double*, const unsigned*,
- int*);
- void dtpsv_(const char*, const char*, const char*, const unsigned*,
- double*, double*, const unsigned*);
- void dtrtrs_(const char*, const char*, const char*, const unsigned*,
- const unsigned*, const double*, const unsigned*,
- double*, const unsigned*, int*);
- void dtrsm_(const char*, const char*, const char*, const char*,
- const unsigned*, const unsigned*, const double*, const double*,
- const unsigned*, double*, const unsigned*);
- void dtpmv_(const char*, const char*, const char*, const unsigned*,
- const double*, double*, const unsigned*);
- void dlacpy_(const char*, const unsigned*, const unsigned*, const double*,
- const unsigned*, double*, const unsigned*);
- void dlaset_(const char*, const unsigned*, const unsigned*, const double*,
- const double*, double*, const unsigned*);
- void dtrmm_(const char *, const char *, const char *, const char *,
- unsigned *, unsigned *, const double *, double *, unsigned *,
- double *, unsigned *);
- void dswap_(const unsigned*, double*, const unsigned*, double*,
- const unsigned*);
-
- /******************************************************************/
- //single precision real
- /******************************************************************/
- unsigned isamax_(const unsigned*, const float*, const unsigned*);
- void scopy_(const unsigned*, const float*, const unsigned*,
- float*, const unsigned*);
- void saxpy_(const unsigned*, const float*, const float*,
- const unsigned*, float*, const unsigned*);
- void sscal_(const unsigned*, const float*, const float*,
- const unsigned*);
- float sdot_(const unsigned*, const float*, const unsigned*,
- const float*, const unsigned*);
- float snrm2_(const unsigned*, const float* const, const unsigned*);
-
- void sgtsv_(const unsigned*, const unsigned*, const float*,
- const float*, const float*, const float*, const unsigned*,
- const int*);
- void sgemm_(const char*, const char*, const unsigned*, const unsigned*,
- const unsigned*, const float*, const float*, const unsigned*,
- const float*, const unsigned*, const float*, float*,
- const unsigned*);
- void sger_(const unsigned*, const unsigned*, const float*, const float*,
- const unsigned*, float*, const unsigned*, const float*,
- const unsigned*);
- void sgemv_(const char*, const unsigned*, const unsigned*, const float*,
- const float*, const unsigned*, const float*, const unsigned*,
- const float*, float*, const unsigned*);
- void sorgqr_(const unsigned*, const unsigned*, const unsigned*,
- float*, const unsigned*, float*, float*,
- const unsigned*, int*);
- void sormqr_(const char*, const char*, const unsigned*, const unsigned*,
- const unsigned*, float*, const unsigned*, float*,
- float*, const unsigned*, float*, const unsigned*,
- int*);
- void ssyev_(const char*, const char*, const unsigned*, float*,
- const unsigned*, float*, float*, const unsigned*, int*);
- void stpsv_(const char*, const char*, const char*, const unsigned*,
- float*, float*, const unsigned*);
- void sgeqrf_(const unsigned*, const unsigned*, float*, const unsigned*,
- float*, float*, const unsigned*, int*);
- void sgesvd_(const char*, const char*, const unsigned*, const unsigned*,
- float*, const unsigned*, float*, float*, const unsigned*,
- float*, const unsigned*, float*, const unsigned*, int*);
- void sgetrf_(const unsigned*, const unsigned*, float*, const unsigned*,
- unsigned*, int*);
- void sgetri_(const unsigned*, float*, const unsigned*, unsigned*, float*,
- const unsigned*, int*);
- void sspmv_(const char*, const unsigned*, const float*,
- const float*, const float*, const unsigned*,
- const float*, float*, const unsigned*);
- void strsm_(const char*, const char*, const char*, const char*,
- const unsigned*, const unsigned*, const float*, const float*,
- const unsigned*, float*, const unsigned*);
- void ssptrf_(const char*, const unsigned*, const float*, int*, int*);
- void ssptri_(const char*, const unsigned*, const float*, int*, float*, int*);
- void spotrf_(const char*, const unsigned*, const float*, const unsigned*,
- int*);
- void spotri_(const char*, const unsigned*, const float*, const unsigned*,
- int*);
- void spptrf_(const char*, const unsigned*, const float*, int*);
- void spptri_(const char*, const unsigned*, const float*, int*);
- void stptrs_(const char*, const char*, const char*, const unsigned*,
- const unsigned*, float*, float*, const unsigned*, int*);
- void strtrs_(const char*, const char*, const char*, const unsigned*,
- const unsigned*, const float*, const unsigned*,
- float*, const unsigned*, int*);
- void stpmv_(const char*, const char*, const char*, const unsigned*,
- const float*, float*, const unsigned*);
- void sswap_(const unsigned*, float*, const unsigned*, float*,
- const unsigned*);
-}
-
-
-namespace blas
-{
- inline void swap(const unsigned n, double* x, const unsigned incx,
- double* y, const unsigned incy )
- {
- dswap_(&n, x, &incx, y, &incy);
- }
-
- inline void swap(const unsigned n, float* x, const unsigned incx,
- float* y, const unsigned incy )
- {
- sswap_(&n, x, &incx, y, &incy);
- }
-
- inline void laset(const unsigned m, const unsigned n, const double a,
- const double b, double* A, unsigned ldA)
- {
- dlaset_(JOB_STR, &m, &n, &a, &b, A, &ldA);
- }
- inline void lasetu(const unsigned m, const unsigned n, const double a,
- const double b, double* A, unsigned ldA)
- {
- dlaset_(JOB_STR+5, &m, &n, &a, &b, A, &ldA);
- }
- inline void lasetl(const unsigned m, const unsigned n, const double a,
- const double b, double* A, unsigned ldA)
- {
- dlaset_(JOB_STR+6, &m, &n, &a, &b, A, &ldA);
- }
-
- inline unsigned maxi(const unsigned n, double* const v)
- {
- return idamax_(&n, v, &N_ONE);
- }
- inline unsigned maxi(const unsigned n, float* const v)
- {
- return isamax_(&n, v, &N_ONE);
- }
-
- inline void load(const unsigned n, double e, double* const v)
- {
- for (unsigned i=0; i<n; ++i) v[i] = e;
- }
- inline void load(const unsigned n, float e, float* const v)
- {
- for (unsigned i=0; i<n; ++i) v[i] = e;
- }
-
- inline void setzero(const unsigned n, unsigned* const v)
- {
- for (unsigned i=0; i<n; ++i) v[i] = 0;
- }
- inline void setzero(const unsigned n, double* const v)
- {
- load(n, D_ZERO, v);
- }
- inline void setzero(const unsigned n, float* const v)
- {
- load(n, S_ZERO, v);
- }
- inline double nrm2(const unsigned n, const double* const v)
- {
- return dnrm2_(&n, v, &N_ONE);
- }
- inline float nrm2(const unsigned n, const float* const v)
- {
- return snrm2_(&n, v, &N_ONE);
- }
-
- inline void copy(const unsigned n, const double*const orig, double* dest)
- {
- dcopy_(&n, orig, &N_ONE, dest, &N_ONE);
- }
- inline void copy(const unsigned n, const float*const orig, float* dest)
- {
- scopy_(&n, orig, &N_ONE, dest, &N_ONE);
- }
-
- inline void lacpy(const unsigned m, const unsigned n, double* A,
- const unsigned ldA, double* B, const unsigned ldB)
- {
- dlacpy_(JOB_STR, &m, &n, A, &ldA, B, &ldB);
- }
- inline void lacpyu(const unsigned m, const unsigned n, double* A,
- const unsigned ldA, double* B, const unsigned ldB)
- {
- dlacpy_(JOB_STR+5, &m, &n, A, &ldA, B, &ldB);
- }
-
-
- inline void copy(const unsigned n, double* orig, const unsigned inco,
- double* dest, const unsigned incd)
- {
- dcopy_(&n, orig, &inco, dest, &incd);
- }
- inline void copy(const unsigned n, float* orig, const unsigned inco,
- float* dest, const unsigned incd)
- {
- scopy_(&n, orig, &inco, dest, &incd);
- }
-
- // Conv.Copy double2float
- inline void copy(const unsigned n, double* orig, float* dest)
- {
- for (unsigned i=0; i<n; i++) dest[i] = static_cast<float>(orig[i]);
- }
-
- // Conv.Copy float2double
- inline void copy(const unsigned n, float* orig, double* dest)
- {
- for (unsigned i=0; i<n; i++) dest[i] = static_cast<double>(orig[i]);
- }
-
- // Scalar product conj(x)*y
- inline double scpr(const unsigned n, const double* const v1,
- const double* const v2)
- {
- return ddot_(&n, v1, &N_ONE, v2, &N_ONE);
- }
- inline float scpr(const unsigned n, const float* const v1,
- const float* const v2)
- {
- return sdot_(&n, v1, &N_ONE, v2, &N_ONE);
- }
-
- inline double sqrsum(const unsigned n, double* const v)
- {
- return ddot_(&n, v, &N_ONE, v, &N_ONE);
- }
-
- inline float sqrsum(const unsigned n, float* const v)
- {
- return sdot_(&n, v, &N_ONE, v, &N_ONE);
- }
-
- inline void add(const unsigned n, double* const x, double* const y)
- {
- daxpy_(&n, &D_ONE, x, &N_ONE, y, &N_ONE);
- }
- inline void add(const unsigned n, float* const x, float* const y)
- {
- saxpy_(&n, &S_ONE, x, &N_ONE, y, &N_ONE);
- }
-
- inline void axpy(const unsigned n, const double d, const double* const x,
- double* const y)
- {
- daxpy_(&n, &d, x, &N_ONE, y, &N_ONE);
- }
- inline void axpy(const unsigned n, const float d, const float* const x,
- float* const y)
- {
- saxpy_(&n, &d, x, &N_ONE, y, &N_ONE);
- }
-
- inline void scal(const unsigned n, const double d, double* const x, unsigned incx)
- { dscal_(&n, &d, x, &incx); }
- inline void scal(const unsigned n, const float d, float* const x, unsigned incx)
- { sscal_(&n, &d, x, &incx); }
-
- template<class T>
- inline void scal(const unsigned n, const T d, T* const x)
- { scal(n, d, x, N_ONE); }
-
- template<class T> inline void normalize(unsigned n, T* x)
- {
- T s = 1.0/blas::nrm2(n, x);
- blas::scal(n, s, x);
- }
-
- template<class T> inline void mkOrth(unsigned n, const T* v1, T* v2)
- {
- T s = -blas::scpr(n, v1, v2);
- blas::axpy(n, s, v1, v2);
- }
-
- template<class T>
- inline void scal(const unsigned n, const T d, T* const x, T* const y)
- {
- for (unsigned i=0; i<n; ++i) y[i] = d * x[i];
- }
-
-
-
-
- // y = d Ax
- inline void gemv(const unsigned m, const unsigned n, double d,
- const double* A, double *x, double *y)
- {
- dgemv_(JOB_STR, &m, &n, &d, A, &m, x, &N_ONE, &D_ZERO, y, &N_ONE);
- }
- inline void gemv(const unsigned m, const unsigned n, float d, const float* A,
- float *x, float *y)
- {
- sgemv_(JOB_STR, &m, &n, &d, A, &m, x, &N_ONE, &S_ZERO, y, &N_ONE);
- }
-
- // y += d Ax
- inline void gemva(const unsigned m, const unsigned n, double d, const double* A,
- const double *x, double *y)
- {
- dgemv_(JOB_STR, &m, &n, &d, A, &m, x, &N_ONE, &D_ONE, y, &N_ONE);
- }
- inline void gemva(const unsigned m, const unsigned n, float d, const float* A,
- const float *x, float *y)
- {
- sgemv_(JOB_STR, &m, &n, &d, A, &m, x, &N_ONE, &S_ONE, y, &N_ONE);
- }
-
- // y = d A^H x
- inline void gemhv(const unsigned m, const unsigned n, double d, const double* A,
- const double *x, double *y)
- {
- dgemv_(JOB_STR+1, &m, &n, &d, A, &m, x, &N_ONE, &D_ZERO, y, &N_ONE);
- }
- inline void gemhv(const unsigned m, const unsigned n, float d, const float* A,
- const float *x, float *y)
- {
- sgemv_(JOB_STR+1, &m, &n, &d, A, &m, x, &N_ONE, &S_ZERO, y, &N_ONE);
- }
-
- // y += d A^H x
- inline void gemhva(const unsigned m, const unsigned n, double d,
- const double* A, unsigned ldA, const double *x, unsigned incx,
- double *y, unsigned incy)
- {
- dgemv_(JOB_STR+1, &m, &n, &d, A, &ldA, x, &incx, &D_ONE, y, &incy);
- }
- inline void gemhva(const unsigned m, const unsigned n, float d,
- const float* A, unsigned ldA, const float *x, unsigned incx,
- float *y, unsigned incy)
- {
- sgemv_(JOB_STR+1, &m, &n, &d, A, &ldA, x, &incx, &S_ONE, y, &incy);
- }
-
- inline void gemhva(const unsigned m, const unsigned n, double d, const double* A,
- const double *x, double *y)
- { gemhva(m, n, d, A, m, x, N_ONE, y, N_ONE); }
- inline void gemhva(const unsigned m, const unsigned n, float d, const float* A,
- const float *x, float *y)
- { gemhva(m, n, d, A, m, x, N_ONE, y, N_ONE); }
-
- // y += d A x (A symm. dense in packed format)
- inline void gemsva(const unsigned n, double d, double* A,
- double *x, double *y)
- {
- dspmv_(JOB_STR+5, &n, &d, A, x, &N_ONE, &D_ONE, y, &N_ONE);
- }
- inline void gemsva(const unsigned n, float d, float* A,
- float *x, float *y)
- {
- sspmv_(JOB_STR+5, &n, &d, A, x, &N_ONE, &S_ONE, y, &N_ONE);
- }
-
- // sovles Ax=B, A is a triangluar Matrix
- inline void gtsv(const unsigned* n, const double* DiagLower,
- const double* Diag, const double* DiagUpper,
- const double* B, const int* INFO)
- {
- dgtsv_(n, &N_ONE, DiagLower, Diag, DiagUpper, B, n, INFO);
- }
- inline void gtsv(const unsigned* n, const float* DiagLower,
- const float* Diag, const float* DiagUpper,
- const float* B, const int* INFO)
- {
- sgtsv_(n, &N_ONE, DiagLower, Diag, DiagUpper, B, n, INFO);
- }
-
- // C = d A B, A is m x p, B is p x n
- inline void gemm(const unsigned m, const unsigned p, const unsigned n,
- const double d, const double* const A, const unsigned ldA,
- const double* const B, const unsigned ldB,
- double* C, const unsigned ldC)
- {
- dgemm_(JOB_STR, JOB_STR, &m, &n, &p, &d, A, &ldA, B, &ldB,
- &D_ZERO, C, &ldC);
- }
- inline void gemm(const unsigned m, const unsigned p, const unsigned n,
- const float d, const float* const A, const unsigned ldA,
- const float* const B, const unsigned ldB,
- float* C, const unsigned ldC)
- {
- sgemm_(JOB_STR, JOB_STR, &m, &n, &p, &d, A, &ldA, B, &ldB,
- &S_ZERO, C, &ldC);
- }
-
- // C += d A B, A is m x p, B is p x n
- inline void gemma(const unsigned m, const unsigned p, const unsigned n,
- const double d, const double* const A, const unsigned ldA,
- const double* const B, const unsigned ldB,
- double* C, const unsigned ldC)
- {
- dgemm_(JOB_STR, JOB_STR, &m, &n, &p, &d, A, &ldA, B, &ldB,
- &D_ONE, C, &ldC);
- }
- inline void gemma(const unsigned m, const unsigned p, const unsigned n,
- const float d, const float* const A, const unsigned ldA,
- const float* const B, const unsigned ldB,
- float* C, const unsigned ldC)
- {
- sgemm_(JOB_STR, JOB_STR, &m, &n, &p, &d, A, &ldA, B, &ldB,
- &S_ONE, C, &ldC);
- }
-
- // C = d A^H B, A is m x p, B is m x n
- inline void gemhm(const unsigned m, const unsigned p, const unsigned n,
- const double d, const double* A, const unsigned ldA,
- const double *B, const unsigned ldB,
- double* C, const unsigned ldC)
- {
- dgemm_(JOB_STR+1, JOB_STR, &p, &n, &m, &d, A, &ldA, B, &ldB,
- &D_ZERO, C, &ldC);
- }
- inline void gemhm(const unsigned m, const unsigned p, const unsigned n,
- const float d, const float* const A, const unsigned ldA,
- const float* const B, const unsigned ldB,
- float* C, const unsigned ldC)
- {
- sgemm_(JOB_STR+1, JOB_STR, &p, &n, &m, &d, A, &ldA, B, &ldB,
- &S_ZERO, C, &ldC);
- }
-
- // C += d A^H B, A is m x p, B is m x n
- inline void gemhma(unsigned m, unsigned p, unsigned n, double d,
- const double* const A, const unsigned ldA, const double* const B,
- const unsigned ldB, double* C, unsigned ldC)
- {
- dgemm_(JOB_STR+1, JOB_STR, &p, &n, &m, &d, A, &ldA, B, &ldB,
- &D_ONE, C, &ldC);
- }
- inline void gemhma(unsigned m, unsigned p, unsigned n, float d,
- const float* const A, const unsigned ldA, const float* const B,
- const unsigned ldB, float* C, unsigned ldC)
- {
- sgemm_(JOB_STR+1, JOB_STR, &p, &n, &m, &d, A, &ldA, B, &ldB,
- &S_ONE, C, &ldC);
- }
-
- // C = d A B^H, A is m x p, B is n x p
- inline void gemmh(const unsigned m, const unsigned p, const unsigned n,
- const double d, const double* const A, const unsigned ldA,
- const double* const B, const unsigned ldB,
- double* C, const unsigned ldC)
- {
- dgemm_(JOB_STR, JOB_STR+1, &m, &n, &p, &d, A, &ldA, B, &ldB,
- &D_ZERO, C, &ldC);
- }
- inline void gemmh(const unsigned m, const unsigned p, const unsigned n,
- const float d, const float* const A, const unsigned ldA,
- const float* const B, const unsigned ldB,
- float* C, const unsigned ldC)
- {
- sgemm_(JOB_STR, JOB_STR+1, &m, &n, &p, &d, A, &ldA, B, &ldB,
- &S_ZERO, C, &ldC);
- }
-
- // C += d A B^H, A is m x p, B is n x p
- inline void gemmha(const unsigned m, const unsigned p, const unsigned n,
- const double d, const double* const A, const unsigned ldA,
- const double* const B, const unsigned ldB,
- double* C, const unsigned ldC)
- {
- dgemm_(JOB_STR, JOB_STR+1, &m, &n, &p, &d, A, &ldA, B, &ldB,
- &D_ONE, C, &ldC);
- }
- inline void gemmha(const unsigned m, const unsigned p, const unsigned n,
- const float d, const float* const A, const unsigned ldA,
- const float* const B, const unsigned ldB,
- float* C, const unsigned ldC)
- {
- sgemm_(JOB_STR, JOB_STR+1, &m, &n, &p, &d, A, &ldA, B, &ldB,
- &S_ONE, C, &ldC);
- }
- inline void gemmha(const unsigned m, const unsigned p, const unsigned n,
- const double* const A, const unsigned ldA,
- const double* const B, const unsigned ldB,
- double* C, const unsigned ldC)
- {
- dgemm_(JOB_STR, JOB_STR+1, &m, &n, &p, &D_ONE, A, &ldA, B, &ldB,
- &D_ONE, C, &ldC);
- }
- inline void gemmha(const unsigned m, const unsigned p, const unsigned n,
- const float* const A, const unsigned ldA,
- const float* const B, const unsigned ldB,
- float* C, const unsigned ldC)
- {
- sgemm_(JOB_STR, JOB_STR+1, &m, &n, &p, &S_ONE, A, &ldA, B, &ldB,
- &S_ONE, C, &ldC);
- }
-
- // C = d A^H B^H, A is p x m, B is n x p
- inline void gemhmh(const unsigned m, const unsigned p, const unsigned n,
- const double d, const double* const A, const unsigned ldA,
- const double* const B, const unsigned ldB,
- double* C, const unsigned ldC)
- {
- dgemm_(JOB_STR+1, JOB_STR+1, &m, &n, &p, &d, A, &ldA, B, &ldB,
- &D_ZERO, C, &ldC);
- }
- inline void gemhmh(const unsigned m, const unsigned p, const unsigned n,
- const float d, const float* const A, const unsigned ldA,
- const float* const B, const unsigned ldB,
- float* C, const unsigned ldC)
- {
- sgemm_(JOB_STR+1, JOB_STR+1, &m, &n, &p, &d, A, &ldA, B, &ldB,
- &S_ZERO, C, &ldC);
- }
-
- //C += d*AB, A is mxm (packed upper half is stored), B is mxn and regular matrix
- inline void sygemma(const unsigned m, const unsigned n,
- const double* const A, const double* const B,
- const double d, double* const C)
- {
- for(unsigned i=0;i<m;i++){
- for(unsigned j=0;j<n;j++){
- for(unsigned k=i;k<m;k++){
- if(i==k){
- C[j*m+i] += d*A[i+k*(k+1)/2]*B[k+j*m];
- }else{
- C[j*m+i] += d*A[i+k*(k+1)/2]*B[k+j*m];
- C[j*m+k] += d*A[i+k*(k+1)/2]*B[i+j*m];
- }
- }
- }
- }
- }
- inline void sygemma(const unsigned m, const unsigned n,
- const float* const A, const float* const B,
- const float d, float* const C)
- {
- for(unsigned i=0;i<m;i++){
- for(unsigned j=0;j<n;j++){
- for(unsigned k=i;k<m;k++){
- if(i==k){
- C[j*m+i] += d*A[i+k*(k+1)/2]*B[k+j*m];
- }else{
- C[j*m+i] += d*A[i+k*(k+1)/2]*B[k+j*m];
- C[j*m+k] += d*A[i+k*(k+1)/2]*B[i+j*m];
- }
- }
- }
- }
- }
-
- //C += d*AB, A is mxn and regular matrix, B is nxn (packed upper half is stored)
- inline void gesymma(const unsigned m, const unsigned n,
- const double* const A, const double* const B,
- const double d, double* const C)
- {
- for(unsigned i=0;i<m;i++){
- for(unsigned j=0;j<n;j++){
- for(unsigned k=j;k<n;k++){
- if(j==k)
- C[j*m+i] += d*A[i+k*m]*B[k+j*(j+1)/2];
- else{
- C[j*m+i] += d*A[i+k*m]*B[j+k*(k+1)/2];
- C[k*m+i] += d*A[i+j*m]*B[j+k*(k+1)/2];
- }
- }
- }
- }
- }
- inline void gesymma(const unsigned m, const unsigned n,
- const float* const A, const float* const B,
- const float d, float* const C)
- {
- for(unsigned i=0;i<m;i++){
- for(unsigned j=0;j<n;j++){
- for(unsigned k=j;k<n;k++){
- if(j==k)
- C[j*m+i] += d*A[i+k*m]*B[k+j*(j+1)/2];
- else{
- C[j*m+i] += d*A[i+k*m]*B[j+k*(k+1)/2];
- C[k*m+i] += d*A[i+j*m]*B[j+k*(k+1)/2];
- }
- }
- }
- }
- }
-
- // C += d A^H A, C is a symm. matrix (packed upper half is stored), A is mxn
- inline void symhm(const unsigned m, const unsigned n, const double* const A,
- const double d, double* C)
- {
- for (unsigned j=0; j<n; ++j) {
- for (unsigned i=0; i<=j; ++i) {
- double sum = 0.0;
- for (unsigned k=0; k<m; ++k) sum += A[k+i*m] * A[k+j*m];
- C[i+j*(j+1)/2] += d * sum;
- }
- }
- }
- inline void symhm(const unsigned m, const unsigned n, const float* const A,
- const float d, float* C)
- {
- for (unsigned j=0; j<n; ++j) {
- for (unsigned i=0; i<=j; ++i) {
- float sum = 0.0;
- for (unsigned k=0; k<m; ++k) sum += A[k+i*m] * A[k+j*m];
- C[i+j*(j+1)/2] += d * sum;
- }
- }
- }
-
- // C += d A A^H, C is a symm. matrix (packed upper half is stored), A is mxn
- inline void symmh(unsigned m, unsigned n, double* A, double d, double* C)
- {
- for (unsigned k=0; k<n; ++k) {
- for (unsigned j=0; j<=n; ++j) {
- double e = d * A[j+k*m];
- for (unsigned i=0; i<j; ++i) C[i+j*(j+1)/2] += e * A[i+k*m];
- }
- }
- }
- inline void symmh(unsigned m, unsigned n, float* A, float d, float* C)
- {
- for (unsigned k=0; k<n; ++k) {
- for (unsigned j=0; j<n; ++j) {
- float e = d * A[j+k*m];
- for (unsigned i=0; i<=j; ++i) C[i+j*(j+1)/2] += e * A[i+k*m];
- }
- }
- }
-
- // Singular Value Decomposition
- inline int gesvdS(unsigned m, unsigned n, double* A, double* S,
- double* U, unsigned ldU, double* VT, unsigned ldVT,
- unsigned nwk, double* wk)
- {
- int INF;
- dgesvd_(JOB_STR+3, JOB_STR+3, &m, &n, A, &m, S, U, &ldU, VT, &ldVT,
- wk, &nwk, &INF);
- return INF;
- }
- inline int gesvd(unsigned m, unsigned n, double* A, double* S,
- double* VT, unsigned ldVT, unsigned nwk, double* wk)
- {
- int INF;
- dgesvd_(JOB_STR+2, JOB_STR+3, &m, &n, A, &m, S, A, &m, VT, &ldVT,
- wk, &nwk, &INF);
- return INF;
- }
- inline int gesvd(unsigned m, unsigned n, float* A, float* S,
- float* VT, unsigned ldVT, unsigned nwk, float* wk)
- {
- int INF;
- sgesvd_(JOB_STR+2, JOB_STR+3, &m, &n, A, &m, S, A, &m, VT, &ldVT,
- wk, &nwk, &INF);
- return INF;
- }
-
- inline int gesvd(unsigned m, unsigned n, double* A, double* S,
- double* U, unsigned ldU, double* VT, unsigned ldVT,
- unsigned nwk, double* wk)
- {
- int INF;
- dgesvd_(JOB_STR+9, JOB_STR+9, &m, &n, A, &m, S, U, &ldU, VT, &ldVT,
- wk, &nwk, &INF);
- return INF;
- }
-
- // compute singular values
- inline int svals(unsigned m, unsigned n, double* A, double* S,
- unsigned nwk, double* wk)
- {
- int INF;
- dgesvd_(JOB_STR, JOB_STR, &m, &n, A, &m, S, A, &m, A, &n, wk, &nwk, &INF);
- return INF;
- }
- inline int svals(unsigned m, unsigned n, float* A, float* S,
- unsigned nwk, float* wk)
- {
- int INF;
- sgesvd_(JOB_STR, JOB_STR, &m, &n, A, &m, S, A, &m, A, &n, wk, &nwk, &INF);
- return INF;
- }
-
- // triangular factorisation
- inline int getrf(const unsigned n, double* A, unsigned* ipiv)
- {
- int INF;
- dgetrf_(&n, &n, A, &n, ipiv, &INF);
- return INF;
- }
- inline int getrf(const unsigned n, float* A, unsigned* ipiv)
- {
- int INF;
- sgetrf_(&n, &n, A, &n, ipiv, &INF);
- return INF;
- }
-
- // upper triangular packed MV
- inline void utrpv(const unsigned n, double* A, double* x)
- {
- dtpmv_(JOB_STR+5, JOB_STR, JOB_STR, &n, A, x, &N_ONE);
- }
- inline void utrpv(const unsigned n, float* A, float* x)
- {
- stpmv_(JOB_STR+5, JOB_STR, JOB_STR, &n, A, x, &N_ONE);
- }
-
- // lower triangular packed transpose MV
- inline void ltrphv(const unsigned n, double* A, double* x)
- {
- dtpmv_(JOB_STR+6, JOB_STR+1, JOB_STR+5, &n, A, x, &N_ONE);
- }
- inline void ltrphv(const unsigned n, float* A, float* x)
- {
- stpmv_(JOB_STR+6, JOB_STR+1, JOB_STR+5, &n, A, x, &N_ONE);
- }
-
- // QR factorisation
- inline int geqrf(const unsigned m, const unsigned n, double* A,
- double* tau, unsigned nwk, double* wk)
- {
- int INF;
- dgeqrf_(&m, &n, A, &m, tau, wk, &nwk, &INF);
- return INF;
- }
- inline int geqrf(const unsigned m, const unsigned n, float* A,
- float* tau, unsigned nwk, float* wk)
- {
- int INF;
- sgeqrf_(&m, &n, A, &m, tau, wk, &nwk, &INF);
- return INF;
- }
- inline int geqrf(const unsigned m, const unsigned n, double* A,
- const unsigned ldA, double* tau, unsigned nwk, double* wk)
- {
- int INF;
- dgeqrf_(&m, &n, A, &ldA, tau, wk, &nwk, &INF);
- return INF;
- }
- inline int geqrf(const unsigned m, const unsigned n, float* A,
- const unsigned ldA, float* tau, unsigned nwk, float* wk)
- {
- int INF;
- sgeqrf_(&m, &n, A, &ldA, tau, wk, &nwk, &INF);
- return INF;
- }
-
- // Multiply a general Matrix with the Q-Matrix (QR factorization), Q C
- inline int ormqr(const unsigned m, const unsigned n, const unsigned p,
- double* A, double* tau, double* C,
- unsigned nwk, double* wk)
- {
- int INF;
- dormqr_(JOB_STR+6, JOB_STR, &m, &n, &p, A, &m, tau, C, &m, wk, &nwk, &INF);
- return INF;
- }
- inline int ormqr(const unsigned m, const unsigned n, const unsigned p,
- float* A, float* tau, float* C,
- unsigned nwk, float* wk)
- {
- int INF;
- sormqr_(JOB_STR+6, JOB_STR, &m, &n, &p, A, &m, tau, C, &m, wk, &nwk, &INF);
- return INF;
- }
- inline int ormqr(const unsigned m, const unsigned n, const unsigned p,
- double* A, const unsigned ldA, double* tau, double* C,
- const unsigned ldC, unsigned nwk, double* wk)
- {
- int INF;
- dormqr_(JOB_STR+6, JOB_STR, &m, &n, &p, A, &ldA, tau, C, &ldC, wk, &nwk,
- &INF);
- return INF;
- }
- inline int ormqr(const unsigned m, const unsigned n, const unsigned p,
- float* A, const unsigned ldA, float* tau, float* C,
- const unsigned ldC, unsigned nwk, float* wk)
- {
- int INF;
- sormqr_(JOB_STR+6, JOB_STR, &m, &n, &p, A, &ldA, tau, C, &ldC, wk, &nwk,
- &INF);
- return INF;
- }
-
- // Q^H C
- inline int ormqrh(const unsigned m, const unsigned n, const unsigned p,
- double* A, const unsigned ldA, double* tau, double* C,
- const unsigned ldC, unsigned nwk, double* wk)
- {
- int INF;
- dormqr_(JOB_STR+6, JOB_STR+1, &m, &n, &p, A, &ldA, tau, C, &ldC, wk, &nwk,
- &INF);
- return INF;
- }
- inline int ormqrh(const unsigned m, const unsigned n, const unsigned p,
- float* A, const unsigned ldA, float* tau, float* C,
- const unsigned ldC, unsigned nwk, float* wk)
- {
- int INF;
- sormqr_(JOB_STR+6, JOB_STR+1, &m, &n, &p, A, &ldA, tau, C, &ldC, wk, &nwk,
- &INF);
- return INF;
- }
-
- inline int ormqrh(const unsigned m, const unsigned n, const unsigned p,
- double* A, const unsigned ldA, double* tau, double* C,
- unsigned nwk, double* wk)
- {
- int INF;
- dormqr_(JOB_STR+6, JOB_STR+1, &m, &n, &p, A, &ldA, tau, C, &m, wk, &nwk,
- &INF);
- return INF;
- }
- inline int ormqrh(const unsigned m, const unsigned n, const unsigned p,
- float* A, const unsigned ldA, float* tau, float* C,
- unsigned nwk, float* wk)
- {
- int INF;
- sormqr_(JOB_STR+6, JOB_STR+1, &m, &n, &p, A, &ldA, tau, C, &m, wk, &nwk,
- &INF);
- return INF;
- }
-
- inline int ormqrh(const unsigned m, const unsigned n, const unsigned p,
- double* A, double* tau, double* C,
- unsigned nwk, double* wk)
- {
- int INF;
- dormqr_(JOB_STR+6, JOB_STR+1, &m, &n, &p, A, &m, tau, C, &m, wk, &nwk, &INF);
- return INF;
- }
- inline int ormqrh(const unsigned m, const unsigned n, const unsigned p,
- float* A, float* tau, float* C,
- unsigned nwk, float* wk)
- {
- int INF;
- sormqr_(JOB_STR+6, JOB_STR+1, &m, &n, &p, A, &m, tau, C, &m, wk, &nwk, &INF);
- return INF;
- }
- inline int morqr(const unsigned m, const unsigned n, const unsigned p,
- double* A, const unsigned ldA, double* tau, double* C,
- const unsigned ldC, unsigned nwk, double* wk)
- {
- int INF;
- dormqr_(JOB_STR+8, JOB_STR, &m, &n, &p, A, &ldA, tau, C, &ldC, wk, &nwk,
- &INF);
- return INF;
- }
- inline int morqr(const unsigned m, const unsigned n, const unsigned p,
- float* A, const unsigned ldA, float* tau, float* C,
- const unsigned ldC, unsigned nwk, float* wk)
- {
- int INF;
- sormqr_(JOB_STR+8, JOB_STR, &m, &n, &p, A, &ldA, tau, C, &ldC, wk, &nwk,
- &INF);
- return INF;
- }
-
- inline void ger(unsigned M, unsigned N, double d, double* X, unsigned INCX,
- double* y, unsigned INCY, double* A, unsigned LDA)
- {
- dger_(&M, &N, &d, X, &INCX, y, &INCY, A, &LDA);
- }
-
- inline void ger(unsigned M, unsigned N, float d, float* X, unsigned INCX,
- float* y, unsigned INCY, float* A, unsigned LDA)
- {
- sger_(&M, &N, &d, X, &INCX, y, &INCY, A, &LDA);
- }
-
- // return Q-Matrix (QR factorization) in A
- inline int orgqr(const unsigned m, const unsigned n, double* A, double* tau,
- unsigned nwk, double* wk)
- {
- int INF;
- dorgqr_(&m, &n, &n, A, &m, tau, wk, &nwk, &INF);
- return INF;
- }
- inline int orgqr(const unsigned m, const unsigned n, float* A, float* tau,
- unsigned nwk, float* wk)
- {
- int INF;
- sorgqr_(&m, &n, &n, A, &m, tau, wk, &nwk, &INF);
- return INF;
- }
-
- // B=A^H, A in R^(m,n)
- inline void transpose(unsigned m, unsigned n, double* A, double* B)
- {
- for (unsigned i=0; i<m; ++i) dcopy_(&n, A+i, &m, B+i*n, &N_ONE);
- }
- inline void transpose(unsigned m, unsigned n, float* A, float* B)
- {
- for (unsigned i=0; i<m; ++i) scopy_(&n, A+i, &m, B+i*n, &N_ONE);
- }
-
- // product of an upper triangular matrix U and a matrix A, A:=U A
- inline void utrgemm(unsigned m, unsigned n, double* U, unsigned ldU,
- double* A, unsigned ldA)
- {
- dtrmm_(JOB_STR+6, JOB_STR+5, JOB_STR, JOB_STR, &m, &n, &D_ONE, U, &ldU,
- A, &ldA);
- }
-
-
- // A += d U U^H, where U is upper triangular in packed storage
- // only the upper triangular part of A is computed
- inline void putrmmh(double d, unsigned n, double* U, double* A)
- {
- for (unsigned j=0; j<n; ++j) {
- for (unsigned l=j; l<n; ++l) {
- unsigned idl = l*(l+1)/2;
- double e = d * U[j+idl];
- for (unsigned i=0; i<=j; ++i) A[i+j*(j+1)/2] += e * U[i+idl];
- }
- }
- }
-
- inline void putrmmh(float d, unsigned n, float* U, float* A)
- {
- for (unsigned j=0; j<n; ++j) {
- for (unsigned l=j; l<n; ++l) {
- unsigned idl = l*(l+1)/2;
- float e = d * U[j+idl];
- for (unsigned i=0; i<=j; ++i) A[i+j*(j+1)/2] += e * U[i+idl];
- }
- }
- }
-
- inline void fill0_ltr(unsigned n, double* A)
- {
- for (unsigned j=0; j<n; ++j) {
- *A++ = static_cast<double>(j);
- for (unsigned i=j+1; i<n; ++i) *A++ = D_ZERO;
- }
- }
- inline void fill0_ltr(unsigned n, float* A)
- {
- for (unsigned j=0; j<n; ++j) {
- *A++ = static_cast<float>(j);
- for (unsigned i=j+1; i<n; ++i) *A++ = S_ZERO;
- }
- }
-
- // fill nxn matrix A with identity
- inline void fillId(unsigned n, double *A)
- {
- for (unsigned j=0; j<n; ++j) {
- unsigned i = 0;
- for (; i<j; ++i) *A++ = D_ZERO;
- *A++ = D_ONE;
- ++i;
- for (; i<n; ++i) *A++ = D_ZERO;
- }
- }
- inline void fillId(unsigned n, float *A)
- {
- for (unsigned j=0; j<n; ++j) {
- unsigned i = 0;
- for (; i<j; ++i) *A++ = S_ZERO;
- *A++ = S_ONE;
- ++i;
- for (; i<n; ++i) *A++ = S_ZERO;
- }
- }
-
- // fill nxn upper triang. matrix A with identity (packed storage)
- inline void fillId_utr(unsigned n, double *A)
- {
- for (unsigned j=0; j<n; ++j) {
- for (unsigned i=0; i<j; ++i) *A++ = D_ZERO;
- *A++ = D_ONE;
- }
- }
- inline void fillId_utr(unsigned n, float *A)
- {
- for (unsigned j=0; j<n; ++j) {
- for (unsigned i=0; i<j; ++i) *A++ = S_ZERO;
- *A++ = S_ONE;
- }
- }
-
- // fill nxn normalized lower triang. matrix A with identity (packed storage)
- inline void fillId_ltr(unsigned n, double *A)
- {
- for (unsigned i=0; i<n; ++i) {
- for (unsigned j=0; j<i; ++j) *A++ = D_ZERO;
- // for pivoting, a ltr is assumed to have ones on the diagonal
- *A++ = static_cast<double>(i);
- }
- }
- inline void fillId_ltr(unsigned n, float *A)
- {
- for (unsigned i=0; i<n; ++i) {
- for (unsigned j=0; j<i; ++j) *A++ = S_ZERO;
- // for pivoting, a ltr is assumed to have ones on the diagonal
- *A++ = static_cast<float>(i);
- }
- }
-
-}
-
-
-namespace lapack
-{
-
- // general inversion
- inline void geinv(const unsigned n, double* const A)
- {
- unsigned* const ipiv = new unsigned[n];
- assert(ipiv!=NULL);
-
- const unsigned lwork = 4*n;
- double* const work = new double[lwork];
- assert(work!=NULL);
-
- int INFO = blas::getrf(n, A, ipiv);
- assert(INFO==0);
-
- dgetri_(&n, A, &n, ipiv, work, &lwork, &INFO);
- assert(INFO==0);
-
- delete [] work;
- delete [] ipiv;
- }
- inline void geinv(const unsigned n, float* const A)
- {
- unsigned* const ipiv = new unsigned[n];
- assert(ipiv!=NULL);
-
- const unsigned lwork = 4*n;
- float* const work = new float[lwork];
- assert(work!=NULL);
-
- int INFO = blas::getrf(n, A, ipiv);
- assert(INFO==0);
-
- sgetri_(&n, A, &n, ipiv, work, &lwork, &INFO);
- assert(INFO==0);
-
- delete [] work;
- delete [] ipiv;
- }
-
- inline void geinv_sym(const unsigned n, double* const A)
- {
- int* const ipiv = new int[n];
- assert(ipiv!=NULL);
-
- const unsigned lwork = 4*n;
- double* const work = new double[lwork];
- assert(work!=NULL);
-
- int INFO;
- dsptrf_(JOB_STR+5, &n, A, ipiv, &INFO);
- assert(INFO==0);
-
- dsptri_(JOB_STR+5, &n, A, ipiv, work, &INFO);
- assert(INFO==0);
-
- delete [] work;
- delete [] ipiv;
- }
- inline void geinv_sym(const unsigned n, float* const A)
- {
- int* const ipiv = new int[n];
- assert(ipiv!=NULL);
-
- const unsigned lwork = 4*n;
- float* const work = new float[lwork];
- assert(work!=NULL);
-
- int INFO;
- ssptrf_(JOB_STR+5, &n, A, ipiv, &INFO);
- assert(INFO==0);
-
- ssptri_(JOB_STR+5, &n, A, ipiv, work, &INFO);
- assert(INFO==0);
-
- delete [] work;
- delete [] ipiv;
- }
-
- // packed triangular factorisation of positive definite matrix
- inline int pptrf(const unsigned n, double* A)
- {
- int inf;
- dpptrf_(JOB_STR+5, &n, A, &inf);
- return inf;
- }
- inline int pptrf(const unsigned n, float* A)
- {
- int inf;
- spptrf_(JOB_STR+5, &n, A, &inf);
- return inf;
- }
-
- // packed triangular factorization of hermitian matrix
- inline int sptrf(const unsigned n, double* A, int* ipiv)
- {
- int inf;
- dsptrf_(JOB_STR+5, &n, A, ipiv, &inf);
- return inf;
- }
- inline int sptrf(const unsigned n, float* A, int* ipiv)
- {
- int inf;
- ssptrf_(JOB_STR+5, &n, A, ipiv, &inf);
- return inf;
- }
-
- // lower triangular solve
- inline void ltrs(const unsigned n, double* A,
- const unsigned p, double* B, const unsigned ldB)
- {
- // dtptrs_(JOB_STR+6, JOB_STR, JOB_STR+5, &n, &p, A, B, &ldB, &inf);
- for (unsigned i=0; i<p; ++i)
- dtpsv_(JOB_STR+6, JOB_STR, JOB_STR+5, &n, A, B+i*ldB, &N_ONE);
- }
- inline void ltrs(const unsigned n, float* A,
- const unsigned p, float* B, const unsigned ldB)
- {
- // stptrs_(JOB_STR+6, JOB_STR, JOB_STR+5, &n, &p, A, B, &ldB, &inf);
- for (unsigned i=0; i<p; ++i)
- stpsv_(JOB_STR+6, JOB_STR, JOB_STR+5, &n, A, B+i*ldB, &N_ONE);
- }
-
- // lower triangular transpose solve
- inline void ltrhs(const unsigned n, double* A,
- const unsigned p, double* B, const unsigned ldB)
- {
- // dtptrs_(JOB_STR+6, JOB_STR+1, JOB_STR+5, &n, &p, A, B, &ldB, &inf);
- for (unsigned i=0; i<p; ++i)
- dtpsv_(JOB_STR+6, JOB_STR+1, JOB_STR+5, &n, A, B+i*ldB, &N_ONE);
- }
- inline void ltrhs(const unsigned n, float* A,
- const unsigned p, float* B, const unsigned ldB)
- {
- // stptrs_(JOB_STR+6, JOB_STR+1, JOB_STR+5, &n, &p, A, B, &ldB, &inf);
- for (unsigned i=0; i<p; ++i)
- stpsv_(JOB_STR+6, JOB_STR+1, JOB_STR+5, &n, A, B+i*ldB, &N_ONE);
- }
-
- // unit upper triangular solve (with L and R stored in one matrix)
- // XR=B, R is pxp, B is nxp
- inline void utrcs(const unsigned p, const double* LR, const unsigned ldLR,
- const unsigned n, double* X, const unsigned ldX)
- {
- dtrsm_(JOB_STR+8, JOB_STR+5, JOB_STR, JOB_STR+5, &n, &p, &D_ONE,
- LR, &ldLR, X, &ldX);
- }
-
- inline void utrcs(const unsigned p, const float* LR, const unsigned ldLR,
- const unsigned n, float* X, const unsigned ldX)
- {
- strsm_(JOB_STR+8, JOB_STR+5, JOB_STR, JOB_STR+5, &n, &p, &S_ONE,
- LR, &ldLR, X, &ldX);
- }
-
- // unit upper triangular solve (with L and R stored in one matrix)
- // RX=B, R is nxn, B is nxp
- inline void utlcs(const unsigned n, float* LR, const unsigned ldLR,
- const unsigned p, float* X, const unsigned ldX)
- {
- strsm_(JOB_STR+6, JOB_STR+5, JOB_STR, JOB_STR+5, &n, &p, &S_ONE,
- LR, &ldLR, X, &ldX);
- }
- inline void utlcs(const unsigned n, double* LR, const unsigned ldLR,
- const unsigned p, double* X, const unsigned ldX)
- {
- dtrsm_(JOB_STR+6, JOB_STR+5, JOB_STR, JOB_STR+5, &n, &p, &D_ONE,
- LR, &ldLR, X, &ldX);
- }
-
- // unit lower triangular solve (with L and R stored in one matrix)
- // XL=B, L is pxp, B is nxp
- inline void ltrcs(const unsigned p, float* LR, const unsigned ldLR,
- const unsigned n, float* X, const unsigned ldX)
- {
- strsm_(JOB_STR+8, JOB_STR+6, JOB_STR, JOB_STR+5, &n, &p, &S_ONE,
- LR, &ldLR, X, &ldX);
- }
- inline void ltrcs(const unsigned p, double* LR, const unsigned ldLR,
- const unsigned n, double* X, const unsigned ldX)
- {
- dtrsm_(JOB_STR+8, JOB_STR+6, JOB_STR, JOB_STR+5, &n, &p, &D_ONE,
- LR, &ldLR, X, &ldX);
- }
-
-
- // unit lower triangular transposed solve (with L and R stored in one matrix)
- // XL^T=B, L is pxp, B is nxp
- inline void ltrhcs(const unsigned p, double* LR, const unsigned ldLR,
- const unsigned n, double* X, const unsigned ldX)
- {
- dtrsm_(JOB_STR+8, JOB_STR+6, JOB_STR+1, JOB_STR+5, &n, &p, &D_ONE,
- LR, &ldLR, X, &ldX);
- }
-
- // unit lower triangular solve (with L and R stored in one matrix)
- // LX=B, L is nxn, B is nxp
- inline void ltlcs(const unsigned n, float* LR, const unsigned ldLR,
- const unsigned p, float* X, const unsigned ldX)
- {
- strsm_(JOB_STR+6, JOB_STR+6, JOB_STR, JOB_STR+5, &n, &p, &S_ONE,
- LR, &ldLR, X, &ldX);
- }
- inline void ltlcs(const unsigned n, double* LR, const unsigned ldLR,
- const unsigned p, double* X, const unsigned ldX)
- {
- dtrsm_(JOB_STR+6, JOB_STR+6, JOB_STR, JOB_STR+5, &n, &p, &D_ONE,
- LR, &ldLR, X, &ldX);
- }
-
- // upper triangular solve
- inline void utrs(const unsigned n, double* A,
- const unsigned p, double* B, const unsigned ldB)
- {
- // dtptrs_(JOB_STR+5, JOB_STR, JOB_STR, &n, &p, A, B, &ldB, &inf);
- for (unsigned i=0; i<p; ++i)
- dtpsv_(JOB_STR+5, JOB_STR, JOB_STR, &n, A, B+i*ldB, &N_ONE);
- }
- inline void utrs(const unsigned n, float* A,
- const unsigned p, float* B, const unsigned ldB)
- {
- //stptrs_(JOB_STR+5, JOB_STR, JOB_STR, &n, &p, A, B, &ldB, &inf);
- for (unsigned i=0; i<p; ++i)
- stpsv_(JOB_STR+5, JOB_STR, JOB_STR, &n, A, B+i*ldB, &N_ONE);
- }
-
- // upper triangluar transpose solve
- inline void utrhs(const unsigned n, double* A,
- const unsigned p, double* B, const unsigned ldB)
- {
- //dtptrs_(JOB_STR+5, JOB_STR+1, JOB_STR, &n, &p, A, B, &ldB, &inf);
- for (unsigned i=0; i<p; ++i)
- dtpsv_(JOB_STR+5, JOB_STR+1, JOB_STR, &n, A, B+i*ldB, &N_ONE);
- }
- inline void utrhs(const unsigned n, float* A,
- const unsigned p, float* B, const unsigned ldB)
- {
- //stptrs_(JOB_STR+5, JOB_STR+1, JOB_STR, &n, &p, A, B, &ldB, &inf);
- for (unsigned i=0; i<p; ++i)
- stpsv_(JOB_STR+5, JOB_STR+1, JOB_STR, &n, A, B+i*ldB, &N_ONE);
- }
-}
-
-#endif
diff --git a/MathLib/LinAlg/Solvers/solver.h b/MathLib/LinAlg/Solvers/solver.h
deleted file mode 100644
index c6717555ca8d754427302a3fbe03c829c99ded84..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/solver.h
+++ /dev/null
@@ -1,21 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-09-27
- * \brief Solver includes.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef SOLVER_H_
-#define SOLVER_H_
-
-#include "CG.h"
-#include "BiCGStab.h"
-
-#endif /* SOLVER_H_ */
diff --git a/MathLib/LinAlg/Sparse/CRSMatrix.h b/MathLib/LinAlg/Sparse/CRSMatrix.h
deleted file mode 100644
index e87e9acbbeb853318b71345bafb1407d70be440f..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/CRSMatrix.h
+++ /dev/null
@@ -1,456 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-09-20
- * \brief Definition of the CRSMatrix class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef CRSMATRIX_H
-#define CRSMATRIX_H
-
-#include <cassert>
-#include <fstream>
-#include <iosfwd>
-#include <string>
-
-// MathLib
-#include "SparseMatrixBase.h"
-#include "MatrixSparsityPattern.h"
-#include "sparse.h"
-#include "amuxCRS.h"
-#include "../Preconditioner/generateDiagPrecond.h"
-
-namespace MathLib {
-
-template<typename FP_TYPE, typename IDX_TYPE>
-class CRSMatrix: public SparseMatrixBase<FP_TYPE, IDX_TYPE>
-{
-public:
- typedef FP_TYPE FP_T;
-
-public:
- explicit CRSMatrix(std::string const &fname) :
- SparseMatrixBase<FP_TYPE, IDX_TYPE>(),
- _row_ptr(NULL), _col_idx(NULL), _data(NULL)
- {
- std::ifstream in(fname.c_str(), std::ios::in | std::ios::binary);
- if (in) {
- CS_read(in, SparseMatrixBase<FP_TYPE, IDX_TYPE>::_n_rows, _row_ptr, _col_idx, _data);
- SparseMatrixBase<FP_TYPE, IDX_TYPE>::_n_cols = SparseMatrixBase<FP_TYPE, IDX_TYPE>::_n_rows;
- in.close();
- } else {
- std::cout << "cannot open " << fname << std::endl;
- }
- }
-
- explicit CRSMatrix(IDX_TYPE n, IDX_TYPE *iA, IDX_TYPE *jA, FP_TYPE* A) :
- SparseMatrixBase<FP_TYPE, IDX_TYPE>(n,n),
- _row_ptr(iA), _col_idx(jA), _data(A)
- {}
-
- explicit CRSMatrix(MatrixSparsityPattern const& mat_sparsity_pattern) :
- SparseMatrixBase<FP_TYPE, IDX_TYPE>(mat_sparsity_pattern.getNRows(),
- mat_sparsity_pattern.getNRows()),
- _row_ptr(nullptr), _col_idx(nullptr), _data(nullptr)
- {
- // reserve memory for _row_ptr
- _row_ptr = new IDX_TYPE [this->_n_rows + 1];
- // initialize _row_ptr
- _row_ptr[0] = 0;
- for (std::size_t row(0); row < this->_n_rows; row++) {
- _row_ptr[row + 1] = _row_ptr[row]
- + std::distance(mat_sparsity_pattern.getRowBeginIterator(row),
- mat_sparsity_pattern.getRowEndIterator(row));
- }
-
- std::size_t const nnz = _row_ptr[this->_n_rows];
- // reserve memory for _col_idx
- _col_idx = new IDX_TYPE [nnz];
- // fill _col_idx
- for (std::size_t row(0); row < this->_n_rows; row++) {
- std::copy(mat_sparsity_pattern.getRowBeginIterator(row),
- mat_sparsity_pattern.getRowEndIterator(row),
- &_col_idx[_row_ptr[row]]);
- }
-
- // reserve memory for _data
- _data = new FP_TYPE [nnz];
- setZero();
- }
-
- /// Reset data entries to zero.
- virtual void setZero()
- {
- std::fill_n(_data, _row_ptr[this->_n_rows], 0);
- }
-
- virtual ~CRSMatrix()
- {
- delete [] _row_ptr;
- delete [] _col_idx;
- delete [] _data;
- }
-
- virtual void amux(FP_TYPE const d,
- FP_TYPE const* const __restrict__ x,
- FP_TYPE* __restrict__ y) const
- {
- amuxCRS<FP_TYPE, IDX_TYPE>(d, this->getNRows(), _row_ptr, _col_idx, _data, x, y);
- }
-
- virtual void precondApply(FP_TYPE* /*x*/) const
- {}
-
- /**
- * get the number of non-zero entries
- * @return number of non-zero entries
- */
- IDX_TYPE getNNZ() const { return _row_ptr[this->_n_rows]; }
-
- /**
- * This method inserts/overwrites a non-zero matrix entry.
- * Precondition: the entry have to be in the sparsity pattern!
- * @param row the row number
- * @param col the column number
- * @param val the value that should be set at pos row,col
- * @return true, if the entry is contained in the sparsity pattern, else false
- */
- bool setValue(IDX_TYPE row, IDX_TYPE col, FP_TYPE val)
- {
- assert(0 <= row && row < this->_n_rows);
-
- // linear search - for matrices with many entries per row binary search is much faster
- const IDX_TYPE idx_end (_row_ptr[row+1]);
- IDX_TYPE j(_row_ptr[row]), k;
-
- while (j<idx_end && (k=_col_idx[j]) <= col) {
- if (k == col) {
- _data[j] = val;
- return true;
- }
- j++;
- }
- return false;
- }
-
- /**
- * This method adds value val to an existing matrix entry at position row,col.
- * Precondition: the entry have to be in the sparsity pattern!
- * @param row the row number
- * @param col the column number
- * @param val the value that should be set at pos row,col
- * @return true, if the entry is contained in the sparsity pattern, else false
- */
- bool addValue(IDX_TYPE row, IDX_TYPE col, FP_TYPE val)
- {
- assert(0 <= row && row < this->_n_rows);
-
- // linear search - for matrices with many entries per row binary search is much faster
- const IDX_TYPE idx_end (_row_ptr[row+1]);
- IDX_TYPE j(_row_ptr[row]), k;
-
- while (j<idx_end && (k=_col_idx[j]) <= col) {
- if (k == col) {
- #pragma omp atomic
- _data[j] += val;
- return true;
- }
- j++;
- }
- return false;
- }
-
- /**
- * This is an access operator to a non-zero matrix entry. If the value of
- * a non-existing matrix entry is requested it will be 0.0 returned.
- * @param row the row number
- * @param col the column number
- * @return The corresponding matrix entry or 0.0.
- */
- FP_TYPE getValue(IDX_TYPE row, IDX_TYPE col)
- {
- assert(0 <= row && row < this->_n_rows);
-
- // linear search - for matrices with many entries per row binary search is much faster
- const IDX_TYPE idx_end (_row_ptr[row+1]);
- IDX_TYPE j(_row_ptr[row]), k;
-
- while (j<idx_end && (k=_col_idx[j]) <= col) {
- if (k == col) {
- return _data[j];
- }
- j++;
- }
- return 0.0;
- }
-
- /**
- * This is the constant access operator to a non-zero matrix entry.
- * Precondition: the entries have to be in the sparsity pattern!
- * @param row the row number
- * @param col the column number
- * @return The corresponding matrix entry.
- */
- FP_TYPE operator() (IDX_TYPE row, IDX_TYPE col) const
- {
- assert(0 <= row && row < this->_n_rows);
-
- // linear search - for matrices with many entries per row binary search is much faster
- const IDX_TYPE idx_end (_row_ptr[row+1]);
- IDX_TYPE j(_row_ptr[row]), k;
-
- while (j<idx_end && (k=_col_idx[j]) <= col) {
- if (k == col) {
- return _data[j];
- }
- j++;
- }
- return 0.0;
- }
-
- /**
- * get const access to the row pointer array of CRS matrix
- * @return the index array _row_ptr
- */
- IDX_TYPE const* getRowPtrArray() const { return _row_ptr; }
-
- /**
- * get const access to the column index array of CRS matrix
- * @return the index array _col_idx
- */
- IDX_TYPE const* getColIdxArray ()const { return _col_idx; }
-
- /**
- * get the matrix entries within an array of CRS matrix
- */
- FP_TYPE const* getEntryArray() const { return _data; }
-
- /**
- * erase rows and columns from sparse matrix
- * @param n_rows_cols number of rows / columns to remove
- * @param rows_cols sorted list of rows/columns that should be removed
- */
- void eraseEntries(IDX_TYPE n_rows_cols, IDX_TYPE const* const rows_cols)
- {
- //*** remove the rows
- removeRows(n_rows_cols, rows_cols);
- //*** transpose
- transpose();
- //*** remove columns in original means removing rows in the transposed
- removeRows(n_rows_cols, rows_cols);
- //*** transpose again
- transpose();
- }
-
- /**
- * get the j-th column of the sparse matrix
- * @param j the column number that should be returned
- * @param column_entries the column entries (have to be allocated
- */
- void getColumn(IDX_TYPE j, FP_TYPE* column_entries) const
- {
- for (IDX_TYPE k(0); k < this->_n_rows; k++) {
- const IDX_TYPE end_row(_row_ptr[k+1]);
- IDX_TYPE i(_row_ptr[k+1]);
- while (i<end_row && _col_idx[i] != j) {
- i++;
- }
- if (i==end_row) {
- column_entries[k] = 0.0;
- } else {
- column_entries[k] = _data[i];
- }
- }
- }
-
- CRSMatrix<FP_TYPE, IDX_TYPE>* getTranspose() const
- {
- CRSMatrix<FP_TYPE, IDX_TYPE>* transposed_mat(new CRSMatrix<FP_TYPE, IDX_TYPE>(*this));
- transposed_mat->transpose();
- return transposed_mat;
- }
-
-protected:
- CRSMatrix(CRSMatrix const& rhs) :
- SparseMatrixBase<FP_TYPE, IDX_TYPE> (rhs.getNRows(), rhs.getNCols()),
- _row_ptr(new IDX_TYPE[rhs.getNRows() + 1]), _col_idx(new IDX_TYPE[rhs.getNNZ()]),
- _data(new FP_TYPE[rhs.getNNZ()])
- {
- // copy the data
- IDX_TYPE const* row_ptr(rhs.getRowPtrArray());
- for (IDX_TYPE k(0); k <= this->_n_rows; k++) {
- _row_ptr[k] = row_ptr[k];
- }
-
- IDX_TYPE nnz(rhs.getNNZ());
- IDX_TYPE const*const col_idx(rhs.getColIdxArray());
- for (IDX_TYPE k(0); k<nnz; k++) {
- _col_idx[k] = col_idx[k];
- }
-
- FP_TYPE const*const data(rhs.getEntryArray());
- for (IDX_TYPE k(0); k<nnz; k++) {
- _data[k] = data[k];
- }
- }
-
- void removeRows (IDX_TYPE n_rows_cols, IDX_TYPE const*const rows)
- {
- //*** determine the number of new rows and the number of entries without the rows
- const IDX_TYPE n_new_rows(this->_n_rows - n_rows_cols);
- IDX_TYPE *row_ptr_new(new IDX_TYPE[n_new_rows+1]);
- row_ptr_new[0] = 0;
- IDX_TYPE row_cnt (1), erase_row_cnt(0);
- for (unsigned k(0); k < this->_n_rows; k++) {
- if (erase_row_cnt < n_rows_cols) {
- if (k != rows[erase_row_cnt]) {
- row_ptr_new[row_cnt] = _row_ptr[k+1] - _row_ptr[k];
- row_cnt++;
- } else {
- erase_row_cnt++;
- }
- } else {
- row_ptr_new[row_cnt] = _row_ptr[k+1] - _row_ptr[k];
- row_cnt++;
- }
- }
-
- //*** sum up the entries
- for (IDX_TYPE k(0); k<n_new_rows; k++) {
- row_ptr_new[k+1] = row_ptr_new[k+1] + row_ptr_new[k];
- }
-
- //*** create new memory for col_idx and data
- IDX_TYPE nnz_new(row_ptr_new[n_new_rows]);
- IDX_TYPE *col_idx_new (new IDX_TYPE[nnz_new]);
- FP_TYPE *data_new (new FP_TYPE[nnz_new]);
-
- //*** copy the entries
- // initialization
- IDX_TYPE *row_ptr_new_tmp(new IDX_TYPE[n_new_rows+1]);
- for (unsigned k(0); k<=n_new_rows; k++) {
- row_ptr_new_tmp[k] = row_ptr_new[k];
- }
- erase_row_cnt = 0;
- row_cnt = 0;
- // copy column index and data entries
- for (IDX_TYPE k(0); k < this->_n_rows; k++) {
- if (erase_row_cnt < n_rows_cols) {
- if (k != rows[erase_row_cnt]) {
- const IDX_TYPE end (_row_ptr[k+1]);
- // walk through row
- for (IDX_TYPE j(_row_ptr[k]); j<end; j++) {
- col_idx_new[row_ptr_new_tmp[row_cnt]] = _col_idx[j];
- data_new[row_ptr_new_tmp[row_cnt]] = _data[j];
- row_ptr_new_tmp[row_cnt]++;
- }
- row_cnt++;
- } else {
- erase_row_cnt++;
- }
- } else {
- const IDX_TYPE end (_row_ptr[k+1]);
- // walk through row
- for (IDX_TYPE j(_row_ptr[k]); j<end; j++) {
- col_idx_new[row_ptr_new_tmp[row_cnt]] = _col_idx[j];
- data_new[row_ptr_new_tmp[row_cnt]] = _data[j];
- row_ptr_new_tmp[row_cnt]++;
- }
- row_cnt++;
- }
- }
-
- this->_n_rows -= n_rows_cols;
- std::swap (row_ptr_new, _row_ptr);
- std::swap (col_idx_new, _col_idx);
- std::swap (data_new, _data);
-
- delete [] row_ptr_new_tmp;
- delete [] row_ptr_new;
- delete [] col_idx_new;
- delete [] data_new;
- }
-
- void transpose ()
- {
- // create a helper array row_ptr_nnz
- IDX_TYPE *row_ptr_nnz(new IDX_TYPE[this->_n_cols+1]);
- for (IDX_TYPE k(0); k <= this->_n_cols; k++) {
- row_ptr_nnz[k] = 0;
- }
-
- // count entries per row in the transposed matrix
- IDX_TYPE nnz(_row_ptr[this->_n_rows]);
- for (IDX_TYPE k(0); k < nnz; k++) {
- row_ptr_nnz[_col_idx[k]]++;
- }
-
- // create row_ptr_trans
- IDX_TYPE *row_ptr_trans(new IDX_TYPE[this->_n_cols + 1]);
- row_ptr_trans[0] = 0;
- for (IDX_TYPE k(0); k < this->_n_cols; k++) {
- row_ptr_trans[k+1] = row_ptr_trans[k] + row_ptr_nnz[k];
- }
-
- // make a copy of row_ptr_trans
- for (IDX_TYPE k(0); k <= this->_n_cols; k++) {
- row_ptr_nnz[k] = row_ptr_trans[k];
- }
-
- // create arrays col_idx_trans and data_trans
- assert(nnz == row_ptr_trans[this->_n_cols]);
- IDX_TYPE *col_idx_trans(new IDX_TYPE[nnz]);
- FP_TYPE *data_trans(new FP_TYPE[nnz]);
-
- // fill arrays col_idx_trans and data_trans
- for (IDX_TYPE i(0); i < this->_n_rows; i++) {
- const IDX_TYPE row_end(_row_ptr[i + 1]);
- for (IDX_TYPE j(_row_ptr[i]); j < row_end; j++) {
- const IDX_TYPE k(_col_idx[j]);
- col_idx_trans[row_ptr_nnz[k]] = i;
- data_trans[row_ptr_nnz[k]] = _data[j];
- row_ptr_nnz[k]++;
- }
- }
-
- std::swap(this->_n_rows, this->_n_cols);
- std::swap(row_ptr_trans, _row_ptr);
- std::swap(col_idx_trans, _col_idx);
- std::swap(data_trans, _data);
-
- delete[] row_ptr_nnz;
- delete[] row_ptr_trans;
- delete[] col_idx_trans;
- delete[] data_trans;
- }
-
-#ifndef NDEBUG
- void printMat() const
- {
- for (IDX_TYPE k(0); k < this->_n_rows; k++) {
- std::cout << k << ": " << std::flush;
- const IDX_TYPE row_end(_row_ptr[k+1]);
- for (IDX_TYPE j(_row_ptr[k]); j<row_end; j++) {
- std::cout << _col_idx[j] << " " << std::flush;
- }
- std::cout << std::endl;
- }
- }
-#endif
-
- IDX_TYPE *_row_ptr;
- IDX_TYPE *_col_idx;
- FP_TYPE* _data;
-};
-
-} // end namespace MathLib
-
-#endif
-
diff --git a/MathLib/LinAlg/Sparse/CRSMatrixDiagPrecond.h b/MathLib/LinAlg/Sparse/CRSMatrixDiagPrecond.h
deleted file mode 100644
index cc00ec9cf2ddb0976739c2be73cc3cadc8821817..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/CRSMatrixDiagPrecond.h
+++ /dev/null
@@ -1,87 +0,0 @@
-#ifndef CRSMATRIXDIAGPRECOND_H
-#define CRSMATRIXDIAGPRECOND_H
-
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
-#include "CRSMatrix.h"
-#include "sparse.h"
-#include "../Preconditioner/generateDiagPrecond.h"
-
-namespace MathLib {
-
-/**
- * Class CRSMatrixDiagPrecond represents a matrix in compressed row storage
- * format associated with a diagonal preconditioner.
- *
- * The user can either read the matrix from a file or give corresponding arrays
- * to an alternative constructor. In both cases the user have to calculate the
- * preconditioner explicit via calcPrecond() method!
- */
-class CRSMatrixDiagPrecond : public CRSMatrix<double, unsigned>
-{
-public:
- /**
- * Constructor takes a file name. The file is read in binary format
- * by the constructor of the base class (template) CRSMatrix.
- * Further details you can see in function CS_read().
- *
- * The user have to calculate the preconditioner explicit via calcPrecond() method!
- *
- * @param fname the name of the file that contains the matrix in
- * binary compressed row storage format
- */
- CRSMatrixDiagPrecond(std::string const &fname) :
- CRSMatrix<double, unsigned> (fname), _inv_diag(NULL)
- {}
-
- /**
- * Constructs a matrix object from given data.
- *
- * The user have to calculate the preconditioner explicit via calcPrecond() method!
- * @param n number of rows / columns of the matrix
- * @param iA row pointer of matrix in compressed row storage format
- * @param jA column index of matrix in compressed row storage format
- * @param A data entries of matrix in compressed row storage format
- */
- CRSMatrixDiagPrecond(unsigned n, unsigned *iA, unsigned *jA, double* A) :
- CRSMatrix<double, unsigned> (n, iA, jA, A), _inv_diag(NULL)
- {}
-
- void calcPrecond()
- {
- delete [] _inv_diag;
- _inv_diag = new double[_n_rows];
-
- if (!generateDiagPrecond(_n_rows, _row_ptr, _col_idx, _data, _inv_diag)) {
- std::cout << "Could not create diagonal preconditioner" << std::endl;
- }
-// if (!generateDiagPrecondRowSum(_n_rows, _row_ptr, _data, _inv_diag)) {
-// std::cout << "Could not create diagonal preconditioner" << std::endl;
-// }
-// if (!generateDiagPrecondRowMax(_n_rows, _row_ptr, _data, _inv_diag)) {
-// std::cout << "Could not create diagonal preconditioner" << std::endl;
-// }
-
- }
-
- void precondApply(double* x) const
- {
- for (unsigned k=0; k<_n_rows; ++k) {
- x[k] = _inv_diag[k]*x[k];
- }
- }
-
- ~CRSMatrixDiagPrecond()
- {
- delete [] _inv_diag;
- }
-private:
- double *_inv_diag;
-};
-
-}
-
-#endif
-
diff --git a/MathLib/LinAlg/Sparse/CRSMatrixOpenMP.h b/MathLib/LinAlg/Sparse/CRSMatrixOpenMP.h
deleted file mode 100644
index 66f4a04aaf1e5577bbaf01473564f039d35102fd..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/CRSMatrixOpenMP.h
+++ /dev/null
@@ -1,53 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-08-08
- * \brief Definition of the CRSMatrixOpenMP class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef CRSMATRIXOPENMP_H_
-#define CRSMATRIXOPENMP_H_
-
-#ifdef _OPENMP
-#include <string>
-
-#include "CRSMatrix.h"
-#include "amuxCRS.h"
-
-namespace MathLib {
-
-template<typename FP_TYPE, typename IDX_TYPE> class CRSMatrixOpenMP : public CRSMatrix<FP_TYPE, IDX_TYPE> {
-public:
- CRSMatrixOpenMP(std::string const &fname) :
- CRSMatrix<FP_TYPE, IDX_TYPE>(fname)
- {}
-
- CRSMatrixOpenMP(unsigned n, IDX_TYPE *iA, IDX_TYPE *jA, FP_TYPE* A) :
- CRSMatrix<FP_TYPE, IDX_TYPE>(n, iA, jA, A)
- {}
-
- CRSMatrixOpenMP(unsigned n1) :
- CRSMatrix<FP_TYPE, IDX_TYPE>(n1)
- {}
-
- virtual ~CRSMatrixOpenMP()
- {}
- virtual void amux(FP_TYPE const d, FP_TYPE const* const __restrict__ x,
- FP_TYPE* __restrict__ y) const
- {
- amuxCRSParallelOpenMP(d, this->_n_rows, this->_row_ptr, this->_col_idx, this->_data, x, y);
- }
-};
-
-} // end namespace MathLib
-
-#endif // _OPENMP
-
-#endif /* CRSMATRIXOPENMP_H_ */
diff --git a/MathLib/LinAlg/Sparse/CRSMatrixPThreads.h b/MathLib/LinAlg/Sparse/CRSMatrixPThreads.h
deleted file mode 100644
index 377914bbcf11ff9e15afc8c8808d2c8d80c9ba6e..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/CRSMatrixPThreads.h
+++ /dev/null
@@ -1,108 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-08-02
- * \brief Definition of the CRSMatrixPThreads class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef CRSMATRIXPTHREADS_H
-#define CRSMATRIXPTHREADS_H
-
-#include <string>
-
-#include "SparseMatrixBase.h"
-#include "sparse.h"
-#include "CRSMatrix.h"
-#include "amuxCRS.h"
-
-namespace MathLib {
-
-template<class T> class CRSMatrixPThreads : public CRSMatrix<T,unsigned>
-{
-public:
- CRSMatrixPThreads(std::string const &fname, unsigned num_of_threads) :
- CRSMatrix<T,unsigned>(fname), _n_threads (num_of_threads),
- _workload_intervals(new unsigned[num_of_threads+1])
- {
- calcWorkload();
- }
-
- CRSMatrixPThreads(unsigned n, unsigned *iA, unsigned *jA, T* A, unsigned num_of_threads) :
- CRSMatrix<T,unsigned>(n, iA, jA, A), _n_threads (num_of_threads),
- _workload_intervals(new unsigned[num_of_threads+1])
- {
- calcWorkload();
- }
-
- CRSMatrixPThreads(unsigned n1) :
- CRSMatrix<T,unsigned>(n1), _n_threads (1),
- _workload_intervals(new unsigned[_n_threads+1])
- {
- calcWorkload();
- }
-
- virtual ~CRSMatrixPThreads()
- {
- delete [] _workload_intervals;
- }
-
- virtual void amux(T d, T const * const x, T *y) const
- {
- amuxCRSParallelPThreads(d, SparseMatrixBase<T, unsigned>::_n_rows,
- CRSMatrix<T, unsigned>::_row_ptr, CRSMatrix<T, unsigned>::_col_idx,
- CRSMatrix<T, unsigned>::_data, x, y, _n_threads, _workload_intervals);
- }
-
-protected:
- void calcWorkload()
- {
- _workload_intervals[0] = 0;
- _workload_intervals[_n_threads] = SparseMatrixBase<T, unsigned>::_n_rows;
-
- const unsigned work_per_core (this->getNNZ()/_n_threads);
- for (unsigned k(1); k<_n_threads; k++) {
- unsigned upper_bound_kth_core(k * work_per_core);
- // search in _row_ptr array for the appropriate index
- unsigned beg (_workload_intervals[k-1]);
- unsigned end (_workload_intervals[_n_threads]);
- bool found (false);
- while (beg < end && !found) {
- unsigned m ((end+beg)/2);
-
- if (upper_bound_kth_core == this->_row_ptr[m]) {
- _workload_intervals[k] = m;
- found = true;
- } else {
- if (upper_bound_kth_core < this->_row_ptr[m]) {
- end = m;
- } else {
- beg = m+1;
- }
- }
- }
- if (!found)
- _workload_intervals[k] = beg;
- }
-
- for (unsigned k(0); k<_n_threads; k++) {
- std::cout << "proc " << k << ": [" << _workload_intervals[k] << "," << _workload_intervals[k+1] << ") - "
- << _workload_intervals[k+1] - _workload_intervals[k] << " rows and "
- << this->_row_ptr[_workload_intervals[k+1]] - this->_row_ptr[_workload_intervals[k]] << " entries" << std::endl;
- }
- }
-
- const unsigned _n_threads;
- unsigned *_workload_intervals;
-};
-
-} // end namespace MathLib
-
-#endif
-
diff --git a/MathLib/LinAlg/Sparse/CRSSymMatrix.h b/MathLib/LinAlg/Sparse/CRSSymMatrix.h
deleted file mode 100644
index fafb9e6d3ed2980494ade06d2bd17514f63d9ced..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/CRSSymMatrix.h
+++ /dev/null
@@ -1,74 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-07-22
- * \brief Definition of the CRSSymMatrix class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef CRSSYMMATRIX_H_
-#define CRSSYMMATRIX_H_
-
-#include "CRSMatrix.h"
-
-template<class T> class CRSSymMatrix : public CRSMatrix<T>
-{
-public:
- typedef T FP_T;
-
-public:
- CRSSymMatrix(std::string const &fname)
- : CRSMatrix<T> (fname)
- {
- unsigned nnz (0);
-
- // count number of non-zeros in the upper triangular part
- for (unsigned i = 0; i < SparseMatrixBase<T>::_n_rows; i++) {
- unsigned idx = CRSMatrix<T>::_row_ptr[i+1];
- for (unsigned j = CRSMatrix<T>::_row_ptr[i]; j < idx; j++)
- if (CRSMatrix<T>::_col_idx[j] >= i)
- ++nnz;
- }
-
- double *A_new (new double[nnz]);
- unsigned *jA_new (new unsigned[nnz]);
- unsigned *iA_new (new unsigned[SparseMatrixBase<T>::_n_rows+1]);
-
- iA_new[0] = nnz = 0;
-
- for (unsigned i = 0; i < SparseMatrixBase<T>::_n_rows; i++) {
- const unsigned idx (CRSMatrix<T>::_row_ptr[i+1]);
- for (unsigned j = CRSMatrix<T>::_row_ptr[i]; j < idx; j++) {
- if (CRSMatrix<T>::_col_idx[j] >= i) {
- A_new[nnz] = CRSMatrix<T>::_data[j];
- jA_new[nnz++] = CRSMatrix<T>::_col_idx[j];
- }
- }
- iA_new[i+1] = nnz;
- }
-
- std::swap(CRSMatrix<T>::_row_ptr, iA_new);
- std::swap(CRSMatrix<T>::_col_idx, jA_new);
- std::swap(CRSMatrix<T>::_data, A_new);
-
- delete[] iA_new;
- delete[] jA_new;
- delete[] A_new;
- }
-
- virtual ~CRSSymMatrix() {}
-
- void amux(T d, T const * const x, T *y) const
- {
- amuxCRSSym (d, SparseMatrixBase<T>::_n_rows, CRSMatrix<T>::_row_ptr, CRSMatrix<T>::_col_idx, CRSMatrix<T>::_data, x, y);
- }
-
-};
-
-#endif /* CRSSYMMATRIX_H_ */
diff --git a/MathLib/LinAlg/Sparse/CRSTools-impl.h b/MathLib/LinAlg/Sparse/CRSTools-impl.h
deleted file mode 100644
index 5316970ab3bc09caf4a8059a8b9080f1220185a1..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/CRSTools-impl.h
+++ /dev/null
@@ -1,70 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <memory>
-#include "logog/include/logog.hpp"
-
-#include "CRSTools.h"
-
-#include "CRSMatrix.h"
-
-namespace MathLib
-{
-
-template <typename VEC_T, typename FP_TYPE>
-void applyKnownSolution(CRSMatrix<FP_TYPE, typename VEC_T::IndexType>*& mat,
- VEC_T &rhs, std::vector<typename VEC_T::IndexType> const& rows,
- std::vector<FP_TYPE> const& vals)
-{
- std::unique_ptr<MathLib::CRSMatrix<FP_TYPE, typename VEC_T::IndexType>> mat_t(mat->getTranspose());
-
- // row pointer of the transposed matrix
- typename VEC_T::IndexType const*const iAt(mat_t->getRowPtrArray());
- // column indices of the transposed matrix
- typename VEC_T::IndexType const*const jAt(mat_t->getColIdxArray());
- // entries of the transposed matrix
- double * At(const_cast<double *>(mat_t->getEntryArray()));
-
- // b_i -= A(i,k)*val, i!=k => b_i -= A(k,i)^T * val
- // set A^T(k,i) = 0, i!=k (i.e. set column entries of original matrix to
- // zero)
- for (std::size_t r(0); r<rows.size(); ++r) {
- auto const row = rows[r];
- auto const val = vals[r];
- for (typename VEC_T::IndexType j(iAt[row]); j<iAt[row+1]; ++j) {
- if (jAt[j] == row) // skip diagonal entry
- continue;
- rhs.add(jAt[j], -At[j] * val);
- At[j] = 0.0;
- }
- }
-
- delete mat;
- mat = mat_t->getTranspose();
- typename VEC_T::IndexType const*const iA(mat->getRowPtrArray()); // row ptrs
- typename VEC_T::IndexType const*const jA(mat->getColIdxArray()); // col idx
- double * entries(const_cast<double*>(mat->getEntryArray()));
-
- // set row entries, except the diagonal entry, to zero
- for (std::size_t r(0); r<rows.size(); ++r) {
- auto const row = rows[r];
- for (typename VEC_T::IndexType j = iA[row]; j < iA[row + 1]; ++j)
- {
- if (jA[j] == row) {
- entries[j] = 1.0; // A(row,row) = 1.0
- // rhs[row] = A(row,row) * vals[r]
- rhs.set(row, vals[r]);
- } else
- entries[j] = 0.0;
- }
- }
-}
-
-} // MathLib
-
diff --git a/MathLib/LinAlg/Sparse/CRSTools.h b/MathLib/LinAlg/Sparse/CRSTools.h
deleted file mode 100644
index 737e79e720d0a667897a7c03af5b43d15ffbcd2d..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/CRSTools.h
+++ /dev/null
@@ -1,46 +0,0 @@
-/**
- * \brief Integration of Dirichlet boundary conditions into the system of
- * linear equations.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef CRSTOOLS_H_
-#define CRSTOOLS_H_
-
-#include <vector>
-
-namespace MathLib
-{
-
-template<typename FP_TYPE, typename IDX_TYPE> class CRSMatrix;
-template<typename IDX_TYPE> class DenseVector;
-
-/**
- * Integrate Dirichlet boundary conditions into a system of linear equations.
- *
- * This function introduces the constants into the system by setting
- * appropriated row and column entries of the matrix to zero (except the
- * diagonal entries) and modifying values within the right hand side vector.
- *
- * @param mat Coefficient matrix
- * @param rhs RHS vector
- * @param rows a vector of known solution entry IDs
- * @param vals a vector of known solutions
- */
-template <typename VEC_T, typename FP_TYPE = double>
-void applyKnownSolution(CRSMatrix<FP_TYPE, typename VEC_T::IndexType>*& mat,
- VEC_T &rhs, std::vector<std::size_t> const& rows,
- std::vector<FP_TYPE> const& vals);
-
-} // MathLib
-
-#include "CRSTools-impl.h"
-
-#endif // CRSTOOLS_H_
-
diff --git a/MathLib/LinAlg/Sparse/CRSTranspose.h b/MathLib/LinAlg/Sparse/CRSTranspose.h
deleted file mode 100644
index b7591e620ef2041f96ac95f1a884d3621c31ec87..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/CRSTranspose.h
+++ /dev/null
@@ -1,51 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-09-27
- * \brief Definition of the CRSTranspose class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef CRSTRANSPOSE_H_
-#define CRSTRANSPOSE_H_
-
-void CS_transp(unsigned n, unsigned *iA, unsigned* jA, double* A,
- unsigned *iB, unsigned *jB, double* B)
-{
- unsigned nnz = iA[n];
- unsigned *inz(new unsigned[n]);
-
- for (unsigned i=0; i<n; ++i) inz[i] = 0;
-
- // compute number of entries of each column in A
- for (unsigned l=0; l<nnz; ++l) inz[jA[l]]++;
-
- // create iB
- iB[0] = nnz = 0;
- for (unsigned i=0; i<n; ++i) {
- nnz += inz[i];
- inz[i] = 0;
- iB[i+1] = nnz;
- }
-
- // create arrays jB, B
- unsigned l = iA[0];
- for (unsigned i=0; i<n; ++i) {
- while (l<iA[i+1]) {
- unsigned j = jA[l];
- unsigned k = iB[j] + inz[j]++;
- B[k] = A[l++];
- jB[k] = i;
- }
- }
-
- delete [] inz;
-}
-
-#endif /* CRSTRANSPOSE_H_ */
diff --git a/MathLib/LinAlg/Sparse/MatrixSparsityPattern.cpp b/MathLib/LinAlg/Sparse/MatrixSparsityPattern.cpp
deleted file mode 100644
index 3ab04e23b6762d889a56c6b60fa6dc69f38b938a..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/MatrixSparsityPattern.cpp
+++ /dev/null
@@ -1,46 +0,0 @@
-/**
- * @file MatrixSparsityPattern.cpp
- * @author Thomas Fischer
- * @date Apr 12, 2013
- *
- * @copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/LICENSE.txt
- */
-
-#include "LinAlg/Sparse/MatrixSparsityPattern.h"
-
-namespace MathLib
-{
-
-MatrixSparsityPattern::MatrixSparsityPattern(std::size_t const n_rows) :
- _pattern(n_rows)
-{}
-
-MatrixSparsityPattern::~MatrixSparsityPattern()
-{}
-
-std::size_t MatrixSparsityPattern::getNRows() const
-{
- return _pattern.size();
-}
-
-MatrixSparsityPattern::ConstRowIterator MatrixSparsityPattern::getRowBeginIterator(
- std::size_t const row) const
-{
- return _pattern[row].cbegin();
-}
-MatrixSparsityPattern::ConstRowIterator MatrixSparsityPattern::getRowEndIterator(
- std::size_t const row) const
-{
- return _pattern[row].cend();
-}
-
-void MatrixSparsityPattern::insert(std::size_t const row, std::size_t const col)
-{
- _pattern[row].insert(col);
-}
-
-} // end namespace MathLib
diff --git a/MathLib/LinAlg/Sparse/MatrixSparsityPattern.h b/MathLib/LinAlg/Sparse/MatrixSparsityPattern.h
deleted file mode 100644
index 853487f057ae48d05214a371c94fb10e3f407f3d..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/MatrixSparsityPattern.h
+++ /dev/null
@@ -1,60 +0,0 @@
-/**
- * @file MatrixSparsityPattern.h
- * @author Thomas Fischer
- * @date Apr 12, 2013
- *
- * @copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/LICENSE.txt
- */
-
-#ifndef MATRIXSPARSITYPATTERN_H_
-#define MATRIXSPARSITYPATTERN_H_
-
-#include <set>
-#include <vector>
-
-#include "BaseLib/CodingTools.h"
-
-namespace MathLib
-{
-
-/// \brief Class for representation of matrix sparsity pattern, required for
-/// creation of sparse matrices.
-///
-/// \details Current implementation requires only number of matrix's rows,
-/// allowing thus non-rectangular sparsity patterns. The class is based on
-/// std::set container which automatically check for multiple insertions.
-class MatrixSparsityPattern
-{
-public:
- /// Constant iterator over sorted entries of a row.
- typedef std::set<std::size_t>::const_iterator ConstRowIterator;
-
- explicit MatrixSparsityPattern(std::size_t const n_rows);
- virtual ~MatrixSparsityPattern();
-
- /// Returns number of sparsity pattern rows.
- std::size_t getNRows() const;
-
- /// Constant iterator over sorted entries of a row.
- ConstRowIterator getRowBeginIterator(std::size_t const row) const;
- /// Constant iterator over sorted entries of a row.
- ConstRowIterator getRowEndIterator(std::size_t const row) const;
-
- /// Inserts an entry in the sparsity pattern.
- /// \param row The row index the entry should be inserted to. The row index must be less or equal to the value returned by getNRows().
- /// \param col The column index. A new entry will be created if needed.
- void insert(std::size_t const row, std::size_t const col);
-
-private:
- DISALLOW_COPY_AND_ASSIGN(MatrixSparsityPattern);
-
- std::vector<std::set<std::size_t> > _pattern;
-};
-
-} // end namespace MathLib
-
-#endif // MATRIXSPARSITYPATTERN_H_
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.cpp b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.cpp
deleted file mode 100644
index feb93e94c7d1f4a48ef54066380440301556f1fc..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.cpp
+++ /dev/null
@@ -1,211 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-02
- * \brief Implementation of the AdjMat class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "AdjMat.h"
-
-#include <algorithm>
-#include <utility>
-
-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)
- std::sort(jAn+iAn[i], jAn+iAn[i+1]);
- return new AdjMat(nsize, iAn, jAn, nullptr);
-}
-
-/**
- * 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)
- */
-static
-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];
-
- std::swap(jA, jAn);
- std::swap(iA, iAn);
-
- delete[] jAn;
- delete[] con;
- delete[] co;
- delete[] iAn;
-}
-
-static
-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++;
- }
- }
-
- std::swap(jA, jAn);
- std::swap(iA, iAn);
-
- delete[] jAn;
- delete[] iAn;
- delete[] cnt;
-}
-
-void AdjMat::makeSymmetric()
-{
- // store upper triangular mat values
- genAdjMat(this->_n_rows, _row_ptr, _col_idx);
- // mirror the upper triangular part into lower
- genFullAdjMat(this->_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
deleted file mode 100644
index 116154af20ff2b1c8a5ec5e1e157eb55b3399062..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.h
+++ /dev/null
@@ -1,70 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-02
- * \brief Definition of the AdjMat class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef ADJMAT_H_
-#define ADJMAT_H_
-
-#include "MathLib/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()
- {}
-
- /**
- *
- */
- 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:
- DISALLOW_COPY_AND_ASSIGN(AdjMat);
-};
-
-}
-
-#endif /* ADJMAT_H_ */
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReordered.cpp b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReordered.cpp
deleted file mode 100644
index 6fcfde26e93acb3d4640d91da6e67bf2bb52dfb4..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReordered.cpp
+++ /dev/null
@@ -1,79 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-03
- * \brief Implementation of the CRSMatrixReordered class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "CRSMatrixReordered.h"
-
-#include <utility>
-
-#include "BaseLib/quicksort.h"
-
-namespace MathLib {
-
-CRSMatrixReordered::CRSMatrixReordered(std::string const &fname) :
- CRSMatrix<double,unsigned>(fname)
-{}
-
-CRSMatrixReordered::CRSMatrixReordered(unsigned n, unsigned *iA, unsigned *jA, double* A) :
- CRSMatrix<double, unsigned> (n, iA, jA, A)
-{}
-
-CRSMatrixReordered::~CRSMatrixReordered()
-{}
-
-void CRSMatrixReordered::reorderMatrix(unsigned const*const op_perm, unsigned const*const po_perm)
-{
- unsigned i; // row and col idx in permuted matrix
- unsigned j, idx; // pointer in jA
-
- const unsigned size(getNRows());
-
- unsigned *pos(new unsigned[size + 1]);
- for (i = 0; i < size; i++) {
- const unsigned original_row(op_perm[i]);
- pos[i] = _row_ptr[original_row+1] - _row_ptr[original_row];
- }
- pos[size] = 0;
-
- unsigned *iAn(new unsigned[size + 1]);
- iAn[0] = 0;
- for (i = 0; i < size; i++)
- iAn[i + 1] = iAn[i] + pos[i];
- for (i = 0; i < size; i++)
- pos[i] = iAn[i];
-
- unsigned *jAn(new unsigned[iAn[size]]);
- double *An(new double[iAn[size]]);
- for (i = 0; i < size; i++) {
- const unsigned original_row(op_perm[i]);
- idx = _row_ptr[original_row+1];
- for (j = _row_ptr[original_row]; j < idx; j++) {
- jAn[pos[i]] = po_perm[_col_idx[j]];
- An[pos[i]++] = _data[j];
- }
- }
-
- delete[] pos;
- for (i = 0; i < size; ++i)
- BaseLib::quicksort(jAn, static_cast<std::size_t>(iAn[i]), static_cast<std::size_t>(iAn[i + 1]), An);
-
- std::swap(iAn, _row_ptr);
- std::swap(jAn, _col_idx);
- std::swap(An, _data);
-
- delete [] iAn;
- delete [] jAn;
- delete [] An;
-}
-
-} // end namespace MathLib
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReordered.h b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReordered.h
deleted file mode 100644
index 1acf92378beaf8e753cea901c7e5f5b089d371c8..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReordered.h
+++ /dev/null
@@ -1,35 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-03
- * \brief Definition of the CRSMatrixReordered class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef CRSMATRIXREORDERED_H_
-#define CRSMATRIXREORDERED_H_
-
-#include <string>
-
-#include "MathLib/LinAlg/Sparse/CRSMatrix.h"
-
-namespace MathLib {
-
-class CRSMatrixReordered: public MathLib::CRSMatrix<double,unsigned>
-{
-public:
- CRSMatrixReordered(std::string const &fname);
- CRSMatrixReordered(unsigned n, unsigned *iA, unsigned *jA, double* A);
- virtual ~CRSMatrixReordered();
- void reorderMatrix(unsigned const*const op_perm, unsigned const*const po_perm);
-};
-
-}
-
-#endif /* CRSMATRIXREORDERED_H_ */
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReorderedOpenMP.cpp b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReorderedOpenMP.cpp
deleted file mode 100644
index 3c1460603e784a0eacadf8de504f96ce273ce41d..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReorderedOpenMP.cpp
+++ /dev/null
@@ -1,36 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-12
- * \brief Implementation of the CRSMatrixReorderedOpenMP class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifdef _OPENMP
-#include "LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReorderedOpenMP.h"
-#include "LinAlg/Sparse/amuxCRS.h"
-
-namespace MathLib {
-
-CRSMatrixReorderedOpenMP::CRSMatrixReorderedOpenMP(unsigned n, unsigned *iA, unsigned *jA, double* A) :
- CRSMatrixReordered(n, iA, jA, A)
-{
-}
-
-CRSMatrixReorderedOpenMP::~CRSMatrixReorderedOpenMP()
-{}
-
-void CRSMatrixReorderedOpenMP::amux(double d, double const * const x, double *y) const
-{
- amuxCRSParallelOpenMP(d, this->_n_rows, this->_row_ptr, this->_col_idx, this->_data, x, y);
-}
-
-} // end namespace MathLib
-
-#endif // _OPENMP
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReorderedOpenMP.h b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReorderedOpenMP.h
deleted file mode 100644
index fc0b86537314d286fb190f3e9752250d721030a8..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReorderedOpenMP.h
+++ /dev/null
@@ -1,34 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-20
- * \brief Definition of the CRSMatrixReorderedOpenMP class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef CRSMATRIXREORDEREDOPENMP_H_
-#define CRSMATRIXREORDEREDOPENMP_H_
-
-#ifdef _OPENMP
-#include "CRSMatrixReordered.h"
-
-namespace MathLib {
-
-class CRSMatrixReorderedOpenMP : public CRSMatrixReordered {
-public:
- CRSMatrixReorderedOpenMP(unsigned n, unsigned *iA, unsigned *jA, double* A);
- virtual ~CRSMatrixReorderedOpenMP();
-
- virtual void amux(double d, double const * const __restrict__ x, double *__restrict__ y) const;
-};
-
-}
-#endif // _OPENMP
-
-#endif /* CRSMATRIXREORDEREDOPENMP_H_ */
diff --git a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.cpp b/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.cpp
deleted file mode 100644
index e1148e85785e49c1328e26bdca04120fc113ad2f..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.cpp
+++ /dev/null
@@ -1,215 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-02
- * \brief Implementation of the Cluster class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "metis.h"
-
-// BaseLib
-
-#include "LinAlg/Sparse/NestedDissectionPermutation/Cluster.h"
-//#include "blas.h"
-#include "Cluster.h"
-#include "Separator.h"
-#include "AdjMat.h"
-
-namespace MathLib {
-
-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) {
-
- idx_t n_rows(static_cast<idx_t>(_l_adj_mat->getNRows()));
-
- idx_t *xadj(new idx_t[n_rows+1]);
- unsigned const*const original_row_ptr(_l_adj_mat->getRowPtrArray());
- for(idx_t k(0); k<=n_rows; k++) {
- xadj[k] = original_row_ptr[k];
- }
-
- unsigned nnz(_l_adj_mat->getNNZ());
- idx_t *adjncy(new idx_t[nnz]);
- unsigned const*const original_adjncy(_l_adj_mat->getColIdxArray());
- for(unsigned k(0); k<nnz; k++) {
- adjncy[k] = original_adjncy[k];
- }
-// unsigned nparts = 2;
- idx_t options[METIS_NOPTIONS]; // for METIS
- METIS_SetDefaultOptions(options);
-// options[METIS OPTION PTYPE] = METIS PTYPE RB;
-// options[METIS OPTION OBJTYPE] = METIS OBJTYPE CUT;
-// options[METIS OPTION CTYPE] = METIS CTYPE SHEM;
-// options[] = ;
-// options[] = ;
-// options[] = ;
-
-// unsigned sepsize(0); // for METIS
- idx_t *vwgt(new idx_t[n_rows + 1]);
-// const unsigned nnz(xadj[n_rows]);
-// unsigned *adjwgt(new unsigned[nnz]);
- for (idx_t 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_ComputeVertexSeparator(&n_rows, xadj, adjncy, vwgt, &options,
-// &sepsize, part);
- idx_t *loc_op_perm(new idx_t[n_rows]);
- idx_t *loc_po_perm(new idx_t[n_rows]);
- for (idx_t k(0); k<n_rows; k++) {
- loc_op_perm[k] = _g_op_perm[k];
- }
- for (idx_t k(0); k<n_rows; k++) {
- loc_po_perm[k] = _g_po_perm[k];
- }
- METIS_NodeND(&n_rows, xadj, adjncy, vwgt, options, loc_op_perm, loc_po_perm);
- for (idx_t k(0); k<n_rows; k++) {
- _g_op_perm[k] = loc_op_perm[k];
- }
- for (idx_t k(0); k<n_rows; k++) {
- _g_po_perm[k] = loc_po_perm[k];
- }
- delete [] loc_op_perm;
- delete [] loc_po_perm;
- delete [] vwgt;
- delete [] adjncy;
- delete [] xadj;
-// // 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
- std::swap(l_op_perm[beg], l_op_perm[end]);
- std::swap(l_po_perm[l_op_perm[beg]], l_po_perm[l_op_perm[end]]);
- std::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
- std::swap(l_op_perm[beg], l_op_perm[end]);
- std::swap(l_po_perm[l_op_perm[beg]], l_po_perm[l_op_perm[end]]);
- std::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
deleted file mode 100644
index 9a2be5da64b7fe0293363c059c8ed15dfaf080e5..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.h
+++ /dev/null
@@ -1,85 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-02
- * \brief Definition of the Cluster class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#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
- */
- 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
deleted file mode 100644
index c671588c389d76a53f03fd30e9190a5bae5948cc..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/ClusterBase.cpp
+++ /dev/null
@@ -1,72 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-02
- * \brief Implementation of the ClusterBase class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-//#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), _g_adj_mat(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();
-
- // make a copy of the local row_ptr array
- unsigned const* l_row_ptr(_l_adj_mat->getRowPtrArray());
- unsigned *g_row_ptr(new unsigned[n + 1]);
- for (unsigned k = 0; k <= n; ++k)
- g_row_ptr[k] = l_row_ptr[k];
- // make a copy of the local col_idx array
- unsigned const* l_col_idx(_l_adj_mat->getColIdxArray());
- const unsigned g_nnz(g_row_ptr[n]);
- unsigned *g_col_idx(new unsigned[g_nnz]);
- for (unsigned k = 0; k < g_nnz; ++k)
- g_col_idx[k] = l_col_idx[k];
- // generate global matrix from local matrix
- // (only in the root of cluster tree)
- _g_adj_mat = new AdjMat(n, g_row_ptr, g_col_idx);
-
-
-}
-
-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
deleted file mode 100644
index 7ed1c36943fecb00b2d43abd6a1b65feb3f7e139..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/ClusterBase.h
+++ /dev/null
@@ -1,112 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-02
- * \brief Definition of the ClusterBase class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#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
- */
- 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;
-
-#ifndef NDEBUG
- AdjMat const* getGlobalAdjMat() const { return _g_adj_mat; }
-#endif
-
-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
deleted file mode 100644
index 640bdf8a66a8e79124fc877c86495eed6db3f447..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Separator.cpp
+++ /dev/null
@@ -1,52 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-02-01
- * \brief Implementation of the Separator class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-// BaseLib
-
-#include <algorithm>
-#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
- std::swap(l_op_perm [beg], l_op_perm[end]);
- std::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
deleted file mode 100644
index 39957e8582f3f5c7f44e625f3dc27520f92b6929..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/NestedDissectionPermutation/Separator.h
+++ /dev/null
@@ -1,64 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-02
- * \brief Definition of the Separator class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#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
deleted file mode 100644
index 69d3c3a432aafdb898bce0b856dbae10995406c0..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/SparseMatrixBase.h
+++ /dev/null
@@ -1,48 +0,0 @@
-#ifndef SPARSEMATRIXBASE_H
-#define SPARSEMATRIXBASE_H
-
-namespace MathLib {
-
-template<typename FP_TYPE, typename IDX_TYPE> class SparseMatrixBase
-{
-public:
- SparseMatrixBase(IDX_TYPE n1, IDX_TYPE n2) :
- _n_rows(n1), _n_cols(n2)
- {}
- SparseMatrixBase() :
- _n_rows(static_cast<IDX_TYPE>(0)), _n_cols(static_cast<IDX_TYPE>(0))
- {}
- /**
- * y = d * A * x
- * @param d scalar factor
- * @param x vector to multiply with
- * @param y result vector
- */
- virtual void amux(FP_TYPE const d, FP_TYPE const* const __restrict__ x,
- FP_TYPE* __restrict__ y) const = 0;
- virtual ~SparseMatrixBase() {}
- /**
- * get the number of rows
- * @return the number of rows
- */
- IDX_TYPE getNRows () const { return _n_rows; }
- /**
- * get the number of columns
- * @return the number of columns
- */
- IDX_TYPE getNCols () const { return _n_cols; }
-
-protected:
- /**
- * the number of rows
- */
- IDX_TYPE _n_rows;
- /**
- * the number of columns
- */
- IDX_TYPE _n_cols;
-};
-
-} // end namespace MathLib
-
-#endif
diff --git a/MathLib/LinAlg/Sparse/Sparsity.h b/MathLib/LinAlg/Sparse/Sparsity.h
deleted file mode 100644
index f1273f2c53257099e28ebc4bb154184336e0abb0..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/Sparsity.h
+++ /dev/null
@@ -1,31 +0,0 @@
-/**
- * \file
- * \author Norihiro Watanabe
- * \date 2012-06-25
- * \brief Definition of the Sparsity class.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef SPARSITY_H_
-#define SPARSITY_H_
-
-#include <vector>
-#include <set>
-
-namespace MathLib
-{
-
-/**
- * \brief Row-major sparse pattern
- */
-typedef std::vector<std::set<std::size_t> > RowMajorSparsity;
-
-} //MathLib
-
-#endif // SPARSITY_H_
diff --git a/MathLib/LinAlg/Sparse/amuxCRS.cpp b/MathLib/LinAlg/Sparse/amuxCRS.cpp
deleted file mode 100644
index 702f9e6fc735d851fa7983d8122100b5eef940b4..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/amuxCRS.cpp
+++ /dev/null
@@ -1,177 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-09-20
- * \brief Implementation of amuxCRS functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "amuxCRS.h"
-#include <cstddef>
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
-#ifdef HAVE_PTHREADS
-#include <pthread.h>
-#endif
-
-namespace MathLib {
-
-struct MatMultThreadParam {
- MatMultThreadParam (double scalar_factor, unsigned beg_row, unsigned end_row,
- unsigned const * const iA, unsigned const * const jA,
- double const * const A, double const * const x, double* y) :
- _a (scalar_factor), _beg_row(beg_row), _end_row(end_row),
- _row_ptr(iA), _col_idx(jA), _data(A), _x(x), _y(y)
- {}
-
- double _a;
- unsigned _beg_row;
- unsigned _end_row;
- unsigned const * const _row_ptr;
- unsigned const * const _col_idx;
- double const * const _data;
- double const * const _x;
- double * _y;
-};
-
-extern "C" {
-void* amuxCRSpthread (void* ptr)
-{
- MatMultThreadParam *thread_param(static_cast<MatMultThreadParam*>(ptr));
- const double a(thread_param->_a);
- const unsigned beg_row(thread_param->_beg_row);
- const unsigned end_row(thread_param->_end_row);
- unsigned const * const iA(thread_param->_row_ptr);
- unsigned const * const jA(thread_param->_col_idx);
- double const * const A(thread_param->_data);
- double const * const x(thread_param->_x);
- double* y(thread_param->_y);
-
- for (unsigned i(beg_row); i<end_row; i++) {
- y[i] = A[iA[i]] * x[jA[iA[i]]];
- const unsigned end (iA[i+1]);
- for (unsigned j(iA[i]+1); j<end; j++) {
- y[i] += A[j] * x[jA[j]];
- }
- y[i] *= a;
- }
- return NULL;
-}
-} // end extern "C"
-
-void amuxCRSParallelPThreads (double a,
- unsigned n, unsigned const * const iA, unsigned const * const jA,
- double const * const A, double const * const x, double* y,
- unsigned num_of_pthreads)
-{
-#ifdef HAVE_PTHREADS
- // fill thread data objects
- MatMultThreadParam** thread_param_array (new MatMultThreadParam*[num_of_pthreads]);
- double step_size (static_cast<double>(n)/(num_of_pthreads));
- for (unsigned k(0); k<num_of_pthreads; k++) {
- const unsigned beg (static_cast<unsigned>(k*step_size));
- const unsigned end (static_cast<unsigned>((k+1)*step_size));
- thread_param_array[k] = new MatMultThreadParam (a, beg, end, iA, jA, A, x, y);
- }
-
- // allocate thread_array and return value array
- pthread_t *thread_array (new pthread_t[num_of_pthreads]);
- int *ret_vals (new int[num_of_pthreads]);
-
- // create threads
- for (unsigned k(0); k<num_of_pthreads; k++) {
- ret_vals[k] = pthread_create( &(thread_array[k]), NULL, amuxCRSpthread, thread_param_array[k]);
- }
-
- // join threads
- for (unsigned k(0); k<num_of_pthreads; k++) {
- pthread_join (thread_array[k], NULL);
- }
-
- delete [] ret_vals;
- for (unsigned k(0); k<num_of_pthreads; k++)
- delete thread_param_array[k];
- delete [] thread_param_array;
- delete [] thread_array;
-#else
- (void)num_of_pthreads;
- amuxCRS (a, n, iA, jA, A, x, y);
-#endif
-}
-
-void amuxCRSParallelPThreads (double a,
- unsigned n, unsigned const * const iA, unsigned const * const jA,
- double const * const A, double const * const x, double* y,
- unsigned num_of_pthreads, unsigned const*const workload_intervals)
-{
- (void) n; // Unused if HAVE_PTHREADS is not defined.
-
-#ifdef HAVE_PTHREADS
- // fill thread data objects
- MatMultThreadParam** thread_param_array (new MatMultThreadParam*[num_of_pthreads]);
- for (unsigned k(0); k<num_of_pthreads; k++) {
- thread_param_array[k] = new MatMultThreadParam (a, workload_intervals[k], workload_intervals[k+1], iA, jA, A, x, y);
- }
-
- // allocate thread_array and return value array
- pthread_t *thread_array (new pthread_t[num_of_pthreads]);
- int *ret_vals (new int[num_of_pthreads]);
-
- // create threads
- for (unsigned k(0); k<num_of_pthreads; k++) {
- ret_vals[k] = pthread_create( &(thread_array[k]), NULL, amuxCRSpthread, thread_param_array[k]);
- }
-
- // join threads
- for (unsigned k(0); k<num_of_pthreads; k++) {
- pthread_join (thread_array[k], NULL);
- }
-
- delete [] ret_vals;
- for (unsigned k(0); k<num_of_pthreads; k++)
- delete thread_param_array[k];
- delete [] thread_param_array;
- delete [] thread_array;
-#else
- (void)num_of_pthreads;
- amuxCRS (a, n, iA, jA, A, x, y);
-#endif
-}
-
-
-void amuxCRSSym (double a,
- unsigned n, unsigned const * const iA, unsigned const * const jA,
- double const * const A, double const * const x, double* y)
-{
- for (unsigned i(0); i<n; i++) {
- y[i] = 0.0;
- }
-
- for (unsigned i(0); i<n; i++) {
- unsigned j (iA[i]);
- // handle diagonal
- if (jA[j] == i) {
- y[i] += A[j] * x[jA[j]];
- j++;
- }
- const unsigned end (iA[i+1]);
- for (; j<end; j++) {
- y[i] += A[j] * x[jA[j]];
- y[jA[j]] += A[j] * x[i];
- }
- }
-
- for (unsigned i(0); i<n; i++) {
- y[i] *= a;
- }
-}
-
-} // end namespace MathLib
diff --git a/MathLib/LinAlg/Sparse/amuxCRS.h b/MathLib/LinAlg/Sparse/amuxCRS.h
deleted file mode 100644
index 94c26d308e9687c30656216b8fc6ee3caf98dd6d..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/amuxCRS.h
+++ /dev/null
@@ -1,81 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-09-20
- * \brief Definition of amuxCRS functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef AMUXCRS_H
-#define AMUXCRS_H
-
-namespace MathLib {
-
-template<typename FP_TYPE, typename IDX_TYPE>
-void amuxCRS(FP_TYPE a, IDX_TYPE n, IDX_TYPE const * const iA, IDX_TYPE const * const jA,
- FP_TYPE const * const A, FP_TYPE const * const x, FP_TYPE* y)
-{
- for (IDX_TYPE i(0); i < n; i++) {
- const IDX_TYPE end(iA[i + 1]);
- y[i] = A[iA[i]] * x[jA[iA[i]]];
- for (IDX_TYPE j(iA[i]+1); j < end; j++) {
- y[i] += A[j] * x[jA[j]];
- }
- y[i] *= a;
- }
-}
-
-void amuxCRSParallelPThreads (double a,
- unsigned n, unsigned const * const iA, unsigned const * const jA,
- double const * const A, double const * const x, double* y,
- unsigned num_of_pthreads);
-
-void amuxCRSParallelPThreads (double a,
- unsigned n, unsigned const * const iA, unsigned const * const jA,
- double const * const A, double const * const x, double* y,
- unsigned num_of_pthreads, unsigned const*const workload_intervals);
-
-#ifdef _OPENMP
-template<typename FP_TYPE, typename IDX_TYPE>
-void amuxCRSParallelOpenMP (FP_TYPE a, unsigned n,
- IDX_TYPE const * const __restrict__ iA,
- IDX_TYPE const * const __restrict__ jA, FP_TYPE const * const A,
- FP_TYPE const * const __restrict__ x, FP_TYPE* __restrict__ y)
-{
- OPENMP_LOOP_TYPE i;
- IDX_TYPE j;
- FP_TYPE t;
- {
-#pragma omp parallel for private(i, j, t)
-#ifdef WIN32
-#pragma warning ( push )
-#pragma warning ( disable: 4018 )
-#endif
- for (i = 0; i < n; i++) {
- const IDX_TYPE end(iA[i + 1]);
- t = A[iA[i]] * x[jA[iA[i]]];
- for (j = iA[i]+1; j < end; j++) {
- t += A[j] * x[jA[j]];
- }
- y[i] = t * a;
- }
- }
-#ifdef WIN32
-#pragma warning ( pop )
-#endif
-}
-#endif
-
-void amuxCRSSym (double a,
- unsigned n, unsigned const * const iA, unsigned const * const jA,
- double const * const A, double const * const x, double* y);
-
-} // end namespace MathLib
-
-#endif
diff --git a/MathLib/LinAlg/Sparse/sparse.h b/MathLib/LinAlg/Sparse/sparse.h
deleted file mode 100644
index 839890eca812be9ed2fdf8a4e27b3cefb9d9b048..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Sparse/sparse.h
+++ /dev/null
@@ -1,65 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date no date
- * \brief Definition of sparse matrix IO functions.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef SPARSE_H
-#define SPARSE_H
-
-#include <iostream>
-#include <cassert>
-
-//extern void CS_write(char*, unsigned, unsigned const*, unsigned const*, double const*);
-//extern void CS_read(char*, unsigned&, unsigned*&, unsigned*&, double*&);
-
-template<class T> void CS_write(std::ostream &os, unsigned n, unsigned const* iA, unsigned const* jA, T const* A)
-{
- os.write(reinterpret_cast<char*>(&n), sizeof(unsigned));
- os.write(reinterpret_cast<char*>(const_cast<unsigned*>(iA)), (n + 1) * sizeof(unsigned));
- os.write(reinterpret_cast<char*>(const_cast<unsigned*>(jA)), iA[n] * sizeof(unsigned));
- os.write(reinterpret_cast<char const*>(A), iA[n] * sizeof(T));
-}
-
-template<class T> void CS_read(std::istream &is, unsigned &n, unsigned* &iA, unsigned* &jA, T* &A)
-{
- is.read(reinterpret_cast<char*>(&n), sizeof(unsigned));
- if (iA != NULL) {
- delete[] iA;
- delete[] jA;
- delete[] A;
- }
- iA = new unsigned[n + 1];
- assert(iA != NULL);
- is.read(reinterpret_cast<char*>(iA), (n + 1) * sizeof(unsigned));
-
- jA = new unsigned[iA[n]];
- assert(jA != NULL);
- is.read(reinterpret_cast<char*>(jA), iA[n] * sizeof(unsigned));
-
- A = new T[iA[n]];
- assert(A != NULL);
- is.read(reinterpret_cast<char*>(A), iA[n] * sizeof(T));
-
-#ifndef NDEBUG
- // do simple checks
- if (iA[0] != 0) std::cerr << "\nCRS matrix: array iA doesn't start with 0\n";
-
- unsigned i = 0;
- while (i < iA[n] && jA[i] < n)
- ++i;
- if (i < iA[n]) std::cerr << "\nCRS matrix: the " << i
- << "th entry of jA has the value " << jA[i] << ", which is out of bounds.\n";
-#endif
-}
-
-#endif
-
diff --git a/ProcessLib/NumericsConfig.h b/ProcessLib/NumericsConfig.h
index 696f51484bddfdcc961711b48712274d2277d163..dcdba10ec651473f5082158303a86f85461b699a 100644
--- a/ProcessLib/NumericsConfig.h
+++ b/ProcessLib/NumericsConfig.h
@@ -37,20 +37,6 @@
using LinearSolverType = MathLib::EigenLisLinearSolver;
}
-#elif defined(USE_LIS)
-
- #include "MathLib/LinAlg/Lis/LisMatrix.h"
- #include "MathLib/LinAlg/Lis/LisVector.h"
- #include "MathLib/LinAlg/Lis/LisLinearSolver.h"
-
-namespace detail
-{
- using GlobalVectorType = MathLib::LisVector;
- using GlobalMatrixType = MathLib::LisMatrix;
-
- using LinearSolverType = MathLib::LisLinearSolver;
-}
-
#elif defined(USE_PETSC)
#include "MathLib/LinAlg/PETSc/PETScVector.h"
#include "MathLib/LinAlg/PETSc/PETScMatrix.h"
diff --git a/SimpleTests/MatrixTests/CMakeLists.txt b/SimpleTests/MatrixTests/CMakeLists.txt
index 445aa4d05c7ed75887ec4857e03048966adac296..f5db03f53f0d95ecfcd3fd16312c192f92720a89 100644
--- a/SimpleTests/MatrixTests/CMakeLists.txt
+++ b/SimpleTests/MatrixTests/CMakeLists.txt
@@ -1,81 +1,4 @@
-if(CMAKE_USE_PTHREADS_INIT)
- set(HAVE_PTHREADS TRUE)
-endif()
-
# Create the executable
-add_executable(MatMult
- MatMult.cpp
- ${SOURCES}
- ${HEADERS}
-)
-set_target_properties(MatMult PROPERTIES FOLDER SimpleTests)
-target_link_libraries(MatMult logog)
-
-if(HAVE_PTHREADS)
- add_executable(MatVecMultPthreads
- MatVecMultPthreads.cpp
- ${SOURCES}
- ${HEADERS}
- )
- set_target_properties(MatVecMultPthreads PROPERTIES FOLDER SimpleTests)
- target_link_libraries(MatVecMultPthreads
- pthread
- BaseLib
- MathLib
- logog
- )
-endif()
-
-set_target_properties(MatMult PROPERTIES FOLDER SimpleTests)
-target_link_libraries(MatMult
- logog
- BaseLib
- MathLib
-)
-
-add_executable(MatTestRemoveRowsCols
- MatTestRemoveRowsCols.cpp
- ${SOURCES}
- ${HEADERS}
-)
-set_target_properties(MatTestRemoveRowsCols PROPERTIES FOLDER SimpleTests)
-target_link_libraries(MatTestRemoveRowsCols
- BaseLib
- MathLib
-)
-
-if(METIS_FOUND)
- add_executable(MatVecMultNDPerm
- MatVecMultNDPerm.cpp
- ${SOURCES}
- ${HEADERS}
- )
- set_target_properties(MatVecMultNDPerm PROPERTIES FOLDER SimpleTests)
-
- target_link_libraries(MatVecMultNDPerm
- BaseLib
- MathLib
- logog
- ${METIS_LIBRARIES}
- )
-
- if(OPENMP_FOUND)
- add_executable(MatVecMultNDPermOpenMP
- MatVecMultNDPermOpenMP.cpp
- ${SOURCES}
- ${HEADERS}
- )
- set_target_properties(MatVecMultNDPermOpenMP PROPERTIES FOLDER SimpleTests)
-
- target_link_libraries(MatVecMultNDPermOpenMP
- BaseLib
- MathLib
- logog
- ${METIS_LIBRARIES}
- )
- endif()
-endif()
-
add_executable(DenseGaussEliminationChecker
DenseGaussEliminationChecker.cpp
${SOURCES}
@@ -87,4 +10,3 @@ target_link_libraries(DenseGaussEliminationChecker
BaseLib
MathLib
)
-
diff --git a/SimpleTests/MatrixTests/MatMult.cpp b/SimpleTests/MatrixTests/MatMult.cpp
deleted file mode 100644
index 26dfd4780094e730edf14435d25358071c556eb1..0000000000000000000000000000000000000000
--- a/SimpleTests/MatrixTests/MatMult.cpp
+++ /dev/null
@@ -1,192 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-03
- * \brief Implementation of tests.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <fstream>
-#include <iostream>
-#include <cmath>
-#include <limits>
-#include <cstdlib>
-
-#ifdef UNIX
-#include <sys/unistd.h>
-#endif
-
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
-#include <logog/include/logog.hpp>
-#include <logog/include/formatter.hpp>
-#include <tclap/CmdLine.h>
-
-#include "BaseLib/BuildInfo.h"
-#include "BaseLib/CPUTime.h"
-#include "BaseLib/RunTime.h"
-
-#include "MathLib/LinAlg/Sparse/CRSMatrix.h"
-#include "MathLib/LinAlg/Sparse/CRSMatrixOpenMP.h"
-#include "MathLib/LinAlg/Sparse/CRSMatrixPThreads.h"
-
-/**
- * new formatter for logog
- */
-class FormatterCustom : public logog::FormatterGCC
-{
- virtual TOPIC_FLAGS GetTopicFlags( const logog::Topic &topic )
- {
- return ( Formatter::GetTopicFlags( topic ) &
- ~( TOPIC_FILE_NAME_FLAG | TOPIC_LINE_NUMBER_FLAG ));
- }
-};
-
-int main(int argc, char *argv[])
-{
- LOGOG_INITIALIZE();
-
- TCLAP::CmdLine cmd("Simple matrix vector multiplication test", ' ', "0.1");
-
- // Define a value argument and add it to the command line.
- // A value arg defines a flag and a type of value that it expects,
- // such as "-m matrix".
- TCLAP::ValueArg<std::string> matrix_arg("m", "matrix", "input matrix file", true, "", "string");
-
- // Add the argument mesh_arg to the CmdLine object. The CmdLine object
- // uses this Arg to parse the command line.
- cmd.add( matrix_arg );
-
- TCLAP::ValueArg<unsigned> n_cores_arg("p", "number-cores", "number of cores to use", false, 1, "number");
- cmd.add( n_cores_arg );
-
- TCLAP::ValueArg<unsigned> n_mults_arg("n", "number-of-multiplications", "number of multiplications to perform", true, 10, "number");
- cmd.add( n_mults_arg );
-
- TCLAP::ValueArg<std::string> output_arg("o", "output", "output file", false, "", "string");
- cmd.add( output_arg );
-
- TCLAP::ValueArg<unsigned> verbosity_arg("v", "verbose", "level of verbosity [0 very low information, 1 much information]", false, 0, "string");
- cmd.add( verbosity_arg );
-
- cmd.parse( argc, argv );
-
- // read the number of multiplication to execute
- unsigned n_mults (n_mults_arg.getValue());
- std::string fname_mat (matrix_arg.getValue());
-
- FormatterCustom *custom_format (new FormatterCustom);
- logog::Cout *logogCout(new logog::Cout);
- logogCout->SetFormatter(*custom_format);
-
- logog::LogFile *logog_file(NULL);
- if (! output_arg.getValue().empty()) {
- logog_file = new logog::LogFile(output_arg.getValue().c_str());
- logog_file->SetFormatter( *custom_format );
- }
-
- // read number of threads
- unsigned n_threads (n_cores_arg.getValue());
-
- INFO("%s was build with compiler %s",
- argv[0],
- BaseLib::BuildInfo::cmake_cxx_compiler.c_str());
-#ifdef NDEBUG
- INFO("CXX_FLAGS: %s %s",
- BaseLib::BuildInfo::cmake_cxx_flags.c_str(),
- BaseLib::BuildInfo::cmake_cxx_flags_release.c_str());
-#else
- INFO("CXX_FLAGS: %s %s",
- BaseLib::BuildInfo::cmake_cxx_flags.c_str(),
- BaseLib::BuildInfo::cmake_cxx_flags_debug.c_str());
-#endif
-
-#ifdef UNIX
- const int max_host_name_len (255);
- char *hostname(new char[max_host_name_len]);
- if (gethostname(hostname, max_host_name_len) == 0)
- INFO("hostname: %s", hostname);
- delete [] host_name_len;
-#endif
-
- // *** reading matrix in crs format from file
- std::ifstream in(fname_mat.c_str(), std::ios::in | std::ios::binary);
- double *A(NULL);
- unsigned *iA(NULL), *jA(NULL), n;
- if (in) {
- INFO("reading matrix from %s ...", fname_mat.c_str());
- BaseLib::RunTime timer;
- timer.start();
- CS_read(in, n, iA, jA, A);
- INFO("\t- took %e s", timer.elapsed());
- } else {
- INFO("error reading matrix from %s", fname_mat.c_str());
- return -1;
- }
- unsigned nnz(iA[n]);
- INFO("\tParameters read: n=%d, nnz=%d", n, nnz);
-
-#ifdef _OPENMP
- omp_set_num_threads(n_threads);
- unsigned *mat_entries_per_core(new unsigned[n_threads]);
- for (unsigned k(0); k<n_threads; k++) {
- mat_entries_per_core[k] = 0;
- }
-
- OPENMP_LOOP_TYPE i;
- {
-#pragma omp parallel for
- for (i = 0; i < n; i++) {
- mat_entries_per_core[omp_get_thread_num()] += iA[i + 1] - iA[i];
- }
- }
-
- INFO("*** work per core ***");
- for (unsigned k(0); k<n_threads; k++) {
- INFO("\t%d\t%d", k, mat_entries_per_core[k]);
- }
-#endif
-
-#ifdef _OPENMP
- omp_set_num_threads(n_threads);
- MathLib::CRSMatrixOpenMP<double, unsigned> mat (n, iA, jA, A);
-#else
- MathLib::CRSMatrix<double, unsigned> mat (n, iA, jA, A);
-#endif
-
- double *x(new double[n]);
- double *y(new double[n]);
-
- for (unsigned k(0); k<n; ++k)
- x[k] = 1.0;
-
- INFO("*** %d matrix vector multiplications (MVM) with Toms amuxCRS (%d threads) ...", n_mults, n_threads);
- BaseLib::RunTime run_timer;
- BaseLib::CPUTime cpu_timer;
- run_timer.start();
- cpu_timer.start();
- for (std::size_t k(0); k<n_mults; k++) {
- mat.amux (1.0, x, y);
- }
-
- INFO("\t[MVM] - took %e sec cpu time, %e sec run time", cpu_timer.elapsed(), run_timer.elapsed());
-
- delete [] x;
- delete [] y;
-
- delete custom_format;
- delete logogCout;
- delete logog_file;
- LOGOG_SHUTDOWN();
-
- return 0;
-}
-
diff --git a/SimpleTests/MatrixTests/MatTestRemoveRowsCols.cpp b/SimpleTests/MatrixTests/MatTestRemoveRowsCols.cpp
deleted file mode 100644
index 6bdcb27d5fa18b362502bf2d89f6c83ae9f5fd24..0000000000000000000000000000000000000000
--- a/SimpleTests/MatrixTests/MatTestRemoveRowsCols.cpp
+++ /dev/null
@@ -1,79 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-11-08
- * \brief Test for removing entries in matrices.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- */
-
-#include <fstream>
-#include <iostream>
-
-#include "BaseLib/RunTime.h"
-#include "BaseLib/CPUTime.h"
-
-#include "MathLib/LinAlg/Sparse/CRSMatrix.h"
-
-int main(int argc, char *argv[])
-{
- if (argc < 3) {
- std::cout << "Usage: " << argv[0] << " input-matrix output-matrix" << std::endl;
- return 1;
- }
-
- std::string fname_mat (argv[1]);
- bool verbose (true);
-
- // *** reading matrix in crs format from file
- std::ifstream in(fname_mat.c_str(), std::ios::in | std::ios::binary);
- double *A(NULL);
- unsigned *iA(NULL), *jA(NULL), n = 0;
- if (in) {
- if (verbose) {
- std::cout << "reading matrix from " << fname_mat << " ... " << std::flush;
- }
- BaseLib::RunTime timer;
- timer.start();
- CS_read(in, n, iA, jA, A);
- in.close();
- if (verbose) {
- std::cout << "ok, " << timer.elapsed() << " s" << std::endl;
- }
- } else {
- std::cout << "error reading matrix from " << fname_mat << std::endl;
- }
- unsigned nnz(iA[n]);
- if (verbose) {
- std::cout << "Parameters read: n=" << n << ", nnz=" << nnz << std::endl;
- }
-
- MathLib::CRSMatrix<double, unsigned> *mat (new MathLib::CRSMatrix<double, unsigned>(n, iA, jA, A));
-
- const unsigned n_rows_cols_to_erase(300);
- unsigned *rows_cols_to_erase(new unsigned[n_rows_cols_to_erase]);
-
- for (unsigned k(0); k<n_rows_cols_to_erase; k++) {
- rows_cols_to_erase[k] = (k+1)*2;
- }
-
- BaseLib::RunTime timer;
- std::cout << "erasing " << n_rows_cols_to_erase << " rows and columns ... " << std::flush;
- timer.start();
- mat->eraseEntries(n_rows_cols_to_erase, rows_cols_to_erase);
- std::cout << "ok, " << timer.elapsed() << " s" << std::endl;
- delete[] rows_cols_to_erase;
-
- fname_mat = argv[2];
- std::ofstream out (fname_mat.c_str(), std::ios::binary);
- CS_write (out, mat->getNRows(), mat->getRowPtrArray(), mat->getColIdxArray(), mat->getEntryArray());
- out.close();
-
- std::cout << "wrote " << fname_mat << " with " << mat->getNRows() << " rows and " << mat->getRowPtrArray()[mat->getNRows()] << " entries" << std::endl;
-
- delete mat;
-}
diff --git a/SimpleTests/MatrixTests/MatVecMultNDPerm.cpp b/SimpleTests/MatrixTests/MatVecMultNDPerm.cpp
deleted file mode 100644
index e64af5a77e8e6160ef82b5af803bddb66aa09247..0000000000000000000000000000000000000000
--- a/SimpleTests/MatrixTests/MatVecMultNDPerm.cpp
+++ /dev/null
@@ -1,188 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-03
- * \brief Implementation of tests.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <cstdlib>
-
-#ifdef UNIX
-#include <sys/unistd.h>
-#endif
-
-#include <tclap/CmdLine.h>
-#include <logog/include/logog.hpp>
-#include <logog/include/formatter.hpp>
-
-#include "BaseLib/BuildInfo.h"
-#include "BaseLib/RunTime.h"
-#include "BaseLib/CPUTime.h"
-#include "BaseLib/LogogSimpleFormatter.h"
-
-#include "MathLib/LinAlg/Sparse/sparse.h"
-#include "MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.h"
-#include "MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReordered.h"
-#include "MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.h"
-
-int main(int argc, char *argv[])
-{
- LOGOG_INITIALIZE();
-
- TCLAP::CmdLine cmd("The purpose of this program is the speed test of sparse matrix vector multiplication (MVM), where the matrix is stored in CRS format. Before executing the MVM a nested dissection reordering is performed.", ' ', "0.1");
-
- // Define a value argument and add it to the command line.
- // A value arg defines a flag and a type of value that it expects,
- // such as "-m matrix".
- TCLAP::ValueArg<std::string> matrix_arg("m","matrix","input matrix file in CRS format",true,"","file name of the matrix in CRS format");
-
- // Add the argument matrix_arg to the CmdLine object. The CmdLine object
- // uses this Arg to parse the command line.
- cmd.add( matrix_arg );
-
-// TCLAP::ValueArg<unsigned> n_cores_arg("n", "number-cores", "number of cores to use", true, "1", "number");
-// cmd.add( n_cores_arg );
-
- TCLAP::ValueArg<unsigned> n_mults_arg("n", "number-of-multiplications", "number of multiplications to perform", true, 10, "number of multiplications");
- cmd.add( n_mults_arg );
-
- TCLAP::ValueArg<std::string> output_arg("o", "output", "output file", false, "", "string");
- cmd.add( output_arg );
-
- TCLAP::ValueArg<bool> verbosity_arg("v", "verbose", "level of verbosity [0 very low information, 1 much information]", false, 0, "string");
- cmd.add( verbosity_arg );
-
- cmd.parse( argc, argv );
-
- // read the number of multiplication to execute
- unsigned n_mults (n_mults_arg.getValue());
- std::string fname_mat (matrix_arg.getValue());
- bool verbose (verbosity_arg.getValue());
-
- BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
- logog::Cout *logogCout(new logog::Cout);
- logogCout->SetFormatter(*custom_format);
-
- INFO("%s was build with compiler %s",
- argv[0],
- BaseLib::BuildInfo::cmake_cxx_compiler.c_str());
-#ifdef NDEBUG
- INFO("CXX_FLAGS: %s %s",
- BaseLib::BuildInfo::cmake_cxx_flags.c_str(),
- BaseLib::BuildInfo::cmake_cxx_flags_release.c_str());
-#else
- INFO("CXX_FLAGS: %s %s",
- BaseLib::BuildInfo::cmake_cxx_flags.c_str(),
- BaseLib::BuildInfo::cmake_cxx_flags_debug.c_str());
-#endif
-
-#ifdef UNIX
- const std::size_t length(256);
- char *hostname(new char[length]);
- gethostname (hostname, length);
- INFO("hostname: %s", hostname);
- delete [] hostname;
-#endif
-
- // *** reading matrix in crs format from file
- std::ifstream in(fname_mat.c_str(), std::ios::in | std::ios::binary);
- double *A(NULL);
- unsigned *iA(NULL), *jA(NULL), n;
- if (in) {
- if (verbose) {
- INFO("reading matrix from %s ...", fname_mat.c_str());
- }
- BaseLib::RunTime timer;
- timer.start();
- CS_read(in, n, iA, jA, A);
- if (verbose) {
- INFO("\t- took %e s", timer.elapsed());
- }
- } else {
- ERR("error reading matrix from %s", fname_mat.c_str());
- return -1;
- }
- unsigned nnz(iA[n]);
- if (verbose) {
- INFO("\tParameters read: n=%d, nnz=%d", n, nnz);
- }
-
- MathLib::CRSMatrixReordered mat(n, iA, jA, A);
-
- double *x(new double[n]);
- double *y(new double[n]);
-
- for (unsigned k(0); k<n; ++k)
- x[k] = 1.0;
-
- // create time measurement objects
- BaseLib::RunTime run_timer;
- BaseLib::CPUTime cpu_timer;
-
- // calculate the nested dissection reordering
- if (verbose) {
- INFO("*** calculating nested dissection (ND) permutation of matrix ...");
- }
- run_timer.start();
- cpu_timer.start();
- MathLib::Cluster cluster_tree(n, iA, jA);
- unsigned *op_perm(new unsigned[n]);
- unsigned *po_perm(new unsigned[n]);
- for (unsigned k(0); k<n; k++)
- op_perm[k] = po_perm[k] = k;
- cluster_tree.createClusterTree(op_perm, po_perm, 1000);
- if (verbose) {
- INFO("\t[ND] - took %e sec \t%e sec", cpu_timer.elapsed(), run_timer.elapsed());
- }
-
- // applying the nested dissection reordering
- if (verbose) {
- INFO("\t[ND] applying nested dissection permutation to FEM matrix ... ");
- }
- run_timer.start();
- cpu_timer.start();
- mat.reorderMatrix(op_perm, po_perm);
- if (verbose) {
- INFO("\t[ND]: - took %e sec\t%e sec", cpu_timer.elapsed(), run_timer.elapsed());
- }
-
-#ifndef NDEBUG
-// std::string fname_mat_out(fname_mat.substr(0,fname_mat.length()-4)+"-reordered.bin");
-// std::ofstream os (fname_mat_out.c_str(), std::ios::binary);
-// if (os) {
-// std::cout << "writing matrix to " << fname_mat_out << " ... " << std::flush;
-// CS_write(os, n, mat.getRowPtrArray(), mat.getColIdxArray(), mat.getEntryArray());
-// std::cout << "done" << std::endl;
-// }
-#endif
-
- if (verbose) {
- INFO("*** %d matrix vector multiplications (MVM) with Toms amuxCRS ... ", n_mults);
- }
- run_timer.start();
- cpu_timer.start();
- for (std::size_t k(0); k<n_mults; k++) {
- mat.amux (1.0, x, y);
- }
-
- if (verbose) {
- INFO("\t[MVM] - took %e sec\t %e sec", cpu_timer.elapsed(), run_timer.elapsed());
- }
-
- delete [] x;
- delete [] y;
-
- delete custom_format;
- delete logogCout;
- LOGOG_SHUTDOWN();
-
-
- return 0;
-}
diff --git a/SimpleTests/MatrixTests/MatVecMultNDPermOpenMP.cpp b/SimpleTests/MatrixTests/MatVecMultNDPermOpenMP.cpp
deleted file mode 100644
index a9d8fdf312acedb744221ffa5d321cef1ce26190..0000000000000000000000000000000000000000
--- a/SimpleTests/MatrixTests/MatVecMultNDPermOpenMP.cpp
+++ /dev/null
@@ -1,197 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2012-01-20
- * \brief Implementation of tests.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <cstdlib>
-
-#ifdef UNIX
-#include <sys/unistd.h>
-#endif
-
-
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
-#include <tclap/CmdLine.h>
-#include <logog/include/logog.hpp>
-#include <logog/include/formatter.hpp>
-
-#include "BaseLib/CPUTime.h"
-#include "BaseLib/RunTime.h"
-#include "BaseLib/LogogSimpleFormatter.h"
-#include "BaseLib/BuildInfo.h"
-
-#include "MathLib/LinAlg/Sparse/NestedDissectionPermutation/AdjMat.h"
-#include "MathLib/LinAlg/Sparse/NestedDissectionPermutation/CRSMatrixReorderedOpenMP.h"
-#include "MathLib/LinAlg/Sparse/NestedDissectionPermutation/Cluster.h"
-
-int main(int argc, char *argv[])
-{
-#ifndef _OPENMP
- static_assert(false, "This code must be compiled with _OPENMP macro enabled.");
-#endif
-
- LOGOG_INITIALIZE();
-
- TCLAP::CmdLine cmd("The purpose of this program is the speed test of sparse matrix vector multiplication (MVM) employing OpenMP technique, where the matrix is stored in CRS format. Before executing the MVM a nested dissection reordering is performed.", ' ', "0.1");
-
- // Define a value argument and add it to the command line.
- // A value arg defines a flag and a type of value that it expects,
- // such as "-m matrix".
- TCLAP::ValueArg<std::string> matrix_arg("m","matrix","input matrix file in CRS format",true,"","file name of the matrix in CRS format");
-
- // Add the argument matrix_arg to the CmdLine object. The CmdLine object
- // uses this Arg to parse the command line.
- cmd.add( matrix_arg );
-
- TCLAP::ValueArg<unsigned> n_cores_arg("p", "number-cores", "number of cores to use", true, 1, "number of cores");
- cmd.add( n_cores_arg );
-
- TCLAP::ValueArg<unsigned> n_mults_arg("n", "number-of-multiplications", "number of multiplications to perform", true, 10, "number of multiplications");
- cmd.add( n_mults_arg );
-
- TCLAP::ValueArg<std::string> output_arg("o", "output", "output file", false, "", "string");
- cmd.add( output_arg );
-
- TCLAP::ValueArg<bool> verbosity_arg("v", "verbose", "level of verbosity [0 very low information, 1 much information]", false, 0, "string");
- cmd.add( verbosity_arg );
-
- cmd.parse( argc, argv );
-
- // read the number of multiplication to execute
- unsigned n_mults (n_mults_arg.getValue());
- std::string fname_mat (matrix_arg.getValue());
- bool verbose (verbosity_arg.getValue());
-
- BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
- logog::Cout *logogCout(new logog::Cout);
- logogCout->SetFormatter(*custom_format);
-
- // read number of threads
- unsigned n_threads (n_cores_arg.getValue());
-
- INFO("%s was build with compiler %s",
- argv[0],
- BaseLib::BuildInfo::cmake_cxx_compiler.c_str());
-#ifdef NDEBUG
- INFO("CXX_FLAGS: %s %s",
- BaseLib::BuildInfo::cmake_cxx_flags.c_str(),
- BaseLib::BuildInfo::cmake_cxx_flags_release.c_str());
-#else
- INFO("CXX_FLAGS: %s %s",
- BaseLib::BuildInfo::cmake_cxx_flags.c_str(),
- BaseLib::BuildInfo::cmake_cxx_flags_debug.c_str());
-#endif
-
-#ifdef UNIX
- const std::size_t length(256);
- char *hostname(new char[length]);
- gethostname (hostname, length);
- INFO("hostname: %s", hostname);
- delete [] hostname;
-#endif
-
- // *** reading matrix in crs format from file
- std::ifstream in(fname_mat.c_str(), std::ios::in | std::ios::binary);
- double *A(NULL);
- unsigned *iA(NULL), *jA(NULL), n;
- if (in) {
- if (verbose) {
- INFO("reading matrix from %s ...", fname_mat.c_str());
- }
- BaseLib::RunTime timer;
- timer.start();
- CS_read(in, n, iA, jA, A);
- if (verbose) {
- INFO("\t- took %e s", timer.elapsed());
- }
- } else {
- ERR("error reading matrix from %s", fname_mat.c_str());
- return -1;
- }
- unsigned nnz(iA[n]);
- if (verbose) {
- INFO("\tParameters read: n=%d, nnz=%d", n, nnz);
- }
-
-#ifdef _OPENMP
- omp_set_num_threads(n_threads);
- MathLib::CRSMatrixReorderedOpenMP mat(n, iA, jA, A);
-#else
- delete [] iA;
- delete [] jA;
- delete [] A;
- ERROR("program is not using OpenMP");
- return -1;
-#endif
- double *x(new double[n]);
- double *y(new double[n]);
-
- for (unsigned k(0); k<n; ++k)
- x[k] = 1.0;
-
- // create time measurement objects
- BaseLib::RunTime run_timer;
- BaseLib::CPUTime cpu_timer;
-
- // calculate the nested dissection reordering
- if (verbose) {
- INFO("*** calculating nested dissection (ND) permutation of matrix ...");
- }
- run_timer.start();
- cpu_timer.start();
- MathLib::Cluster cluster_tree(n, iA, jA);
- unsigned *op_perm(new unsigned[n]);
- unsigned *po_perm(new unsigned[n]);
- for (unsigned k(0); k<n; k++)
- op_perm[k] = po_perm[k] = k;
- cluster_tree.createClusterTree(op_perm, po_perm, 1000);
- if (verbose) {
- INFO("\t[ND] - took %e sec \t%e sec", cpu_timer.elapsed(), run_timer.elapsed());
- }
-
- // applying the nested dissection reordering
- if (verbose) {
- INFO("\t[ND] applying nested dissection permutation to FEM matrix ... ");
- }
- run_timer.start();
- cpu_timer.start();
- mat.reorderMatrix(op_perm, po_perm);
- if (verbose) {
- INFO("\t[ND]: - took %e sec\t%e sec", cpu_timer.elapsed(), run_timer.elapsed());
- }
-
- if (verbose) {
- INFO("*** %d matrix vector multiplications (MVM) with Toms amuxCRS (%d threads)... ", n_mults, n_threads);
- }
-
- run_timer.start();
- cpu_timer.start();
- for (std::size_t k(0); k<n_mults; k++) {
- mat.amux (1.0, x, y);
- }
-
- if (verbose) {
- INFO("\t[MVM] - took %e sec cpu time, %e sec run time", cpu_timer.elapsed(), run_timer.elapsed());
- }
-
- delete [] x;
- delete [] y;
-
- delete custom_format;
- delete logogCout;
- LOGOG_SHUTDOWN();
-
- return 0;
-}
diff --git a/SimpleTests/MatrixTests/MatVecMultPthreads.cpp b/SimpleTests/MatrixTests/MatVecMultPthreads.cpp
deleted file mode 100644
index aa50412972bc1576c8bc09a3b10c627be5b3c03d..0000000000000000000000000000000000000000
--- a/SimpleTests/MatrixTests/MatVecMultPthreads.cpp
+++ /dev/null
@@ -1,152 +0,0 @@
-/**
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.net/LICENSE.txt
- *
- * \file MatVecMultPthreads.cpp
- *
- * Created on Jul 3, 2012 by Thomas Fischer
- */
-
-#include <fstream>
-#include <iostream>
-#include <cmath>
-#include <limits>
-#include <cstdlib>
-
-#ifdef UNIX
-#include <sys/unistd.h>
-#endif
-
-#ifdef HAVE_PTHREADS
-#include <pthread.h>
-#endif
-
-#include <logog/include/logog.hpp>
-#include <logog/include/formatter.hpp>
-#include <tclap/CmdLine.h>
-
-#include "MathLib/LinAlg/Sparse/CRSMatrix.h"
-#include "MathLib/LinAlg/Sparse/CRSMatrixOpenMP.h"
-#include "MathLib/LinAlg/Sparse/CRSMatrixPThreads.h"
-
-#include "BaseLib/BuildInfo.h"
-#include "BaseLib/CPUTime.h"
-#include "BaseLib/RunTime.h"
-#include "BaseLib/LogogSimpleFormatter.h"
-
-int main(int argc, char *argv[])
-{
- LOGOG_INITIALIZE();
-
- TCLAP::CmdLine cmd("Simple matrix vector multiplication test employing pthreads", ' ', "0.1");
-
- // Define a value argument and add it to the command line.
- // A value arg defines a flag and a type of value that it expects,
- // such as "-m matrix".
- TCLAP::ValueArg<std::string> matrix_arg("m", "matrix", "input matrix file", true, "", "string");
-
- // Add the argument mesh_arg to the CmdLine object. The CmdLine object
- // uses this Arg to parse the command line.
- cmd.add( matrix_arg );
-
- TCLAP::ValueArg<unsigned> n_cores_arg("p", "number-cores", "number of cores to use", false, 1, "number");
- cmd.add( n_cores_arg );
-
- TCLAP::ValueArg<unsigned> n_mults_arg("n", "number-of-multiplications", "number of multiplications to perform", true, 10, "number");
- cmd.add( n_mults_arg );
-
- TCLAP::ValueArg<std::string> output_arg("o", "output", "output file", false, "", "string");
- cmd.add( output_arg );
-
- TCLAP::ValueArg<bool> verbosity_arg("v", "verbose", "level of verbosity [0 very low information, 1 much information]", false, 0, "string");
- cmd.add( verbosity_arg );
-
- cmd.parse( argc, argv );
-
- std::string fname_mat (matrix_arg.getValue());
-
- BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
- logog::Cout *logogCout(new logog::Cout);
- logogCout->SetFormatter(*custom_format);
-
- logog::LogFile *logog_file(nullptr);
- if (! output_arg.getValue().empty()) {
- logog_file = new logog::LogFile(output_arg.getValue().c_str());
- logog_file->SetFormatter( *custom_format );
- }
-
-INFO("%s was build with compiler %s",
- argv[0],
- BaseLib::BuildInfo::cmake_cxx_compiler.c_str());
-#ifdef NDEBUG
- INFO("CXX_FLAGS: %s %s",
- BaseLib::BuildInfo::cmake_cxx_flags.c_str(),
- BaseLib::BuildInfo::cmake_cxx_flags_release.c_str());
-#else
- INFO("CXX_FLAGS: %s %s",
- BaseLib::BuildInfo::cmake_cxx_flags.c_str(),
- BaseLib::BuildInfo::cmake_cxx_flags_debug.c_str());
-#endif
-
-#ifdef UNIX
- const int max_host_name_len (255);
- char *hostname(new char[max_host_name_len]);
- if (gethostname(hostname, max_host_name_len) == 0)
- INFO("hostname: %s", hostname);
- delete [] host_name_len;
-#endif
-
- // *** reading matrix in crs format from file
- std::ifstream in(fname_mat.c_str(), std::ios::in | std::ios::binary);
- double *A(NULL);
- unsigned *iA(NULL), *jA(NULL), n;
- if (in) {
- INFO("reading matrix from %s ...", fname_mat.c_str());
- BaseLib::RunTime timer;
- timer.start();
- CS_read(in, n, iA, jA, A);
- INFO("\t- took %e s", timer.elapsed());
- } else {
- INFO("error reading matrix from %s", fname_mat.c_str());
- return -1;
- }
- unsigned nnz(iA[n]);
- INFO("\tParameters read: n=%d, nnz=%d", n, nnz);
-
-#ifdef HAVE_PTHREADS
- unsigned n_threads(n_cores_arg.getValue());
- MathLib::CRSMatrixPThreads<double> mat (n, iA, jA, A, n_threads);
-
- double *x(new double[n]);
- double *y(new double[n]);
-
- for (unsigned k(0); k<n; ++k)
- x[k] = 1.0;
-
- // read the number of multiplication to execute
- unsigned n_mults (n_mults_arg.getValue());
-
- INFO("*** %d matrix vector multiplications (MVM) with Toms amuxCRS (%d threads) ...", n_mults, n_threads);
- BaseLib::RunTime run_timer;
- BaseLib::CPUTime cpu_timer;
- run_timer.start();
- cpu_timer.start();
- for (std::size_t k(0); k<n_mults; k++) {
- mat.amux (1.0, x, y);
- }
-
- INFO("\t[MVM] - took %e sec cpu time, %e sec run time", cpu_timer.elapsed(), run_timer.elapsed());
-
- delete [] x;
- delete [] y;
-#endif
-
- delete custom_format;
- delete logogCout;
- delete logog_file;
- LOGOG_SHUTDOWN();
-
- return 0;
-}
diff --git a/SimpleTests/SolverTests/BiCGStabDiagPrecond.cpp b/SimpleTests/SolverTests/BiCGStabDiagPrecond.cpp
deleted file mode 100644
index 1765f45eaceec739fb0ecf4674461b6eaa9bae8d..0000000000000000000000000000000000000000
--- a/SimpleTests/SolverTests/BiCGStabDiagPrecond.cpp
+++ /dev/null
@@ -1,83 +0,0 @@
-#include <fstream>
-#include <iostream>
-#include <cmath>
-#include <cstdlib>
-
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
-#include "BaseLib/RunTime.h"
-#include "BaseLib/CPUTime.h"
-
-#include "MathLib/LinAlg/Solvers/BiCGStab.h"
-#include "MathLib/LinAlg/Sparse/CRSMatrixDiagPrecond.h"
-#include "MathLib/vector_io.h"
-
-int main(int argc, char *argv[])
-{
- if (argc != 4) {
- std::cout << "Usage: " << argv[0] << " matrix rhs number-of-threads" << std::endl;
- return -1;
- }
-
- // read number of threads
-// unsigned num_omp_threads (1);
-// num_omp_threads = atoi (argv[3]);
-
- // *** reading matrix in crs format from file
- std::string fname(argv[1]);
- MathLib::CRSMatrixDiagPrecond *mat (new MathLib::CRSMatrixDiagPrecond(fname));
-
- unsigned n (mat->getNRows());
- bool verbose (true);
- if (verbose)
- std::cout << "Parameters read: n=" << n << std::endl;
-
- double *x(new double[n]);
- double *b(new double[n]);
-
- // *** init start vector x
- for (std::size_t k(0); k<n; k++) {
- x[k] = 1.0;
- }
- // *** read rhs
- fname = argv[2];
- std::ifstream in(fname.c_str());
- if (in) {
- read (in, n, b);
- in.close();
- } else {
- std::cout << "problem reading rhs - initializing b with 1.0" << std::endl;
- for (std::size_t k(0); k<n; k++) {
- b[k] = 1.0;
- }
- }
-
-
- if (verbose)
- std::cout << "solving system with BiCGStab method (diagonal preconditioner) ... " << std::flush;
-
- double eps (1.0e-6);
- unsigned steps (4000);
- BaseLib::RunTime run_timer;
- BaseLib::CPUTime cpu_timer;
- run_timer.start();
- cpu_timer.start();
-
- MathLib::BiCGStab ((*mat), b, x, eps, steps);
-
- if (verbose) {
- std::cout << " in " << steps << " iterations" << std::endl;
- std::cout << "\t(residuum is " << eps << ") took " << cpu_timer.elapsed() << " sec time and " << run_timer.elapsed() << " sec" << std::endl;
- } else {
- std::cout << cpu_timer.elapsed() << std::endl;
- }
-
- delete mat;
- delete [] x;
- delete [] b;
-
- return 0;
-}
-
diff --git a/SimpleTests/SolverTests/CMakeLists.txt b/SimpleTests/SolverTests/CMakeLists.txt
deleted file mode 100644
index e25b8f744c2edb5b8d19fa0836b1f6090c8a7e5e..0000000000000000000000000000000000000000
--- a/SimpleTests/SolverTests/CMakeLists.txt
+++ /dev/null
@@ -1,59 +0,0 @@
-add_executable(ConjugateGradientUnpreconditioned
- ConjugateGradientUnpreconditioned.cpp
- ${SOURCES}
- ${HEADERS}
-)
-set_target_properties(ConjugateGradientUnpreconditioned PROPERTIES FOLDER SimpleTests)
-
-add_executable(ConjugateGradientDiagPrecond
- ConjugateGradientDiagonalPreconditioned.cpp
- ${SOURCES}
- ${HEADERS}
-)
-set_target_properties(ConjugateGradientDiagPrecond PROPERTIES FOLDER SimpleTests)
-
-add_executable(BiCGStabDiagPrecond
- BiCGStabDiagPrecond.cpp
- ${SOURCES}
- ${HEADERS}
-)
-set_target_properties(BiCGStabDiagPrecond PROPERTIES FOLDER SimpleTests)
-
-add_executable(GMResDiagPrecond
- GMResDiagPrecond.cpp
- ${SOURCES}
- ${HEADERS}
-)
-set_target_properties(GMResDiagPrecond PROPERTIES FOLDER SimpleTests)
-
-target_link_libraries(ConjugateGradientUnpreconditioned
- ${BLAS_LIBRARIES}
- ${LAPACK_LIBRARIES}
- ${ADDITIONAL_LIBS}
- MathLib
- BaseLib
-)
-
-target_link_libraries(ConjugateGradientDiagPrecond
- ${BLAS_LIBRARIES}
- ${LAPACK_LIBRARIES}
- ${ADDITIONAL_LIBS}
- MathLib
- BaseLib
-)
-
-target_link_libraries(BiCGStabDiagPrecond
- ${BLAS_LIBRARIES}
- ${LAPACK_LIBRARIES}
- ${ADDITIONAL_LIBS}
- MathLib
- BaseLib
-)
-
-target_link_libraries(GMResDiagPrecond
- ${BLAS_LIBRARIES}
- ${LAPACK_LIBRARIES}
- ${ADDITIONAL_LIBS}
- MathLib
- BaseLib
-)
diff --git a/SimpleTests/SolverTests/ConjugateGradientDiagonalPreconditioned.cpp b/SimpleTests/SolverTests/ConjugateGradientDiagonalPreconditioned.cpp
deleted file mode 100644
index 046074c7a81f7f3382ba03ea4a2122327b8fb9fe..0000000000000000000000000000000000000000
--- a/SimpleTests/SolverTests/ConjugateGradientDiagonalPreconditioned.cpp
+++ /dev/null
@@ -1,91 +0,0 @@
-#include <fstream>
-#include <iostream>
-#include <cmath>
-#include <cstdlib>
-
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
-#include "BaseLib/RunTime.h"
-#include "BaseLib/CPUTime.h"
-
-#include "MathLib/LinAlg/Solvers/CG.h"
-#include "MathLib/LinAlg/Sparse/CRSMatrixDiagPrecond.h"
-#include "MathLib/vector_io.h"
-
-int main(int argc, char *argv[])
-{
- if (argc != 4) {
- std::cout << "Usage: " << argv[0] << " matrix rhs number-of-threads" << std::endl;
- return -1;
- }
-
- // read number of threads
- unsigned num_omp_threads (1);
- num_omp_threads = atoi (argv[3]);
-
- // *** reading matrix in crs format from file
- std::string fname(argv[1]);
- MathLib::CRSMatrixDiagPrecond *mat (new MathLib::CRSMatrixDiagPrecond(fname));
-
- unsigned n (mat->getNRows());
- bool verbose (true);
- if (verbose)
- std::cout << "Parameters read: n=" << n << std::endl;
-
- double *x(new double[n]);
- double *b(new double[n]);
-
- // *** init start vector x
- for (std::size_t k(0); k<n; k++) {
- x[k] = 0.0;
- }
- // *** read rhs
- fname = argv[2];
- std::ifstream in(fname.c_str());
- if (in) {
- read (in, n, b);
- in.close();
- } else {
- std::cout << "problem reading rhs - initializing b with 1.0" << std::endl;
- for (std::size_t k(0); k<n; k++) {
- b[k] = 1.0;
- }
- }
-
-
- if (verbose)
- std::cout << "solving system with PCG method (diagonal preconditioner) ... " << std::flush;
-
- double eps (1.0e-6);
- unsigned steps (4000);
- BaseLib::RunTime run_timer;
- BaseLib::CPUTime cpu_timer;
- run_timer.start();
- cpu_timer.start();
-
- if (num_omp_threads == 1) {
- MathLib::CG(mat, b, x, eps, steps);
- } else {
- #ifdef _OPENMP
- MathLib::CGParallel(mat, b, x, eps, steps);
- #else
- std::cout << "OpenMP is not switched on" << std::endl;
- #endif
- }
-
- if (verbose) {
- std::cout << " in " << steps << " iterations" << std::endl;
- std::cout << "\t(residuum is " << eps << ") took " << cpu_timer.elapsed() << " sec time and " << run_timer.elapsed() << " sec" << std::endl;
- } else {
- std::cout << cpu_timer.elapsed() << std::endl;
- }
-
- delete mat;
- delete [] x;
- delete [] b;
-
- return 0;
-}
-
diff --git a/SimpleTests/SolverTests/ConjugateGradientUnpreconditioned.cpp b/SimpleTests/SolverTests/ConjugateGradientUnpreconditioned.cpp
deleted file mode 100644
index 2f722b58c3fbeea2b65149007f4c301e16f39f48..0000000000000000000000000000000000000000
--- a/SimpleTests/SolverTests/ConjugateGradientUnpreconditioned.cpp
+++ /dev/null
@@ -1,71 +0,0 @@
-#include <fstream>
-#include <iostream>
-
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
-#include "MathLib/LinAlg/Solvers/CG.h"
-#include "MathLib/LinAlg/Sparse/CRSMatrix.h"
-#include "MathLib/vector_io.h"
-
-int main(int argc, char *argv[])
-{
- (void) argc;
- (void) argv;
-
- // *** reading matrix in crs format from file
- std::string fname("/work/Thomas Fischer/data/testmat.bin");
-// std::ifstream in(fname.c_str(), std::ios::binary);
- MathLib::CRSMatrix<double, unsigned> *mat (new MathLib::CRSMatrix<double, unsigned>(fname));
-/* double *A(NULL);
- unsigned *iA(NULL), *jA(NULL), n;
- if (in) {
- std::cout << "reading matrix from " << fname << " ... " << std::flush;
- CS_read(in, n, iA, jA, A);
- in.close();
- std::cout << " ok" << std::endl;
- }
- unsigned nnz(iA[n]);
-*/
- unsigned n (mat->getNRows());
- std::cout << "Parameters read: n=" << n << std::endl;
-
- double *x(new double[n]);
- double *b(new double[n]);
-
- // *** init start vector x
- for (std::size_t k(0); k<n; k++) {
- x[k] = 1.0;
- }
- // *** read rhs
- fname = "/work/Thomas Fischer/data/rhs.dat";
- std::ifstream in(fname.c_str());
- if (in) {
- read (in, n, b);
- in.close();
- } else {
- std::cout << "problem reading rhs - initializing b with 1.0" << std::endl;
- for (std::size_t k(0); k<n; k++) {
- b[k] = 1.0;
- }
- }
-
- std::cout << "solving system with CG method ... " << std::flush;
- time_t start_time, end_time;
- time(&start_time);
-// double cg_time (cputime(0.0));
- double eps (1.0e-6);
- unsigned steps (4000);
- CG (mat, b, x, eps, steps);
-// cg_time = cputime(cg_time);
- time(&end_time);
- std::cout << " in " << steps << " iterations (residuum is " << eps << ") took " << /*cg_time <<*/ " sec time and " << (end_time-start_time) << " sec" << std::endl;
-
- delete mat;
- delete [] x;
- delete [] b;
-
- return 0;
-}
-
diff --git a/SimpleTests/SolverTests/GMResDiagPrecond.cpp b/SimpleTests/SolverTests/GMResDiagPrecond.cpp
deleted file mode 100644
index 1680db2b5f3a9806113908829c37c493a27e8f33..0000000000000000000000000000000000000000
--- a/SimpleTests/SolverTests/GMResDiagPrecond.cpp
+++ /dev/null
@@ -1,88 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date 2011-10-05
- * \brief Test for GMRes preconditioner.
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <iostream>
-#include <cstdlib>
-
-#include "BaseLib/CPUTime.h"
-#include "BaseLib/RunTime.h"
-#include "MathLib/LinAlg/Solvers/GMRes.h"
-#include "MathLib/LinAlg/Sparse/CRSMatrixDiagPrecond.h"
-#include "MathLib/vector_io.h"
-
-int main(int argc, char *argv[])
-{
- if (argc != 3) {
- std::cout << "Usage: " << argv[0] << " matrix rhs" << std::endl;
- return -1;
- }
-
- // *** reading matrix in crs format from file
- std::string fname(argv[1]);
- MathLib::CRSMatrixDiagPrecond *mat(new MathLib::CRSMatrixDiagPrecond(fname));
-
- unsigned n(mat->getNRows());
- bool verbose(true);
- if (verbose) std::cout << "Parameters read: n=" << n << std::endl;
-
- double *x(new double[n]);
- double *b(new double[n]);
-
- // *** init start vector x
- for (std::size_t k(0); k < n; k++) {
- x[k] = 1.0;
- }
- // *** read rhs
- fname = argv[2];
- std::ifstream in(fname.c_str());
- if (in) {
- read(in, n, b);
- in.close();
- } else {
- std::cout << "problem reading rhs - initializing b with 1.0"
- << std::endl;
- for (std::size_t k(0); k < n; k++) {
- b[k] = 1.0;
- }
- }
-
- if (verbose)
- std::cout << "solving system with GMRes method (diagonal preconditioner) ... "
- << std::flush;
-
- double eps(1.0e-6);
- unsigned steps(4000);
- BaseLib::RunTime run_timer;
- BaseLib::CPUTime cpu_timer;
- run_timer.start();
- cpu_timer.start();
-
- MathLib::GMRes((*mat), b, x, eps, 30, steps);
-
- if (verbose) {
- std::cout << " in " << steps << " iterations" << std::endl;
- std::cout << "\t(residuum is " << eps << ") took "
- << cpu_timer.elapsed() << " sec time and "
- << run_timer.elapsed() << " sec" << std::endl;
- } else {
- std::cout << cpu_timer.elapsed() << std::endl;
- }
-
- delete mat;
- delete[] x;
- delete[] b;
-
- return 0;
-}
-
diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp
index 78c52800379d93aaa77d15499b631466a5738bcb..af7cbfd8123bb3e878d0ed92b0dd882f543133fe 100644
--- a/Tests/MathLib/TestLinearSolver.cpp
+++ b/Tests/MathLib/TestLinearSolver.cpp
@@ -240,24 +240,6 @@ TEST(Math, CheckInterface_EigenLis)
}
#endif
-#ifdef USE_LIS
-TEST(Math, CheckInterface_Lis)
-{
- // set solver options using Boost property tree
- boost::property_tree::ptree t_root;
- boost::property_tree::ptree t_solver;
- t_root.put("lis", "-i cg -p none -tol 1e-15 -maxiter 1000");
- BaseLib::ConfigTree conf(t_root, "",
- BaseLib::ConfigTree::onerror, BaseLib::ConfigTree::onwarning);
-
- using IntType = MathLib::LisMatrix::IndexType;
-
- MathLib::LisMatrix A(Example1<IntType>::dim_eqs);
- checkLinearSolverInterface<MathLib::LisMatrix, MathLib::LisVector,
- MathLib::LisLinearSolver, IntType>(A, conf);
-}
-#endif
-
#ifdef USE_PETSC
TEST(MPITest_Math, CheckInterface_PETSc_Linear_Solver_basic)
{