Comments 6
У вас шаблонная функция, а епсилон фиксированная? А вообще для этого есть std::numeric_limits<T>::epsilon
abs(*(alfa.m_data + j * alfa.m_columns + j)) < 1.0e-36
Что за магическое число?
while (norm < 1.0e-17);
Уверены, что тут никогда не будет деления на 0?
sum / *(alpha.m_data + i * alpha.m_columns + i);А что произойдет, если я передам T = std::string?
template <typename T>
T MATRIX<T>::normI()
{
T norm = 0;Для копирования столбцов по-хорошему нужно перегружать оператор присваивания, а не писать отдельную функцию CopyColumn. Еще и VECTOR по неконстантной ссылке.
template <typename T>
void MATRIX<T>::CopyColumn(VECTOR<T>& v, int j)
{
Да и зачем вам свой matrix класс, если кроме обращения элементов тут почти ничего нет, не все операторы перегружены, безопасности исключений тут нет. Использовали бы тогда уж std::vector<std::vector>> там хотя-бы более безопасно работать.
Или на питоне, на порядки же проще да и не напишите вы на С++ быстрее numpy.
Ну и вы проверяли на численную стабильность? Обычно для этого используют LUP разложение. А у вас вроде такого нет.
А вообще смысл статьи? Нового ничего нет. Если учебный проект, то сделан он халтурно, как в плане кода, так и в плане математики. Особенности вычислений с плавающей точкой не учтены от слова совсем.
Если порассуждать про какие-либо принципиально новые методы вычислений то может обратить внимание на численно-символические. Составим простой скрипт на Maxima или любой другой CAS:
Q:makelist(makelist(a[j,i],i,1,4),j,1,4);
A: apply ('matrix,Q);
B:A^^-1;
horner(num(B[1,1]));
horner(denom(B[1,1]));Сформировать например, матрицу 4х4, взять ей обратную и посмотреть что получается с числителем и знаменателем, например, для первого элемента применив схему Горнера.
Разумеется, что это следует из определителей матриц-миноров, однако, можно выявить индексную арифметику для некоего цифрового фильтра с заданными коэффициентами и операцией умножения с накоплением. В этом случае знаменатель считается один раз (он общий для всех элементов обратной матрицы). А вот числители вычисляются для каждого элемента примерно так:
Умножение с накоплением идёт по внутренним скобкам как коэффициент умноженный на отсчёт. Иными словами мы представляем матрицу в виде массивов отсчётов как это бы делали для цифрового фильтра - умножение на коэффициент. То есть фактически всё распадается на работу с индексами и формирование промежуточных массивов с результатами. И если расположить элементы как показано то получается что-то вроде "бабочки" для БПФ. Наверняка может быть матрицы с размерностью кратной степени 2 имеют такой сверхбыстрый алгоритм.
Ага, сиквел к прошлой статье подъехал. Зачем это на Хабре? Линал запрогали и выложили... Никаких новых идей, как будто лаба во вычматам.
Так, теперь код-ревью. Иду сверху вниз по файлам.
MatrixVectorTemplate.hpp - тут вижу using namespace std;. Вас за такое на части разорвать могут. Кстати, инклюдить надо не math.h, а cmath.
Дальше, у вас умерла кодировка в файле => русский язык выглядит, как тарабарщина.
Зачем метод класса inline int size()? Инлайнить компилятор будет и без ваших подсказок. А почему не помечено словом const? Хотите, чтобы у константных векторов нельзя было спросить размер?
Почему указатель m_data инициализируете числовой константой 0? У вас отняли nullptr? Ну хотя бы NULL, но всё равно плохо.
А почему не используете списки инициализации в конструкторах?
Из вашей реализации копи-конструктора я вижу, что у вас отвратительные знания основ ООП в C++, т.к. вы зачем-то проверяете, нет ли у вас там аллоцированного массива и предварительно удаляете, если есть.
В реализации оператора присваивания отсутствует защита от присваивания себе.
Почему-то вы считаете, что умножать векторы разной длины - норм (возвращаем ноль), а вот складывать - не норм (делаем throw "Âåêòîðû èìåþò ðàçíóþ ðàçìåðíîñòü";).
Почему бросаете не что-то, наследующеся от std::exception? Почему в блоке после выброса исключения какой-то ещё код?
Кстати, почему вы понапихали френдов в код? Обычно бинарные операторы, тип результата которых совпадает с типами обоих аргументов, реализуют через копию + составное присваивание...
А ещё странно, что вы где-то используете this для доступа к полям, а где-то - нет. Нужен общий стиль. C-style касты понапиханы... В целом, очень много непонятных названий, особенно всякие однобуквенные переменные. Дефайны в коде есть для констант - тоже плохо... Магические числа в рандомных местах (всякие эпсилоны в сравнениях).
std::rand используете вместо качественных генераторов, а в функции randomDouble при неудачном рнг делите на ноль. Кстати, непонятно, почему дабл в названии есть - это можно инстанцировать и для целых типов.
Уфф... Язык знаете плохо, общая задумка слабая, совсем непонятно, зачем работаете с голой памятью вместо уже имеющихся в языке классов.
а теперь сравните скорость вашего велосипеда с каким-нибудь настоящим BLAS
/// <summary>
/// Äåñòðóêòîð: îñâîáîæäåíèå ïàìÿòè, çàíèìàåìîé ïîä ìàòðèöó
/// </summary>
template <typename T>
MATRIX<T>::~MATRIX()
{
if (!m_data) delete[] m_data;
}
🤭
Матрицы и векторы: вычисление обратной матрицы