Solver ( Author: }:+()___ [Smile] )

solver.cpp

// Author: }:+()___ [Smile]
#include <memory>
#include <cassert>


class Solver
{
  int mat_size, max_coeff, cur_row, cur_coeff;
  std::unique_ptr<int []> row_start, col_index;
  std::unique_ptr<double []> mat, residue, result;
  std::unique_ptr<double []> diag, scale, tmp, tmp1, bcg_p, bcg_q;
  double qqi;
  double prescale()
  {
    memset(diag.get(), 0, mat_size * sizeof(double));
    memset(tmp1.get(), 0, mat_size * sizeof(double));
    memset(bcg_p.get(), 0, mat_size * sizeof(double));
    memset(bcg_q.get(), 0, mat_size * sizeof(double));

    for(int i = 0; i < mat_size; i++)
      for(int j = row_start[i]; j < row_start[i + 1]; j++)
        if(col_index[j] > i)tmp1[i] -= mat[j];
        else if(col_index[j] == i)diag[i] += mat[j];
    for(int i = 0; i < mat_size; i++)
    {
      tmp[i] = diag[i];
      for(int j = row_start[i]; j < row_start[i + 1]; j++)
        if(col_index[j] < i)tmp[i] += mat[j] * tmp[col_index[j]] * tmp1[col_index[j]];
      scale[i] = sqrt(tmp[i] = 1 / tmp[i]);
    }

    for(int i = 0; i < mat_size; i++)
    {
      diag[i] = tmp[i] * diag[i] - 2;
      for(int j = row_start[i]; j < row_start[i + 1]; j++)
        if(col_index[j] != i)mat[j] *= scale[i] * scale[col_index[j]];
    }
    for(int i = mat_size - 1; i != -1; i--)
    {
      result[i] = tmp1[i] = result[i] / scale[i];
      for(int j = row_start[i]; j < row_start[i + 1]; j++)
        if(col_index[j] > i)result[i] += mat[j] * tmp1[col_index[j]];
    }

    double rr = 0;
    for(int i = 0; i < mat_size; i++)
    {
      tmp[i] = residue[i] * scale[i] - diag[i] * tmp1[i] - result[i];
      for(int j = row_start[i]; j < row_start[i + 1]; j++)
        if(col_index[j] < i)tmp[i] -= mat[j] * tmp[col_index[j]];
      residue[i] = tmp[i] - tmp1[i];  rr += residue[i] * residue[i];
    }
    qqi = 0;  return rr;
  }
  double step()
  {
    for(int i = mat_size - 1; i != -1; i--)
    {
      tmp[i] = residue[i];
      for(int j = row_start[i]; j < row_start[i + 1]; j++)
        if(col_index[j] > i)tmp[i] -= mat[j] * tmp[col_index[j]];
    }

    double tq = 0;
    for(int i = 0; i < mat_size; i++)
    {
      tmp1[i] = residue[i] + diag[i] * tmp[i];
      for(int j = row_start[i]; j < row_start[i + 1]; j++)
        if(col_index[j] < i)tmp1[i] -= mat[j] * tmp1[col_index[j]];
      tmp[i] += tmp1[i];  tq += tmp[i] * bcg_q[i];
    }

    double bk = tq * qqi, rq = 0, qq = 0;
    for(int i = 0; i < mat_size; i++)
    {
      bcg_q[i] = tmp[i] - bk * bcg_q[i];  bcg_p[i] = residue[i] - bk * bcg_p[i];
      rq += residue[i] * bcg_q[i];  qq += bcg_q[i] * bcg_q[i];
    }

    qqi = 1 / qq;  double ak = rq * qqi, rr = 0;
    for(int i = 0; i < mat_size; i++)
    {
      result[i] += ak * bcg_p[i];  residue[i] -= ak * bcg_q[i];  rr += residue[i] * residue[i];
    }
    return rr;
  }
  void postscale()
  {
    for(int i = mat_size - 1; i != -1; i--)
    {
      tmp1[i] = result[i];
      for(int j = row_start[i]; j < row_start[i + 1]; j++)
        if(col_index[j] > i)tmp1[i] -= mat[j] * tmp1[col_index[j]];
      result[i] = scale[i] * tmp1[i];
    }
  }
public:
  Solver(int n, int nc):mat_size(n),max_coeff(nc),cur_row(0),cur_coeff(0),
    row_start(new int[n + 1]), col_index(new int[nc]), mat(new double[nc]), residue(new double[n]), result(new double[n]),
    diag(new double[n]), scale(new double[n]), tmp(new double[n]), tmp1(new double[n]), bcg_p(new double[n]), bcg_q(new double[n])
  {
    row_start[0] = 0;
  }

  void reset(){cur_row = cur_coeff = 0;}

  void next_row()
  {
    assert(cur_row < mat_size);
    row_start[++cur_row] = cur_coeff;
  }

  void add_coeff(int col, double val)
  {
    assert(cur_coeff < max_coeff);
    assert(col >= 0 && col < mat_size);
    col_index[cur_coeff] = col;  mat[cur_coeff++] = val;
  }

  double *right_part(){return residue.get();}

  double &right_part(int row)
  {
    assert(row >= 0 && row < mat_size);
    return residue[row];
  }

  void set_right_part(int row, double val)
  {
    assert(row >= 0 && row < mat_size);
    residue[row] = val;
  }

  double *initial(){return result.get();}

  double &initial(int row)
  {
    assert(row >= 0 && row < mat_size);
    return result[row];
  }

  void set_initial(int row, double val)
  {
    assert(row >= 0 && row < mat_size);
    result[row] = val;
  }

  const double *solve(double eps, int max_iter)
  {
    assert(cur_row == mat_size);  cur_row = cur_coeff = 0;
    double err = prescale();  if(!(err >= 0))return nullptr;  // bad matrix

    eps *= eps;  eps *= err;
    for(int i = 0; i < max_iter; i++)
    {
      if(err < eps)
      {
        postscale();  return result.get();
      }
      err = step();
    }
    return nullptr;  // bad convergence
  }
  int matrix_size()const{return mat_size;}
};

bool test_solver()
{
  const int n = 16;
  Solver solver(n, 3 * n);
  for(int i = 0; i < n; i++)
  {
    solver.add_coeff(i, 2);
    if(i)solver.add_coeff(i - 1, -1);
    if(i + 1 < n)solver.add_coeff(i + 1, -1);

    solver.set_right_part(i, 1);
    solver.set_initial(i, 0);
    solver.next_row();
  }

  const double *res = solver.solve(1e-6, 2);
  if(!res)return false;

  vector<string> out;
  for(int i = 0; i < n; i++)qap_add_back(out)=FToS(res[i]);
  auto resstr=join(out," ")+"\n";

  for(int i = 0; i < n; i++)
  {
    double test = 2 * res[i];
    if(i)test -= res[i - 1];
    if(i + 1 < n)test -= res[i + 1];
    resstr+=FToS(test)+" ";
  }
  int gg=1;
  return true;
}