Commits (75)
......@@ -75,15 +75,15 @@ trigger-tutorial:
strategy: depend
trigger-mlmc:
trigger-mluq:
stage: downstream
only:
variables:
- ($RUN_PROJECTS == "all" || $RUN_PROJECTS == "mlmc")
- ($RUN_PROJECTS == "all" || $RUN_PROJECTS == "mluq")
variables:
UPSTREAM_COMMIT: ${CI_COMMIT_SHA}
trigger:
project: mpp/mlmc
project: mpp/mluq
branch: feature
strategy: depend
......@@ -127,6 +127,19 @@ trigger-spacetime:
strategy: depend
trigger-dgwave:
stage: downstream
only:
variables:
- ($RUN_PROJECTS == "all" || $RUN_PROJECTS == "dgwave")
variables:
UPSTREAM_COMMIT: ${CI_COMMIT_SHA}
trigger:
project: mpp/dgwave
branch: master
strategy: depend
deploy-mpp:
stage: deploy
variables:
......
## What is M++?
A c++ library providing Finite Element Methods (FEM) to solve Partial Differential Equations (PDE).
It uses MPI to distribute the computations on parallel processes and nodes.
Associated projects:
- https://git.scc.kit.edu/mpp/dgwave
- https://git.scc.kit.edu/mpp/tutorial
- https://git.scc.kit.edu/mpp/spacetimeproject
- https://git.scc.kit.edu/mpp/navierstokes
- https://git.scc.kit.edu/mpp/cardmech
- https://git.scc.kit.edu/mpp/mluq
### Prerequisites
* M++ uses CMake (https://cmake.org/download/) as building tool.
* The packages BLAS (Basic Linear Algebra Subroutines
https://www.netlib.org/blas/) and LAPACK (Linear Algebra PACKage https://www.netlib.org/lapack/) are also required.
* To distribute your computations on parallel processes you need to have Open MPI > 4.0
(Massage Passing Interface https://www.open-mpi.org/).
### Installing M++
* To clone the project with the necessary submodules run:
```git clone --recurse-submodules https://git.scc.kit.edu/mpp/mpp.git```
* Alternatively, clone it and update the submodules by hand:
```git clone https://git.scc.kit.edu/mpp/mpp.git```
```cd mpp```
```git submodule update --init --recursive```
* To install M++, run in the ```mpp``` root directory:
```mkdir build```
```cd build```
```cmake ..```
```make -j```
* Alternatively, simply run:
```./init.sh```
### Run tests
* Unit tests can be executed via:
```cd build```
```ctest```
```python3 mppyrun.py --mpi_tests=1 --mute=1```
* Examples of how to use the library are given by the associated projects.
E.g. have a look at https://git.scc.kit.edu/mpp/tutorial.
git submodule update --init --recursive
mkdir build
cd build
cmake ..
make -j
#ifndef INTERVALARITHMETIC_H
#define INTERVALARITHMETIC_H
#include <string>
void SetIAOutputPrecision(int width = -1, int digits = -1);
std::string OutputDown(double real);
......
......@@ -3,18 +3,11 @@
#include "Assertion.hpp"
#include "Parallel.hpp"
#include "Config.hpp"
#include "Logger.hpp"
void Assertion::assertion(bool b, const char *s, const char *f, int L) {
if (b) return;
std::string msg = "\n\033[1;35mAssert: " + std::string(s) + "\033[0m\n";
if (ParallelProgrammingModel::isInitialized())
msg += " on proc " + std::to_string(PPM->proc()) + "\n";
msg += " in " + std::string(f) + ":" + std::to_string(L) + "\n";
msg += " on line " + std::to_string(L) + "\n\n";
std::cerr << msg;
std::cerr << std::flush;
assert(0);
assertion(b, std::string(s), f, L);
}
void Assertion::assertion(bool b, std::string s, const char *f, int L) {
......@@ -56,13 +49,7 @@ void Assertion::exit(std::string s, const char *f, int L) {
}
void Assertion::warning(const char *s, const char *f, int L) {
if (!PPM->master()) return;
std::string msg = "\n\033[1;33mWarning: " + std::string(s) + "\033[0m\n"
+ " in " + std::string(f) + ":" + std::to_string(L) + "\n";
+ " on line " + std::to_string(L) +
"\n\n";
std::cout << msg;
std::cout << std::flush;
warning(std::string(s), f, L);
}
void Assertion::warning(std::string s, const char *f, int L) {
......@@ -73,21 +60,19 @@ void Assertion::warning(std::string s, const char *f, int L) {
"\n\n";
std::cout << msg;
std::cout << std::flush;
MLogging &mlog = dynamic_cast<MLogging &>(mout);
mlog.AddWarningMsg("Warning: " + s + "\n"
+ " in " + std::string(f) + ":"+std::to_string(L)+ "\n"
+ " on line " + std::to_string(L) +
"\n\n");
}
void Assertion::warningOnProc(int proc, const char *s, const char *f, int L) {
if (PPM->proc() != proc) return;
std::string msg = "\n\033[1;33mWarning: " + std::string(s) + "\033[0m\n";
if (ParallelProgrammingModel::isInitialized())
msg += " on proc " + std::to_string(PPM->proc()) + "\n";
msg += " in " + std::string(f)+ ":" + std::to_string(L) + "\n";
msg += " on line " + std::to_string(L) + "\n\n";
std::cout << msg;
std::cout << std::flush;
warningOnProc(proc, std::string(s), f, L);
}
void Assertion::warningOnProc(int proc, std::string s, const char *f, int L) {
if (PPM->proc() != proc) return;
if (PPM->proc() == proc) {
std::string msg = "\n\033[1;33mWarning: " + s + "\033[0m\n";
if (ParallelProgrammingModel::isInitialized())
msg += " on proc " + std::to_string(PPM->proc()) + "\n";
......@@ -95,16 +80,20 @@ void Assertion::warningOnProc(int proc, std::string s, const char *f, int L) {
msg += " on line " + std::to_string(L) + "\n\n";
std::cout << msg;
std::cout << std::flush;
}
if (PPM->Proc() == 0) { // second if since proc can also be 0!
MLogging &mlog = dynamic_cast<MLogging &>(mout);
std::string msg = "Warning: " + s + "\n";
if (ParallelProgrammingModel::isInitialized())
msg += " on proc " + std::to_string(PPM->proc()) + "\n";
msg += " in " + std::string(f) + ":" + std::to_string(L) + "\n";
msg += " on line " + std::to_string(L) + "\n\n";
mlog.AddWarningMsg(msg);
}
}
void Assertion::warningOnProc(const char *s, const char *f, int L) {
std::string msg = "\n\033[1;33mWarning: " + std::string(s) + "\033[0m\n";
if (ParallelProgrammingModel::isInitialized())
msg += " on proc " + std::to_string(PPM->proc()) + "\n";
msg += " in " + std::string(f) + ":" + std::to_string(L) + "\n";
msg += " on line " + std::to_string(L) + "\n\n";
std::cout << msg;
std::cout << std::flush;
warningOnProc(std::string(s), f, L);
}
void Assertion::warningOnProc(std::string s, const char *f, int L) {
......@@ -115,4 +104,50 @@ void Assertion::warningOnProc(std::string s, const char *f, int L) {
msg += " on line " + std::to_string(L) + "\n\n";
std::cout << msg;
std::cout << std::flush;
if (PPM->Proc() == 0) {
MLogging &mlog = dynamic_cast<MLogging &>(mout);
std::string msg = "Warning: " + s + "\n";
if (ParallelProgrammingModel::isInitialized())
msg += " on proc " + std::to_string(PPM->proc()) + "\n";
msg += " in " + std::string(f) + ":" + std::to_string(L) + "\n";
msg += " on line " + std::to_string(L) + "\n\n";
mlog.AddWarningMsg("Warning: " + s + "\n"
+ " in " + std::string(f) + ":" + std::to_string(L) + "\n"
+ " on line " + std::to_string(L) +
"\n\n");
}
}
void Assertion::synchronizeErrors(bool failed, std::string s, const char *f, int L) {
std::string msg = "";
if (failed) {
exit(s.c_str(), f, L);
msg += "Error: " + s + "\n";
if (ParallelProgrammingModel::isInitialized())
msg += " on proc " + std::to_string(PPM->proc()) + "\n";
msg += " in " + std::string(f) + ":" + std::to_string(L) + "\n";
msg += " on line " + std::to_string(L) + "\n\n";
}
failed = PPM->Or(failed);
if (failed) {
const char *c_msg = msg.c_str();
ExchangeBuffer exBuffer;
for (int i = 0; i < msg.length(); ++i)
exBuffer.Send(0) << c_msg[i];
exBuffer.Communicate();
if (PPM->proc() == 0) {
MLogging &mlog = dynamic_cast<MLogging &>(mout);
for (short p = 0; p < PPM->size(); ++p) {
std::string msg_p;
while (exBuffer.Receive(p).size() < exBuffer.ReceiveSize(p)) {
char c;
exBuffer.Receive(p) >> c;
msg_p += c;
}
mlog.AddErrorMsg(msg_p);
}
}
std::exit(1);
}
}
......@@ -4,6 +4,7 @@
#include <string>
#include <typeinfo>
#include <cassert>
#include <exception>
#include "mpi.h"
......@@ -38,8 +39,37 @@ public:
static void warningOnProc(const char *s, const char *f, int L);
static void warningOnProc(std::string s, const char *f, int L);
static void synchronizeErrors(bool failed, std::string = "", const char *f = "", int L = -1);
};
struct MppException : public std::exception {
std::string error;
const char *file;
int line;
MppException(std::string error, const char *file, int line)
: error(error), file(file), line(line) {}
const char *what() const throw() { return error.c_str(); }
};
#define THROW(s) throw MppException(s, __FILE__, __LINE__);
#define TRY try
#define CATCH(s) \
catch (const MppException &exception) { \
Assertion::synchronizeErrors(true, \
exception.error, \
exception.file, \
exception.line); \
} catch (...) { \
Assertion::synchronizeErrors(true, s, __FILE__, __LINE__); \
} \
Assertion::synchronizeErrors(false);
#define Exit(s) { Assertion::exit(s, __FILE__, __LINE__); exit(1); }
#define Warning(s) { Assertion::warning(s, __FILE__, __LINE__); }
......
......@@ -238,6 +238,18 @@ void MLogging::printFromPLogging(const string &line_p, int proc) {
*this << std::endl;
}
void MLogging::AddErrorMsg(const std::string &msg) {
if (msg.size() == 0) return;
error_msg += msg;
error_cnt++;
}
void MLogging::AddWarningMsg(const std::string &msg) {
if (msg.size() == 0) return;
warning_msg += msg;
warning_cnt++;
}
Logging &MLogging::operator<<(std::ostream &(*f)(std::ostream &)) {
if (f == (std::basic_ostream<char> &(*)(std::basic_ostream<char> &)) &std::endl) {
if (onMaster) {
......@@ -462,15 +474,40 @@ Logging &MLogging::operator<<(const IACInterval &s) {
#endif
std::string getMsg(int error_cnt, int warning_cnt) {
if (error_cnt != 0) {
if (warning_cnt != 0) {
return " with " + std::to_string(error_cnt) + " error(s) and " + std::to_string(warning_cnt)
+ " warning(s):\n";
} else {
return " with " + std::to_string(error_cnt) + " error(s):\n";
}
} else {
if (warning_cnt != 0) {
return " with " + std::to_string(warning_cnt)
+ " warning(s):\n";
} else {
return "";
}
}
}
MLogging::~MLogging() {
if (!onMaster) return;
Date End;
std::string msg = getMsg(error_cnt, warning_cnt);
if (error_cnt != 0)
msg += error_msg;
if (warning_cnt != 0)
msg += warning_msg;
if (screenEnabled)
std::cout << "end program after " << End - startTime
<< " on " << numProcs << " procs at " << End << std::endl;
<< " on " << numProcs << " procs at " << End << msg << std::endl;
if (fileEnabled)
out << "end program after " << End - startTime
<< " on " << numProcs << " procs at " << End << std::endl;
<< " on " << numProcs << " procs at " << End << msg << std::endl;
}
PLogging::PLogging() {
......
......@@ -24,6 +24,11 @@ class MLogging : public Logging {
std::stack<Date> timeStamps;
std::stack<std::string> blockStamps;
std::string error_msg = "";
int error_cnt = 0;
std::string warning_msg = "";
int warning_cnt = 0;
//This information is stored here to avoid dependencies at program end
int numProcs;
bool onMaster;
......@@ -90,6 +95,10 @@ public:
void EndBlockWithoutIndent(bool mute = false) override;
void printFromPLogging(const std::string &line_p, int proc) override;
void AddErrorMsg(const std::string &msg);
void AddWarningMsg(const std::string &msg);
};
class PLogging : public Logging {
......
......@@ -7,6 +7,7 @@ set (math_src
basicalgebra/HermCMatrix.cpp
basicalgebra/CMatrix.cpp
basicalgebra/SparseRMatrix.cpp
basicsolver/BasicSolver.cpp
basicsolver/BasicNonLinearSolver.cpp
basicsolver/BasicSparseSolver.cpp
spectrum/Spectrum.cpp
......
......@@ -8,6 +8,7 @@
template<typename T = double, int dim = SpaceDimension>
class TensorT {
protected:
/// Stores values of Tensor in Rows
std::array<VectorFieldT<T, dim>, dim> t{};
public:
constexpr TensorT() {
......@@ -348,6 +349,26 @@ public:
return n;
}
VectorFieldT<T, dim> Row(int i) const {
if constexpr(dim == 1)
return VectorFieldT<T, dim> (t[i][0]);
if constexpr(dim == 2)
return VectorFieldT<T, dim> (t[i][0], t[i][1]);
if constexpr(dim == 3)
return VectorFieldT<T, dim> (t[i][0], t[i][1], t[i][2]);
else Exit("Not definet for dimensions higher than 3.")
}
VectorFieldT<T, dim> Col(int i) const {
if constexpr(dim == 1)
return VectorFieldT<T, dim> (t[0][i]);
if constexpr(dim == 2)
return VectorFieldT<T, dim> (t[0][i], t[1][i]);
if constexpr(dim == 3)
return VectorFieldT<T, dim> (t[0][i], t[1][i], t[2][i]);
else Exit("Not definet for dimensions higher than 3.")
}
friend bool operator==(const TensorT<T, dim> &R, const TensorT<T, dim> &S) {
bool equal = true;
for (int i = 0; i < dim; ++i)
......@@ -740,6 +761,23 @@ constexpr inline TensorT<T, sDim> Product(const VectorFieldT<T, sDim> &V,
else Exit("Dyadic Product is not defined for dimensions higher than 3")
}
template<typename T, int sDim>
constexpr TensorT<T, sDim> Cross(const TensorT<T, sDim> R, const TensorT<T, sDim> S){
if constexpr(sDim != 3) Exit("Cross Product not defined for dimensions other than 3")
return TensorT<T, sDim>(
R[1][1] * S[2][2] - R[1][2] * S[2][1] + R[2][2] * S[1][1] - R[2][1] * S[1][2],
R[1][2] * S[2][0] - R[1][0] * S[2][2] + R[2][0] * S[1][2] - R[2][2] * S[1][0],
R[1][0] * S[2][1] - R[1][1] * S[2][0] + R[2][1] * S[1][0] - R[2][0] * S[1][1],
R[0][2] * S[2][1] - R[0][1] * S[2][2] + R[2][1] * S[0][2] - R[2][2] * S[0][1],
R[2][2] * S[0][0] - R[2][0] * S[0][2] + R[0][0] * S[2][2] - R[0][2] * S[2][0],
R[2][0] * S[0][1] - R[2][1] * S[0][0] + R[0][1] * S[2][0] - R[0][0] * S[2][1],
R[0][1] * S[1][2] - R[0][2] * S[1][1] + R[1][2] * S[0][1] - R[1][1] * S[0][2],
R[0][2] * S[1][0] - R[0][0] * S[1][2] + R[1][0] * S[0][2] - R[1][2] * S[0][0],
R[0][0] * S[1][1] - R[0][1] * S[1][0] + R[1][1] * S[0][0] - R[1][0] * S[0][1]
);
}
template<typename T, int dim>
constexpr inline TensorT<T, dim> diag(const VectorFieldT<T, dim> &D) {
TensorT<T, dim> S;
......
......@@ -34,6 +34,40 @@ AntisymRMatrixT<REAL> &AntisymRMatrixT<REAL>::transpose() {
return *this;
}
template<>
AntisymRMatrixT<double> &AntisymRMatrixT<double>::Invert() {
if (dim % 2 == 1)
Exit("Antisymmetric matrix not invertible for dim= " + std::to_string(dim))
double a[dim * dim];
for (int i = 0; i < dim; ++i)
for (int j = 0; j < dim; ++j)
a[i * dim + j] = (*this)(i, j);
int I[dim];
int info;
int lwork = 2 * dim;
double work[2 * lwork];
dgetrf_(&dim, &dim, a, &dim, I, &info);
if (info < 0) {
Exit("Error in LAPACK: " + std::to_string(-info) + "-th argument had an illegal value")
} else if (info >0) {
Exit("Error in LAPACK: A(" + std::to_string(-info) + "," + std::to_string(-info) + " is exactly zero.")
}
dgetri_(&dim, a, &dim, I, work, &lwork, &info);
if (info < 0) {
Exit("Error in LAPACK: " + std::to_string(-info) + "-th argument had an illegal value")
} else if (info >0) {
Exit("Error in LAPACK: A(" + std::to_string(-info) + "," + std::to_string(-info) + " is exactly zero.")
}
for (int i = 1; i < dim; ++i)
for (int j = 0; j < i; ++j)
(*this)(a[i * dim + j], i, j);
return *this;
}
template<typename REAL>
void AntisymRMatrixT<REAL>::Accumulate() {
PPM->Sum(z);
......
......@@ -3,6 +3,14 @@
#include "RVector.hpp"
//==================================================================================================
//BLAS/LAPACK routines
//==================================================================================================
//LU factorization
extern "C" void dgetrf_(int *M, int *N, void *A, int *LDA, int *IPIV, int *INFO);
//Invert matrix using LU factorization
extern "C" void dgetri_(int *N, void *A, int *LDA, int *IPIV, void *WORK, int *LWORK, int *INFO);
//==================================================================================================
template<typename REAL=double>
class AntisymRMatrixT {
......@@ -113,6 +121,8 @@ public:
AntisymRMatrixT &transpose();
AntisymRMatrixT &Invert();
void Accumulate();
Saver &save(Saver &saver) const;
......@@ -182,6 +192,12 @@ inline AntisymRMatrixT<REAL> transpose(const AntisymRMatrixT<REAL> &A) {
return B.transpose();
}
template<typename REAL>
inline AntisymRMatrixT<REAL> invert(const AntisymRMatrixT<REAL> &a) {
AntisymRMatrixT<REAL> b(a);
return b.Invert();
}
template<typename REAL>
inline Saver &operator<<(Saver &saver, const AntisymRMatrixT<REAL> &A) {
return A.save(saver);
......
......@@ -14,42 +14,42 @@ extern "C" void zgetrf_(int *M, int *N, void *A, int *LDA, int *IPIV, int *INFO)
extern "C" void zgetri_(int *N, void *A, int *LDA, int *IPIV, void *WORK, int *LWORK, int *INFO);
//==================================================================================================
template<typename REAL, typename COMPLEX>
CVectorT<REAL, COMPLEX> CMatrixT<REAL, COMPLEX>::diag() {
CVectorT<REAL, COMPLEX> d(std::min(Nrows, Ncols));
template<typename REAL>
CVectorT<REAL> CMatrixT<REAL>::diag() {
CVectorT<REAL> d(std::min(Nrows, Ncols));
for (int i = 0; i < std::min(Nrows, Ncols); ++i)
d[i] = z[i * Ncols + i];
return d;
}
template<typename REAL, typename COMPLEX>
void CMatrixT<REAL, COMPLEX>::Accumulate() {
template<typename REAL>
void CMatrixT<REAL>::Accumulate() {
PPM->Sum(z);
}
template<typename REAL, typename COMPLEX>
CMatrixT<REAL, COMPLEX> &CMatrixT<REAL, COMPLEX>::conj() {
template<typename REAL>
CMatrixT<REAL> &CMatrixT<REAL>::conj() {
for (int i = 0; i < z.size(); ++i) z[i] = baConj(z[i]);
return *this;
}
template<typename REAL, typename COMPLEX>
RMatrixT<REAL> CMatrixT<REAL, COMPLEX>::real() const {
template<typename REAL>
RMatrixT<REAL> CMatrixT<REAL>::real() const {
RMatrixT<REAL> re(Nrows, Ncols);
for (int i = 0; i < z.size(); ++i) re.Data()[i] = baReal(z[i]);
return re;
}
template<typename REAL, typename COMPLEX>
RMatrixT<REAL> CMatrixT<REAL, COMPLEX>::imag() const {
template<typename REAL>
RMatrixT<REAL> CMatrixT<REAL>::imag() const {
RMatrixT<REAL> im(Nrows, Ncols);
for (int i = 0; i < z.size(); ++i) im.Data()[i] = baImag(z[i]);
return im;
}
template<typename REAL, typename COMPLEX>
CMatrixT<REAL, COMPLEX> &CMatrixT<REAL, COMPLEX>::transpose() {
CMatrixT<REAL, COMPLEX> b(*this);
template<typename REAL>
CMatrixT<REAL> &CMatrixT<REAL>::transpose() {
CMatrixT<REAL> b(*this);
resize(b.Ncols, b.Nrows);
for (int i = 0; i < Nrows; ++i)
for (int j = 0; j < Ncols; ++j)
......@@ -57,30 +57,30 @@ CMatrixT<REAL, COMPLEX> &CMatrixT<REAL, COMPLEX>::transpose() {
return *this;
}
template<typename REAL, typename COMPLEX>
COMPLEX CMatrixT<REAL, COMPLEX>::Mean() const {
template<typename REAL>
COMPLEX_TYPE<REAL> CMatrixT<REAL>::Mean() const {
COMPLEX mean{};
for (auto const &value_vec: z)
mean += value_vec;
return mean / z.size();
}
template<typename REAL, typename COMPLEX>
REAL CMatrixT<REAL, COMPLEX>::Variance() const {
template<typename REAL>
REAL CMatrixT<REAL>::Variance() const {
REAL var_real = real().Variance();
REAL var_imag = imag().Variance();
return var_real + var_imag;
}
template<typename REAL, typename COMPLEX>
CMatrixT<REAL, COMPLEX> &CMatrixT<REAL, COMPLEX>::adjoint() {
template<typename REAL>
CMatrixT<REAL> &CMatrixT<REAL>::adjoint() {
transpose();
conj();
return *this;
}
template<typename REAL, typename COMPLEX>
CMatrixT<REAL, COMPLEX> &CMatrixT<REAL, COMPLEX>::Identity() {
template<typename REAL>
CMatrixT<REAL> &CMatrixT<REAL>::Identity() {
for (int i = 0; i < Nrows; i++)
for (int j = 0; j < Ncols; j++)
z[i * Ncols + j] = COMPLEX(i == j);
......@@ -88,7 +88,7 @@ CMatrixT<REAL, COMPLEX> &CMatrixT<REAL, COMPLEX>::Identity() {
}
template<>
CMatrixT<double, std::complex<double>> &CMatrixT<double, std::complex<double>>::Invert() {
CMatrixT<double> &CMatrixT<double>::Invert() {
if (Nrows != Ncols)
Exit("Nonquadratic matix in invert")
......@@ -112,15 +112,15 @@ CMatrixT<double, std::complex<double>> &CMatrixT<double, std::complex<double>>::
return *this;
}
template<typename REAL, typename COMPLEX>
Saver &CMatrixT<REAL, COMPLEX>::save(Saver &saver) const {
template<typename REAL>
Saver &CMatrixT<REAL>::save(Saver &saver) const {
saver << Nrows << Ncols;
for (int i = 0; i < z.size(); ++i) saver << z[i];
return saver;
}
template<typename REAL, typename COMPLEX>
Loader &CMatrixT<REAL, COMPLEX>::load(Loader &loader) {
template<typename REAL>
Loader &CMatrixT<REAL>::load(Loader &loader) {
int r, c;
loader >> r >> c;
resize(r, c);
......@@ -128,20 +128,20 @@ Loader &CMatrixT<REAL, COMPLEX>::load(Loader &loader) {
return loader;
}
template<typename REAL, typename COMPLEX>
void CMatrixT<REAL, COMPLEX>::resize(int rows, int cols) {
template<typename REAL>
void CMatrixT<REAL>::resize(int rows, int cols) {
z = std::vector<COMPLEX>(rows * cols, COMPLEX{});
Nrows = rows;
Ncols = cols;
}
template
class CMatrixT<double, std::complex<double>>;
class CMatrixT<double>;
#ifdef BUILD_IA
template
class CMatrixT<IAInterval, IACInterval>;
class CMatrixT<IAInterval>;
CMatrix mid(const IACMatrix &IA_A) {
CMatrix A(IA_A.rows(), IA_A.cols());
......
This diff is collapsed.
......@@ -4,55 +4,55 @@
#include "CFunctions.hpp"
template<typename REAL, typename COMPLEX>
CVectorT<REAL, COMPLEX>::CVectorT(int length) : z(length, 0.0) {}
template<typename REAL>
CVectorT<REAL>::CVectorT(int length) : z(length, 0.0) {}
template<typename REAL, typename COMPLEX>
CVectorT<REAL, COMPLEX> &CVectorT<REAL, COMPLEX>::conj() {
template<typename REAL>
CVectorT<REAL> &CVectorT<REAL>::conj() {
for (int i = 0; i < z.size(); ++i) z[i] = baConj(z[i]);
return *this;
}
template<typename REAL, typename COMPLEX>
RVectorT<REAL> CVectorT<REAL, COMPLEX>::real() const {
template<typename REAL>
RVectorT<REAL> CVectorT<REAL>::real() const {
RVectorT<REAL> re(z.size());
for (int i = 0; i < z.size(); ++i) re[i] = baReal(z[i]);
return re;
}
template<typename REAL, typename COMPLEX>
RVectorT<REAL> CVectorT<REAL, COMPLEX>::imag() const {
template<typename REAL>
RVectorT<REAL> CVectorT<REAL>::imag() const {
RVectorT<REAL> im(z.size());
for (int i = 0; i < z.size(); ++i) im[i] = baImag(z[i]);
return im;
}
template<typename REAL, typename COMPLEX>
REAL CVectorT<REAL, COMPLEX>::normSqr() const {
template<typename REAL>
REAL CVectorT<REAL>::normSqr() const {
REAL sum = REAL{};
for (int i = 0; i < z.size(); ++i) sum += baAbsSqr(z[i]);
return sum;
}
template<typename REAL, typename COMPLEX>
REAL CVectorT<REAL, COMPLEX>::norm() const {
template<typename REAL>
REAL CVectorT<REAL>::norm() const {
return sqrt(normSqr());
}
template<typename REAL, typename COMPLEX>
void CVectorT<REAL, COMPLEX>::Accumulate() {
template<typename REAL>
void CVectorT<REAL>::Accumulate() {
PPM->Sum(z);
}
template<typename REAL, typename COMPLEX>
Saver &CVectorT<REAL, COMPLEX>::save(Saver &saver) const {
template<typename REAL>
Saver &CVectorT<REAL>::save(Saver &saver) const {
saver << size();
for (int i = 0; i < z.size(); ++i) saver << z[i];
return saver;
}
template<typename REAL, typename COMPLEX>
Loader &CVectorT<REAL, COMPLEX>::load(Loader &loader) {
template<typename REAL>
Loader &CVectorT<REAL>::load(Loader &loader) {
int N;
loader >> N;
resize(N);
......@@ -60,28 +60,28 @@ Loader &CVectorT<REAL, COMPLEX>::load(Loader &loader) {
return loader;
}
template<typename REAL, typename COMPLEX>
COMPLEX CVectorT<REAL, COMPLEX>::Mean() const {
COMPLEX mean{};
template<typename REAL>
COMPLEX_TYPE<REAL> CVectorT<REAL>::Mean() const {
COMPLEX_TYPE<REAL> mean{};
for (auto const &value_vec: z)
mean += value_vec;
return mean / z.size();
}
template<typename REAL, typename COMPLEX>
REAL CVectorT<REAL, COMPLEX>::Variance() const {
template<typename REAL>
REAL CVectorT<REAL>::Variance() const {
REAL var_real = real().Variance();
REAL var_imag = imag().Variance();
return var_real + var_imag;
}
template
class CVectorT<double, std::complex<double>>;
class CVectorT<double>;
#ifdef BUILD_IA
template
class CVectorT<IAInterval, IACInterval>;
class CVectorT<IAInterval>;
CVector mid(const IACVector &IA_v) {
CVector v(IA_v.size());
......
......@@ -3,369 +3,382 @@
#include "RVector.hpp"
#include <complex>
#ifdef BUILD_IA
#include "IACInterval.hpp"
#endif
template<typename T>
using COMPLEX_TYPE
= typename std::conditional<std::is_same_v<T, double>, std::complex<double>,
typename std::conditional<std::is_same_v<T, IAInterval>, IACInterval, void>::type>::type;
#include <complex>
#else
template<typename T>
using COMPLEX_TYPE
= typename std::conditional<std::is_same_v<T, double>, std::complex<double>, void>::type;
#endif
namespace std {
inline std::complex<double> operator*(const std::complex<double> &a, int b) {
inline std::complex<double> operator*(const std::complex<double> &a, int b) {
return a * double(b);
}
}
inline std::complex<double> operator*(int a, const std::complex<double> &b) {
inline std::complex<double> operator*(int a, const std::complex<double> &b) {
return double(a) * b;
}
}
inline std::complex<double> operator/(const std::complex<double> &a, int b) {
inline std::complex<double> operator/(const std::complex<double> &a, int b) {
return a / double(b);
}
}
}
template<typename REAL=double, typename COMPLEX=std::complex<double>>
template<typename REAL=double>
class CVectorT {
using COMPLEX = COMPLEX_TYPE<REAL>;
protected:
std::vector<COMPLEX> z;
std::vector<COMPLEX> z;
public:
CVectorT() = default;
explicit CVectorT(int length);
template<typename REAL1, typename COMPLEX1>
explicit CVectorT(const CVectorT<REAL1, COMPLEX1> &v) : z(v.size()) {
for (int i = 0; i < z.size(); ++i) z[i] = v[i];
}
template<typename REAL1>
explicit CVectorT(const RVectorT<REAL1> &v) : z(v.size()) {
for (int i = 0; i < z.size(); ++i) z[i] = v[i];
}
template<typename COMPLEX1>
CVectorT(const std::vector<COMPLEX1> &v) : z(v.begin(), v.end()) {}
template<typename T>
constexpr CVectorT(const T &a, int length) : z(length) {
for (int i = 0; i < z.size(); ++i) z[i] = a;
}
template<typename T>
CVectorT &operator=(const T &a) {
for (int i = 0; i < z.size(); ++i) z[i] = a;
return *this;
}
template<typename REAL1, typename COMPLEX1>
CVectorT &operator=(const CVectorT<REAL1, COMPLEX1> &v) {
z.resize(v.size());
for (int i = 0; i < z.size(); ++i) z[i] = v[i];
return *this;
}
CVectorT() = default;
explicit CVectorT(int length);
template<typename REAL1>
explicit CVectorT(const CVectorT<REAL1> &v) : z(v.size()) {
for (int i = 0; i < z.size(); ++i) z[i] = v[i];
}
template<typename REAL1>
explicit CVectorT(const RVectorT<REAL1> &v) : z(v.size()) {
for (int i = 0; i < z.size(); ++i) z[i] = v[i];
}
template<typename COMPLEX1>
CVectorT(const std::vector<COMPLEX1> &v) : z(v.begin(), v.end()) {}
template<typename T>
constexpr CVectorT(const T &a, int length) : z(length) {
for (int i = 0; i < z.size(); ++i) z[i] = a;
}
template<typename T>
CVectorT &operator=(const T &a) {
for (int i = 0; i < z.size(); ++i) z[i] = a;
return *this;
}
template<typename REAL1>
CVectorT &operator=(const CVectorT<REAL1> &v) {
z.resize(v.size());
for (int i = 0; i < z.size(); ++i) z[i] = v[i];
return *this;
}
template<typename REAL1>
CVectorT &operator=(const RVectorT<REAL1> &v) {
z.resize(v.size());
for (int i = 0; i < z.size(); ++i) z[i] = v[i];
return *this;
}
template<typename REAL1>
CVectorT &operator=(const RVectorT<REAL1> &v) {
z.resize(v.size());
for (int i = 0; i < z.size(); ++i) z[i] = v[i];
return *this;
}
template<typename T>
CVectorT &operator+=(const T &a) {
for (int i = 0; i < z.size(); ++i) z[i] += a;
return *this;
}
template<typename T>
CVectorT &operator+=(const T &a) {
for (int i = 0; i < z.size(); ++i) z[i] += a;
return *this;
}
template<typename REAL1>
CVectorT &operator+=(const CVectorT<REAL1> &v) {
if (z.size() != v.size()) Exit("Size does not fit!");
for (int i = 0; i < z.size(); ++i) z[i] += v[i];
return *this;
}
template<typename REAL1, typename COMPLEX1>
CVectorT &operator+=(const CVectorT<REAL1, COMPLEX1> &v) {
if (z.size() != v.size()) Exit("Size does not fit!");
for (int i = 0; i < z.size(); ++i) z[i] += v[i];
return *this;
}
template<typename REAL1>
CVectorT &operator+=(const RVectorT<REAL1> &v) {
if (z.size() != v.size()) Exit("Size does not fit!");
for (int i = 0; i < z.size(); ++i) z[i] += v[i];
return *this;
}
template<typename REAL1>
CVectorT &operator+=(const RVectorT<REAL1> &v) {
if (z.size() != v.size()) Exit("Size does not fit!");
for (int i = 0; i < z.size(); ++i) z[i] += v[i];
return *this;
}
template<typename T>
CVectorT &operator-=(const T &a) {
for (int i = 0; i < z.size(); ++i) z[i] -= a;
return *this;
}
template<typename T>
CVectorT &operator-=(const T &a) {
for (int i = 0; i < z.size(); ++i) z[i] -= a;
return *this;
}
template<typename REAL1>
CVectorT &operator-=(const CVectorT<REAL1> &v) {
if (z.size() != v.size()) Exit("Size does not fit!");
for (int i = 0; i < z.size(); ++i) z[i] -= v[i];
return *this;
}
template<typename REAL1, typename COMPLEX1>
CVectorT &operator-=(const CVectorT<REAL1, COMPLEX1> &v) {
if (z.size() != v.size()) Exit("Size does not fit!");
for (int i = 0; i < z.size(); ++i) z[i] -= v[i];
return *this;
}
template<typename REAL1>
CVectorT &operator-=(const RVectorT<REAL1> &v) {
if (z.size() != v.size()) Exit("Size does not fit!");
for (int i = 0; i < z.size(); ++i) z[i] -= v[i];
return *this;
}
template<typename REAL1>
CVectorT &operator-=(const RVectorT<REAL1> &v) {
if (z.size() != v.size()) Exit("Size does not fit!");
for (int i = 0; i < z.size(); ++i) z[i] -= v[i];
return *this;
}
template<typename T>
CVectorT &operator*=(const T &a) {
for (int i = 0; i < z.size(); ++i) z[i] *= a;
return *this;
}
template<typename T>
CVectorT &operator*=(const T &a) {
for (int i = 0; i < z.size(); ++i) z[i] *= a;
return *this;
}
template<typename REAL1>
CVectorT &operator*=(const RVectorT<REAL1> &v) {
for (int i = 0; i < z.size(); ++i) z[i] *= v[i];
return *this;
}
template<typename REAL1>
CVectorT &operator*=(const RVectorT<REAL1> &v) {
for (int i = 0; i < z.size(); ++i) z[i] *= v[i];
return *this;
}
// This a component wise operation
template<typename REAL1>
CVectorT &operator*=(const CVectorT<REAL1> &v) {
for (int i = 0; i < z.size(); ++i) z[i] *= v[i];
return *this;
}
// This a component wise operation
template<typename REAL1, typename COMPLEX1>
CVectorT &operator*=(const CVectorT<REAL1, COMPLEX1> &v) {
for (int i = 0; i < z.size(); ++i) z[i] *= v[i];
return *this;
}
template<typename T>
CVectorT &operator/=(const T &a) {
for (int i = 0; i < z.size(); ++i) z[i] /= a;
return *this;
}
template<typename T>
CVectorT &operator/=(const T &a) {
for (int i = 0; i < z.size(); ++i) z[i] /= a;
return *this;
}
const std::vector<COMPLEX> &asVector() const { return z; }
const std::vector<COMPLEX> &asVector() const { return z; }
std::vector<COMPLEX> &asVector() { return z; }
int Dim() const { return z.size(); }
int Dim() const { return z.size(); }
int size() const { return z.size(); }
int size() const { return z.size(); }
void resize(int N) { z.resize(N); }
void resize(int N) { z.resize(N); }
COMPLEX &operator[](int i) { return z[i]; }
COMPLEX &operator[](int i) { return z[i]; }
const COMPLEX &operator[](int i) const { return z[i]; }
const COMPLEX &operator[](int i) const { return z[i]; }
auto begin() const { return z.begin(); }
auto begin() const { return z.begin(); }
auto begin() { return z.begin(); }
auto begin() { return z.begin(); }
auto rbegin() const { return z.rbegin(); }
auto rbegin() const { return z.rbegin(); }
auto end() const { return z.end(); }
auto end() const { return z.end(); }
auto end() { return z.end(); }
auto end() { return z.end(); }
auto rend() const { return z.rend(); }
auto rend() const { return z.rend(); }
auto push_back(const COMPLEX &a) { return z.push_back(a); }
auto push_back(const COMPLEX &a) { return z.push_back(a); }
void clear() { z.clear(); }
void clear() { z.clear(); }
auto append(const CVectorT &v) {
return z.insert(z.end(), v.begin(), v.end());
}
template<typename COMPLEX1>
auto append(const CVectorT<COMPLEX1> &v) {
return z.insert(z.end(), v.begin(), v.end());
}
CVectorT &conj();
CVectorT &conj();
RVectorT<REAL> real() const;
RVectorT<REAL> real() const;
RVectorT<REAL> imag() const;
RVectorT<REAL> imag() const;
REAL normSqr() const;
REAL normSqr() const;
REAL norm() const;