Вычисления собственных значений и векторов

Конкретно нужен:
Метод Якоби вычисления собственных значений и векторов.

Других не предлагать.

http://prografix.narod.ru/rus_jacobi.html

Сорсец на Buider6, нужен - готовая лаба + метод QR и спектрального разложения?

Посмотри сурс под Borland C++ 3.1

#include <math.h>
#include <conio.h>
#include <stdio.h>
#include <graphics.h>
#include <mem.h>

#define   DEF   4

const     double EPS = 0.00001;
const     N  = 1;
const int n1 = 8;
typedef  double  VECTOR[n1];
typedef  VECTOR  MATRIX[n1];

int Koren = 0;

void  OutMatrix(MATRIX _A)
{
     clrscr();
     for(int i=0; i<n1; i++)
      for(int j=0; j<n1; j++)
     {
         gotoxy( j*8+ 3, i+5);
         cprintf("\%4.2f", _A[i][j]);
         //if(i==j)cprintf("-L");
     }
     getch();
}

double P(VECTOR P1, double lambda)
{
     double s = pow(lambda, (double)DEF );
     for(int i=1; i<=DEF; i++)
       s += P1[i-1]*pow(lambda, (double)(DEF - i ) );

     return s;
}

double _P(VECTOR P1, double lambda, double h)
{
     return ( P(P1, lambda + h) - P(P1, lambda))/h;
}

void  MultMatrix(MATRIX M1, MATRIX &M2)
{
     int     i, j;
     MATRIX  Res;
     for( i=0; i<n1; i++)
     for( j=0; j<n1; j++)Res[i][j] = 0;

     for( i=0; i<n1; i++)
     for( j=0; j<n1; j++)
      {
          for(int k=0; k<n1; k++)
           Res[i][j] += M1[k][i]*M2[j][k];
      }
     for( i=0; i<n1; i++)
     for( j=0; j<n1; j++)M2[i][j] = Res[i][j];
}

void  Fix_Leverrier(MATRIX _A, VECTOR xP)
{
     VECTOR   xS;
     MATRIX   xA;
     double    lsd;
     int    i, j;
     for( i=0; i<n1; i++)
     for( j=0; j<n1; j++)
       xA[i][j] = _A[i][j];

     for(int k=0; k<n1; k++)
      {
       xS[k] = 0;
       for( i=0; i<n1; i++)
           xS[k] += xA[i][i];

       MultMatrix(_A, xA);
//        OutMatrix(xA);
//        getch();
      }
     xP[0] = xS[0];
     for( k=1; k<DEF; k++)
      {
        xP[k] = 0;
        for( i=0; i<k; i++)
           xP[k] += xP[i]*xS[k-i-1];

        xP[k] = (xS[k] - xP[k])/(double)(k+1);
      }

     clrscr();
     cprintf("\r\nP = x^\%d", DEF);
     for( k=0; k<DEF; k++)
     {
      if(xP[k]>0)
       cprintf("+");

      cprintf("\%6.2f*x^\%d", xP[k], DEF-k-1);
     }
     getch();
}

void  Newton(double alfa, double beta, VECTOR P1, VECTOR &xA)
{
     double  xk, xk1, a = alfa, h = 0.1, tmp;

     do{
       if( P(P1, a)*P(P1, xk = a += h) < 0)
       {
         do{
            xk1 = xk - P(P1, xk)/_P(P1, xk, 0.00001);
            tmp = (xk - xk1)/2;
            xk  = xk1;
         }while( fabs(tmp) < EPS);
         xA[Koren++] = (xk + xk1)/2;
       }
     }while(a < beta);
}

void  MFill(MATRIX c, MATRIX a)
{
     for(int i=0; i<n1; i++)
     for(int j=0; j<n1; j++)
          c[i][j] = a[i][j];
}
void  VFill(VECTOR &c, VECTOR a)
{
     for(int i=0; i<n1; i++)
          c[i] = a[i];
}
void  MVGauss(MATRIX a,VECTOR &b)
{
     int  i, j, l;
     MATRIX  A;
     VECTOR  B;

     MFill(A, a);
     VFill(B, b);

     for(l=0; l<n1; l++)
      {
      MFill(a, A);
      VFill(b, B);

      for(j=l; j<n1; j++)
         for(i=l+1; i<n1; i++)
         {
           A[i][j] = a[i][j] - (a[i][l] * a[l][j] ) / a[l][l];
           B[i]    = b[i]    - (b[l] * a[i][l] )  / a[l][l];
         }
      }
      MFill(a, A);
      VFill(b, B);
}


void  Polynom(double alfa, double beta, VECTOR Z, VECTOR xTRM)
{
     int    Gd=0, Gm=0;
     double  h = 0.0005;
     double  x  = alfa, y, MULT = 20;
     double  MX = 250, MY = 200;
     char   CH[20];

     initgraph(&Gd, &Gm, "c:\\language\\bc\\bgi");
     setcolor(7);
     line(0,MY,640,MY);
     line(MX,0,MX,480);

     setcolor(10);
     for(int i=0; i<Koren; i++)
      {
       line(xTRM[i]*MULT + MX, 0, xTRM[i]*MULT + MX,480);
       sprintf(CH, "\%6.4f", xTRM[i] );
       outtextxy(xTRM[i]*MULT + MX + 3, 20, CH);
      }
     do
      {
         y = P(Z, x += h);
         putpixel(x*MULT + MX, y/4 + MY, 15);
      }while(x < beta);
     getch();
     closegraph();
}

void  main( void )
{
     double    alfa = -15.0, beta = 20.0;
     int   i, j;
     VECTOR xP, xZ, xG;
     MATRIX _A;

     memset(&xP,0,sizeof(VECTOR));
     memset(&xZ,0,sizeof(VECTOR));
     memset(&xG,0,sizeof(VECTOR));

     for( i=0; i<n1; i++)
      for( j=0; j<n1; j++)
      if(i==j)_A[i][j] = i/(double)4.0 + 0.1*(double)N;
      else  _A[i][j] = (i+j)/(double)8.0;
     OutMatrix(_A);

     Fix_Leverrier(_A, xP);
     Newton(alfa, beta, xP, xZ);
     Polynom(alfa, beta, xP, xZ);

     printf("\n");
     for( i=0; i<Koren; i++)
      printf("\r\nX[\%d] = \%.10f", i, xZ[i]);

     getch();
     double MaxxZ = xZ[0];
     for( i=0; i<Koren; i++)
     if(MaxxZ < fabs(xZ[i])) MaxxZ = xZ[i];

     for( i=0; i<n1; i++)
      for( j=0; j<n1; j++)
      if(i==j)_A[i][j] = i/(double)4.0 + 0.1*(double)N - MaxxZ;
      else  _A[i][j] = (i+j)/(double)8.0;

//      OutMatrix(_A);

     MVGauss(_A, xG);
     printf("\n");
     for( i=0; i<n1; i++)
      printf("\r\nX[\%d] = \%.10f", i, xG[i]);


     getch();
}

Описан метод Леверьер (точный метод).
Здесь нужно изменить только ввод матрицы.
Так как всегда считается только для одной матрицы.