Я исчисляю до 2k x 2k на процессоре, используя многопоточность, используя собственную библиотеку, поэтому теперь это в 3,5-3,65 раза быстрее (зависит от размера матрицы), чем при использовании одного потока. Я использовал процессор Intel Xeon 3,5 ГГц E5-1620 v3 и 16 ГБ оперативной памяти. (К сожалению, я удалил старую версию, чтобы добавить точные значения, но если есть какой-то приоритет, я мог бы переписать ПО)
Это мой алгоритм обратной матрицы, с которым я сравнивал. (Правильно, как я сказал, сделав много тестов снова превосходным результатом):
/*
Uses 2D arrays.
Main routines are:
init_2Dvector() that initializes any 2d vector (can be uchar, char, int, float or double)
multiply_2Dvector()
inverse()
*/
#include<iostream>
#include <vector>
#include <stdlib.h>
#include <time.h>
using namespace std;
/*
void print_2Dvector(vector<vector<double> >& vec)
{
size_t xmax, ymax;
ymax = vec.size();
xmax = vec[0].size();
int x, y;
for (y = 0; y < ymax; y++)
{
for (x = 0; x < xmax; x++)
cout << vec[y][x] << " \t";
cout << endl;
}
}*/
void print_2Dvector(vector<vector<double> >& vec,char *format="%lg \t")
{
size_t xmax, ymax;
ymax = vec.size();
xmax = vec[0].size();
int x, y;
for (y = 0; y < ymax; y++)
{
{
for (x = 0; x < xmax; x++)
printf(format, vec[y][x]);
}
cout << endl;
}
}
//Resizes to y_dim,x_dim any kind of 2d array:
template<typename T>
void init_2Dvector(vector<vector<T> >& vec, size_t y_dim, size_t x_dim)
{
vec.resize(y_dim);
for (size_t i = 0; i < vec.size(); i++)
vec[i].resize(x_dim);
}
//Returns vec1*vec2. vec1 and 2 are not touch
vector< vector<double> > multiply_2Dvector(vector< vector<double> > vec1, vector< vector<double> > vec2)
{
size_t xmax, ymax;
ymax = vec1.size();
xmax = vec1[0].size();
vector< vector<double> > vec3;
if ((ymax != vec2[0].size()) || (xmax != vec2.size()))
{
cout << "ERROR on dim2_multiply() dimensions of vec2 not corresponding with vec1 ones" << endl; getchar(); return{};//returns a null
}
init_2Dvector(vec3, ymax, ymax);
cout << "dimensions of vec3=" << vec3.size() << " x " << vec3[0].size() << endl;
double xx;
for (int y = 0; y < ymax; y++)
for (int x = 0; x < ymax; x++)
{
xx = 0.0;
for (int t = 0; t < xmax; t++)
xx += vec1[y][t] * vec2[t][x];
vec3[y][x] = xx;
}
return vec3;//ok
}
//returns inverse of x2, x2 is not modified
vector< vector<double> > inverse(vector< vector<double> > x)
{
if (x.size() != x[0].size())
{
cout << "ERROR on inverse() not square array" << endl; getchar(); return{};//returns a null
}
size_t dim = x.size();
int i, j, ord;
vector< vector<double> > y(dim,vector<double>(dim));//output
//init_2Dvector(y, dim, dim);
//1. Unity array y:
for (i = 0; i < dim; i++)
{
y[i][i] = 1.0;
for (j = i+1; j < dim; j++)
{
y[i][j]= y[j][i] = 0.0;
}
}
double diagon, coef;
double *ptrx, *ptry, *ptrx2, *ptry2;
for (ord = 0; ord<dim; ord++)
{
//2 Hacemos diagonal de x =1
int i2;
if (fabs(x[ord][ord])<1e-15) //Si el elemento diagonal es 0 sumamos una columna que no sea 0 el elemento correspondiente
{
for (i2 = ord + 1; i2<dim; i2++)
{
if (fabs(x[i2][ord])>1e-15) break;
}
if (i2 >= dim)
return{};//error, returns null
for (i = 0; i<dim; i++)//sumo la linea que no es 0 el de la misma fila de ord
{
x[ord][i] += x[i2][i];
y[ord][i] += y[i2][i];
}
}
diagon = 1.0/x[ord][ord];
ptry = &y[ord][0];
ptrx = &x[ord][0];
for (i = 0; i < dim; i++)
{
*ptry++ *= diagon;
*ptrx++ *= diagon;
}
//Hacemos '0' la columna ord salvo elemento diagonal:
for (i = 0; i<dim; i++)//Empezamos por primera fila
{
if (i == ord) continue;
coef = x[i][ord];//elemento ha hacer 0
if (fabs(coef)<1e-15) continue; //si es cero se evita
ptry = &y[i][0];
ptry2 = &y[ord][0];
ptrx = &x[i][0];
ptrx2 = &x[ord][0];
for (j = 0; j < dim; j++)
{
*ptry++ = *ptry - coef * (*ptry2++);//1ª matriz
*ptrx++ = *ptrx - coef * (*ptrx2++);//2ª matriz
}
}
}//end ord
return y;
}
void test_5_inverse()
{
vector< vector<double> > vec1 = {
{0,-5,0,7,33,11,-1},
{72,0,-11,7,9,33,5 },
{-13,31,-5,15,29,30,24 },
{-24,9,8,-23,31,-12,4 },
{-3,-22,4,-24,-5,27,-10 },
{-10,-21,-16,-32,-11,20,14 },
{5,30,13,-32,29,-13,-13 }
};
vector< vector<double> > vec2;
vec2 = inverse(vec1);
vector< vector<double> > vec3;
vec3 = multiply_2Dvector(vec1, vec2);
cout << "initial array (must be unmodified):" << endl;
print_2Dvector(vec1);
cout << "Must be diagon array:" << endl;
print_2Dvector(vec3," %8.3lf");
cout << endl;
}
void test_6_inverse(int dim)
{
vector< vector<double> > vec1(dim, vector<double>(dim));
for (int i=0;i<dim;i++)
for (int j = 0; j < dim; j++)
{
vec1[i][j] = (-1.0 + 2.0*rand() / RAND_MAX) * 10000;
}
vector< vector<double> > vec2;
double ini, end;
ini = (double)clock();
vec2 = inverse(vec1);
end = (double)clock();
cout << "Time inverse =" << (end - ini) / CLOCKS_PER_SEC << endl;
vector< vector<double> > vec3;
vec3 = multiply_2Dvector(vec1, vec2);
cout << "initial array (must be unmodified):" << endl;
//print_2Dvector(vec1);
cout << "Must be diagon array:" << endl;
//print_2Dvector(vec3, " %8.3lf");
cout << endl;
}
int main()
{
vector< vector<double> > vec1;
init_2Dvector(vec1, 10, 5); //size_t ymax = vec1.size(),xmax = vec1[0].size();
//test_2_dimension_vector();
//test_one_dimension_vector();
test_5_inverse();
test_6_inverse(1300);
cout << endl << "=== END ===" << endl; getchar();
return 1;
}
person
mathengineer
schedule
10.07.2018
