1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
|
#ifndef ACTIONGESV_HPP
#define ACTIONGESV_HPP
#include <vector>
#include <iostream>
#include <string>
#include "LinearCongruential.hpp"
extern "C" {
void dgesv_(const int*, const int*, double*, const int*, double*,
double*, const int*, int*);
void dgemv_(const char*, const int*, const int*, const double*,
const double*, const int*, const double*, const int*,
const double*, double*, const int*);
double dnrm2_(const int*, const double*, const int*);
}
template<typename value_t = double>
class ActionGESV {
typedef std::vector<value_t> storage_t;
typedef LinearCongruential<> rangen_t;
public:
ActionGESV(const int& size, const unsigned& seed=0)
: size(size), rg(seed), A(rg.fillVector<>(size*size, 0.)), Acopy(A),
x(rg.fillVector(size, 0.)), b(x), bcopy(b)
{ }
void compute() {
const int ONE = 1;
int info;
std::vector<value_t> ipiv(size);
dgesv_(&size, &ONE, &A[0], &size, &ipiv[0], &x[0], &size, &info);
if (info != 0)
std::cerr << "Info: " << info << "\n";
}
double getResidual() {
const double alpha = -1., beta = 1.;
const int ONE = 1;
dgemv_("N", &size, &size, &alpha, &Acopy[0], &size, &x[0], &ONE, &beta,
&b[0], &ONE);
return dnrm2_(&size, &b[0], &ONE)/dnrm2_(&size, &bcopy[0], &ONE);
}
static std::string fileName() {
return "accuracy_general_solve.dat";
}
private:
const int size;
rangen_t rg;
storage_t A, Acopy, x, b, bcopy;
};
#endif // ACTIONGESV_HPP
|
GitWeb
