ABACUS/src/MATRIX/lubksb_CX.cc

#include "ABACUS.h"
using namespace std;

void ABACUS::lubksb_CX (SQMat_CX& a, Vect_INT& indx, Vect_CX& b)
{
  int i, ii=0, ip, j;
  complex<DP> sum;

  int n = a.size();
  for (i = 0; i < n; i++) {
    ip = indx[i];
    sum = b[ip];
    b[ip] = b[i];
    if (ii != 0)
      for (j = ii-1; j < i; j++) sum -= a[i][j] * b[j];
    else if (sum != 0.0)
      ii = i + 1;
    b[i] = sum;
  }
  for (i = n - 1; i >= 0; i--) {
    sum = b[i];
    for (j = i + 1; j < n; j++) sum -= a[i][j] * b[j];
    b[i] = sum/a[i][i];
  }
}