Класс Matrix.


class Matrix
{
double ** p;
// Запрет оператора присваивания
void operator = ( const Matrix & );
public:
const int nRow; // к-во строк
const int nCol; // к-во столбцов

Matrix ( int r, int c );
Matrix ( int r, int c, const double * const * a );
Matrix ( const Matrix & );
~Matrix ();

void fill ( double );
Matrix & operator *= ( double );
operator double * const * () { return p; }
operator const double * const * () const { return p; }
};

Оператор присваивания сделан недоступным, т.к. члены nRow и nCol - константы и не могут быть изменены. Имеются три конструктора. Первый из них просто выделяет память для элементов матрицы:

Matrix::Matrix ( int r, int c ) :
nRow ( r > 0 && c > 0 ? r : 0 ),
nCol ( r > 0 && c > 0 ? c : 0 )
{
if ( nRow == 0 )
{
p = 0;
return;
}
p = new double*[nRow];
*p = new double[nRow*nCol];
for ( int i = 1; i < nRow; ++i ) p[i] = p[i-1] + nCol;
}

Обратите внимание - память под элементы выделяется непрерывным куском и это будет в дальнейшем использоваться. Второй конструктор выделяет память и инициализирует элементы:

Matrix::Matrix ( int r, int c, const double * const * a ) :
nRow ( r > 0 && c > 0 && a != 0 ? r : 0 ),
nCol ( r > 0 && c > 0 && a != 0 ? c : 0 )
{
if ( nRow == 0 )
{
p = 0;
return;
}
p = new double*[nRow];
*p = new double[nRow*nCol];
for ( int i = 0; i < nRow; ++i )
{
if ( i ) p[i] = p[i-1] + nCol;
double * t = p[i];
const double * s = a[i];
for ( int j = 0; j < nCol; ++j ) t[j] = s[j];
}
}

Третий конструктор - это конструктор копии:

Matrix::Matrix ( const Matrix & m ) : nRow ( m.nRow ), nCol ( m.nCol )
{
if ( nRow == 0 )
{
p = 0;
return;
}
p = new double*[nRow];
*p = new double[nRow*nCol];
for ( int i = 0; i < nRow; ++i )
{
if ( i ) p[i] = p[i-1] + nCol;
double * t = p[i];
const double * s = m[i];
for ( int j = 0; j < nCol; ++j ) t[j] = s[j];
}
}

Деструктор освобождает выделенную память:

Matrix::~Matrix()
{
if ( p )
{
delete[] *p;
delete[] p;
}
}

Функция fill заполняет матрицу указанным значением и здесь используется то, что память для элементов матрицы непрерывна:

void Matrix::fill ( double d )
{
if ( p == 0 ) return;
double * t = *p;
for ( int i = nRow*nCol; --i >= 0; ) *t++ = d;
}

Оператор *= умножает матрицу на число:

Matrix & Matrix::operator *= ( double d )
{
if ( p )
{
double * t = *p;
for ( int i = nRow*nCol; --i >= 0; ) *t++ *= d;
}
return *this;
}

Есть четыре матричные нормы ( 1, 2, Фробениуса и бесконечная ):

double norm1 ( const Matrix & A ) // 1-норма
{
if ( A.nCol == 0 ) return 0.;
double res = 0.;
for ( int j = 0; j < A.nCol; ++j )
{
double t = 0.;
for ( int i = 0; i < A.nRow; ++i ) t += fabs ( A[i][j] );
if ( res < t ) res = t;
}
return res;
}

double norm2 ( const Matrix & A ) // 2-норма
{
if ( A.nCol == 0 ) return 0.;
Matrix T ( A.nCol, A.nRow );
trans ( A, T );
Matrix S ( A.nCol, A.nCol );
multi ( T, A, S );
Matrix V ( A.nCol, A.nCol );
double * d = T[0];
jacobi ( S.nCol, S, d, V );
double res = 0.;
for ( int i = S.nCol; --i >= 0; ) if ( res < d[i] ) res = d[i];
return sqrt(res);
}

