ABACUS/src/MATRIX/elmhes.cc

#include "ABACUS.h"
using namespace std;

void ABACUS::elmhes(SQMat_DP& a)
{
  // Reduction to Hessenberg form by elimination method.  NR p.490

  int i, j, m;
  DP y, x;

  int n = a.size();
  for (m = 1; m < n-1; m++) {
    x = 0.0;
    i = m;
    for (j = m; j < n; j++) {
      if (fabs(a[j][m-1]) > fabs(x)) {
	x = a[j][m-1];
	i = j;
      }
    }
    if (i != m) {
      for (j = m - 1; j < n; j++) SWAP(a[i][j], a[m][j]);
      for (j = 0; j < n; j++) SWAP(a[j][i], a[j][m]);
    }
    if (x != 0.0) {
      for (i = m + 1; i < n; i++) {
	if ((y = a[i][m-1]) != 0.0) {
	  y /= x;
	  a[i][m-1] = y;
	  for (j = m; j < n; j++) a[i][j] -= y*a[m][j];
	  for (j = 0; j < n; j++) a[j][m] += y*a[j][i];
	}
      }
    }
  }
}