// Dieses Programm integriert die Diffusionsgleichung mittels
// ADI(Alternativ Direction Implicit)-Methode und gibt die Ergebnisse in
// Form von PNG-Grafiken aus, aus denen sich z.B. mit Programmen wie
// ffmpeg kleine Filme erzeugen lassen. Die Anzahl der ausgegebenen
// Bilder läßt sich über die Variable "fps" steuern (siehe unten). Für
// die lineare Algebra wird dabei die Bibliothek GLS (Gnu
// Scientific-Library) verwendet, weswegen mit "-ldislin -lgsl
// -lgslcblas" gelinkt werden muss.

#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

// Die folgende Zeile muss evtl. an Eure lokale Dislin-Installation
// angepasst werden.
#include "/usr/local/dislin/dislin.h"

// Eine nützliche Definition zum Runden
#define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))

// Zeitintervall [t_0,t_0]
double t_0 = 0.0;
double t_1 = 100.0;
// Raumbereich [x_0,x_1]x[y_0,y_1]
double x_0 = 0.0;
double x_1 = 1.0;
double y_0 = 0.0;
double y_1 = 1.0;

// Raum- und Zeitschritt. Aus diesen Größen wird weiter unten die
// Anzahl der Schritte berechnet.
double dx = 0.5;
double dy = 0.5;
double dt = 0.005;

// Diffusionskonstanten
double Diff_u = 5.0;
double Diff_v = 12.0;
// Reaktionskonstanten
//double a = 4.0;
//double b = 15.59;
static const double sigma = Diff_u / Diff_v;
static const double a = 2.0 * sqrt(sigma) / (1.0 - sigma);
tatic const double b = 1 + pow(sigma,2);

// fps ist die Anzahl der Bilder, die pro simulierter Sekunde ausgegeben werden sollen
int fps = 2;

void print_matrix(gsl_matrix *A)
{
  int i = 0;
  int j = 0;

  for (i = 0; i<A->size1; i++) {
    for (j = 0; j<A->size2; j++) {
      printf("%.3f ", gsl_matrix_get(A, i, j));
    }
    printf("\n");
  }
}

// Die Funktionen startwert und startableitung werden zu Beginn des
// Programms verwendet, um die Anfangsbedingung festzulegen.
double startwert(){
  return 0.00002 * (rand() % 10001);
}

double f(double u, double v) {
  return a - (b + 1)*u + u*u*v;
}

double g(double u, double v) {
  return b*u - u*u*v;
}

