/********************************************************
*   Calculate the coefficients of the characteristic    *
*   polynomial P(l)=A-l*I of a square symmetric real    *
*         matrix A(i,j) by Lanczos's method             *
* ----------------------------------------------------- *
* Ref.: "Algèbre, Algorithmes et programmes en Pascal   *
*        By Jean-Louis Jardrin, DUNOD Paris, 1988"      *
*        [BIBLI 10].                                    *
* ----------------------------------------------------- *
* SAMPLE RUN:                                           *
* (Find the caracteristic polynomial of matrix:         *
*          10  7   8   7                                *
*      A =  7  5   6   5                                *
*           8  6  10   9                                *
*           7  5   9  10                                *
*                                                       *
* Input matrix size (max 50): 4                         *
*  Line 1:                                              *
*   Element 1: 10                                       *
*   Element 2: 7                                        *
*   Element 3: 8                                        *
*   Element 4: 7                                        *
*  Line 2:                                              *
*   Element 2: 5                                        *
*   Element 3: 6                                        *
*   Element 4: 5                                        *
*  Line 3:                                              *
*   Element 3: 10                                       *
*   Element 4: 9                                        *
*  Line 4:                                              *
*   Element 4: 10                                       *
*                                                       *
* Chracteristic polynomial P(l):                        *
* Coefficient for degree 4: 1.000000                    *
* Coefficient for degree 3: -35.000000                  *
* Coefficient for degree 2: 146.000000                  *
* Coefficient for degree 1: -100.000000                 *
* Coefficient for degree 0: 1.000000                    *
*                                                       *
* (The characteristic polynomial of matrix A is:        *
*   P(l)=l^4-35*l^3+146*l^2-100*l+1 )                   *
*                                                       *
*                   C++ version by J-P Moreau, Paris    *
*  (with dynamically allocated vectors and matrices).   *
*                          (www.jpmoreau.fr)            *
********************************************************/
#include <basis.h>
#include <vmblock.h>

int  it,k,n;
LONG_REAL eps;
LONG_REAL **A, **Tau;
LONG_REAL *D,*E,*P;

void *vmblock = NULL; 


//read symmetric matrix A from screen
void  Read_data(int n) {
  int i,j;
  for (i=1; i<=n; i++) {
    printf("\n Line %d\n", i);
    for (j=i; j<=n; j++) {
      printf("  Element %d: ", j); scanf("%lf",&A[i][j]);
	  A[j][i] = A[i][j];
    }
  }
}

/********************************************************
* The TDSL procedure tridiagonalizes the symmetric real *
* matrix A(n,n) by Lanczos's method. Only the main dia- *
* gonal terms and the supradiagonal terms are calculated*
* ----------------------------------------------------- *
* Inputs:                                               *
*          eps  precision (double)                      *
*          n    size of square A matrix (integer)       *
*          A    square real symmetric matrix (n,n)      *
* Outputs:                                              *
*          it   flag=0 if method failed, else =1        *
*          D    vector of size n+1 storing the n        *
*               elements of the main diagonal.          *
*          E    vector of size n+1 storing the n-1      *
*               elements of the supradiagonal.          *
*          Tau  Utility matrix used by the method.      *
*                                                       *
*  The pointers to a real vector or a real matrix must  *
*  be declared and allocated in the calling program by  *
*  using the functions of module VMBLOCK.CPP.           *
********************************************************/
void TDSL(double eps,int n, LONG_REAL **A, int *it,LONG_REAL *D,
		  LONG_REAL *E, LONG_REAL **Tau)   {
  int i,j,k; LONG_REAL s,u0,v0; LONG_REAL *W;
  void *vmblock1 = NULL;
  //initialize local vector W(n+1)
  vmblock1 = vminit();
  W = (LONG_REAL *) vmalloc(vmblock1, VEKTOR,  n+1, 0);
  Tau[1][1]=1.0; for (i=2; i<=n; i++)  Tau[i][1]=0.0;
  u0=1.0; *it=1; k=1;
  do {
    for (i=1; i<=n; i++) {
      s=0.0;
      for (j=1; j<=n; j++)  s += A[i][j]*Tau[j][k];
      W[i]=s;
    }
    s=0.0; for (i=1; i<=n; i++)  s += Tau[i][k]*W[i];
    D[k]=s/u0;
    if (k!=n) {
      for (i=1; i<=n; i++) {
        if (k==1)
          Tau[i][k+1]=W[i]-D[k]*Tau[i][k];
        else
          Tau[i][k+1]=W[i]-D[k]*Tau[i][k]-E[k-1]*Tau[i][k-1];
		printf(" i=%d k+1=%d Tau=%f\n", i, k+1, Tau[i][k+1]); 
      }
	  v0=0.0; for (i=1; i<=n; i++)  v0 += (Tau[i][k+1]*Tau[i][k+1]);
      if (v0<eps)
        *it=0;
      else {
        E[k]=v0/u0; u0=v0;
      }
    } //if k 
    k++;
  } while(k<=n && *it==1);
  printf("\n *** it=%d\n",*it);
  free(vmblock1);
}

