Обращение матриц на C

Обращение это нахождение обратной?

Если да, то это же то же самое, что и преобразование матрицы к I + одна строчка для совершения тех же действий со второй матрицей, которая изначально I...

эххххх, молодость!


#define __HICOLOR__
#include "sergsoft.cpp"

#define NDim 25
#define MS 4
#define sqr(X) ( (X) * (X) )
#define WAIT\
  printf ("\n Press any key ...");\
  getc (stdin);
dword XX, YY, i, j;
double tmp, tmp1;

double Norma (double *M)
{
  double Res = 0.0, Sum;
  word i, j;

  for ( i = 0; i < MS; i++ )
  {
    Sum = 0.0;
    for ( j = 0; j < MS; j++ )
      Sum += ABS (M[i * MS + j]);
    if ( Sum > Res)
      Res = Sum;
  }

  return Res;
}

void CreateMatrix (double *M, double S)
{
  word i;

  memset (M, 0, sizeof (double) * SQR (MS));
  for ( i = 0; i < MS; i++ )
  {
    M[i * MS + i] = - (double)(2 * i + 1) * MS + (double)i * S - (double)(2 * i * i);
    if ( i != 0 )
      M[i * MS + i - 1] = (double)i * (MS  - double(i) + 1);
    if ( i != (MS - 1) )
      M[i * MS + i + 1] = ((double)i + 1) * (MS + S - (double)i);
  }
}

void PrintMatrix (FILE *F, double *M)
{
  word i;

  for ( i = 0; i < SQR (MS); i++ )
  {
    if ( (i % MS) == 0)
      fprintf (F, "\n\n");
    if ( (double)((signed long)M[i]) == M[i] )
      fprintf (F, "%i ш ", (signed long)M[i]);
    else
      fprintf (F, "%lf ш ", M[i]);
  }
}
 /*
  procedure DECOMP(
                   n  : integer;
             var   A  : floatmatrix ;
             var cond : float   ;
             var ipvt : ivector;
             var work : floatvector
                                   );*/

void Solve (word n, double *a, double *b, signed short *ipvt)
{
  signed short kb, km1, nm1, kp1, i, k, m;
  double t;

  if ( n == 1 )
  {
    b[0] = b[0] / a[0];
    return ;
  }
  nm1 = n - 1;
  for ( k = 1; k <= nm1; k++ )
  {
    kp1 = k + 1;
    m = ipvt[k - 1];
    t = b[m - 1];
    b[m - 1] = b[k - 1];
    b[k - 1] = t;
    for ( i = kp1; i <= n; i++ )
      b[i - 1] = b[i - 1] + a[(i - 1) * MS + (k - 1)] * t;
  }
  for ( kb = 1; kb <= nm1; kb++ )
  {
    km1 = n - kb;
    k = km1 + 1;
    b[k - 1] = b[k - 1] / a[(k - 1) * MS + (k - 1)];
    t = -b[k - 1];
    for ( i = 1; i <= km1; i++ )
      b[i - 1] = b[i - 1] + a[(i - 1) * MS + k - 1] * t;
  }
  b[0] = b[0] / a[0];
}

void Reverce (double *M, signed short *I, double *Res, double *Work)
{
  word i, j;

  for ( i = 0; i < MS; i++ )
  {
    memset (Work, 0, MS * sizeof (double));
    Work[i] = 1.0;
    Solve (MS, M, Work, I);
    for ( j = 0; j < MS; j++ )
      Res[j * MS + i] = Work[j];
  }
}

void Mult (double *M1, double *M2, double *Y)
{
  word i, j, k;

  for ( i = 0; i < MS; i++ )
    for ( j = 0; j < MS; j++ )
    {
      Y[i * MS + j] = 0;
      for ( k = 0; k < MS; k++ )
        (Y[i * MS + j]) += (M1[k * MS + j]) * (M2[i * MS + k]);
    }
}

void SubE (double *M)
{
  word i;

  for ( i = 0; i < MS; i++ )
    M[MS * i + i] -= 1.0;
}

