/*
  СЛАУ. Задача 9
  Написать, протестировать и оценить эффективность функции вычисления
  обратной матрицы методом Фаддеева.
  Исследовать зависимость точности и времени обращения от типа представления
  (float, double) и размерности матриц.
*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <alloc.h>
#include <conio.h>

#define RND (5-(float)rand()/3276.7)

FILE *fp;

// Функция umn_m() - перемножение матриц
void umn_m(float *a, float *b, float *c, int n, int m, int p)
{
  int i, j, k;
  float s;

  for(i=0; i<n; i++)
    for(j=0; j<p; j++)
      {
	for(s=k=0; k<m; k++)
	  s+=*(a+i*m+k)**(b+k*p+j);
	*(c+i*p+j)=s;
      }
}

// Функция spur() - вычисление следа матрицы
float spur(float *a, int n)
{
  int i;
  float sp=0;

  for(i=0; i<n; i++)
    sp+=*(a+i*n+i);
  return sp;
}

// Функция faddeev() - обращение матриц методом Фаддеева
// Обратная матрица получается на месте исходной!
void faddeev(float *a, int n)
{
  int i, j, m;
  float *ai, *b, p;

// Выделяем память под вспомогательныe матрицы ai(n,n) и b(n,n)
  ai=(float*)malloc(n*n*sizeof(float));
  b=(float*)malloc(n*n*sizeof(float));

  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      *(ai+i*n+j)=*(a+i*n+j);
  for(m=1; m<=n; m++)
    {
      p=spur(ai, n)/m;
      if(m==n) break;
      for(i=0; i<n; i++)
	for(j=0; j<n; j++)
	  if(i==j) *(b+i*n+j)=*(ai+i*n+j)-p;
	   else *(b+i*n+j)=*(ai+i*n+j);
      umn_m(a, b, ai, n, n, n);
    }
  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      *(a+i*n+j)=*(b+i*n+j)/p;
  free(ai); free(b);
}

void main()
{
  int i, j, n;
  long *t, t1, t2, dt;
  float *a, *a1, *c, sc, sum, k_norm;

  clrscr();
// Откpываем файл fad_d.txt для записи результатов
  if((fp=fopen("fad_f.txt","a"))==NULL)
    {
      printf("\n Ошибка пpи откpытии файла ");
      exit(0);
    }
  fprintf(fp,"\n Исследование программы обращения матриц "
	     "методом Фаддеева (float) \n");
  printf("\n Исследование программы обращения матриц "
	 "методом Фаддеева (float) \n");
 m1:
  printf("Введите размерность матрицы n<60 : ");
  scanf("%d", &n);

// Выделяем память под матрицы A, A1, C
  a=(float*)malloc(n*n*sizeof(float));
  a1=(float*)malloc(n*n*sizeof(float));
  c=(float*)malloc(n*n*sizeof(float));

// Формируем исходную матрицу
  srand(time(NULL));
  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      *(a1+i*n+j)=*(a+i*n+j)=RND;

// Обращаем матрицу A и фиксируем время обращения
  t1=time(t);
  faddeev(a1,n);
  t2=time(t);
  dt=t2-t1;
// Проверяем: C=A*A1, вычисляем k-норму матрицы С
  umn_m(a,a1,c,n,n,n);
  for(sum=i=0; i<n; i++)
    for(j=0; j<n; j++)
      {
	if(i==j) sc=*(c+i*n+j)-1;
	 else sc=*(c+i*n+j);
	sum+=sc*sc;
      }
  k_norm=sqrt(sum/(n*n));
  fprintf(fp,"\n   n = %2d \t t(c) = %ld \t k_norm = %.2e  ", n, dt, k_norm );
  printf("   n = %2d \t t(c) = %ld \t k_norm = %.2e  ", n, dt, k_norm );
  free(a); free(a1); free(c);
  printf("\n\n Еще? Введи 1 : ");
  scanf("%d", &i);
  if(i==1)
   {
     gotoxy(1, wherey()-1);
     delline();
     goto m1;
   }
  fprintf(fp,"\n\n Исследование окончено \n");
  printf("\n\n Исследование окончено \n");
  fclose(fp); getch();
}
/* Конец программы */

