/********************************************************
* Calculate the coefficients of the characteristic poly-*
* nomial P(l)=A-l*I of a real tridiagonal matrix A(i,j) *
* ----------------------------------------------------- *
* 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:         *
*           4 1 0 0 0 0                                 *
*           2 4 1 0 0 0                                 *
*      A =  0 2 4 1 0 0                                 *
*           0 0 2 4 1 0                                 *
*           0 0 0 2 4 1                                 *
*           0 0 0 0 2 4   )                             *
*                                                       *
* Input matrix size (max 100): 6                        *
* Input 5 elements of subdiagonal:                      *
*   Element 1: 2                                        *
*   Element 2: 2                                        *
*   Element 3: 2                                        *
*   Element 4: 2                                        *
*   Element 5: 2                                        *
* Input 6 elements of main diagonal:                    *
*   Element 1: 4                                        *
*   Element 2: 4                                        *
*   Element 3: 4                                        *
*   Element 4: 4                                        *
*   Element 5: 4                                        *
*   Element 6: 4                                        *
* Input 5 elements of supradiagonal:                    *
*   Element 1: 1                                        *
*   Element 2: 1                                        *
*   Element 3: 1                                        *
*   Element 4: 1                                        *
*   Element 5: 1                                        *
*                                                       *
* Chracteristic polynomial P(l):                        *
* Coefficient for degree 6:      1.0000                 *
* Coefficient for degree 5:    -24.0000                 *
* Coefficient for degree 4:    230.0000                 *
* Coefficient for degree 3:  -1120.0000                 *
* Coefficient for degree 2:   2904.0000                 *
* Coefficient for degree 1:  -3776.0000                 *
* Coefficient for degree 0:   1912.0000                 *
*                                                       *
* (The characteristic polynomial of matrix A is:        *
*   P(l)=l^6-24*l^5+230*l^4-1120*l^3+2904*l^2-3376*l    *
*        +1912   )                                      *
*                                                       *
*                 C++ version by J-P Moreau, Paris      *
*               (with dynamically allocated vectors).   *
*                       (www.jpmoreau.fr)               *
********************************************************/
#include <basis.h>
#include <vmblock.h>

int   k,n;
REAL  *B,*C,*D,*E,*P;

void  *vmblock = NULL;  

/****************************************************
* The procedure PCTR calculates the coefficients of *
* the characteristic polynomial P(l)=A-l*I of a real*
* square tridiagonal matrix A(i,j).                 *
*                                                   *
* 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 must be declared    *
* and allocated in the calling program by using the *
* functions of module VMBLOCK.CPP.                  *
****************************************************/
void PCTR(int n, REAL *B, REAL *D, REAL *P) {
  int i,k; REAL *T;
  void  *vmblock1 = NULL; 
  vmblock1 = vminit();
  T = (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 Read_data(int n) {
  int i;
  printf("\n Input %d elements of subdiagonal:\n",n-1);
  for (i=1; i<n; i++) {
    printf("   Element %d: ", i); scanf("%lf",&C[i]);
  }
  printf("\n Input %d elements of main diagonal:\n",n);
  for (i=1; i<=n; i++) {
    printf("   Element %d: ", i); scanf("%lf",&D[i]);
  }
  printf("\n Input %d elements of supradiagonal:\n",n-1);
  for (i=1; i<n; i++) {
    printf("   Element %d: ", i); scanf("%lf",&E[i]);
  }
}


int main()  {

  //read size of matrix A(i,j)
  printf("\n Input matrix size: "); scanf("%d",&n);

  vmblock = vminit(); 

  //allocate memory for vectors of size n+1
  B       = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);
  C       = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);  
  D       = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);
  E       = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);  
  P       = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0); 

  if (! vmcomplete(vmblock))
  {
    LogError ("No Memory", 0, __FILE__, __LINE__);
    return 1;
  }

  //read elements of matrix A(i,j)
  Read_data(n);
  // prepare B vector
  for (k=1; k<n; k++) B[k]=C[k]*E[k];
  // call routine to calculate coefficients
  PCTR(n,B,D,P);
  // print results
  printf("\n Characteristic polynomial P(l):\n");
  for (k=1; k<=n+1; k++) 
    printf(" Coefficient for degree %d: %12.4f\n", n-k+1, P[k]);
  // exit section
  printf("\n\n");
  free(vmblock);
  return 0;
}

//end of file carpol.cpp