// gcc lop.c -Wmost j.o
#include <assert.h>
#include <math.h>
#include <stdio.h>
typedef double R;
void pm(int n, char * m, R a[n][n]){
   int i; for(i=0; i<n; ++i) {int j; for(j=0; j<n; ++j) printf("%e ", a[i][j]);
     printf(" %s\n", m);}}
#define db 1
#if db
// Eigen test is m(v_i)*vec_i = val_i*vec_i ?
void evt(int n, R m[n][n], R val[n], R vec[n][n])
  {int j=n; while(j--){R u[n]; // check eigenvector j
  {int c=n; while(c--) { // which component of vector
  {int y=n; while(y--) {R s=0; {int q=n; while(q--) s += m[q][c]*vec[j][q];} u[c]=s;}}}}
  {int c=n; while(c--) assert(fabs(val[j]*vec[j][c] - u[c]) < 1.e-6);}}}
#endif

void eigen(int n, R m[n][n], R l[n], R vc[n][n]);
void lev(int n, R C[n][n], R d[n], R eval[n], R g[n][n]){
pm(n, "C", C);
if(1){int j=n; while(j--){int k=n; while(k--) assert(C[j][k]==C[k][j]);}}
R ECE[n][n]; R dE[n];
{int j=n; while(j--) dE[j] = sqrt(d[j]);}
{int j=n; while(j--){int k=n; while(k--) ECE[j][k] = C[j][k]*dE[j]*dE[k];}}
pm(n, "ECE", ECE);
eigen(n, ECE, eval, g); pm(n, "g", g); if(1) evt(n, ECE, eval, g);
{int j; for(j=0; j<n; ++j) printf("%e ", eval[j]); printf("eval\n");}
 {int j=n; while(j--){int c=n; while(c--) g[j][c] *= dE[j];}}}

#if 1
int main(){R c[][2] = {{1, 1}, {1, 2}}; R C[] = {2, 1};
   R vc[2][2]; R vl[2]; lev(2, c, C, vl, vc); printf("Yes!\n");
   {R DC[2][2];
    {int j=2; while(j--){int k=2; while(k--) DC[j][k] = C[j]*c[j][k];}}
    pm(2, "DC", DC); pm(2, "vc", vc);
    evt(2, DC, vl, vc); printf("Yum!\n");} return 0;}
#endif
#if 0
int main(){R c[][3] = {{2., -3.4, 3}, {-3.4, 2.1, 1.2}, {3., 1.2, -1.7}};
  R C[] = {2.1, 1.2, 3.3};
  R vc[3][3]; R vl[3]; lev(3, c, C, vl, vc); printf("Yes!\n");
  {R DC[3][3];
    {int j=3; while(j--){int k=3; while(k--) DC[j][k] = C[j]*c[j][k];}}
    pm(3, "DC", DC); pm(3, "vc", vc);
    evt(3, DC, vl, vc); printf("Yum!\n");}
  return 0;}
#endif