double normF ( const Matrix & A ) // норма Фробениуса
{
if ( A.nCol == 0 ) return 0.;
double res = 0.;
const double * t = A[0];
for ( int i = A.nRow*A.nCol; --i >= 0; )
{
res += *t * *t;
++t;
}
return sqrt(res);
}

double normU ( const Matrix & A ) // бесконечная норма
{
if ( A.nCol == 0 ) return 0.;
double res = 0.;
for ( int i = 0; i < A.nRow; ++i )
{
double t = 0.;
for ( int j = 0; j < A.nCol; ++j ) t += fabs ( A[i][j] );
if ( res < t ) res = t;
}
return res;
}

Описание функции jacobi находится здесь.

Далее следует группа функций у которых есть параметры типа Matrix. Все они проверяют входные данные и возвращают соответствующее логическое значение.

Функция вычисляющая детерминант использует метод Гаусса для построения треугольной матрицы для которой детерминант равен произведению диагональных элементов с учётом знака. Т.к. детерминант определён только для квадратных матриц в начале функции есть соответствующая проверка:

bool determinant ( const Matrix & m, double & d )
{
if ( m.nCol != m.nRow || m.nCol == 0 ) return false;
SLU_Gauss slu ( m.nCol, m );
d = slu.determinant();
return true;
}

Функция copy предназначена вместо оператора присваивания. Здесь также используется то, что память для элементов матрицы непрерывна:

bool copy ( const Matrix & a, Matrix & b )
{
if ( a.nCol != b.nCol || a.nRow != b.nRow || a.nCol == 0 ) return false;
if ( a[0] == b[0] ) return true; // a и b - одна и та же матрица
const double * pa = *a;
double * pb = *b;
for ( int i = b.nRow*b.nCol; --i >= 0; ) *pb++ = *pa++;
return true;
}

Функция trans осуществляет транспонирование матрицы. Вначале проверяется совместны ли размерности матриц a и b. Далее, в зависимости от того одну ли и ту же матрицу представляют a и b, происходит один из двух вариантов транспонирования. Т.е. можно писать trans ( a, a ):

bool trans ( const Matrix & a, Matrix & b )
{
if ( a.nCol != b.nRow || a.nRow != b.nCol || a.nCol == 0 ) return false;
if ( a[0] == b[0] ) // a и b - одна и та же матрица
{
for ( int i = b.nRow; --i > 0; )
{
double * p = b[i];
for ( int j = i; --j >= 0; )
{
double t = p[j];
p[j] = b[j][i];
b[j][i] = t;
}
}
}
else
{
for ( int i = a.nRow; --i >= 0; )
{
const double * p = a[i];
for ( int j = a.nCol; --j >= 0; )
{
b[j][i] = p[j];
}
}
}
return true;
}

Функция inverse порождает обратную матрицу. Матрицы a и b должны быть квадратными и иметь одинаковую размерность. Здесь также возможна запись inverse ( a, a ), и в случае неудачи матрица a не изменится:

bool inverse ( const Matrix & a, Matrix & b ) // a * b = 1
{
if ( a.nCol != b.nCol || a.nRow != b.nRow ) return false;
if ( a.nCol != a.nRow || a.nRow == 0 ) return false;
SLU_Gauss slu ( a.nRow, a );
Matrix m ( b );
double * c = m[0];
for ( int i = m.nRow; --i >= 0; ) c[i] = 0.;
for ( i = m.nRow; --i >= 0; )
{
c[i] = 1.;
if ( ! slu.solve ( c, m[i] ) ) return false;
if ( i ) c[i] = 0.;
}
return trans ( m, b );
}

Следующие три функции осуществляют операции сложения, разности и произведения матриц:

