#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
typedef float * FU ; // Чтоб не писать присвоим FU тип данных float*
typedef double * DU ; // Чтоб не писать присвоим DU тип данных double*

void pechat(FU, int, int, char *);// Функция печати матрицы
void hand(FU,int,int,char *);    // Функция ножного ввода матрицы
void trans(FU,int);     // Функция транспонирования матрицы
void prisv(FU,FU,int,int); // Функция присваивания матрицы 2=1
FU memory(int, int);// Выделение памяти под матрицу размерностью n,m
void jacobi(int n,char eivec,FU a,FU d,FU v,int rot);

     //Это MAIN ..... Велик и могуч наш MAIN!
     // Редкая переменная доберется до его конца в первоначальном виде.
void main(void)
{ int r; FU M1,MD,MV;
  clrscr();  //Установка текстового режима
  printf("\n Input размерность матриц r:");scanf("%d",&r);
  M1=memory(r,r);
  MD=memory(r,1);
  MV=memory(r,r);
  hand(M1,r,r,"Матрица М1");
  jacobi(r,1,M1,MD,MV,60);
 // hand(M2,r,r,"Матрица M2");
  pechat(M1,r,r,"\n M1 :");
  pechat(MD,r,1,"\n MD :");
  pechat(MV,r,r,"\n MV :");
  free(M1); free(MD);
   free(MV);

}
/*--Печать матрицы и текста--На экран--*/
void pechat(FU a,int m,int n,char *str)
{
  int i,j;

  puts(str);
  for(i=0;i<m;i++)
    {
      for(j=0;j<n;j++)
	printf("%6.5f ",*(a+i*n+j));
      printf("\n");
    }
   printf("Ткни че нибудь ");getch();
}
/*--Выделение памяти--*/
FU memory(int n,int m)
{
  FU mas;
  mas=(FU)malloc(m*n*sizeof(float));
  return(mas);
}
/*--Присвоение матрице B значения матрицы A--*/
void prisv(FU a,FU b,int n,int m)
{ int i,j;
  for(i=0;i<n;i++)
   for(j=0;j<m;j++)
    *(b+i*m+j)=*(a+i*m+j);
}
/*--Ручной вход--*/
void hand(FU a,int n,int m,char *help)
{
  int x=1,y=5,i,j;
  float p;
  clrscr();
  gotoxy(1,1);puts(help);printf("Размерность матрицы [%d,%d]",n,m);
  for(i=0;i<n;i++)
    { for(j=0;j<m;j++)
      {gotoxy(30,2);printf("Элемент [%d,%d]",i+1,j+1);
      gotoxy(x,y);x+=13; scanf("%f",&p);
      *(a+i*m+j)=p;
      }
      x=1; y+=2;
    }
}
//Транспортирование матрицы
void trans(FU a,int n)
{ float buf;
  int i,j;
  for(i=0;i<n;i++)
      for(j=i;j<n;j++)
      {	buf=*(a+i*n+j);
	*(a+i*n+j)=*(a+j*n+i);
	*(a+j*n+i)=buf;
      }
}
//Метод Якобиби-96 _ChicagoTEST v1.0
void jacobi(int n,char eivec,FU a,FU d,FU v,int rot)
{ float sm,c,s,t,h,g,tau,theta,tresh;
  int p,q,i,j;FU b,z;
  b=memory(n,1);z=memory(n,1);

  if(eivec)
  { for(p=0;p<n;p++)
     for(q=0;q<n;q++)
     if(p==q) *(v+p*n+q)=1.0;else *(v+p*n+q)=0;
  }
  for(p=0;p<n;p++)
  {b[p]=d[p]=*(a+n*p+p);z[p]=0;}
  rot=0;
  for(i=0;i<150;i++)
  { sm=0;
    for(p=0;p<n-1;p++)
     for(q=p+1;q<n;q++)
     sm+=fabs(*(a+p*n+q));

    if(sm==0) {printf("May day!!! May day!!! I'm down ;(");goto out;}

    if(i<4)  tresh=0.2*sm/(n*n);    else    tresh=0;

    for(p=0;p<n-1;p++)
     for(q=p+1;q<n;q++)
     { g=100.0*fabs(*(a+p*n+q));

       if(i>3&&( (fabs(d[p])+g)==fabs(d[p]))&&( (fabs(d[q])+g)==fabs(d[q]))) *(a+p*n+q)=0;

       else if(fabs(*(a+p*n+q))>tresh)
	    { h=d[q]-d[p];
	      if( ((fabs(h)+g)==fabs(h)) ) t=*(a+p*n+q)/h;
	      else {
		     theta=0.5*h / (*(a+p*n+q));
		     t=1.0/(fabs(theta) + sqrt( 1.0 + theta*theta) );
		     if (theta<0) t=-1.0*t;
		    }
	      c=1.0/sqrt(1.0+t*t);
	      s=t*c;
	      tau=s/(1.0+c);
	       printf("\nRotation N=%d I=%d",rot,i);
	      h=t*(*(a+p*n+q));
	      z[p]=z[p]-h;
	      z[q]=z[q]+h;
	      d[p]=d[p]-h;
	      d[q]=d[q]+h;
	      *(a+p*n+q)=0;
	      for(j=0;j<p-1;j++)
	      { g=*(a+j*n+p);h=*(a+j*n+q);
		*(a+j*n+p)=g-s*(h+g*tau);
		*(a+j*n+q)=h+s*(g-h*tau);
	      }
	      for(j=p+1;j<q-1;j++)
	      { g=*(a+p*n+j);h=*(a+j*n+q);
		*(a+p*n+j)=g-s*(h+g*tau);
		*(a+j*n+p)=h+s*(g-h*tau);
	      }
	      for(j=q+1;j<n;j++)
	      { g=*(a+p*n+j);h=*(a+q*n+j);
		*(a+p*n+j)=g-s*(h+g*tau);
		*(a+q*n+j)=h+s*(g-h*tau);
	      }
	      for(j=0;j<n;j++)
	      { g=*(v+j*n+p);h=*(v+j*n+q);
		*(v+j*n+p)=g-s*(h+g*tau);
		*(v+j*n+q)=h+s*(g-h*tau);
	      }
	      rot++;

	    }
     }
     for(p=0;p<n;p++)
     {d[p]=b[p]=b[p]+z[p];z[p]=0;}
  }
 out: free(b);free(z);  //Проверено МИН НЕТ!
}

