coefficient_matrix()

Per poter applicare il metodo di Gauss è necessario scrivere il sistema () nella forma $Ax = b$ dove $x$ è il vettore di incognite che corrisponde alla concentrazione al tempo $n+1$ ($C_{n+1}$), $b$ è il vettore dei coefficienti noti ed $A$ è la matrice dei coefficienti. La funzione coefficient_matrix() riempie la matrice $A$, ponendola uguale a $\{2\, Id +D\: (- B')\}$. La () ora assumerà la forma $AC_{n+1} = \{2\, Id + D\: B'\} C_n$.

Per risolvere $A = \{2\, Id +D\: (- B')\}$ il simulatore non esegue operazioni tra matrici (addizione, prodotto tra matrice e vettore, prodotto di una matrice per uno scalare e prodotto tra matrici), se avessi usato questo metodo, la sua complessità sarebbe $\Theta(n^3)$. Ho sfruttato il fatto che gli elementi della matrice $A$ sono per la maggior parte nulli.

La funzione è la seguente:

float **coefficient_matrix(int Dx, int Dy, double Dt,
                           Matrice *D) {
  float **A;
  int k;
  int dim = D -> r * D -> c;

1   if ((A = (float **)malloc(dim * sizeof(float *))) == NULL)
2     exit(1);

3   for (k = 0; k < dim; k++)
4     if ((A[k] = (float *)calloc(dim, sizeof(float))) == NULL) {
5       printf(``the memory is below to run fick.\n'');
6       exit(1);
7     }

8   for (k = 0; k < dim; k++) {
9     A[k][k] = 2*(1 + Lx(k) + Ly(k));   /* 2(1 + Lx + Ly) */

10    if (k - 1 >= 0 && (k-1)%D->c != D->c - 1)
11      A[k][k-1] = -1 * Lx(k);            /* -Lx */

12    if (k + 1 < dim && (k+1)%D->c != 0)
13      A[k][k+1] = -1 * Lx(k);            /* -Lx */

14    if (k - D->c >= 0)
15      A[k][k - D->c] = -1 * Ly(k);       /* -Ly */

16    if (k + D->c < dim)
17      A[k][k+D->c] = -1 * Ly(k);         /* -Ly */
18  }
19  return A;
20}

Dopo la dichiarazione di alcune variabili locali, nelle linee 1-7 si riserva la memoria per la matrice dinamica A. Usando la funzione calloc() la memoria viene automaticamente inizializzata a zero. Dato che la matrice ha al più cinque elementi diversi da zero per ogni riga, il ciclo (linee 8-18) ha il compito di scrivere tali valori. Le istruzioni if all' interno del ciclo garantiscono che l' area in esame non venga considerata cilindrica, cioè eliminano la possibilità che una zona situata al margine sia influenzata dalla concentrazione della zona del margine opposto.

La funzione ha complessità $\Theta(n)$ dove $n$ è il numero di righe della matrice.