Lutest.cc

#include <iostream.h>
#include "lapack.h"
#include "veclib.h"
#include "Lutest.h"
#include "Random.h"

Lutest::Lutest(int n) {
  size = n;

  cout << "lapack example program" << endl;
  cout << "solving ax=b via lu decomposition" << endl;
  cout << "   a             :" << size << " by " << size << endl;

  Random rnd(0L);
  a = rnd.randomMatrix(size, size);
  b = rnd.randomVector(size);
}

 Vector &Lutest::solve(void) {
  // Call the lapack routine dgesv to solve the linear system a * x = b.
  // dgesv replaces a by its lu decomposition and b by the solution.
  // Since we need the original data for checking, we make copies.
  
  Matrix adecomp = a.copy();
  Vector &x       = b.copy();
  
  int  info;
  int *ipiv  = new int[size];
  int one = 1;
  dgesv(&size, &one, adecomp.getArray(), &size, ipiv, 
	x.getArray(), &size, &info);

  cout << "info code returned by dgesv = " << info << endl;

  delete[] ipiv;

  return x;
}

void Lutest::check(Vector &x) {
  // check the result

  double eps = dlamch("epsilon", 7);      // get machine precision

  // get infinity-norms of a and x
  Vector work(size);
  double anorm = dlange("i", &size, &size, a.getArray(), 
			&size, work.getArray(), 1);
  int one = 1;
  double xnorm = dlange("i", &size, &one, x.getArray(), 
			&size, work.getArray(), 1);

  // compute a*x - b and its norm
  // dgemm replaces b by the result
  Vector b1 = b.copy();

  double oned = 1.0;
  double minusoned = -1.0;
  dgemm("n", "n", &size, &one, &size, &oned, a.getArray(), &size, 
	x.getArray(), &size, &minusoned, b1.getArray(), &size, 1, 1);
  double resnorm = dlange("i", &size, &one, b1.getArray(), &size, 
			  work.getArray(), 1);
  
  // check residual error
  double resid = resnorm / (anorm*xnorm*eps*size);
  if (resid < 10.0) {
    cout << "solution is correct! residual:" << endl;
  } else {
    cout << "solution is incorrect! residual:" << endl;
  }
  cout << "  ||a*x - b|| / ( ||x||*||a||*eps*n ) = " << resid << endl;
}