bool plus ( const Matrix & a, const Matrix & b, Matrix & c ) // c = a + b
{
if ( a.nCol != b.nCol || a.nRow != b.nRow ) return false;
if ( a.nCol != c.nCol || a.nRow != c.nRow || a.nRow == 0 ) return false;
const double * pa = *a;
const double * pb = *b;
double * pc = *c;
for ( int i = a.nRow*a.nCol; --i >= 0; ) *pc++ = *pa++ + *pb++;
return true;
}

bool minus ( const Matrix & a, const Matrix & b, Matrix & c ) // c = a - b
{
if ( a.nCol != b.nCol || a.nRow != b.nRow ) return false;
if ( a.nCol != c.nCol || a.nRow != c.nRow || a.nRow == 0 ) return false;
const double * pa = *a;
const double * pb = *b;
double * pc = *c;
for ( int i = a.nRow*a.nCol; --i >= 0; ) *pc++ = *pa++ - *pb++;
return true;
}

bool multi ( const Matrix & a, const Matrix & b, Matrix & c ) // c = a * b
{
if ( a.nCol != b.nRow || a.nRow != c.nRow || b.nCol != c.nCol || a.nRow == 0 ) return false;
if ( a[0] == c[0] || b[0] == c[0] )
{
Matrix m ( c.nRow, c.nCol );
return multi ( a, b, m ) && copy ( m, c );
}
for ( int i = c.nRow; --i >= 0; )
{
const double * p = a[i];
for ( int j = c.nCol; --j >= 0; )
{
double s = 0.;
for ( int k = a.nCol; --k >= 0; ) s += p[k] * b[k][j];
c[i][j] = s;
}
}
return true;
}

Функция svd осуществляет сингулярное разложение матрицы A в виде: A = U*W*V, где U и V - ортогональные матрицы, a W имеет ненулевые элементы только на главной диагонале.

bool svd ( const Matrix & A, Matrix & U, Matrix & W, Matrix & V );

Шаблон Matr предназначен для осуществления арифметических операций с матрицами на уровне операторов. В этом случае контроль размерностей матриц происходит на этапе компиляции:

template < int r, int c >
class Matr : public Matrix
{
public:
Matr () : Matrix ( r, c ) {}
Matr ( const Matr & m ) : Matrix ( m ) {}

Matr & operator = ( const Matr & m )
{
copy ( m, *this );
return *this;
}

Matr & operator += ( const Matr & m )
{
plus ( *this, m, *this );
return *this;
}

Matr & operator -= ( const Matr & m )
{
plus ( *this, m, *this );
return *this;
}

Matr & operator *= ( double d )
{
Matrix::operator *= ( d );
return *this;
}
};

template < int m, int n >
const Matr < m, n > operator + ( const Matr < m, n > & a, const Matr < m, n > & b )
{
Matr < m, n > c;
plus ( a, b, c );
return c;
}

template < int m, int n >
const Matr < m, n > operator - ( const Matr < m, n > & a, const Matr < m, n > & b )
{
Matr < m, n > c;
minus ( a, b, c );
return c;
}

template < int m, int n >
const Matr < n, m > operator * ( const Matr < n, m > & a, double d )
{
Matr < n, m > b = a;
b *= d;
return b;
}

template < int m, int n >
const Matr < n, m > operator * ( double d, const Matr < n, m > & a )
{
Matr < n, m > b = a;
b *= d;
return b;
}

template < int l, int m, int n >
const Matr < l, n > operator * ( const Matr < l, m > & a, const Matr < m, n > & b )
{
Matr < l, n > c;
multi ( a, b, c );
return c;
}

template < int m, int n >
const Matr < n, m > trans ( const Matr < m, n > & a )
{
Matr < n, m > b;
trans ( a, b );
return b;
}

template < int m, int n >
const Matr < n, m > inverse ( const Matr < m, n > & a )
{
Matr < n, m > b;
inverse ( a, b );
return b;
}

Есть проблема с умножением матрицы на матрицу. Мой компилятор VC6 почему-то не находит подходящую функцию при записи a * b, но мне сообщили, что в седьмой версии это проходит.

Исходники находятся в mathem.zip.