<commit_before>/** \addtogroup examples 
  * @{ 
  * \defgroup algebraic_multigrid algebraic_multigrid
  * @{ 
  * \brief Benchmark for smoothed algebraic multgrid
  */
#include <ctf.hpp>
using namespace CTF;

void smooth_jacobi(Matrix<> & A, Vector<> & x, Vector <> & b, int nsmooth){
  Vector<> d(x.len, *x.wrld);
  d["i"] = A["ii"];
  Transform<>([](double & d){ d= fabs(d) > 0.0 ? 1./d : 0.0; })(d["i"]);
  Matrix<> R(A);
  R["ii"] = 0.0;

  //20 iterations of Jacobi, should probably be a parameter or some convergence check instead
  for (int i=0; i<nsmooth; i++){
    x["i"] = -1.0*R["ij"]*x["j"];
    x["i"] += b["i"];
    x["i"] *= d["i"];
  }
}

void vcycle(Matrix<> & A, Vector<> & x, Vector<> & b, Matrix<> * T, int n, int nlevel, int nsmooth){
  //do smoothing using Jacobi
  Vector<> r(b);
  r["i"] -= A["ij"]*x["j"];
  double rnorm0 = r.norm2();
  smooth_jacobi(A,x,b,nsmooth);
  r["i"] = b["i"];
  r["i"] -= A["ij"]*x["j"];
  double rnorm = r.norm2();
  if (nlevel == 0){
    if (A.wrld->rank == 0) printf("At level %d (coarsest level), residual norm was %1.2E initially\n",nlevel,rnorm0);
    if (A.wrld->rank == 0) printf("At level %d (coarsest level), residual norm was %1.2E after smooth\n",nlevel,rnorm);
    return; 
  }
  Vector<> x2(x);
  Vector<> r2(x);

  smooth_jacobi(A,x2,b,nsmooth);
  r2["i"] = b["i"];
  r2["i"] -= A["ij"]*x2["j"];
  double rnorm_alt = r2.norm2();
  

  int m = T[0].lens[1];

//  Vector<> d(x.len, *x.wrld);
//  d["i"] = A["ii"];
//  P["ik"] -= d["i"]*A["ij"]*T[0]["jk"];
  //smooth the restriction/interpolation operator P = (I-omega*diag(A)^{-1}*A)T
  Matrix<> P(T[0].lens[0], T[0].lens[1], SP, *T[0].wrld);
  Matrix<> D(n,n,SP,*A.wrld);
  D["ii"] = A["ii"];
  double omega=.1;
  Transform<>([=](double & d){ d= omega/d; })(D["ii"]);
  P["ik"] = A["ij"]*T[0]["jk"];
  P["ik"] =  D["il"]*P["lk"];
  P["ij"] += T[0]["ij"];

  //restrict residual vector
  Vector<> PTr(m, *x.wrld);
  PTr["i"] = P["ji"]*r["j"];
 
  //coarses initial guess should be zeros
  Vector<> zx(m, *b.wrld);
  
  int atr = 0;
  if (A.is_sparse){ 
    atr = atr | SP;
  }
  //atr = atr | A.symm;
  Matrix<> AP(n, m, atr, *x.wrld);
  Matrix<> PTAP(m, m, atr, *x.wrld);
 
  //restrict A via triple matrix product, should probably be done outside v-cycle
  AP["lj"] = A["lk"]*P["kj"];
  PTAP["ij"] = P["li"]*AP["lj"];
 
  //recurse into coarser level
  vcycle(PTAP, zx, PTr, T+1, m, nlevel-1, nsmooth);
  
//  zx.print();

  //interpolate solution to residual equation at coraser level back
  x["i"] += P["ij"]*zx["j"]; 
 
  //smooth new solution
  r["i"] = b["i"];
  r["i"] -= A["ij"]*x["j"];
  double rnorm2 = r.norm2();
  smooth_jacobi(A,x,b,nsmooth);
  r["i"] = b["i"];
  r["i"] -= A["ij"]*x["j"];
  double rnorm3 = r.norm2();
    if (A.wrld->rank == 0) printf("At level %d, residual norm was %1.2E initially\n",nlevel,rnorm0);
  if (A.wrld->rank == 0) printf("At level %d, residual norm was %1.2E after initial smooth\n",nlevel,rnorm);
  if (A.wrld->rank == 0) printf("At level %d, residual norm was %1.2E after coarse recursion\n",nlevel,rnorm2);
  if (A.wrld->rank == 0) printf("At level %d, residual norm was %1.2E after final smooth\n",nlevel,rnorm3);
  if (A.wrld->rank == 0) printf("At level %d, residual norm was %1.2E if we skipped the coarsening\n",nlevel,rnorm_alt);
}

/**
 * \brief computes Multigrid for a 3D regular discretization
 */
