LUsolve

La funzione LUsolve viene invocata solo quando si utilizza il metodo implicito, essa trova il vettore $C_{n+1}$ memorizzandolo in v, combinando le sostituzioni in avanti e a ritroso (par. ). Si ricava $y$ usando la sostituzione in avanti nelle linee 4-9, e poi trova $v$ usando la sostituzione all' indietro nelle linee 15-21. Nella sostituzione all' indietro le linee 15-18 si occupano di porre il valore della concentrazione a $255$ in caso questo lo superi, altrimenti nelle linee 19-20 si approssima il valore ottenuto dalla sostituzione a ritroso e si memorizza in $v$.

void LUsolve(float *LU, float *b, int dim, Uint8 *v) {
  int i, j;
  float *y, *vett, sum;

1   if (((y = (float *)malloc(dim*sizeof(float))) == NULL) ||
2       ((vett = (float *)malloc(dim*sizeof(float))) == NULL))
3     exit(1);

4   for (i = 0; i < dim; i++) {
5     sum = 0;
6     for(j = 0; j < i; j++)
7       sum += LU[indd(i,j)] * y[j];
8     y[i] = b[i] - sum;
9   }

10  for (i = dim - 1; i >= 0; i--) {
11    sum = 0;
12    for(j = i + 1; j < dim; j++)
13      sum += LU[indd(i,j)] * vett[j];
14    vett[i] = ((y[i] - sum) / LU[indd(i,i)]);

15    if (((y[i] - sum) / LU[indd(i,i)]) > 255) {
16      printf("%f\n", (y[i] - sum) / LU[indd(i,i)] );
17      vett[i] = 255;
18    }
19    else vett[i] = ((y[i] - sum) / LU[indd(i,i)]) + 0.5;
20    v[i] = (Uint8)vett[i];
21  }
22  free(y);
23  free(vett);
24}

Il tempo di esecuzione è $\Theta(n^2 + n^2) = \Theta(n^2)$.