RichardsFastSlow/src/math/algo.cpp

#include "algo.hpp"

void Thomas(size_t n,double* a,double* b,double* c,double* d,double* x){
  double last=c[0]=c[0]/b[0];
  for(size_t i=1;i<n-1;++i){
    last=c[i]=c[i]/(b[i]-a[i-1]*last);
  }
  last=d[0]=d[0]/b[0];
  for(size_t i=1;i<n;++i){
    last=d[i]=(d[i]-a[i-1]*last)/(b[i]-a[i-1]*c[i-1]);
  }
  last=x[n-1]=d[n-1];
  for(int i=n-1;i>0;--i){
    last=x[i-1]=d[i-1]-c[i-1]*last;
  }
}