int algebraic_multigrid(int     n,
                        double  sp_frac,
                        int     nlvl,
                        int     ndiv,<commit_msg>got rid of extra indent<commit_after>/** \addtogroup examples 
  * @{ 
  * \defgroup algebraic_multigrid algebraic_multigrid
  * @{ 
  * \brief Benchmark for smoothed algebraic multgrid
  */
#include <ctf.hpp>
using namespace CTF;

void smooth_jacobi(Matrix<> & A, Vector<> & x, Vector <> & b, int nsmooth){
  Vector<> d(x.len, *x.wrld);
  d["i"] = A["ii"];
  Transform<>([](double & d){ d= fabs(d) > 0.0 ? 1./d : 0.0; })(d["i"]);
  Matrix<> R(A);
  R["ii"] = 0.0;

  //20 iterations of Jacobi, should probably be a parameter or some convergence check instead
  for (int i=0; i<nsmooth; i++){
    x["i"] = -1.0*R["ij"]*x["j"];
    x["i"] += b["i"];
    x["i"] *= d["i"];
  }
}

void vcycle(Matrix<> & A, Vector<> & x, Vector<> & b, Matrix<> * T, int n, int nlevel, int nsmooth){
  //do smoothing using Jacobi
  Vector<> r(b);
  r["i"] -= A["ij"]*x["j"];
  double rnorm0 = r.norm2();
  smooth_jacobi(A,x,b,nsmooth);
  r["i"] = b["i"];
  r["i"] -= A["ij"]*x["j"];
  double rnorm = r.norm2();
  if (nlevel == 0){
    if (A.wrld->rank == 0) printf("At level %d (coarsest level), residual norm was %1.2E initially\n",nlevel,rnorm0);
    if (A.wrld->rank == 0) printf("At level %d (coarsest level), residual norm was %1.2E after smooth\n",nlevel,rnorm);
    return; 
  }
  Vector<> x2(x);
  Vector<> r2(x);

  smooth_jacobi(A,x2,b,nsmooth);
  r2["i"] = b["i"];
  r2["i"] -= A["ij"]*x2["j"];
  double rnorm_alt = r2.norm2();
  

  int m = T[0].lens[1];

//  Vector<> d(x.len, *x.wrld);
//  d["i"] = A["ii"];
//  P["ik"] -= d["i"]*A["ij"]*T[0]["jk"];
  //smooth the restriction/interpolation operator P = (I-omega*diag(A)^{-1}*A)T
  Matrix<> P(T[0].lens[0], T[0].lens[1], SP, *T[0].wrld);
  Matrix<> D(n,n,SP,*A.wrld);
  D["ii"] = A["ii"];
  double omega=.1;
  Transform<>([=](double & d){ d= omega/d; })(D["ii"]);
  P["ik"] = A["ij"]*T[0]["jk"];
  P["ik"] =  D["il"]*P["lk"];
  P["ij"] += T[0]["ij"];

  //restrict residual vector
  Vector<> PTr(m, *x.wrld);
  PTr["i"] = P["ji"]*r["j"];
 
  //coarses initial guess should be zeros
  Vector<> zx(m, *b.wrld);
  
  int atr = 0;
  if (A.is_sparse){ 
    atr = atr | SP;
  }
  //atr = atr | A.symm;
  Matrix<> AP(n, m, atr, *x.wrld);
  Matrix<> PTAP(m, m, atr, *x.wrld);
 
  //restrict A via triple matrix product, should probably be done outside v-cycle
  AP["lj"] = A["lk"]*P["kj"];
  PTAP["ij"] = P["li"]*AP["lj"];
 
  //recurse into coarser level
  vcycle(PTAP, zx, PTr, T+1, m, nlevel-1, nsmooth);
  
//  zx.print();

  //interpolate solution to residual equation at coraser level back
  x["i"] += P["ij"]*zx["j"]; 
 
  //smooth new solution
  r["i"] = b["i"];
  r["i"] -= A["ij"]*x["j"];
  double rnorm2 = r.norm2();
  smooth_jacobi(A,x,b,nsmooth);
  r["i"] = b["i"];
  r["i"] -= A["ij"]*x["j"];
  double rnorm3 = r.norm2();
  if (A.wrld->rank == 0) printf("At level %d, residual norm was %1.2E initially\n",nlevel,rnorm0);
  if (A.wrld->rank == 0) printf("At level %d, residual norm was %1.2E after initial smooth\n",nlevel,rnorm);
  if (A.wrld->rank == 0) printf("At level %d, residual norm was %1.2E after coarse recursion\n",nlevel,rnorm2);
  if (A.wrld->rank == 0) printf("At level %d, residual norm was %1.2E after final smooth\n",nlevel,rnorm3);
  if (A.wrld->rank == 0) printf("At level %d, residual norm was %1.2E if we skipped the coarsening\n",nlevel,rnorm_alt);
}

/**
 * \brief computes Multigrid for a 3D regular discretization
 */
int algebraic_multigrid(int     n,
                        double  sp_frac,
                        int     nlvl,
                        int     ndiv,
{'fraction_non_alphanumeric': 0.2365074879909579, 'fraction_numerical': 0.01949703306018649, 'mean_word_length': 3.231320980274955, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 63, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6696980', 'n_tokens_mistral': 3099, 'n_tokens_neox': 2954, 'n_words': 723}