/****************************************************
* The procedure PCTR calculates the coefficients of *
* the characteristic polynomial P(l)=A-l*I of a real*
* square tridiagonal matrix A(i,j) by Souriau's     *
* method.                                           *
*                                                   *
* Note: the roots of P(l) are the eigenvalues of    *
*       matrix A(i,j).                              *
* ------------------------------------------------- *
* INPUTS:                                           *
*          n: zize of matrix (integer)              *
*          B: vector of codiagonal terms products   *
*             (see main program).                   *
*          D: vector of main diagonal terms         *
* OUTPUT:                                           *
*          P: vector of coefficients of P(l)        *
* ------------------------------------------------- *
*  The pointers to a real vector or a real matrix   *
*  must be declared and allocated in the calling    *
*  program by using the functions of module         *
*  VMBLOCK.CPP.                                     *
*                                                   *
*  For a demonstration of function PCTR alone, see  *
*  program carpol.cpp.                              *
****************************************************/
void PCTR(int n, LONG_REAL *B, LONG_REAL *D, LONG_REAL *P) {
  int i,k; LONG_REAL *T;
  void  *vmblock1 = NULL; 
  vmblock1 = vminit();
  T = (LONG_REAL *) vmalloc(vmblock1, VEKTOR, n, 0);
  P[2]=D[1]; T[1]=1.0; P[1]=-T[1];
  for (k=2; k<=n; k++) {
    P[k+1]=D[k]*P[k]-B[k-1]*T[k-1];
    for (i=k; i>2; i--) {
      T[i]=P[i];
      P[i]=-T[i]+D[k]*P[i-1]-B[k-1]*T[i-2];
    }
    T[2]=P[2]; P[2]=-T[2]+D[k]*P[1];
    T[1]=P[1]; P[1]=-T[1];
  }
  free(vmblock1);
}


void main() {

  //read size n of matrix A(n,n)
  printf("\n Input matrix size: "); scanf("%d", &n);
  //allocate matrices and vectors
  vmblock = vminit();
  A   = (LONG_REAL **)vmalloc(vmblock, MATRIX,  n+1, n+1);
  Tau = (LONG_REAL **)vmalloc(vmblock, MATRIX,  n+1, n+1);
  D   = (LONG_REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);
  E   = (LONG_REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);
  P   = (LONG_REAL *) vmalloc(vmblock, VEKTOR,  n+2, 0);
  //read A matrix from screen
  Read_data(n);
  //fix precision
  eps=1e-10;
  //tridiagonalize symmetric A matrix by Lanczos's method
  TDSL(eps,n,A,&it,D,E,Tau);
  //if success, use method for a tridiagonal matrix
  if (it==1)  PCTR(n,E,D,P);
  //print results
  switch(it) {
   case 0: 
	 printf("\n Method failed.\n"); break;
   case 1:
     printf("\n Characteristic polynomial P(l):\n");
     for (k=1; k<n+2; k++)
       printf(" Coefficient for degree %d: %f\n", n-k+1, P[k]);
  }
  //exit section
  printf("\n\n");
  free(vmblock);

}

// end of file carpol3.cpp