Если да, то это же то же самое, что и преобразование матрицы к I + одна строчка для совершения тех же действий со второй матрицей, которая изначально I...
И пожалуйста простите ананимусу глупые вопросы :) ...
1. Я конечно понимаю, что код программы есть лучший комментарий,
но иногда хотелось бы знать, что делает та или иная часть программы,
а заодно и метод по которому она это делает.
Все ж таки код пишут люди и для людей (машине все равно что молотить:)
2. Хотелось бы также узнать, какие преимушества имеет Ваш код в сравнении с другими пакетами (ну например Lapack)
3. И наконец: какая лицензия на все это богатство.
Берешь любую книжку по чисметам, хотя бы того же Бахвалова или Рябенького и делаешь. Или в Инете готовый алгоритм берешь, а вообще говоря для таких вещей gnu fortran 95 есть - под него библиотек до фига и с С/С++ дружит.
Я всё-таки сам буду писать. С библиотеками не буду писать, я не профессиональный программист, чтобы загаживать программу всякой фигнёй и подключать ненужные библиотеки. Код выдирать тоже не буду, я не попсовик! ;)
double mypow(double a, int n)
{
int i;
double ret=1;
for(i = 0; i < n; i++) ret *= a;
return ret;
}
double delete_from_array(double *source, int m, int i, int j, double *dest);
double determinant(double *a, int m);
double cofactor(int m, int i, int j, double *a);
void inverse(double *source, int m, double *dest);
void transpose(double *matrix, int m);
main()
{
inverse(m1, 4, m2);
int i;
for(i = 0; i < 16; i++) printf("%.3G\n", m2[i]);
}
double delete_from_array(double *source, int m, int i, int j, double *dest)
{
int r, s, t;
t = 0;
for(s = 0; s < m; s++)
for(r = 0; r < m; r++)
{
if(s != j && r != i)
{
dest[t] = source[m*s + r];
t++;
}
}
}
double determinant(double *a, int m)
{
if (m == 2) return (a[0]*a[3] - a[2] * a[1]);
if (m == 1) return a[0];
double d = 0.0;
int i;
for(i = 0; i < m; i++) d += a[i*m] * cofactor(m, 0, i, a);
return d;
}
double cofactor(int m, int i, int j, double *a)
{
double d = mypow(-1.0, i+j+2);
double t[(m-1)*(m-1)];
delete_from_array(a, m, i, j, t);
d *= determinant(t, m-1);
return d;
}
void transpose(double *matrix, int m)
{
int i, j;
double temp;
for(j = 0; j < m; j++)
for(i = j+1; i < m; i++)
{
temp = matrix[j*m + i];
matrix[j*m + i] = matrix[i*m + j];
matrix[i*m + j] = temp;
}
}
void inverse(double *source, int m, double *dest)
{
int i;
double C[m*m];
double det = 0.0;
for(i = 0; i < m*m; i++) C[i] = cofactor(m, i%m, i/m, source);
for(i = 0; i < m; i++) det += C[i*m] * source[i*m];
transpose(C, m);
for(i = 0; i < m*m; i++) dest[i] = C[i]/det;
}