SaraP
459 lines
21 KiB
C++
459 lines
21 KiB
C++
////////////////////////////////////////////////////////////////////////////////
|
|
/**
|
|
\file \brief \ru Квадратная матрица чисел N x N.
|
|
\en Square matrix of numbers N x N. \~
|
|
*/
|
|
////////////////////////////////////////////////////////////////////////////////
|
|
|
|
#ifndef __MB_MATRIXNN_H
|
|
#define __MB_MATRIXNN_H
|
|
|
|
|
|
#include <math_define.h>
|
|
#include <alg_indicator.h>
|
|
#include <mb_property_title.h>
|
|
#include <mb_enum.h>
|
|
|
|
|
|
//------------------------------------------------------------------------------
|
|
/** \brief \ru Квадратная матрица чисел N x N.
|
|
\en Square matrix of numbers N x N. \~
|
|
\ingroup Base_Items
|
|
*/
|
|
// ---
|
|
class MATH_CLASS MatrixNN
|
|
{
|
|
private :
|
|
double ** parr; ///< \ru Указатель на первый элемент матрицы. \en A pointer to the first element of the matrix.
|
|
size_t n; ///< \ru Размерность матрицы. \en A dimension of the matrix.
|
|
|
|
protected:
|
|
/// \ru Конструктор. \en Constructor.
|
|
MatrixNN() : parr( nullptr ), n( 0 ) {}
|
|
/// \ru Конструктор по заданной размерности. \en The constructor by a given dimension.
|
|
MatrixNN( size_t dim ) : parr( nullptr ), n( 0 ) { SetSize( dim ); }
|
|
public:
|
|
/// \ru Конструктор ограниченной размерности. \en The constructor of restricted dimension.
|
|
MatrixNN ( const uint16 & dim ) : parr( nullptr ), n( 0 ) { SetSize( dim ); }
|
|
/// \ru Конструктор копирования. \en The copy constructor.
|
|
explicit MatrixNN ( const MatrixNN & );
|
|
/// \ru Деструктор. \en Destructor.
|
|
~MatrixNN() { SetSize( 0 ); }
|
|
|
|
public:
|
|
/// \ru Конструктор по заданной размерности. \en The constructor by a given dimension.
|
|
static MatrixNN * Create( size_t m );
|
|
|
|
public: // Общие методы матриц (двумерных массивов)
|
|
size_t Lines () const { return n; } ///< \ru Дать количество строк матрицы. \en Give the number of matrix rows.
|
|
size_t Columns() const { return n; } ///< \ru Дать количество столбцов матрицы. \en Give the number of matrix columns.
|
|
size_t Count () const { return (n*n); } ///< \ru Количество элементов матрицы. \en Give the number of matrix elements.
|
|
c3d::IndicesPair GetSize() const { return c3d::IndicesPair( n, n ); } ///< \ru Дать размер матрицы. \en Give the size of the matrix.
|
|
bool SetSize( c3d::IndicesPair sz ) { size_t m = std_min( sz.first, sz.second ); return SetSize( m ); } ///< \ru Установить размер. \en Set size.
|
|
bool SetSize( size_t lsz, size_t csz, bool save_vals = false ) { size_t m = std_min( lsz, csz ); return SetSize( m, save_vals ); } ///< \ru Установить размер. \en Set size.
|
|
bool SetSize( size_t dim, bool save_vals = false ); ///< \ru Установить размер. \en Set size.
|
|
|
|
/// \ru Получить элемент матрицы (i,j). \en Get an element of the matrix (i,j).
|
|
const double & GetElem( size_t i, size_t j ) const { C3D_ASSERT( !!parr && i < n && j < n ); return parr[i][j]; }
|
|
/// \ru Установить элемент матрицы (i,j). \en Set an element of the matrix (i,j).
|
|
void SetElem( size_t i, size_t j, double v ) { C3D_ASSERT( !!parr && i < n && j < n ); parr[i][j] = v; }
|
|
/// \ru Добавить элемент матрицы (i,j). \en Add an element of the matrix (i,j).
|
|
void AddElem( size_t i, size_t j, double v ) { C3D_ASSERT( !!parr && i < n && j < n ); parr[i][j] += v; }
|
|
/// \ru Получить элемент матрицы (i,j). \en Get an element of the matrix (i,j).
|
|
const double & operator() ( size_t i, size_t j ) const { C3D_ASSERT( !!parr && i < n && j < n ); return parr[i][j]; }
|
|
/// \ru Обнулить матрицу. \en Set the matrix to null.
|
|
MatrixNN & SetZero();
|
|
/// \ru Инициализировать элементами другой матрицы. \en Initialize by elements of another matrix.
|
|
bool Init( const MatrixNN & );
|
|
/// \ru Оператор присваивания. \en The assignment operator.
|
|
MatrixNN & operator = ( const MatrixNN & mtr ) { Init( mtr ); return *this; }
|
|
/// \ru Поменять местами строки. \en Swap lines.
|
|
bool SwapLines( size_t ln1, size_t ln2 );
|
|
|
|
public:
|
|
/// \ru Установить элемент матрицы (i,j). \en Set an element of the matrix (i,j).
|
|
double & operator() ( size_t i, size_t j ) { C3D_ASSERT( !!parr && i < n && j < n ); return parr[i][j]; }
|
|
/// \ru Выдать адрес начала строки матрицы. \en Get an address of the matrix row start .
|
|
const double * GetLine( size_t i ) const { C3D_ASSERT( !!parr && i < n ); return parr[i]; }
|
|
/// \ru Выдать адрес начала строки матрицы. \en Get an address of the matrix row start .
|
|
double * SetLine( size_t i ) { C3D_ASSERT( !!parr && i < n ); return parr[i]; } // SKIP_SA
|
|
/// \ru Инициировать элемент. \en Initiate an element.
|
|
void Init( size_t i, size_t j, double v ) { C3D_ASSERT( !!parr && i < n && j < n ); parr[i][j] = v; }
|
|
/// \ru Установить строку. \en Set a row.
|
|
void SetLine( size_t i, double * p ) { C3D_ASSERT( !!parr && i < n ); parr[i] = p; }
|
|
/// \ru Выдать адрес матрицы. \en Get an address of the matrix.
|
|
double ** SetParr() { return parr; }
|
|
|
|
/// \ru Сделать матрицу единичной. \en Set the matrix to unit.
|
|
void SetSingle();
|
|
/// \ru Увеличить размерность, добавив строку и столбец в конец. \en Increase the dimension by adding a row and coloumn to the end.
|
|
bool Add();
|
|
/// \ru Удалить строку и столбец. \en Remove a row and column.
|
|
void Delete( size_t i );
|
|
|
|
private:
|
|
// \ru Инициализировать матрицу и установить размерность. \en Initialize a matrix and set dimension.
|
|
void SetParrAndDimension( double ** p, size_t i ) { parr = p; n = i; }
|
|
}; // MatrixNN
|
|
|
|
|
|
//------------------------------------------------------------------------------
|
|
/** \brief \ru Решение системы линейных уравнений методом исключения Гаусса.
|
|
\en System of linear equations is solved by the Gauss method. \~
|
|
\details \ru Решение системы линейных уравнений методом исключения Гаусса. \n
|
|
\en System of linear equations is solved by the Gauss method. \n \~
|
|
\param[in] a - \ru Матрица коэффициентов при неизвестных
|
|
\en Coefficient matrix \~
|
|
\param[in] b - \ru Массив правых частей, в него же помещается результат решения
|
|
\en Array of right parts, on output it contains the result of the solution \~
|
|
\param[in] epsilon - \ru Погрешность решения
|
|
\en Tolerance of solution \~
|
|
\param[in] baseProgBar - \ru Индикатор процесса решения
|
|
\en Progress indicator of solution \~
|
|
\return \ru Код ошибки: если nr_Success (+1), то система решена, если nr_Special, то нет решений или система вырождена.
|
|
\en Error code: if nr_Success (+1), then the system is solved, if nr_Special, there is no solution or the system is degenerate. \~
|
|
\ingroup Base_Items
|
|
*/
|
|
// ---
|
|
template <class Type, class MatrixNN>
|
|
MbeNewtonResult TypedGaussEquation ( MatrixNN & a, Type * b, double epsilon, ProgressBarWrapper * baseProgBar = nullptr )
|
|
{
|
|
ProgressBarWrapper * progBar = nullptr;
|
|
if ( baseProgBar != nullptr ) {
|
|
StrData<MbeProgBarId_PointsSurface> strData( pbarId_Solve_LinearEquationsSystem );
|
|
progBar = &baseProgBar->CreateChildAddRef( strData );
|
|
}
|
|
|
|
ptrdiff_t count = (ptrdiff_t)std_min( a.Lines(), a.Columns() );
|
|
C3D_ASSERT( a.Lines() == a.Columns() );
|
|
|
|
if ( count < 1 ) {
|
|
::FinishProgressBar( progBar );
|
|
::ReleaseItem( progBar );
|
|
return nr_Special;
|
|
}
|
|
|
|
ptrdiff_t i, j, k;
|
|
double m, tmp;
|
|
Type parr;
|
|
|
|
ptrdiff_t halfCount = count / 2;
|
|
|
|
for ( k = 0; k < count - 1; k++ ) {
|
|
// \ru Переставить уравнения так, чтобы a[k][k] != 0 \en Swap equations so that a[k][k] != 0
|
|
ptrdiff_t l = k;
|
|
tmp = ::fabs( a(k, k) );
|
|
for ( i = k + 1; i < count; i++ ) {
|
|
if ( ::fabs( a(i, k) ) > tmp ) {
|
|
l = i;
|
|
tmp = ::fabs( a(i, k) );
|
|
}
|
|
}
|
|
|
|
if ( k == halfCount )
|
|
::SetProgressBarValue( progBar, 25 );
|
|
if ( ::StopProgressBar( progBar ) ) { // \ru Остановка по запросу \en Stop by request
|
|
::ReleaseItem( progBar );
|
|
return nr_Special;
|
|
}
|
|
|
|
if ( l != k ) {
|
|
// (!) Рекомендуется использовать MatrixNN::SwapLines(k,l)
|
|
for ( j = k; j < count; j++ )
|
|
{
|
|
tmp = a(k, j);
|
|
a.SetElem( k, j, a(l, j) );
|
|
a.SetElem( l, j, tmp );
|
|
}
|
|
parr = b[k];
|
|
b[k] = b[l];
|
|
b[l] = parr;
|
|
}
|
|
|
|
if ( ::StopProgressBar( progBar ) ) { // \ru Остановка по запросу \en Stop by request
|
|
::ReleaseItem( progBar );
|
|
return nr_Special;
|
|
}
|
|
|
|
if ( ::fabs( a(k, k) ) < epsilon ) {
|
|
::FinishProgressBar( progBar );
|
|
::ReleaseItem( progBar );
|
|
return nr_Special; // \ru Система не имеет решений \en System doesn't have solutions
|
|
}
|
|
|
|
tmp = 1 / a(k, k);
|
|
|
|
for ( i = k + 1; i < count; i++ )
|
|
{
|
|
if ( a(i, k) != 0.0 ) {
|
|
m = a(i, k) * tmp;
|
|
a.SetElem( i, k, 0.0 );
|
|
for ( j = k + 1; j < count; j++ )
|
|
a.SetElem( i, j, a(i, j) - a(k, j) * m );
|
|
b[i] -= b[k] * m; // SKIP_SA
|
|
}
|
|
}
|
|
|
|
if ( ::StopProgressBar( progBar ) ) { // \ru Остановка по запросу \en Stop by request
|
|
::ReleaseItem( progBar );
|
|
return nr_Special;
|
|
}
|
|
}
|
|
|
|
::SetProgressBarValue( progBar, 50 );
|
|
|
|
if ( ::fabs( a(k, k) ) < epsilon ) {
|
|
::FinishProgressBar( progBar );
|
|
::ReleaseItem( progBar );
|
|
return nr_Special; // \ru Система не имеет решений \en System doesn't have solutions
|
|
}
|
|
|
|
// \ru Обратная подстановка \en Back-substitution
|
|
i = count - 1;
|
|
b[i] *= 1 / a(i, i);
|
|
|
|
for ( i = count - 2; i >= 0; i-- ) {
|
|
j = i + 1;
|
|
parr = b[j] * a(i, j);
|
|
for ( j = i + 2; j < count; j++ )
|
|
parr += b[j] * a(i, j);
|
|
m = 1.0 / a(i, i);
|
|
b[i] = ( b[i] - parr ) * m;
|
|
|
|
if ( i == halfCount )
|
|
::SetProgressBarValue( progBar, 75 );
|
|
if ( ::StopProgressBar( progBar ) ) { // \ru Остановка по запросу \en Stop by request
|
|
::ReleaseItem( progBar );
|
|
return nr_Special;
|
|
}
|
|
}
|
|
|
|
::SetProgressBarValue( progBar, 100 );
|
|
::FinishProgressBar( progBar );
|
|
::ReleaseItem( progBar );
|
|
return nr_Success;
|
|
} // GaussEquation
|
|
|
|
|
|
//------------------------------------------------------------------------------
|
|
/** \brief \ru Решение системы линейных уравнений методом исключения Гаусса.
|
|
\en System of linear equations is solved by the Gauss method. \~
|
|
\details \ru Решение системы линейных уравнений методом исключения Гаусса. \n
|
|
\en System of linear equations is solved by the Gauss method. \n \~
|
|
\param[in] a - \ru Матрица коэффициентов при неизвестных
|
|
\en Coefficient matrix \~
|
|
\param[in] b - \ru Массив правых частей, в него же помещается результат решения
|
|
\en Array of right parts, on output it contains the result of the solution \~
|
|
\param[in] epsilon - \ru Погрешность решения
|
|
\en Tolerance of solution \~
|
|
\param[in] baseProgBar - \ru Индикатор процесса решения
|
|
\en Progress indicator of solution \~
|
|
\return \ru Код ошибки: если nr_Success (+1), то система решена, если nr_Special, то нет решений или система вырождена.
|
|
\en Error code: if nr_Success (+1), then the system is solved, if nr_Special, there is no solution or the system is degenerate. \~
|
|
\ingroup Base_Items
|
|
*/
|
|
// ---
|
|
template <class Type, class SparseMatrix>
|
|
MbeNewtonResult TypedGaussEquationWithBandMatrix( SparseMatrix & a, Type * b, double epsilon, ProgressBarWrapper * baseProgBar = nullptr )
|
|
{
|
|
ProgressBarWrapper * progBar = nullptr;
|
|
if ( baseProgBar != nullptr ) {
|
|
StrData<MbeProgBarId_PointsSurface> strData( pbarId_Solve_LinearEquationsSystem );
|
|
progBar = &baseProgBar->CreateChildAddRef( strData );
|
|
}
|
|
|
|
ptrdiff_t count = (ptrdiff_t)std_min( a.Lines(), a.Columns() );
|
|
C3D_ASSERT( a.Lines() == a.Columns() );
|
|
|
|
if ( count < 1 ) {
|
|
::FinishProgressBar( progBar );
|
|
::ReleaseItem( progBar );
|
|
return nr_Special;
|
|
}
|
|
|
|
ptrdiff_t i, j, k;
|
|
|
|
ptrdiff_t halfCount = count / 2;
|
|
|
|
std::vector<size_t> nzIndices;
|
|
c3d::IndicesPair searchRange;
|
|
|
|
a.NzUpdate();
|
|
|
|
for ( k = 0; k < count - 1; k++ ) {
|
|
// \ru Переставить уравнения так, чтобы a[k][k] != 0 \en Swap equations so that a[k][k] != 0
|
|
ptrdiff_t l = k;
|
|
double tmp = ::fabs( a.GetElem( k, k ) );
|
|
for ( i = k + 1; i < count; i++ ) {
|
|
double faik = ::fabs( a.GetElem( i, k ) );
|
|
if ( faik > tmp ) {
|
|
l = i;
|
|
tmp = faik;
|
|
}
|
|
}
|
|
|
|
if ( k == halfCount )
|
|
::SetProgressBarValue( progBar, 25 );
|
|
if ( ::StopProgressBar( progBar ) ) { // \ru Остановка по запросу \en Stop by request
|
|
::ReleaseItem( progBar );
|
|
return nr_Special;
|
|
}
|
|
|
|
if ( l != k ) {
|
|
a.SwapLines( k, l );
|
|
std::swap( b[k], b[l] );
|
|
}
|
|
|
|
if ( ::StopProgressBar( progBar ) ) { // \ru Остановка по запросу \en Stop by request
|
|
::ReleaseItem( progBar );
|
|
return nr_Special;
|
|
}
|
|
|
|
tmp = a.GetElem(k, k);
|
|
|
|
if ( ::fabs( tmp ) < epsilon ) {
|
|
::FinishProgressBar( progBar );
|
|
::ReleaseItem( progBar );
|
|
return nr_Special; // \ru Система не имеет решений \en System doesn't have solutions
|
|
}
|
|
|
|
tmp = 1.0 / tmp;
|
|
|
|
for ( i = k + 1; i < count; i++ ) {
|
|
double aik = a.GetElem( i, k );
|
|
if ( ::fabs(aik) > NULL_EPSILON ) {
|
|
double m = aik * tmp;
|
|
a.SetElem( i, k, 0.0 );
|
|
aik = 0.0;
|
|
|
|
nzIndices.clear();
|
|
if ( k + 1 < count ) {
|
|
// [k+1,count)
|
|
searchRange.first = (k + 1);
|
|
searchRange.second = (count - 1);
|
|
a.NzIndices( k, searchRange, nzIndices );
|
|
}
|
|
size_t nzCnt = nzIndices.size();
|
|
|
|
for ( size_t nzInd = 0; nzInd < nzCnt; ++nzInd ) {
|
|
j = nzIndices[nzInd];
|
|
double ae = a.GetElem( k, j ) * m;
|
|
if ( ::fabs(ae) > NULL_EPSILON ) {
|
|
ae = a.GetElem( i, j ) - ae;
|
|
a.SetElem( i, j, ae );
|
|
}
|
|
}
|
|
b[i] -= b[k] * m;
|
|
}
|
|
}
|
|
|
|
if ( ::StopProgressBar( progBar ) ) { // \ru Остановка по запросу \en Stop by request
|
|
::ReleaseItem( progBar );
|
|
return nr_Special;
|
|
}
|
|
}
|
|
|
|
a.NzUpdate();
|
|
|
|
::SetProgressBarValue( progBar, 50 );
|
|
|
|
if ( ::fabs( a.GetElem(k, k) ) < epsilon ) {
|
|
::FinishProgressBar( progBar );
|
|
::ReleaseItem( progBar );
|
|
return nr_Special; // \ru Система не имеет решений \en System doesn't have solutions
|
|
}
|
|
|
|
// \ru Обратная подстановка \en Back-substitution
|
|
i = count - 1;
|
|
b[i] *= 1.0 / a.GetElem(i, i);
|
|
|
|
Type be;
|
|
for ( i = count - 2; i >= 0; i-- ) {
|
|
j = i + 1;
|
|
be = b[j] * a.GetElem( i, j );
|
|
|
|
nzIndices.clear();
|
|
if ( i + 2 < count ) {
|
|
// [i+2,count)
|
|
searchRange.first = (i + 2);
|
|
searchRange.second = (count - 1);
|
|
a.NzIndices( i, searchRange, nzIndices );
|
|
}
|
|
size_t nzCnt = nzIndices.size();
|
|
|
|
for ( size_t nzInd = 0; nzInd < nzCnt; ++nzInd ) {
|
|
j = nzIndices[nzInd];
|
|
be += b[j] * a.GetElem( i, j );
|
|
}
|
|
double m = 1.0 / a.GetElem( i, i );
|
|
b[i] = ( b[i] - be ) * m;
|
|
|
|
if ( i == halfCount )
|
|
::SetProgressBarValue( progBar, 75 );
|
|
if ( ::StopProgressBar( progBar ) ) { // \ru Остановка по запросу \en Stop by request
|
|
::ReleaseItem( progBar );
|
|
return nr_Special;
|
|
}
|
|
}
|
|
|
|
::SetProgressBarValue( progBar, 100 );
|
|
::FinishProgressBar( progBar );
|
|
::ReleaseItem( progBar );
|
|
return nr_Success;
|
|
} // GaussEquation
|
|
|
|
|
|
//------------------------------------------------------------------------------
|
|
/**
|
|
\brief \ru Решение системы линейных уравнений с трехдиагональной матрицей методом прогонки.
|
|
\en System of linear equations with a tridiagonal matrix is solved by the sweep method. \~
|
|
\details \ru Решение системы линейных уравнений с трехдиагональной матрицей методом прогонки.
|
|
Используется для построения NURBS-копии незамкнутого кубического сплайна.
|
|
\en System of linear equations with a tridiagonal matrix is solved by the sweep method.
|
|
Used to build a NURBS-copy of a non-closed cubic spline. \~
|
|
\param[in] n - \ru число неизвестных
|
|
\en The number of unknown variables \~
|
|
\param[in] a - \ru Главная диагональ трехдиагональной матрицы, массив double размерности n
|
|
\en Main diagonal of a tridiagonal matrix is an array of doubles of size "n" \~
|
|
\param[in] b - \ru Верхняя диагональ, размерность n-1
|
|
\en Upper diagonal, dimension is n-1 \~
|
|
\param[in] c - \ru <b>c</b> - нижняя диагональ, размерность n-1
|
|
\en <b>c</b> - lower diagonal, dimension is n-1 \~
|
|
\param[in] r - \ru Вектор правой части, массив точек или векторов размерности n; должна быть определена операция умножения на double справа
|
|
\en Vector of the right part is an array of points or vectors of dimension n; the multiplication operation by the double value on the right must be defined \~
|
|
\param[in] solution - \ru Массив решений (точек, векторов), размерности n
|
|
\en Array of solutions (points, vectors) of dimension n \~
|
|
\param[in] epsZero - \ru Погрешность нуля
|
|
\en Tolerance of zero \~
|
|
\return \ru Код ошибки: если nr_Success (+1), то система решена, если nr_Special, то нет решений или система вырождена.
|
|
\en Error code: if nr_Success (+1), then the system is solved, if nr_Special, there is no solution or the system is degenerate. \~
|
|
\ingroup Base_Items
|
|
*/
|
|
// ---
|
|
template <class ArrayDouble, class ArrayType>
|
|
MbeNewtonResult TypedTridiagonalSolve ( const size_t n, ArrayDouble & a, ArrayDouble & b, ArrayDouble & c,
|
|
ArrayType & r, ArrayType & solution, double epsZero )
|
|
{
|
|
if ( ::fabs ( a[0] ) < epsZero )
|
|
return nr_Special; // \ru Прогонка не работает \en Sweep does not work
|
|
// \ru Прямой ход прогонки \en Forward step of sweep
|
|
b[0] /= - a[0];
|
|
r[0] /= a[0];
|
|
for ( size_t i = 1; i < n; ++i ) {
|
|
double dom = a[i] + b[i - 1] * c[i - 1];
|
|
if ( ::fabs ( dom ) < epsZero )
|
|
return nr_Special; // \ru Прогонка не работает \en Sweep does not work
|
|
if ( i < n - 1 ) // \ru Размерность b равна n-1, поэтому не можем вычислять b[n-1], да и не нужно оно. \en B dimension is equal to n-1, therefore we can not calculate b [n-1], and it isn't necessary.
|
|
b[i] /= - dom;
|
|
r[i] = ( r[i] - r[i - 1] * c[i - 1] ) / dom;
|
|
}
|
|
// \ru Обратный ход \en Backward substitution
|
|
solution[n - 1] = r[n - 1];
|
|
for ( ptrdiff_t i = n - 2; i >= 0; --i )
|
|
solution[i] = r[i] + solution[i + 1] * b[i];
|
|
return nr_Success;
|
|
}
|
|
|
|
|
|
#endif // __MB_MATRIXNN_H
|