#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; }