ABACUS/src/RICHARDSON/Richardson.cc

/**********************************************************

This software is part of J.-S. Caux's ABACUS library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  src/RICHARDSON/Rischardson.cc

Purpose:  Richardson model

***********************************************************/


#include "ABACUS.h"

using namespace std;
using namespace ABACUS;

namespace ABACUS {

  bool Solve_Richardson_Quad_S (DP g_int, const Vect_DP epsilon, const Vect_DP sumoneovereps,
				const Vect_CX S, const Vect_CX sumSovereps, int j, complex<DP>& Sjleft, complex<DP>& Sjright)
  {
    // This solves the equation
    // S^2_j - \sum_{i \neq j}^N \frac{S_j - S_i}{\epsilon_j - \epsilon_i} - (1/g) S_j = 0.
    // for S_j, and puts the result in Sj.
    // If the result is real, true is returned, otherwise false.

    //DP bquad = -1.0/g_int - sumoneovereps[j];
    //DP cquad = sumSovereps[j];
    /*
      DP discr = 0.0;
      bool retval = false;
      //if ((discr = bquad * bquad - 4.0 * cquad) >= 0.0) {
      if ((discr = (1.0/g_int + sumoneovereps[j]) * (1.0/g_int + sumoneovereps[j]) - 4.0 * sumSovereps[j]) >= 0.0) {
      discr = sqrt(discr);
      Sjleft = 0.5 * (-(1.0/g_int + sumoneovereps[j]) - discr);
      Sjright = 0.5 * (-(1.0/g_int + sumoneovereps[j]) + discr);
      retval = true;
      }

      return(retval);
    */

    complex<DP> discr = sqrt((1.0/g_int + sumoneovereps[j]) * (1.0/g_int + sumoneovereps[j]) - 4.0 * sumSovereps[j]);
    Sjleft = 0.5 * (1.0/g_int + sumoneovereps[j] - discr);
    Sjright = 0.5 * (1.0/g_int + sumoneovereps[j] + discr);
    return(true);
  }

  bool Solve_Richardson_Quad_S (DP g_int, const Vect_DP epsilon, Vect_CX& S, Vect<bool> leftroot, DP req_prec)
  {
    // This solves the set of equations
    // S^2_j - \sum_{i \neq j}^N \frac{S_j - S_i}{\epsilon_j - \epsilon_i} - (1/g) S_j = 0.

    int Nlevels = epsilon.size();
    Vect_DP absLHSRE (Nlevels); // absolute value of left-hand side of Richardson eqns for S_j written as LHSRE == 0.

    Vect_DP sumoneovereps (0.0, Nlevels);
    Vect_CX sumSovereps (Nlevels);

    for (int j = 0; j < Nlevels; ++j) {
      sumoneovereps[j] = 0.0;
      for (int i = 0; i < Nlevels; ++i) if (i != j) sumoneovereps[j] += 1.0/(epsilon[j] - epsilon[i]);
      sumSovereps[j] = 0.0;
      for (int i = 0; i < Nlevels; ++i) if (i != j) sumSovereps[j] += S[i]/(epsilon[j] - epsilon[i]);
    }

    DP maxabsLHSRE = 0.0;
    int jmax = 0;
    complex<DP> Sleft, Sright;

    DP damping = 1.0;

    do {

      // Calculate the deviations, identify the largest:
      maxabsLHSRE = 0.0;
      for (int j = 0; j < Nlevels; ++j) {
	absLHSRE[j] = abs(S[j] * S[j] - S[j] * sumoneovereps[j] + sumSovereps[j] - S[j]/g_int);
	if (absLHSRE[j] > maxabsLHSRE) {
	  maxabsLHSRE = absLHSRE[j];
	  jmax = j;
	}
      }

      cout << "identified jmax = " << jmax << " with maxabsLHSRE = " << maxabsLHSRE << endl;

      cout << "S[jmax]: from " << S[jmax];

      // Re-solve for jmax:
      if (Solve_Richardson_Quad_S (g_int, epsilon, sumoneovereps, S, sumSovereps, jmax, Sleft, Sright)) {
	if (abs(S[jmax] - Sleft) < abs(S[jmax] - Sright)) S[jmax] = damping * Sleft + (1.0 - damping) * S[jmax];
	else S[jmax] = damping * Sright + (1.0 - damping) * S[jmax];
      }
      else {
	ABACUSerror("Complex jmax root.");
      }
      cout << " to " << S[jmax] << endl;

      // Re-solve also for a random j, given by
      int jrand = rand() % Nlevels;
      cout << "identified jrand = " << jrand << endl;
      cout << "S[jrand]: from " << S[jrand];

      if (Solve_Richardson_Quad_S (g_int, epsilon, sumoneovereps, S, sumSovereps, jrand, Sleft, Sright)) {
	if (abs(S[jrand] - Sleft) < abs(S[jrand] - Sright)) S[jrand] = damping * Sleft + (1.0 - damping) * S[jrand];
	else S[jrand] = damping * Sright + (1.0 - damping) * S[jrand];
      }
      else {
	ABACUSerror("Complex jrand root.");
      }
      cout << " to " << S[jrand] << endl;

      // Recalculate all sums:
      for (int j = 0; j < Nlevels; ++j) {
	sumSovereps[j] = 0.0;
	for (int i = 0; i < Nlevels; ++i) if (i != j) sumSovereps[j] += S[i]/(epsilon[j] - epsilon[i]);
      }

    } while (maxabsLHSRE > req_prec);

    return(true);
  }


} // namespace ABACUS



int main ()
{
  DP g_int = 0.1;

  int Nlevels = 10;

  Vect_DP epsilon (Nlevels);
  for (int i = 0; i < Nlevels; ++i) epsilon[i] = -0.5 * (Nlevels - 1.0) + i;

  DP req_prec = 1.0e-10;

  Vect_CX S (2.0, Nlevels);

  // Initial conditions:
  /*
  // Sa252 g = 0.78
  S[0] = -0.6142;
  S[1] = -0.7344;
  S[2] = -0.9177;
  S[3] = -1.2414;
  S[4] = -2.06103;
  S[5] = 3.06103;
  S[6] = 2.2414;
  S[7] = 1.9178;
  S[8] = 1.7344;
  S[9] = 1.6142;
  //for (int j = 0; j < Nlevels; ++j) S[j] /= g_int;
  */


  Vect<bool> leftroot (true, Nlevels);

  //for (int j = Nlevels/2; j < Nlevels; ++j) leftroot[j] = true;

  leftroot[0] = 0;
  leftroot[1] = 1;
  leftroot[2] = 0;
  leftroot[3] = 1;
  leftroot[4] = 0;
  leftroot[5] = 1;
  leftroot[6] = 1;
  leftroot[7] = 1;
  leftroot[8] = 1;
  leftroot[9] = 1;

  //for (int j = 0; j < Nlevels; ++j) leftroot[j] = !leftroot[j];

  cout << Solve_Richardson_Quad_S (g_int, epsilon, S, leftroot, req_prec) << endl;
  cout << leftroot << endl;
  for (int j = 0; j < Nlevels; ++j) S[j] *= g_int;
  cout << S << endl;

  return(0);
}