// clang a.c -Wmost ../eigen/eigen.o ../../code/invert.o
// gcc -std=c99 a.c -Wmost ../eigen/eigen.o ../../code/invert.o
#include <stdlib.h>
#include <math.h>
#include <stdio.h>

#define n 2
typedef double R;
void eigen(int k, R m[k][k], R l[k], R vc[k][k]);
int slv(int i, int m, double A[i][i+m]);
#define F(j) for(int j=0; j<n; ++j)

void pv(R v[n]) {F(j) printf("%12.8f ", v[j]); printf("\n");}
void pm(R x[n][n]) {F(j) pv(x[j]);}

void sym(int m, R x[n][n]) {F(j) F(k)
    if(x[j][k] != x[k][j]) exit(printf("%d xx\n", m));}

typedef struct{R lambda[n]; R vec[n][n];} PY;  // output
PY le(R c[n][n], R d[n]){sym(1, c);
  R e[n]; F(j) e[j] = sqrt(d[j]); pv(e);
  R ece[n][n]; F(j) F(k) ece[j][k] = e[j]*c[j][k]*e[k];
  pm(c); sym(2, ece);
  PY a; R x[n][n]; eigen(n, ece, a.lambda, x); pm(x);
  F(j) F(k) a.vec[j][k] = e[k]*x[j][k];
  return a;}

int main(){
  R c[][n] = {{-2, 2}, {2, -3}}; R d[] = {0.1, 1./3}; R init[] = {8, 14};
 //  R c[][n] = {{-3, 4, 6}, {4, -1, -5}, {6, -5, -2}};
 //  R d[] = {5, 3.3, 2}; R init[] = {1, 1, 2};
  pm(c); PY py = le(c, d); pm(py.vec); pv(py.lambda);
  R sl[n][n+1]; F(j) {F(k) sl[j][k] = py.vec[k][j]; sl[j][n] = init[j];}
  slv(n, 1, sl);
  for(int t=0; t<21; ++t) {R time = t/10.;
    R xp[n]; F(j) xp[j] = sl[j][n]*exp(py.lambda[j]*time);
    printf("t = %5.1f", time);
    F(j) {R s = 0; F(k) s += py.vec[k][j]*xp[k]; printf("%12.8f ", s);}
    printf("\n");}
  return 0;}