Hello,
C++20 generic class Matrix has been updated, with linear system solution.
Previous implementation details have been explained in previous posts: C++ Matrix generic class, C++ concepts and C++20 Matrix generic class updated.
Move operations [update 25 Jan 2024]
The class template has been updated (the source code is here: https://github.com/Fred68/MatrixT/) with the move constructor and move assignment operator:
Matrix(Matrix&& m) : _row{m._row}, _col{m._col}, dat{m.dat}, _iterators{0}
{
m.dat = nullptr;
m._row = m._col = 0;
m._iterators = 0;
#ifdef _DEBUG
cout << "Matrix(Matrix&& m)" << endl;
#endif
} Matrix<DATA> &operator=(Matrix<DATA> &&m)
{
if(this != &m)
{
_col = m._col;
_row = m._row;
_iterators = m-_iterators;
dat = m.dat;
m.dat = nullptr;
m._row = m._col = 0;
m._iterators = 0;
}
#ifdef _DEBUG
cout << "Matrix::&operator=(Matrix&&)" << endl;
#endif
return *this;
}The original object is discarded, after the move operation, with a generic shallow destructor: for this reason, the allocated memory is released and indexes zeroed.
Iterator
Generic class Matrix<T> has been enriched with an iterator.
Modern C++ iterators, allowing compact cycles like:
for(auto x : values) {}have been put aside.
Instead, a simpler approach has been used.
#ifndef MATRIXITERATORH
#define MATRIXITERATORH
/*************************************/
/* */
/* MatrixIterator.h */
/* */
/*************************************/
#include "Matrix.h"
namespace matrix
{
/* Prototyping */
template<class DATA> class Matrix;
template<class DATA> class MatrixIterator
{
private:
Matrix<DATA> *_m;
int ir, ic;
public:
/* Ctor con argomenti */
MatrixIterator(Matrix<DATA> &m)
{
_m = &m;
ir = ic = 0;
_m->_iterators++;
};
/* Dtor con argomenti */
~MatrixIterator()
{
_m->_iterators--;
}
/* Iteration */
DATA *begin();
DATA *end();
DATA *next();
DATA *peek()
{
return _m->dat + ir * _m->_col + ic;
}
};
#include "MatrixIterator.cpp"
}
#endif
#ifndef MATRIXITERATORCPP
#define MATRIXITERATORCPP
#include "MatrixIterator.h"
using namespace matrix;
template <class DATA> DATA *MatrixIterator<DATA>::begin()
{
ir = ic = 0;
if ((_m->_row == 0) || (_m->_col == 0))
{
throw std::runtime_error(MatrixDef::ERR_ZEROSIZE);
}
return _m->dat;
}
template <class DATA> DATA *MatrixIterator<DATA>::end()
{
ir = _m->_row - 1;
ic = _m->_col - 1;
if ((ir < 0) || (ic < 0))
{
throw std::runtime_error(MatrixDef::ERR_ZEROSIZE);
}
return _m->dat[ir * _m->_col + ic];
}
template <class DATA> DATA *MatrixIterator<DATA>::next()
{
ic++;
if (ic >= _m->_col)
{
ic = 0;
ir ++;
}
if (ir >= _m->_row) // Oltre fine matrice
{
ic = ir = 0; // Azzera l'iteratore
return (DATA*) nullptr;
}
return _m->dat + ir * _m->_col + ic;
}
#endifThe usage of this iterator is:
MatrixIterator<int> it(pr1);
cout << "\npr1=\n" << pr1.to_string(MatrixDef::Cmd::detail) << endl;
for (int *x = it.begin(); x != nullptr; x = it.next())
{
cout << *x << " " << *it.peek() << '\t';
}
The generic iterator is prototyped in Matrix header and declared friend in the class declaration.
Considering a MatrixIterator<T> object contains a pointer to the Matrix<T> object on which it handles the iterations, when the Matrix<T> object is destroyed, the iterator might operate on a deallocated memory address.
For this reason, the generic Matrix class contains the number of created iterators. There is only one MatrixIterator<T> constructor, receiving the reference to a Matrix<T> object, so that no empty iterator can be normally created.
Linear systems
Linear system solver has been encapsulated into a specific class, although only necessary functions could have been added to Matrix<T> class. The reason is that the factorization, the solution and some parameters can be controlled in an easier way.
#ifndef LINEARSYSH
#define LINEARSYSH
#include "Matrix.h"
#include <string>
#include <sstream>
#ifdef _DEBUG
using std::cout;
using std::ostream;
#endif
using std::string;
using std::stringstream;
using std::endl;
namespace matrix
{
//template<class DATA> class Matrix; // Prototyping superfluo
template <class DATA, class MOD> class LinearSys
{
private:
Matrix<DATA> a; // Matrice fattorizzata
Matrix<int> pivot; // Vettore di pivot
const MOD EPS = std::numeric_limits<MOD>::epsilon();
const MOD MIN = std::numeric_limits<MOD>::min();
const MOD MAX = std::numeric_limits<MOD>::max();
MOD _epszero; // Valore minimo considerato > 0
MOD _det; // Determinante
public:
/* Ctor. Vincolo nel costruttore, non serve altrove */
LinearSys() requires RQfloatabs<DATA,MOD>
{
_epszero = EPS;
_det = (MOD) 0;
#ifdef _DEBUG
cout << "LinearSys(): epsilon=" << EPS << ", min=" << MIN << ", max=" << MAX << endl;
#endif
}
/* Soglia sotto la quale un valore si considera nullo */
void set_eps_zero(MOD eps_zero)
{
if(eps_zero > EPS)
{
_epszero = eps_zero;
}
}
void set_eps_zero_ratio(MOD eps_zero)
{
if (eps_zero >= (MOD)1)
{
_epszero = eps_zero * EPS;
}
}
MOD get_eps_zero()
{
return _epszero;
}
/* to_string() */
string to_string()
{
stringstream ss;
ss << "Eps_zero=" << _epszero << endl;
ss << "Det=" << _det << endl;
ss << "PLU=" << a.to_string() << endl;
ss << "Pivot=" << pivot.to_string() << endl;
return ss.str();
}
/* Fattorizzazione LU con pivoting parziale e soluzione. Monegato Metodi e algoritmi... CLUT 2008, pag. 41 e succ.*/
bool factor(Matrix <DATA> &A);
bool solve_check(bool throw_exception = true)
{
bool ok = true;
int n = a.rows();
if (n != a.cols())
{
ok = false;
if(throw_exception) throw std::runtime_error(MatrixDef::ERR_NOTSQUARE);
}
if (n < 1)
{
ok = false;
if (throw_exception) throw std::runtime_error(MatrixDef::ERR_ZEROSIZE);
}
if (abs(_det) < _epszero)
{
ok = false;
if (throw_exception) throw std::runtime_error(MatrixDef::ERR_SINGULAR);
}
if ((pivot.rows() != n - 1) || (pivot.cols() != 1))
{
ok = false;
if (throw_exception) throw std::runtime_error(MatrixDef::ERR_PIVOT);
}
return ok;
}
bool solve(Matrix <DATA> &x, Matrix <DATA> &b);
};
#include "LinearSys.cpp"
}
#endif
#ifndef LINEARSYSCPP
#define LINEARSYSCPP
#include "LinearSys.h"
using namespace matrix;
template <class DATA, class MOD> bool LinearSys<DATA, MOD>::factor(Matrix <DATA> &A)
{
int n;
n = A.rows();
if (n != A.cols())
{
throw std::runtime_error(MatrixDef::ERR_NOTSQUARE);
}
if (n < 1)
{
throw std::runtime_error(MatrixDef::ERR_ZEROSIZE);
}
pivot.dim(n-1, 1); // Vettore colonna n-1 elementi
a = A; // Assegnazione della classe Matrix<>
if (n == 1) // Se matrice 1x1: il determinante è il valore.
{
_det = a(0,0);
if (abs(_det) < _epszero)
{
return false;
}
return true;
}
int k, i, io, j; // Ciclo di calcolo
MOD amax;
DATA tmp;
for (_det = (DATA)1, k = 0; k < n - 1; k++) // Ciclo 1 su tutte le righe k=0...n-2
{
for (amax = (MOD)0, io = i = k; i < n; i++) // Cerca la riga (io), dalla k in poi, con il massimo amax sulla colonna k.
{
if (abs(a(i, k)) >= amax)
{
io = i;
amax = abs(a(i, k));
}
}
pivot(k,0) = io; // Imposta il pivot
if (amax < _epszero) // Se il valore massimo dei coefficienti è (quasi) zero...
{
_det = (MOD) 0; // ...azzera il determinante ed esce.
return false;
}
if (io != k) // Se l'indice non è k...
{
for (j = k; j < n; j++) // ...scambia le righe k e io (solo la parte superiore, l'inferiore è nulla)...
{
tmp = a(k, j);
a(k, j) = a(io, j);
a(io, j) = tmp;
}
_det = -_det; // ...e cambia segno al determinante.
}
for (i = k + 1; i < n; i++) // Ciclo 2, per tutte le righe successive alla riga k attuale
{
a(i, k) = -a(i, k) / a(k, k); // Coefficiente per metodo di Gauss (su riga i, colonna k), sotto diagonale.
for (j = k + 1; j < n; j++) // Combina linearmente le due righe, sopra la diagonale
{
a(i, j) = a(i, j) + a(i, k) * a(k, j);
}
} // Fine ciclo 2
_det = _det * a(k, k);
}
if(abs(a(n - 1, n - 1))<_epszero)
{
_det = (MOD) 0;
return false;
}
_det = _det * a(n - 1, n - 1);
return true;
}
template <class DATA, class MOD> bool LinearSys<DATA, MOD>::solve(Matrix <DATA> &x, Matrix <DATA> &b)
{
int n = a.rows();
solve_check();
if ((b.rows() != n) || (b.cols() != 1)) // Verifica dimensioni vettore termini noti
{
throw std::runtime_error(MatrixDef::ERR_SZMISMATCH);
}
x = b;
if (n == 1)
{
if (abs(a(0, 0)) < _epszero)
{
throw std::runtime_error(MatrixDef::ERR_SINGULAR);
}
x(0,0) = x(0,0) / a(0,0);
return true;
}
int k,j,i;
DATA tmp;
for (k = 0; k < n - 1; k++) // Ciclo su tutte le righe, tranne l'ultima
{
j = pivot(k,0); // Ottiene la riga da scambiare con le riga k
if (j != k) // Le scambia, se necessario
{
tmp = x(j,0);
x(j,0) = x(k,0);
x(k, 0) = tmp;
}
for(i=k+1; i<n; i++)
{
x(i, 0) += a(i, k) * x(k, 0);
}
}
x(n-1, 0) = x(n-1, 0) / a(n-1, n-1);
for (i = n - 2; i >= 0; i--)
{
for (tmp = (DATA)0.0, j = i + 1; j < n; j++)
tmp += a(i, j) * x(j, 0);
x(i, 0) = (x(i, 0) - tmp) / a(i, i);
}
return true;
}
#endif
Note that the:
template <class DATA, class MOD> class LinearSys {}is a generic class with two different data types.
The first one, <DATA>, specifies the data type of the values contained in the Matrix<DATA> working with.
The second one, <MOD>, is the data type of the module (absolute value, function abs) of a <DATA> value. If data are complex<double>, you can (you should….), the absolute value type would be a <double>.
The requirement for the data type is to satisfy the previous requirements (RQfloat) and the existence of a ‘abs()’ function:
template<class T, class A> concept RQfloatabs = RQfloat<T> && RQfloat<A> && requires (T x, A a) { a = abs(x);};The solution is made in two steps: matrix factorization with partial pivoting and solution with back substitutions.
Details and algorithms are clearly explained in Monegato numeric calculation texts.
Compiler and precision
For linear system class, precision depends on floating point data type, according to std::numeric_limits.
In particular, a small number can be considered a zero value if under a threshold, epsilon (the minimum representable difference between 1.0 and the next number), or a multiple, which can be set in a linear system class object.
In C++, a double-precision floating point number is represented, according to IEEE7454, on 64 bit (15 digits). Unfortunately, long double, while it is a different data type, is represented on the same number of bits, in Microsoft C++. Trying to obtain better precision, I compiled the C++ program with Visual Studio 2022 standard compiler and with Intel OneAPI DPC++ Compiler 2023, but not compared the results, yet.
Usage
The usage of the new classes is explained in this test program:
#include "Pippo.h"
#include "Matrix.h"
#include "MatrixIterator.h"
#include "LinearSys.h"
#include <iostream>
//#include <math.h>
//#include <complex.h>
//#include <tgmath.h>
#pragma region USING
// Namespace necessari (evitare: using namespace std;)
using std::cout;
using std::endl;
using std::getchar;
using std::operator<<;
using std::ostream;
using std::complex;
using namespace matrix;
using std::string;
#pragma endregion
// Prototyping
void stop_msg();
static int ramdomint(int i, int j);
static int ramdomint10(int i, int j);
static std::complex<double> cpxrc(int r, int c);
static float ramdomfloat10(int i, int j);
static double ramdomdouble20(int i, int j);
int main()
{
srand((unsigned)time(NULL));
cout << "\n\nMatrix test\n" << endl;
try
{
Matrix<Pippo> P1;
cout << "P1:\n" << P1.to_string() << endl;
Matrix<Pippo> *P2;
P2 = new Matrix<Pippo>(2, 3);
(*P2)(0, 1) = Pippo(1);
(*P2)(1, 2) = Pippo(-5);
Pippo xpip = (*P2)(1, 2);
cout << "xpip:\t" << xpip << endl;
(*P2)(1, 0) = Pippo(-8);
cout << "P2:\n" << P2 << "\tRighe: " << P2->rows() << "\tColonne: " << P2->cols() << endl;
cout << "P2:\n" << P2->to_string(MatrixDef::Cmd::size) << endl;
cout << "P2:\n" << P2->to_string(MatrixDef::Cmd::detail) << endl;
Matrix<Pippo> P2n, P2o;
P2o = P2n = (*P2);
Matrix<Pippo> P2t = !(*P2);
cout << "P2n:\n" << P2n.to_string() << endl;
cout << "P2o:\n" << P2n.to_string() << endl;
cout << "P2t:\n" << P2t.to_string() << endl;
cout << "!P2:\n" << (!(*P2)).to_string() << endl;
Matrix<float> AR[3] = { Matrix<float>(2,2,0.1f), Matrix<float>(2,2,0.2f), Matrix<float>(2,2,0.3f) };
for (int i = 0; i < sizeof(AR) / sizeof(AR[0]); i++)
{
cout << "AR[" << i << "]:\n" << AR[i].to_string() << endl;
}
Matrix<float> F1(2, 2, 0.5);
cout << "F1:\n" << F1.to_string() << endl;
Matrix<float> F2(F1);
cout << "F2:\n" << F2.to_string() << endl;
Matrix<Pippo> P3(*P2);
cout << "P3:\n" << P3.to_string() << endl;
P1.clear();
cout << "P1:\n" << P1.to_string() << endl;
Matrix<Pippo> P4(2, 2, Pippo(3));
cout << "P4:\n" << P4.to_string() << endl;
P4.clear();
cout << "P4:\n" << P4.to_string() << endl;
Matrix<int> J1(2, 2, 100);
cout << "J1:\n" << J1.to_string(',', ';') << endl;
J1.dim(4, 3);
Matrix<int> J2(J1);
cout << "[J2:]\t" << J2.to_string(MatrixDef::Cmd::detail) << endl;
J1(2, 1) = 3;
J1(3, 2) = -99;
J1.dim(10, 10);
cout << "J1:\n" << J1.to_string() << endl;
J1.dim(J2);
cout << "J1:\n" << J1.to_string() << endl;
J1.trim(1, 2, 1, 2);
cout << "J1:\n" << J1.to_string() << endl;
J1.clear();
J1.dim(5, 5);
int i, j;
for (i = 0; i < J1.rows(); i++)
for (j = 0; j < J1.cols(); j++)
J1(i, j) = rand() % 1000;
cout << "J1:\n" << J1.to_string() << endl;
i=4;
cout << "J1row[" << i << "]:\n" << (J1.get_row(i)).to_string() << endl;
i=3;
cout << "J1col[" << i <<"]:\n" << (J1.get_col(i)).to_string() << endl;
cout << "J1sub:\n" << (J1.get_sub(0,2,0,2)).to_string() << endl;
J1.rem_row_col(1, 3);
cout << "J1:\n" << J1.to_string() << endl;
J1.rem_row(2);
cout << "J1:\n" << J1.to_string() << endl;
J1.transpose();
cout << "J1tr:\n" << J1.to_string() << endl;
J1.rem_col(1);
cout << "J1:\n" << J1.to_string() << endl;
#if false
x = (*P2)(1, 3); // Genera eccezione out of bound
#endif
delete P2;
Matrix<int> m(3,4);
for (i = 0; i < m.rows(); i++)
for (j = 0; j < m.cols(); j++)
m(i, j) = rand() % 1000;
pair<int,int> imin = m.min_row_col();
cout << "m=\n" << m.to_string() << endl << "min= " << m.min() << " @ " << '[' << imin.first << ',' << imin.second << ']' << endl;
#if false
imin = P3.min(); // Vincolo non soddisfatto
#endif
Matrix<int> s1(3,4, &ramdomint);
Matrix<int> s2(s1.rows(), s1.cols()/* se +1: errore size mismatch*/, &ramdomint);
cout << "\ns1=\n" << s1.to_string() << endl;
cout << "\ns2=\n" << s2.to_string() << endl;
cout << "\ns1+s2=\n" << (s1 + s2).to_string() << endl;
cout << "\ns1-s2=\n" << (s1 - s2).to_string() << endl;
s1+=s2;
cout << "\ns1+=s2;\ns1=\n" << s1.to_string() << endl;
s1-=s2;
cout << "\ns1-=s2\n;s1=\n" << s1.to_string() << endl;
Matrix<int> pr1(2, 3, ramdomint10);
Matrix<int> pr2(3, 2, ramdomint10);
Matrix<int> pr3 = pr1 * pr2;
Matrix<int> pr4 = pr1 * 2;
int ii = 2;
Matrix<int> pr5 = ii * pr1;
Matrix<int> pr6 = 2 * pr1;
Matrix<long double> ld1(2,2,1.0001l);
cout << "[ld1:]\t" << ld1.to_string(MatrixDef::Cmd::detail) << endl;
Matrix<complex<double>> cp0(2,2,cpxrc);
Matrix<float> a(1,3, ramdomfloat10), bb(1,3, ramdomfloat10);
{
MatrixIterator<int> it(pr1);
cout << "\npr1=\n" << pr1.to_string(MatrixDef::Cmd::detail) << endl;
for (int *x = it.begin(); x != nullptr; x = it.next())
{
cout << *x << " " << *it.peek() << '\t';
}
}
// Float 80 bit non supportato: long double => double in MS C++
// __float128 bit forse supportato
cout << "\npr2=\n" << pr2.to_string() << endl;
cout << "\npr3=\n" << pr3.to_string() << endl;
cout << "\npr1*pr2=\n" << (pr1 * pr2).to_string() << endl;
cout << "\npr4=\n" << pr4.to_string() << endl;
cout << "\npr5=\n" << pr5.to_string() << endl;
cout << "\npr6=\n" << pr6.to_string() << endl;
cout << "\n2*pr1=\n" << (2 * pr1).to_string() << endl;
cout << "\npr1*2=\n" << (pr1 * 2).to_string() << endl;
cout << "\npr6/2=\n" << (pr6 / 2).to_string() << endl;
cout << "sizeof(char)=" << sizeof(char) << endl;
cout << "sizeof(short)=" << sizeof(short) << endl;
cout << "sizeof(int)=" << sizeof(int) << endl;
cout << "sizeof(long int)=" << sizeof(long int) << endl;
cout << "sizeof(float)=" << sizeof(float) << endl;
cout << "sizeof(double)=" << sizeof(double) << endl;
cout << "sizeof(long double)=" << sizeof(long double) << endl;
cout << "\nld1=\n" << ld1.to_string() << endl;
cout << "\ncp0=\n" << cp0.to_string() << endl;
cout << "\na=\n" << a.to_string() << endl;
cout << "\nb=\n" << bb.to_string() << endl;
cout << "\na^bb=\n" << (a^bb) << endl;
Matrix<complex<double>> idm;
idm = Matrix<complex<double>>::Id(3);
cout << "\nId\n" << idm.to_string() << endl;
Matrix<float> mset;
float arrcost[3][2] = {{1.1f,2.0f},{3.0f,4.0f},{-1.0f,-2.0f}};
mset.set(3,2, *arrcost);
cout << "\nmset.to_string:\n" << mset.to_string() << endl;
Matrix<float> mset1;
// mset1.set(2, 2,{{10,20},{-3,-4}});
mset1.set(2, 2, { 10.0f,20.0f,
-3.0f,-4.0f
});
cout << "\nmset1.to_string:\n" << mset1.to_string() << endl;
// mset1.set(2, 2, { 10,20,-3,-4,5}); // Error, sizes do not match
LinearSys<float, float> *ls1 = new LinearSys<float, float>();
LinearSys<double, double> *ls2 = new LinearSys<double, double>();
LinearSys<complex<float>, float> *ls3 = new LinearSys<complex<float>, float>();
LinearSys<complex<double>, double> *ls4 = new LinearSys<complex<double>, double>();
ls4->set_eps_zero_ratio(10);
cout << "\nls1->to_string:\n" << ls1->to_string() << endl;
cout << "\nls2->to_string:\n" << ls2->to_string() << endl;
cout << "\nls3->to_string:\n" << ls3->to_string() << endl;
cout << "\nls4->to_string:\n" << ls4->to_string() << endl;
// LinearSys<float,string> *ls2 = new LinearSys<float, string>(); // Error C7500: nessuna funzione soddisfa i vincoli
//Matrix<double> A(3,3,ramdomdouble20);
Matrix<double> A, xo, b, x;
A.set(3,3, { 1, 2, 3,
-2, 4, -7,
5, 6, 1 });
xo.set(3, 1, { 1, -5, 8});
/*
A.set(4, 4, { 0, 2, 0, -1,
2, -1, 1, -2,
1, 0, -2, 1,
-1, 3, 1, 1
});
xo.set(4, 1, {1, 0, -3, 3});
*/
b = A * xo;
cout << "\nA:\n" << A.to_string() << endl;
cout << "\nxo:\n" << xo.to_string() << endl;
cout << "\nb=A*xo:\n" << b.to_string() << endl;
ls2->factor(A);
cout << "\nls2:Factor()\n" << ls2->to_string() << endl;
/* Verifica con GNU Octave
>> A = [1 2 3; -2 4 -7; 5 6 1]
A = 1 2 3
-2 4 -7
5 6 1
>> [L,U]=lu(A)
L = 0.2000 0.1250 1.0000
-0.4000 1.0000 0
1.0000 0 0
U = 5.0000 6.0000 1.0000
0 6.4000 -6.6000
0 0 3.6250
>> bb=[5 3 10]'
bb = 5
3
10
>> x=A\bb
x = -0.1552
1.6983
0.5862
>> A*x
ans = 5
3
10
*/
if(ls2->solve_check())
{
ls2->solve(x, b);
}
else
{
cout << "\nsolve_check non superato" << endl;
}
cout << "\nx=\n" << x.to_string() << endl;
cout << "\nA*x=\n" << (A*x).to_string() << endl;
cout << "\nx-xo\n" << (x-xo).to_string() << endl;
}
catch (const std::runtime_error ex)
{
cout << ex.what();
stop_msg();
exit(1);
}
stop_msg();
}
void stop_msg()
{
cout << endl << "<enter> per chiudere" << endl;
char tempchar = getchar();
}
static int ramdomint(int i, int j)
{return 1 + rand() % 100;}
static int ramdomint10(int i, int j)
{return (rand() % 10) - (rand() % 5);}
static float ramdomfloat10(int i, int j)
{
return (float) ((rand() % 10) - (rand() % 5));
}
static double ramdomdouble20(int i, int j)
{
return (double)((rand() % 20) - (rand() % 10));
}
static std::complex<double> cpxrc(int r, int c)
{
return std::complex<double>(r*10+ramdomint10(0,0), c*10+ ramdomint10(0, 0));
}
The user class Pippo, used with Matrix<>:
#ifndef PIPPO
#define PIPPO
#define ERR_NO_ASSIGN // Prova vincolo operatore assegnazione
#undef ERR_NO_ASSIGN
#include <iostream>
using std::cout;
using std::endl;
using std::ostream;
class Pippo
{
protected:
int _p = -1;
double _p2;
int _in[3] = {1,2,3};
public:
Pippo()
{
_p = 0;
_p2 = (double)_p / 2.0;
#ifdef _DEBUG
cout << "ctor Pippo()" << endl;
#endif
}
Pippo(int x)
{
_p = x;
#ifdef _DEBUG
cout << "ctor Pippo(int x)" << endl;
#endif
}
int Get() { return _p; }
void Set(int p) { _p = p; }
friend ostream &operator<<(ostream &os, const Pippo &pp);
Pippo &operator=(const Pippo &pp)
#ifdef ERR_NO_ASSIGN
= delete;
#else
{
if (this != &pp)
{
_p = pp._p;
#ifdef _DEBUG
cout << "Pippo operator=()" << endl;
#endif
}
else
{
#ifdef _DEBUG
cout << "Pippo auto assignment in operator=()" << endl;
#endif
}
return *this;
}
#endif
};
#endif
#include "Pippo.h"
ostream& operator<<(ostream& os, const Pippo &pp)
{
os << '{' << pp._p << '}';
return os;
}
The progrm’s output:
Matrix test
P1:
[] R0 x C0
xpip: {-5}
P2:
000002AC9A782A60 Righe: 2 Colonne: 3
P2:
R2 x C3
P2:
R2xC3,iters=0,dtsz=32,empty={0}
P2n:
[{0} {1} {0}
{-8} {0} {-5}
] R2 x C3
P2o:
[{0} {1} {0}
{-8} {0} {-5}
] R2 x C3
P2t:
[{0} {-8}
{1} {0}
{0} {-5}
] R3 x C2
!P2:
[{0} {-8}
{1} {0}
{0} {-5}
] R3 x C2
AR[0]:
[0.1 0.1
0.1 0.1
] R2 x C2
AR[1]:
[0.2 0.2
0.2 0.2
] R2 x C2
AR[2]:
[0.3 0.3
0.3 0.3
] R2 x C2
F1:
[0.5 0.5
0.5 0.5
] R2 x C2
F2:
[0.5 0.5
0.5 0.5
] R2 x C2
P3:
[{0} {1} {0}
{-8} {0} {-5}
] R2 x C3
P1:
[] R0 x C0
P4:
[{3} {3}
{3} {3}
] R2 x C2
P4:
[] R0 x C0
J1:
[100,100;100,100;] R2 x C2
[J2:] R4xC3,iters=0,dtsz=4,empty=7143491
J1:
[100 100 0 0 0 0 0 0 0 0
100 100 0 0 0 0 0 0 0 0
0 3 0 0 0 0 0 0 0 0
0 0 -99 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
] R10 x C10
J1:
[100 100 0
100 100 0
0 3 0
0 0 -99
] R4 x C3
J1:
[100 0
3 0
] R2 x C2
J1:
[329 455 353 584 158
121 121 237 749 390
486 461 637 758 829
852 402 479 98 960
30 36 578 711 852
] R5 x C5
J1row[4]:
[30 36 578 711 852
] R1 x C5
J1col[3]:
[584
749
758
98
711
] R5 x C1
J1sub:
[] R0 x C0
J1:
[329 455 353 158
486 461 637 829
852 402 479 960
30 36 578 852
] R4 x C4
J1:
[329 455 353 158
486 461 637 829
30 36 578 852
] R3 x C4
J1tr:
[329 486 30
455 461 36
353 637 578
158 829 852
] R4 x C3
J1:
[329 30
455 36
353 578
158 852
] R4 x C2
m=
[817 507 359 382
998 278 31 701
571 620 596 482
] R3 x C4
min= 31 @ [1,2]
s1=
[96 40 31 80
61 75 38 65
99 17 94 87
] R3 x C4
s2=
[10 49 39 84
89 49 61 8
36 85 10 40
] R3 x C4
s1+s2=
[106 89 70 164
150 124 99 73
135 102 104 127
] R3 x C4
s1-s2=
[86 -9 -8 -4
-28 26 -23 57
63 -68 84 47
] R3 x C4
s1+=s2;
s1=
[106 89 70 164
150 124 99 73
135 102 104 127
] R3 x C4
s1-=s2
;s1=
[96 40 31 80
61 75 38 65
99 17 94 87
] R3 x C4
[ld1:] R2xC2,iters=0,dtsz=8,empty=6.23035e-307
pr1=
R2xC3,iters=1,dtsz=4,empty=7143491
2 2 -1 -1 -2 -2 9 9 0 0 1 1
pr2=
[2 4
1 4
-1 0
] R3 x C2
pr3=
[5 4
17 36
] R2 x C2
pr1*pr2=
[5 4
17 36
] R2 x C2
pr4=
[4 -2 -4
18 0 2
] R2 x C3
pr5=
[4 -2 -4
18 0 2
] R2 x C3
pr6=
[4 -2 -4
18 0 2
] R2 x C3
2*pr1=
[4 -2 -4
18 0 2
] R2 x C3
pr1*2=
[4 -2 -4
18 0 2
] R2 x C3
pr6/2=
[2 -1 -2
9 0 1
] R2 x C3
sizeof(char)=1
sizeof(short)=2
sizeof(int)=4
sizeof(long int)=4
sizeof(float)=4
sizeof(double)=8
sizeof(long double)=8
ld1=
[1.0001 1.0001
1.0001 1.0001
] R2 x C2
cp0=
[(-2,3) (0,15)
(7,4) (10,15)
] R2 x C2
a=
[5 -2 6
] R1 x C3
b=
[2 7 -1
] R1 x C3
a^bb=
-10
Id
[(1,0) (0,0) (0,0)
(0,0) (1,0) (0,0)
(0,0) (0,0) (1,0)
] R3 x C3
mset.to_string:
[1.1 2
3 4
-1 -2
] R3 x C2
mset1.to_string:
[10 20
-3 -4
] R2 x C2
ls1->to_string:
Eps_zero=1.19209e-07
Det=0
PLU=[] R0 x C0
Pivot=[] R0 x C0
ls2->to_string:
Eps_zero=2.22045e-16
Det=0
PLU=[] R0 x C0
Pivot=[] R0 x C0
ls3->to_string:
Eps_zero=1.19209e-07
Det=0
PLU=[] R0 x C0
Pivot=[] R0 x C0
ls4->to_string:
Eps_zero=2.22045e-15
Det=0
PLU=[] R0 x C0
Pivot=[] R0 x C0
A:
[1 2 3
-2 4 -7
5 6 1
] R3 x C3
xo:
[1
-5
8
] R3 x C1
b=A*xo:
[15
-78
-17
] R3 x C1
ls2:Factor()
Eps_zero=2.22045e-16
Det=-116
PLU=[5 6 1
0.4 6.4 -6.6
-0.2 -0.125 3.625
] R3 x C3
Pivot=[2
1
] R2 x C1
x=
[1
-5
8
] R3 x C1
A*x=
[15
-78
-17
] R3 x C1
x-xo
[0
0
0
] R3 x C1
<enter> per chiudere
Git hub
The complete source has been made available as a public repository on Github,
at https://github.com/Fred68/Matrix.
The source code has been put into a single header file.
See project https://github.com/Fred68/MatrixT/
Best regards !
Update (January 2024)
Hello,
files have been updated !
Everything has been put inside a unique header file. Static constant definitions have been moved inside the class definition:
class MatrixDef
{
public:
enum class Cmd { size, detail };
inline static const string ERR_ALLOC = "Allocation failed";
inline static const string ERR_OUTOFBOUND = "Subscript out of bounds";
inline static const string ERR_WRONGPARAM = "Wrong parameter";
inline static const string ERR_SZMISMATCH = "Sizes do not match";
inline static const string ERR_NOTVECTOR = "Not a vector";
inline static const string ERR_ZEROSIZE = "Zero size";
inline static const string ERR_NOTSQUARE = "Not square";
inline static const string ERR_SINGULAR = "Singular";
inline static const string ERR_PIVOT = "Wrong pivot";
};Some constant definitions:
const MOD EPS = std::numeric_limits<MOD>::epsilon();
const MOD MINV = std::numeric_limits<MOD>::min(); have been removed from the template, because min() didn’t compile, probably due to some ambiguity.
This has been solved by removing the constant (re)definition and using the original std:numeric_limits<…> where needed, letting the correct datatype be checked by the ‘requires‘ statement.
See the ctor, as an example:
LinearSys() requires RQfloatabs<DATA, MOD>
{
_epszero = std::numeric_limits<MOD>::epsilon();
_det = (MOD)0;
#ifdef _DEBUG
cout << "LinearSys(): epsilon=" << _epszero << endl;
#endif
}The updated source code is here: https://github.com/Fred68/MatrixT/
Best regards !