int main () {
  // Die Anzahl der räumlichen und zeitlichen Iterationsschritte wird
  // aus der Größe der entsprechenden Intervalle und den Schrittweiten
  // berechnet.
  int x_steps = round( (x_1 - x_0) / dx); 
  int y_steps = round( (y_1 - y_0) / dy);
  int t_steps = round( (t_1 - t_0) / dt);

  printf("\nx_steps=%d y_steps=%d t_steps=%d\n",x_steps, y_steps, t_steps);

  // Die Variablen fname und titel werden für die Ausgabe mit Dislin
  // benötigt (siehe unten).
  char fname[80];
  char titel[80];

  // Je nach Wahl der globalen Variablen fps und dt wird die Zahl der
  // zeitlichen Iterationen berechnet, die zwischen zwei Ausgaben
  // durchgeführt werden sollen.
  int steps_pro_frame = round(1.0 / fps / dt);
  if ( !steps_pro_frame ) {
    steps_pro_frame = 1;
  }
  
  // Die Zeitintegration erfolgt durch Lösung eines linearen
  // Gleichungssystems. Im Folgenden werden die dafür benötigten Vektoren und
  // eine Matrix definiert.
  int vektorlaenge = (x_steps + 1)*(y_steps + 1);

  gsl_vector *u       = gsl_vector_alloc(vektorlaenge);
  gsl_vector *u_halb  = gsl_vector_alloc(vektorlaenge);
  gsl_vector *v       = gsl_vector_alloc(vektorlaenge);
  gsl_vector *v_halb  = gsl_vector_alloc(vektorlaenge);
  gsl_vector *temp    = gsl_vector_alloc(vektorlaenge);
  gsl_vector *gamma   = gsl_vector_alloc(vektorlaenge);
  gsl_matrix *A_u     = gsl_matrix_alloc(vektorlaenge, vektorlaenge);
  gsl_matrix *B_u     = gsl_matrix_alloc(vektorlaenge, vektorlaenge);
  gsl_matrix *C_u     = gsl_matrix_alloc(vektorlaenge, vektorlaenge);
  gsl_matrix *D_u     = gsl_matrix_alloc(vektorlaenge, vektorlaenge);
  gsl_matrix *A_v     = gsl_matrix_alloc(vektorlaenge, vektorlaenge);
  gsl_matrix *B_v     = gsl_matrix_alloc(vektorlaenge, vektorlaenge);
  gsl_matrix *C_v     = gsl_matrix_alloc(vektorlaenge, vektorlaenge);
  gsl_matrix *D_v     = gsl_matrix_alloc(vektorlaenge, vektorlaenge);
    
  gsl_vector *gamma_u = gsl_vector_alloc(vektorlaenge);
  gsl_vector *gamma_v = gsl_vector_alloc(vektorlaenge);

  // Die folgenden drei Größen werden gsl-intern für die LU-Zerlegung
  // benötigt.
  gsl_permutation *p_A_u = gsl_permutation_alloc(vektorlaenge);
  gsl_permutation *p_C_u = gsl_permutation_alloc(vektorlaenge);
  gsl_permutation *p_A_v = gsl_permutation_alloc(vektorlaenge);
  gsl_permutation *p_C_v = gsl_permutation_alloc(vektorlaenge);
  int sign = 0;

  int i = 0; // i ist der x-Index
  int j = 0; // j ist der y-Index
  int k = 0; // k ist der Zeitindex

  // Division und Potenzierung ist relativ teuer, sollte deshalb nicht
  // öfter als nötig berechnet werden!
  double alpha_u = Diff_u * dt / 2.0 / dx / dx;
  double beta_u  = Diff_u * dt / 2.0 / dy / dy;
  double alpha_v = Diff_v * dt / 2.0 / dx / dx;
  double beta_v  = Diff_v * dt / 2.0 / dy / dy;


  // Das Feld wird entsprechend der Anfangsbedingung mit startwert()
  // initialisiert.
  for (i = 0; i < vektorlaenge; i++) {
      gsl_vector_set(u, i, a + startwert());
      gsl_vector_set(v, i, b / a + startwert());
  }

  // Erstelle die Matrix A_u
  for ( i=0; i<vektorlaenge; i++ ) {
    gsl_matrix_set(A_u, i, i, 1.0+2.0*alpha_u);

    if ( (i % (x_steps + 1) != 0) && (i % (x_steps + 1) != x_steps) ) {
      gsl_matrix_set(A_u, i, i-1, -alpha_u);
      gsl_matrix_set(A_u, i, i+1, -alpha_u);
    } else if ( i % (x_steps + 1) == 0 ) {
      gsl_matrix_set(A_u, i, i+1, -2.0*alpha_u);
    } else {
      gsl_matrix_set(A_u, i, i-1, -2.0*alpha_u);
    } 
  }

  printf("==== Matrix A_u ====\n");
  print_matrix(A_u);

  // Erstelle die Matrix B_u
  for ( i=0; i<vektorlaenge; i++ ) {
    gsl_matrix_set(B_u, i, i, 1.0-2.0*beta_u);

    if ( (i > x_steps ) && (i < vektorlaenge - (x_steps + 1)) ) {
      gsl_matrix_set(B_u, i, i+x_steps+1, beta_u);
      gsl_matrix_set(B_u, i, i-x_steps-1, beta_u);
    } else if ( i <= x_steps ) {
      gsl_matrix_set(B_u, i, i+x_steps+1, 2.0*beta_u);
    } else {
      gsl_matrix_set(B_u, i, i-x_steps-1, 2.0*beta_u);
    } 
  }

  printf("==== Matrix B_u ====\n");
  print_matrix(B_u);

  // Erstelle die Matrix C_u
  for ( i=0; i<vektorlaenge; i++ ) {
    gsl_matrix_set(C_u, i, i, 1.0+2.0*beta_u);

    if ( (i > x_steps ) && (i < vektorlaenge - (x_steps + 1)) ) {
      gsl_matrix_set(C_u, i, i+x_steps+1, -beta_u);
      gsl_matrix_set(C_u, i, i-x_steps-1, -beta_u);
    } else if ( i <= x_steps ) {
      gsl_matrix_set(C_u, i, i+x_steps+1, -2.0*beta_u);
    } else {
      gsl_matrix_set(C_u, i, i-x_steps-1, -2.0*beta_u);
    } 
  }

  printf("==== Matrix C_u ====\n");
  print_matrix(C_u);

  // Erstelle die Matrix D_u
  for ( i=0; i<vektorlaenge; i++ ) {
    gsl_matrix_set(D_u, i, i, 1.0-2.0*alpha_u);

    if ( (i % (x_steps + 1) != 0) && (i % (x_steps + 1) != x_steps) ) {
      gsl_matrix_set(D_u, i, i-1, alpha_u);
      gsl_matrix_set(D_u, i, i+1, alpha_u);
    } else if ( i % (x_steps + 1) == 0 ) {
      gsl_matrix_set(D_u, i, i+1, 2.0*alpha_u);
    } else {
      gsl_matrix_set(D_u, i, i-1, 2.0*alpha_u);
    } 
  }

  printf("==== Matrix D_u ====\n");
  print_matrix(D_u);

  // Erstelle die Matrix A_v
  for ( i=0; i<vektorlaenge; i++ ) {
    gsl_matrix_set(A_v, i, i, 1.0+2.0*alpha_v);

    if ( (i % (x_steps + 1) != 0) && (i % (x_steps + 1) != x_steps) ) {
      gsl_matrix_set(A_v, i, i-1, -alpha_v);
      gsl_matrix_set(A_v, i, i+1, -alpha_v);
    } else if ( i % (x_steps + 1) == 0 ) {
      gsl_matrix_set(A_v, i, i+1, -2.0*alpha_v);
    } else {
      gsl_matrix_set(A_v, i, i-1, -2.0*alpha_v);
    } 
  }

  printf("==== Matrix A_v ====\n");
  print_matrix(A_v);

  // Erstelle die Matrix B_v
  for ( i=0; i<vektorlaenge; i++ ) {
    gsl_matrix_set(B_v, i, i, 1.0-2.0*beta_v);

    if ( (i > x_steps ) && (i < vektorlaenge - (x_steps + 1)) ) {
      gsl_matrix_set(B_v, i, i+x_steps+1, beta_v);
      gsl_matrix_set(B_v, i, i-x_steps-1, beta_v);
    } else if ( i <= x_steps ) {
      gsl_matrix_set(B_v, i, i+x_steps+1, 2.0*beta_v);
    } else {
      gsl_matrix_set(B_v, i, i-x_steps-1, 2.0*beta_v);
    } 
  }

  printf("==== Matrix B_v ====\n");
  print_matrix(B_v);

 // Erstelle die Matrix C_v
  for ( i=0; i<vektorlaenge; i++ ) {
    gsl_matrix_set(C_v, i, i, 1.0+2.0*beta_v);

    if ( (i > x_steps ) && (i < vektorlaenge - (x_steps + 1)) ) {
      gsl_matrix_set(C_v, i, i+x_steps+1, -beta_v);
      gsl_matrix_set(C_v, i, i-x_steps-1, -beta_v);
    } else if ( i <= x_steps ) {
      gsl_matrix_set(C_v, i, i+x_steps+1, -2.0*beta_v);
    } else {
      gsl_matrix_set(C_v, i, i-x_steps-1, -2.0*beta_v);
    } 
  }

  printf("==== Matrix C_v ====\n");
  print_matrix(C_v);

  // Erstelle die Matrix D_v
  for ( i=0; i<vektorlaenge; i++ ) {
    gsl_matrix_set(D_v, i, i, 1.0-2.0*alpha_v);

    if ( (i % (x_steps + 1) != 0) && (i % (x_steps + 1) != x_steps) ) {
      gsl_matrix_set(D_v, i, i-1, alpha_v);
      gsl_matrix_set(D_v, i, i+1, alpha_v);
    } else if ( i % (x_steps + 1) == 0 ) {
      gsl_matrix_set(D_v, i, i+1, 2.0*alpha_v);
    } else {
      gsl_matrix_set(D_v, i, i-1, 2.0*alpha_v);
    } 
  }

  printf("==== Matrix D_v ====\n");
  print_matrix(D_v);

  // Führe zunächst die LU-Zerlegung der Matrizen A und C durch. Dieser
  // Teil des Programms nimmt die meiste Zeit in Anspruch.
  printf("\nLU-Zerlegung von A_u beginnt!\n");
  gsl_linalg_LU_decomp(A_u, p_A_u, &sign);
  printf("\nLU-Zerlegung von C_u beginnt!\n");
  gsl_linalg_LU_decomp(C_u, p_C_u, &sign);
  printf("\nLU-Zerlegung von A_v beginnt!\n");
  gsl_linalg_LU_decomp(A_v, p_A_v, &sign);
  printf("\nLU-Zerlegung von C_v beginnt!\n");
  gsl_linalg_LU_decomp(C_v, p_C_v, &sign);
  printf("\nLU-Zerlegung beendet!\n");

  for ( k = 0; k < t_steps; k++ ) {

    //---{ Ausgabe mit Dislin }-----------------------------------------
    if ( k % steps_pro_frame == 0) {
      scrmod("REVERSE");
      metafl("PNG");
      filmod("DELETE");
      sprintf(fname, "ADI%d.png", round(k / steps_pro_frame) );
      setfil(fname);
      disini();
      texmod("ON");
      
      sprintf(titel, "%dx%d Stuetzstellen, $D_u=%.1f, D_v=%.1f, a=%.2f, b=%.2f, t=%.3f$",x_steps + 1, y_steps + 1, Diff_u, Diff_v, a, b, t_0 + (k + 1) * dt);
      titlin(titel, 4);

      graf3(x_0, x_1, x_0, (x_1 - x_0)/10.0, y_0, y_1, y_0, (y_1 -
            y_0)/10.0, 1.0, 4.5, 0.0, 1.0);
      title();

      double masse = 0.0;
      for (i = 0; i<vektorlaenge; i++) {
        masse += gsl_vector_get(u, i); 
      }
      printf("\n Masse: %f",masse);

      float darst[y_steps + 1][x_steps + 1];
      for (i = 0; i<=y_steps; i++) {
        for (j = 0; j<=x_steps; j++) {
	        darst[i][j] = gsl_vector_get(u, i * (x_steps+1) + j);
	      }
      }

      crvmat((float *) darst, y_steps + 1, x_steps + 1, 25, 25);
      //surmat((float *) darst, y_steps + 1, x_steps + 1, 3, 3);


      disfin();
    }
    //------------------------------------------------------------------

    //---{ Der eigentliche Zeitschritt }-------------------------------
    for (i=0; i<vektorlaenge; i++) {
      gsl_vector_set(gamma_u, i, f(gsl_vector_get(u, i),
            gsl_vector_get(v, i)));
      gsl_vector_set(gamma_v, i, g(gsl_vector_get(u, i),
            gsl_vector_get(v, i)));
    }
    gsl_vector_scale(gamma_u, dt);
    gsl_vector_scale(gamma_v, dt);

    gsl_blas_dgemv(CblasNoTrans, 1.0, B_u, u, 0.0, temp);
    gsl_linalg_LU_solve(A_u, p_A_u, temp, u_halb);

    gsl_blas_dgemv(CblasNoTrans, 1.0, D_u, u_halb, 0.0, temp);
    gsl_linalg_LU_solve(C_u, p_C_u, temp, u);

    gsl_blas_dgemv(CblasNoTrans, 1.0, B_v, v, 0.0, temp);
    gsl_linalg_LU_solve(A_v, p_A_v, temp, v_halb);

    gsl_blas_dgemv(CblasNoTrans, 1.0, D_v, v_halb, 0.0, temp);
    gsl_linalg_LU_solve(C_v, p_C_v, temp, v);

    gsl_vector_add(u, gamma_u);
    gsl_vector_add(v, gamma_v);
    //-----------------------------------------------------------------
  }
  
  // Zuletzt kann der reservierte Speicher wieder freigegeben werden.
  gsl_vector_free(u);
  gsl_vector_free(u_halb);
  gsl_vector_free(v);
  gsl_vector_free(v_halb);
  gsl_vector_free(temp);
  gsl_matrix_free(A_u);
  gsl_matrix_free(B_u);
  gsl_matrix_free(C_u);
  gsl_matrix_free(D_u);
  gsl_matrix_free(A_v);
  gsl_matrix_free(B_v);
  gsl_matrix_free(C_v);
  gsl_matrix_free(D_v);

  return 1;
}