По приведенным описаниям можно составить, к примеру, такой интерфейс для данного класса:
#ifndef __MATRIX_H
#define __MATRIX_H
#ifndef __VECTOR_H
#include "vector.h"
#endif
#ifndef __FSTREAM_H
#include <fstream.h>
#endif
#ifndef __MATH_H
#include <math.h>
#endif
#ifndef __STDIO_H
#include <stdlib.h>
#endif
#ifndef __EXCEPT_H
#include <except.h>
#endif
#ifndef __CSTRING_H
#include <cstring.h>
#endif
#ifndef __COMPLEX_H
#include <complex.h
#endif
/*параметризованный класс для работы с матричными объектами*/
template <class YourOwnFloatType> class matrix
{ //приватные данные
long m,n;/*row, columns - размерность матрицы в строках и столбцах */
vector<YourOwnFloatType> *mtr; /*указатель на вектор данных*/
public://общедоступные члены
matrix(char *);/*загрузка матрицы из файла в формате m n d11 d12 ... dmn */
matrix();//пустая матрица размером 1х1
matrix(long, long);/*пустая матрица заданного размера*/
matrix(long, long, YourOwnFloatType *);/*матрица size1xsize2 из массива */
matrix(matrix &);//конструктор копирования
~matrix();//деструктор
friend matrix<YourOwnFloatType> operator+(matrix <YourOwnFloatType> &, matrix<YourOwnFloatType> &);//сложение
friend matrix<YourOwnFloatType> operator-(matrix <YourOwnFloatType> &, matrix<YourOwnFloatType> &);//вычитание
friend matrix<YourOwnFloatType> operator*(matrix <YourOwnFloatType> &, matrix<YourOwnFloatType> &);//умножение матриц
friend matrix<YourOwnFloatType> operator*( YourOwnFloatType, matrix<YourOwnFloatType> &);//умножение числа на матрицу
friend matrix<YourOwnFloatType> operator*(matrix <YourOwnFloatType> &, YourOwnFloatType );/*умножение матрицы на число */
friend ostream &operator<<(ostream &,matrix <YourOwnFloatType> &);//вывод матрицы в поток
friend istream &operator>>(istream &,matrix <YourOwnFloatType> &);//ввод матрицы из потока
/*первый параметр - квадратная матрица NxN, второй - Nx1 (матричное уравнение) */
friend matrix<YourOwnFloatType> SLAE_Orto(matrix <YourOwnFloatType> &, matrix<YourOwnFloatType> &);//метод ортогонализации
friend matrix<YourOwnFloatType> SLAE_Gauss(matrix <YourOwnFloatType> &, matrix<YourOwnFloatType> &);//метод Гаусса с выбором главного элемента
friend YourOwnFloatType det2(matrix <YourOwnFloatType> &); //аналитический детерминант
friend YourOwnFloatType det(matrix <YourOwnFloatType> &); //численный детерминант
matrix<YourOwnFloatType> minor(long, long);
/*минор - создание матрицы из имеющейся без заданных строки и столбца */
matrix<YourOwnFloatType> operator=(matrix <YourOwnFloatType> &);//присвоение
matrix<YourOwnFloatType> operator*=(matrix <YourOwnFloatType> &x);/*набор сокращённых операций умножения,*/
matrix<YourOwnFloatType> operator+=(matrix <YourOwnFloatType> &x);//сложения
matrix<YourOwnFloatType> operator-=(matrix <YourOwnFloatType> &x); //и вычитания
matrix<YourOwnFloatType> operator^(long);/*степень матрицы как операция */
matrix<YourOwnFloatType> operator^=(long); /*сокращённая степень */
friend matrix<YourOwnFloatType>pow(matrix <YourOwnFloatType> &,long);
//степень матрицы как дружественная функция
matrix<YourOwnFloatType> operator~(); //транспонирование
matrix<YourOwnFloatType> operator!(); /*аналитическое обращение матрицы */
matrix<YourOwnFloatType> operator*(); /*численное обращение матрицы */
YourOwnFloatType operator&();/*численный детерминант унарная операция */
operator YourOwnFloatType();
/*быстрый детерминант как оператор преобразования к типу-параметру*/
matrix<YourOwnFloatType> operator-();/*унарный минус*/
matrix<YourOwnFloatType> operator+();/*унарный плюс*/
friend long operator==(matrix<YourOwnFloatType> &, matrix<YourOwnFloatType> &);/*проверка на равенство*/
friend long operator!=(matrix<YourOwnFloatType> &, matrix<YourOwnFloatType> &);/*проверка на неравенство*/
vector<YourOwnFloatType> &operator[](long a); //индексирование матрицы
long getm() { return m; }//число строк
long getn() { return n; }//число столбцов
};
/*Вспомогательные функции - функция, определяющая знак числа */
inline long sign(long double x)
{
return (x<0) ? -1 : 1;
}
double fabs(complex x) //
{
return abs(x);
}
extern const long double MAX_LONGDOUBLE;
/*
Напомним, что матрица - это тоже вектор, элементами которого являются не числовые объекты, а арифметические вектора. В связи с этим между векторным и матричным классами есть много общего, однако иногда можно наблюдать и существенные отличия. Возьмём, скажем, классический файловый метод хранения данных. Если для вектора достаточно было указать одно служебное число - его длину (размерность), то для матрицы их требуется уже два, причём первое из этих чисел будет определять количество векторов в матрице, а второе - размерность каждого из этих векторов. В файловых операциях с векторами нам поможет определённая в векторном классе операция ввода вектора с некоторого устройства */
template <class YourOwnFloatType> matrix <YourOwnFloatType>::matrix(char *f) /*имя файла с матрицей*/
{
long i;
ifstream fp=f;//пытаемся открыть файл
if(!fp)
throw xmsg("Не могу открыть файл "+string(f)+ "\n");
fp>>m>>n;//размеры матрицы считываем из файла
if(m<=0||n<=0)
throw xmsg("Некорректный размер матрицы \n");
try
{
//здесь работает конструктор без параметров
mtr=new vector<YourOwnFloatType>[m];/*попытка выделения памяти */
}
catch(xalloc)
{
throw xmsg("Не хватает памяти \n");
}
for(i=0;i<m;i++)/*и только здесь вектора расширяется до нужной размерности */
mtr[i]=vector<YourOwnFloatType>(n);
for(i=0;i<m&&fp>>mtr[i];i++);/*при вводе используем перегруженную операцию класса vector*/
}
/*Зная только размеры матрицы, мы можем считать её нулевым элементом данного размера и сконструировать соответствующим образом: */
template <class YourOwnFloatType> matrix <YourOwnFloatType>::matrix(long a, long b): m(a),n(b) //число строк и число столбцов
{
long i;
if(m<=0||n<=0)
throw xmsg("Некорректный размер матрицы \n");
try
{
//здесь работает конструктор без параметров
mtr=new vector<YourOwnFloatType>[m];/*попытка выделения памяти */
}
catch(xalloc)
{
throw xmsg("Не хватает памяти \n");
}
for(i=0;i<m;i++)/*расширяем вектора до нужного размера*/
mtr[i]=vector<YourOwnFloatType>(n);
}
/*Получив о матрице все возможные сведения, имеет смысл, сконструировав её, инициализировать соответствующими элементами: */
template <class YourOwnFloatType>matrix <YourOwnFloatType>::matrix(long a, long b, YourOwnFloatType *mt):m(a),n(b)
{
if(m<=0||n<=0)
throw xmsg("Некорректный размер матрицы \n");
try
{
//здесь работает конструктор без параметров
mtr=new vector<YourOwnFloatType>[m];/*попытка выделения памяти*/
}
catch(xalloc)
{
throw xmsg("Не хватает памяти \n");
}
for(long i=0;i<m;i++)
mtr[i]=vector<YourOwnFloatType>(n);//расширение
for(long i=0,l=0;i<m;i++)
for(long j=0;j<n;j++)
mtr[i][j]=mt[l++];/*сконструировав, копируем данные из массива в матрицу*/
}
/*Обратный случай (когда мы не знаем ни размеров, ни элементов матричных векторов) решается достаточно просто - созданием матрицы из одного единственного вектора длиною в один элемент: */
template <class YourOwnFloatType>
matrix<YourOwnFloatType>::matrix():m(1),n(1)
{
try
{
//здесь работает конструктор без параметров
mtr=new vector<YourOwnFloatType>[m];/*попытка выделения памяти, а обнулится она вектором*/
}
catch(xalloc)
{
throw xmsg("Не хватает памяти \n");
}
}
/*Деструктор: уничтожение матрицы автоматически уничтожает все её векторы: */
template <class YourOwnFloatType>
matrix<YourOwnFloatType>::~matrix()
{
delete []mtr;/*вызов деструкторов для всех векторов*/
}
/* Как и векторы, матрицы тоже бывает необходимым проиндексировать; результатом этого действия, естественно, будет соответствующий вектор: */
template <class YourOwnFloatType> vector <YourOwnFloatType> &matrix<YourOwnFloatType>:: operator[](long a)
{
static vector<YourOwnFloatType> error;
error[0]=MAX_LONGDOUBLE;
if(a>=0&&a<m)
return mtr[a];/*а дальше можно индексировать вектор*/
else
{
cerr<<"Индекс "<<a<<" вне матрицы \n";
return error;
}
}
/*Конструктор копирования. Если в памяти ЭВМ уже есть какая-то матрица, с неё можно снять копию-слепок: */
template <class YourOwnFloatType> matrix <YourOwnFloatType>::matrix(matrix<YourOwnFloatType> &ex): m(ex.m), n(ex.n)
{
try
{
//здесь работает конструктор без параметров
mtr=new vector<YourOwnFloatType>[m];/*попытка выделения памяти*/
}
catch(xalloc)
{
throw xmsg("Не хватает памяти \n");
}
for(long i=0;i<m;i++)
mtr[i]=ex[i];/*очень громоздкое индексирование и присваивание векторов*/
}
/* Сложение матриц одинаковой размерности сводится к сложению соответствующих векторов: */
template <class YourOwnFloatType> matrix <YourOwnFloatType> operator+(matrix <YourOwnFloatType> &f, matrix<YourOwnFloatType> &s)
{
if(f.m!=s.m||f.n!=s.n)
throw xmsg("Размерности слагаемых матриц не совпадают \n");
matrix<YourOwnFloatType> temp(f.m, f.n); /*временная матрица*/
for(long i=0;i<f.m;i++)
temp[i]=f[i]+s[i];//слагаем вектора матриц
return temp;//возвращаем результат
}
/*
Как и для векторов, определим унарный "минус":
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> matrix<YourOwnFloatType>:: operator-()
{
matrix<YourOwnFloatType> temp(m, n);
for(long i=0;i<m;i++)
temp[i]=-(*this)[i];
return temp;//возвращаем результат
}
//унарный плюс
template <class YourOwnFloatType> matrix <YourOwnFloatType> matrix<YourOwnFloatType>:: operator+()
{
return *this; /*возвращаем себя без каких-либо изменений*/
}
/*
В отличие от скалярного произведения векторов, умножение матриц является бинарной алгебраической операцией, однако только для отдельных видов матриц, к тому же эта операция не является коммутативной!
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> operator*(matrix <YourOwnFloatType> &f, matrix<YourOwnFloatType> &s)
{
if(f.n!=s.m)
throw xmsg("Умножение матриц с данными размерами невозможно \n");
matrix<YourOwnFloatType> temp(f.m, s.n);
for(long j=0;j<s.n;j++)
for(long i=0;i<f.m;i++)
for(long k=0;k<f.n;k++)
temp[i][j]+=f[i][k]*s[k][j];
return temp;
}
/*
Операция увеличения матрицы в некоторое число раз:
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> operator*(matrix <YourOwnFloatType> &f, YourOwnFloatType s)
{
matrix<YourOwnFloatType> temp=f;
for(long i=0;i<f.m;i++)
temp[i]=temp[i]*s;/*здесь используем умножение вектора на число*/
return temp;
}
/*
Оператор возведения матрицы в целую степень можно определить так: любая матрица в нулевой степени является единичной, положительная степень определяется через произведение, отрицательная - через обращение матрицы и произведение. Последний случай можно красиво реализовать с использованием рекурсии, как, впрочем, и предыдущий.
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> matrix<YourOwnFloatType>:: operator^(long l)
{
matrix<YourOwnFloatType> temp(getm(),getn());
if(getm()!=getn())
throw xmsg("Для неквадратных матриц степень не определена \n");
if(l==0)
{
for(long i=0;i<getm();i++)
temp[i][i]=1;
return temp;
}
if(l>0)
{
temp=(*this);
for(long i=1;i<l;i++)
temp*=(*this);
return temp;
}
else
return (*(*this))^(-l);
}
//сокращённая операция возведения матрицы в степень
template <class YourOwnFloatType> matrix <YourOwnFloatType> matrix<YourOwnFloatType>:: operator^=(long l)
{
return (*this)=(*this)^l;
}
/*
Степень матрицы в виде функции
*/
template <class YourOwnFloatType> inline matrix <YourOwnFloatType> pow(matrix<YourOwnFloatType> &x, long l)
{
return x^l;
}
/*
Умножение числа на матрицу реализуем через уже имеющуюся операцию умножения матрицы на число для сохранения коммутативности
*/
template <class YourOwnFloatType> inline matrix <YourOwnFloatType> operator*(YourOwnFloatType f, matrix<YourOwnFloatType> &s)
{
return s*f;
}
/*
Вычитание традиционно реализуем через сложение и унарный минус:
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> operator-(matrix <YourOwnFloatType> &f, matrix<YourOwnFloatType> &s)
{
return f+(-s);
}
/*Операция присвоения
При присвоении содержимое матрицы х переписывается в текущую сразу только тогда, когда их размеры совпадают. В противном случае приходится уничтожать текущую матрицу, заново её создавать с новыми размерами и лишь тогда производить повекторное копирование
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> matrix<YourOwnFloatType>::operator= (matrix <YourOwnFloatType> &x)
{
if(m!=x.m||n!=x.n)
{
delete []mtr;//уничтожаем содержимое матрицы
m=x.m, n=x.n;//устанавливаем новые размеры
try
{
mtr=new vector<YourOwnFloatType>[m];/*попытка выделения памяти */
}
catch(xalloc)
{
throw xmsg("Не хватает памяти \n");
}
}
for(long i=0;i<m;i++)
mtr[i]=x[i];
return *this;//возвращение себя
}
/*
Набор сокращённых операций:
умножение
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> matrix<YourOwnFloatType>:: operator*=(matrix<YourOwnFloatType> &x)
{
return (*this)=(*this)*x;
}
/*
сложение
*/
template <class YourOwnFloatType> inline matrix <YourOwnFloatType> matrix<YourOwnFloatType>:: operator+=(matrix<YourOwnFloatType> &x)
{
return (*this)=(*this)+x;
}
/*
вычитание
*/
template <class YourOwnFloatType> inline matrix <YourOwnFloatType> matrix<YourOwnFloatType>::operator-=(matrix<YourOwnFloatType> &x)
{
return (*this)+=-x;/* используем только что сделанное сокращённое сложение*/
}
/*
Очень часто бывает нужна унарная операция транспонирования:
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> matrix<YourOwnFloatType>:: operator~()
{
matrix<YourOwnFloatType> temp(n,m);
for(long j=0;j<n;j++)
for(long i=0;i<m;i++)
temp[j][i]=mtr[i][j];
return temp;//переставлены столбцы и строки
}
/*Сравнение матриц
Сравнивая матрицы, мы сначала учитываем, одинаковой ли они размерности, а потом в случае необходимости сравниваем векторы:
*/
template <class YourOwnFloatType> long operator== (matrix<YourOwnFloatType> &f, matrix<YourOwnFloatType> &s)
{
if(f.m!=s.m||f.n!=s.n)
return 0;
for(long i=0;i<f.m;i++)
if(f[i]!=s[i])
return 0;
return 1;
}
/*
Неравенство
*/
template <class YourOwnFloatType> inline long operator!=(matrix<YourOwnFloatType> &f, matrix<YourOwnFloatType> &s)
{
return !(f==s);//!operator==(f,s);
}
/*Вывод матрицы в поток
Результат иногда хочется посмотреть. Например, так:
*/
template <class YourOwnFloatType> ostream & operator<<(ostream &os,matrix<YourOwnFloatType> &x)
{
for(long i=0;i<x.m;i++)/*построчный вывод векторов матрицы*/
os<<x[i]<<"\n";
return os;
}
/*Или так... (Запись матрицы в файл в градациях серого)*/
template <class YourOwnFloatType>
void writebmp(matrix<YourOwnFloatType> &x, char *name)
{
ofstream f(name, ios::out|ios::binary);/*Битовое изображение */
BITMAPFILEHEADER fh;/*файловый заголовок */
BITMAPINFOHEADER ih;//информационный заголовок
RGBQUAD bmiColors[256];//256 цветов
fh.bfType=0x4d42;//буквы "ВМ"
fh.bfSize=sizeof(fh)+sizeof(ih)+
sizeof(bmiColors)+x.getm()*((x.getn()%4==0)
?x.getn():(x.getn()+4-x.getn()%4));//размер файла
fh.bfReserved1=fh.bfReserved2=0;
fh.bfOffBits=sizeof(fh)+sizeof(ih)+
sizeof(bmiColors);//смещение до данных
f.write((unsigned char*)&fh,sizeof(fh));/*пишем файловый заголовок */
ih.biSize=40;//размер информационного заголовка
ih.biWidth=x.getn();//ширина изображения
ih.biHeight=x.getm();//высота
ih.biPlanes=1;//количество планов изображения
ih.biBitCount=8;//бит на цвет
ih.biCompression=BI_RGB;//без сжатия
ih.biSizeImage=0;//не актуален
ih.biXPelsPerMeter=2963;/*рекомендуемое разрешение по Х*/
ih.biYPelsPerMeter=3158;// по У
ih.biClrUsed=256;//количество используемых цветов
ih.biClrImportant=0;//не актуально
f.write((unsigned char*)&ih,sizeof(ih));/*пишем информационный заголовок*/
for(int i=0;i< (1 << ih.biBitCount);i++)/*цикл по количеству цветов*/
{
/*приравнивая красную, зелёную и голубую составляющие цвета, получаем оттенок серого*/
bmiColors[i].rgbBlue=bmiColors[i].rgbGreen=
bmiColors[i].rgbRed=i;
bmiColors[i].rgbReserved=0;
f.write((unsigned char*)(bmiColors+i),
sizeof(RGBQUAD));//пишем элемент палитры
}
YourOwnFloatType __max=x[0][0], __min=x[0][0];//поиск диапазона значений матрицы
for(i=0;i<x.getm();i++)
for(long j=0;j<x.getn();j++)
{
if(x[i][j]>__max)
__max=x[i][j];
if(x[i][j]<__min)
__min=x[i][j];
}
for(i=x.getm()-1;i>=0;i--)/*строки в ВМР-файле лежат снизу вверх*/
{
for(long j=0;j<x.getn();j++)//записываем строку
{ //формируем номер в палитре оттенков серого
BYTE b=0xff*(x[i][j]-__min)/(__max-__min);
f.write(&b,1);
}
if(x.getn()%4)/*дополняем строку до границы двойного слова*/
for(j=0;j<4-x.getn()%4;j++)
f.write(&j,1);
}
}
/*
Ввод матрицы из потока целиком перекладываем на векторный класс, используя его метод для ввода
*/
template <class YourOwnFloatType> istream & operator>>(istream &is,matrix<YourOwnFloatType> &x)
{
for(long i=0;i<x.m;i++)
is>>x[i];
return is;
}
/*Решение СЛАУ
В матричной форме систему линейных алгебраических уравнений (СЛАУ) можно представить в виде
A*X=B,
где А - матрица коэффициентов при неизвестных,
Х - вектор-столбец этих самих неизвестных, а
В - вектор-столбец свободных членов, взятых с обратными знаками.
Если система имеет решение, то оно будет таким:
А^(-1)*A*X=А^(-1)*B
(А^(-1)*A)*X=А^(-1)*B
Е*X=А^(-1)*B
X=А^(-1)*B,
Существенным моментом является выбор метода решения. Для небольших систем (<200) можно использовать прямые методы (например, метод Гаусса); для реальных расчётов мы рекомендуем метод ортогонализации, в основе которого лежат ортогональные преобразования матриц
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> SLAE_Orto(matrix <YourOwnFloatType> &f, matrix<YourOwnFloatType> &s)
{
matrix<YourOwnFloatType> mtr2(f.m+1,f.m+1), res(f.m,1);
//формируем матрицу из двух для решения СЛАУ
for(long i=0;i<f.m;i++)
for(long j=0;j<f.n;j++)
mtr2[i][j]=f[i][j];
for(long i=0;i<f.m;i++)/*заносим в последнюю строку 0 0 ... 0 1*/
mtr2[i][f.m]=-s[i][0];
mtr2[f.m][f.m]=1;
mtr2[0]=~mtr2[0];//нормируем нулевую строку
/*для улучшения сходимости повторяем этот процесс 3 раза*/
for(long k=0;k<3;k++)
{
for(long l=1;l<f.m+1;l++)
{
for(long i=0;i<l;i++)
{
YourOwnFloatType p=mtr2[l]*mtr2[i];
/*скалярное произведение */
for(long j=0;j<(f.m+1);j++)
mtr2[l][j]-=p*mtr2[i][j];
}
mtr2[l]=~mtr2[l];
}
}
for(long i=0;i<f.m;i++)//переписываем результат
res[i][0]=mtr2[f.m][i]/mtr2[f.m][f.m];
return res;
}
/*
Метод Гаусса с выбором главного элемента
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> SLAE_Gauss (matrix <YourOwnFloatType> &f,matrix<YourOwnFloatType> &s)
{
long i, j, k, num;
matrix<YourOwnFloatType> mtr2(f.m,f.m+1),res(f.m,1);
//формируем матрицу из двух для решения СЛАУ
for(i=0;i<f.m;i++)
for(j=0;j<f.n;j++)
mtr2[i][j]=f[i][j];
for(i=0;i<f.m;i++)
mtr2[i][f.m]=s[i][0];
//выбор главного элемента
for(i=0;i<f.m;i++)
{
YourOwnFloatType max=mtr2[0][i];
for(j=i,num=i;j<f.m;j++)
if(fabs(mtr2[j][i])>max)
max=fabs(mtr2[j][i]),num=j;
if(num!=i)
for(k=0;k<f.m+1;k++)
{
YourOwnFloatType temp=mtr2[num][k];
mtr2[num][k]=mtr2[i][k];
mtr2[i][k]=temp;
}
}
for(i=0;i<f.m;i++)
if(!mtr2[i][i])
throw xmsg("Возможно, матрица вырождена \n");
//Прямой ход Гаусса
for(i=0;i<f.m;i++)
{
double sw=mtr2[i][i];
for(j=0;j<f.m+1;j++)
mtr2[i][j]/=sw;
for(k=i+1;k<f.m;k++)
{
double c=mtr2[k][i];
for(j=0;j<f.m+1;j++)
mtr2[k][j]-=mtr2[i][j]*c;
}
}
//Обратный ход Гаусса
for(i=f.m-2;i>=0;i--)
{
double s=0;
for(j=i+1;j<f.m;j++)
s+=mtr2[i][j]*mtr2[j][f.m];
mtr2[i][f.m]-=s;
}
//переписываем результат
for(i=0;i<f.m;i++)
res[i][0]=mtr2[i][f.m];
return res;
}
/*
Для обращения матрицы и нахождения детерминанта классическими методами нужна функция получения минора матрицы для данного её элемента
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> matrix<YourOwnFloatType>::minor (long a, long b)
{
matrix<YourOwnFloatType> temp(getm()-1,getn()-1);
for(long i1=0,i2=0;i1<getm();i1++)
if(i1!=a)
{
for(long j1=0,j2=0;j1<getn();j1++)
if(j1!=b)
temp[i2][j2]=(*this)[i1][j1],
j2++;
i2++;
}
return temp;
}
/*
Классический детерминант - вычисление разложением по одной из строк
*/
template <class YourOwnFloatType> YourOwnFloatType det2(matrix<YourOwnFloatType> &x)
{
if(x.getm()!=x.getn())
throw xmsg("Детерминант определён только для квадратных матриц \n");
if(x.getm()==1)/*особые случаи - детерминанты 1-го и 2-го порядка*/
return x[0][0];
if(x.getm()==2)
return x[0][0]*x[1][1]-x[0][1]*x[1][0];
YourOwnFloatType result=0;
for(long i=0;i<x.getm();i++)
result+=pow(-1,1+i)*det(x.minor(1,i))*x[1][i];
return result;
}
/*
Медленная (аналитическая) операция обращения квадратной матрицы
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> matrix<YourOwnFloatType>:: operator!()
{
if(getm()!=getn())
throw xmsg("Для неквадратных матриц обращение не определено \n");
matrix<YourOwnFloatType> temp=*this;
YourOwnFloatType _det=det(*this);
if(_det==(YourOwnFloatType)0)
throw xmsg("Обращение особенной матрицы невозможно \n");
for(long i=0;i<getm();i++)
for(long j=0;j<getn();j++)
temp[i][j]=pow(-1,2+i+j)*
det(minor(j,i))/_det;
return temp;
}
/*
Численное нахождение детерминанта методом Гаусса
*/
template <class YourOwnFloatType>
YourOwnFloatType det(matrix<YourOwnFloatType> &y) //быстрый детерминант
{
YourOwnFloatType sw,c,det,max;
long i,j,k,how,num;
matrix<YourOwnFloatType> x=y;
if(x.getm()!=x.getn())
throw xmsg("Детерминант определён только для квадратных матриц \n");
if(x.getm()==1)
return x[0][0];
if(x.getm()==2)
return x[0][0]*x[1][1]-x[0][1]*x[1][0];
for(how=i=0;i<x.getm();i++)
{
max=x[0][i];
for(j=i,num=i;j<x.getm();j++)
if(fabs(x[j][i])>fabs(max))
{
max=fabs(x[j][i]);
num=j;
}
if(num!=i)
{
for(k=0;k<x.getm();k++)
{
YourOwnFloatType temp=x[num][k];
x[num][k]=x[i][k];
x[i][k]=temp;
}
how++;
}
}
for(i=0;i<x.getm();i++)
{
for(sw=x[i][i],j=i+1;j<x.getm();j++)
x[i][j]/=sw;
for(k=i+1;k<x.getm();k++)
for(c=x[k][i],j=0;j<x.getm();j++)
x[k][j]-=x[i][j]*c;
}
for(det=1,i=0;i<x.getm();i++)
det*=x[i][i];
det*=pow(-1,how);
return det;
}
/*
Быстрое (численное) обращение матрицы
*/
template <class YourOwnFloatType> matrix <YourOwnFloatType> matrix<YourOwnFloatType>:: operator*()
{
long i, j, k, l;
YourOwnFloatType v,maxabs,s,tsr;
matrix<YourOwnFloatType> temp(getm(),2*getn());
if(getm()!=getn())
throw xmsg("Для неквадратных матриц обращение не определено \n");
if(det(*this)==(YourOwnFloatType)0)
throw xmsg("Обращение особенной матрицы невозможно \n");
for(i=0;i<getm();i++)
{
for(j=0;j<getn();j++)
temp[i][j]=(*this)[i][j];
for(j=getm();j<2*getm();j++)
temp[i][j]=(j==i+getm()) ? 1 : 0;
}
for(i=0;i<getm();i++)
{
for(maxabs=fabs(temp[i][i]),k=i,l=i+1;
l<getm();l++)
if(fabs(temp[l][i])>fabs(maxabs))
maxabs=fabs(temp[l][i]), k=l;
if(k!=i)
for(j=i;j<2*getm();j++)
v=temp[i][j],
temp[i][j]=temp[k][j], temp[k][j]=v;
}
for(i=0;i<getm();i++)
{
for(s=temp[i][i],j=i+1;j<2*getm();j++)
temp[i][j]/=s;
for(j=i+1;j<getm();j++)
{
for(tsr=temp[j][i],k=i+1;k<2*getm();k++)
temp[j][k]-=temp[i][k]*tsr;
}
}
for(k=getm();k<2*getm();k++)
for(i=getm()-1;i>=0;i--)
{
for(tsr=temp[i][k],j=i+1;j<getm();j++)
tsr-=temp[j][k]*temp[i][j];
temp[i][k]=tsr;
}
matrix<YourOwnFloatType> result=(*this);
for(i=0;i<getm();i++)
for(j=getm();j<2*getm();j++)
result[i][j-getm()]=temp[i][j];
return result;
}
/*
Быстрый детерминант как операция
*/
template <class YourOwnFloatType> inline YourOwnFloatType matrix<YourOwnFloatType>::operator&()
{
return det(*this);
}
/*
Детерминант рассматривается как оператор преобразования матрицы в число
*/
template <class YourOwnFloatType> inline matrix<YourOwnFloatType>::operator YourOwnFloatType()
{
return det(*this);/*используем быстрый детерминант*/
}
#endif
