Результат (см. картинки): Можете тоже проверить. Порядок матрицы совсем не заоблачный. При обращении матрицы использовался метод Жордана с поиском ведущего элемента по всей матрице. Все действительные типы – long double.
Оказывается, при перестановке слагаемых сумма может меняться во много раз!
static long double *xmul;
void mul_matr (long double c[], long double a[], long double b[], int row, int col, int com)
{
for (int i = row; i--; )
for (int j = col; j--; ) {
for (int k = com; k--; )
xmul[k] = a[com * i + k] * b[col * k + j];
c[col * i + j] = sum_matr(xmul, com);
}
}
void mul_matr0 (long double c[], long double a[], long double b[], int row, int col, int com)
{
for (int i = row; i--; )
for (int j = col; j--; ) {
long double tmp = 0;
for (int k = com; k--; )
tmp += a[com * i + k] * b[col * k + j];
c[col * i + j] = tmp;
}
}
long double sum_matr (long double a[], int col)
{
long double ved, koef, tmp;
int im, jm;
while (1) {
// Поиск самой приемлемой пары для суммирования
ved = -1.0;
for (int i = col; --i; )
if (a[i] != 0.0)
for (int j = i; j--; )
if (a[j] != 0.0) {
tmp = a[i] / a[j];
// Двойное предпочтение елементам с разными знаками
koef = tmp < 0.0 ? 2.0 : 1.0;
tmp = fabsl(tmp);
// Отношение д.б. меньшего числа к большему (по модулю)
if (tmp > 1.0) tmp = 1 / tmp;
tmp *= koef;
if(ved < tmp) {
ved = tmp;
im = i;
jm = j;
}
}
if (ved == -1.0) return a[im]; // Одно ненулевое значение
a[im] += a[jm];
a[jm] = 0.0;
}
}