void Decomp (word n, double *a, double &cond, signed short *ipvt, double *work)
{
  signed short nm1, i, j, k, kp1, kb, km1, m;
  double ek, t, anorm, ynorm, znorm;

  ipvt[n - 1] = 1;
  if ( n == 1 )
  {
    cond = 1.0;
    if ( a[0] != 0.0 )
      return ;
    else
    {
      cond = 1.0e+32;
      return;
    }
  }
  nm1 = n - 1;
  anorm = 0.0;
  /*
	for j:=1 to n do begin
		              t:=0.0;
		              for i:=1 to n do  t:=t+abs(a[i,j]);
		              if t > anorm then anorm:=t;
                         end;
   */
  for ( j = 1; j <= n; j++ )
  {
    t = 0.0;
    for ( i = 1; i <= n; i++ )
      t += ABS (a[(i - 1) * MS + (j - 1)]);
    if ( t > anorm )
      anorm = t;
  }
  for ( k = 1;  k <= nm1; k++ )
  {
    kp1 = k + 1;
    m = k;
    for ( i = kp1; i <= n; i++ )
      if ( ABS (a[(i - 1) * MS + k - 1]) > ABS (a[(m - 1) * MS + k - 1]) )
        m = i;
    ipvt[k - 1] = m;
    if ( m != k )
      ipvt[n - 1] = -ipvt[n - 1];
    t = a[(m  - 1) * MS + k - 1];
    a[(m  - 1) * MS + k - 1] = a[(k - 1) * MS + k - 1];
    a[(k - 1) * MS + k - 1] = t;
    if ( t !=  0.0 )
    {
      for ( i = kp1; i <= n; i++ )
        a[(i - 1) * MS + k - 1] = -a[(i - 1) * MS + k - 1] / t;
      for ( j = kp1; j <= n; j++ )
      {
        t = a[(m - 1) * MS + j - 1];
        a[(m - 1) * MS + j - 1] = a[(k - 1) * MS + j - 1];
        a[(k - 1) * MS + j - 1] = t;
        if ( t != 0.0 )
          for ( i = kp1; i <= n; i++ )
            a[(i - 1) * MS + j - 1] = a[(i - 1) * MS + j - 1] + a[(i - 1) * MS + k - 1] * t;
      }
    }
  }
  for ( k = 1; k <= n; k++ )
  {
    t = 0.0;
    if ( k != 1 )
    {
      km1 = k - 1;
      for ( i = 1; i <= km1; i++ )
        t = t + a[(i - 1) * MS + k - 1] * work[i - 1];
    }
    ek = 1.0;
    if ( t < 0.0 )
      ek = -1.0;
    if ( a[(k - 1) * MS + k - 1] == 0.0 )
    {
      cond = 1.0E+32;
      return ;
    }
    work[k - 1] = -(ek + t)  / a[(k - 1) * MS + k - 1];
  }
 for ( kb = 1; k <= nm1; k++ )
 {
    k = n - kb;
    t = 0.0;
    kp1 = k + 1;
    for ( i = kp1; i <= n; i++ )
      t = t + a[(i - 1) * MS + k - 1] * work[k - 1];
    work[k - 1] = t;
    m = ipvt[k - 1];
    if ( m != k )
    {
      t = work[m - 1];
      work[m - 1] = work[k - 1];
      work[k - 1] = t;
    }
 }
  ynorm = 0.0;
  for ( i = 1; i <= n; i++ )
    ynorm += ABS (work[i - 1]);
  Solve (n, a, work, ipvt);
  znorm = 0.0;
  for ( i = 1; i <= n; i++ )
    znorm += ABS (work[i - 1]);
  cond = anorm * znorm / ynorm;
  if ( cond < 1.0 )
    cond = 1.0;
  return ;
}
 

тебе нужна функция Reverce

Спасибо.

>эххххх, молодость!
...больница, любовь

>theserg

Огромное спасибо за код!!

И пожалуйста простите ананимусу глупые вопросы :) ...

1. Я конечно понимаю, что код программы есть лучший комментарий, но иногда хотелось бы знать, что делает та или иная часть программы, а заодно и метод по которому она это делает. Все ж таки код пишут люди и для людей (машине все равно что молотить:)

2. Хотелось бы также узнать, какие преимушества имеет Ваш код в сравнении с другими пакетами (ну например Lapack)

3. И наконец: какая лицензия на все это богатство.

Еще раз спасибо.

Блин, библиотек на С и Фортране по лин.алгебре море. Что, ничего не годится?

это всё перефигаченный пакет Форсайт (сначала с фортрана на паскакаль, потом уже мною с паскакаля на c/c++)

лицензия "бери, ибо добра такого в инете навалом!"

Берешь любую книжку по чисметам, хотя бы того же Бахвалова или Рябенького и делаешь. Или в Инете готовый алгоритм берешь, а вообще говоря для таких вещей gnu fortran 95 есть - под него библиотек до фига и с С/С++ дружит.

ха!

первый нах!

Я всё-таки сам буду писать. С библиотеками не буду писать, я не профессиональный программист, чтобы загаживать программу всякой фигнёй и подключать ненужные библиотеки. Код выдирать тоже не буду, я не попсовик! ;)

Под всякой фигнёй я имел в виду не предложенную мне функцию, а выдернутую из библиотеки.

а ты догадался, что сначала нужно выполнить decomp (LU-разложение), а только потом применять solve?

Я не разбирался ещё, сам писать буду.

Скачай Интеловский MKL и не парься. В нем оптимизированный LAPACK, да еще и с примерами использования.

#include <stdio.h>

double m1[] = {1.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 2.0, 4.0, 9.0, 1.0};
double m2[16];

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;
}