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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
|
//=====================================================
// Copyright (C) 2011 Andrea Arteaga <andyspiros@gmail.com>
//=====================================================
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//
#ifndef LAPACK_LU_HH
#define LAPACK_LU_HH
#include "LinearCongruential.hh"
#include "diff.hh"
vector<int> get_piv(const vector<int>& ipiv)
{
const int size = ipiv.size();
vector<int> ret(size);
for (int i = 0; i < size; ++i)
ret[i] = i;
int ii, jj, tmp;
for (int i = 1; i <= size; ++i) {
ii = i-1;
jj = ipiv[ii]-1;
tmp = ret[ii];
ret[ii] = ret[jj];
ret[jj] = tmp;
}
return ret;
}
double test_LU(const int& N, const unsigned& seed = 0)
{
LinearCongruential lc(seed);
vector<double> A(N*N), Acopy, P(N*N);
/* Fill A and Acopy */
for (int i = 0; i < N*N; ++i)
A[i] = lc.get_01();
Acopy = A;
/* Compute decomposition */
vector<int> ipiv(N);
int info;
dgetrf_(&N, &N, &A[0], &N, &ipiv[0], &info);
vector<int> piv = get_piv(ipiv);
/* Construct P */
for (int r = 0; r < N; ++r) {
int c = piv[r];
P[c+N*r] = 1;
}
/* Test */
const double alpha = 1.;
dtrmm_("R", "L", "N", "U", &N, &N, &alpha, &A[0], &N, &P[0], &N);
dtrmm_("R", "U", "N", "N", &N, &N, &alpha, &A[0], &N, &P[0], &N);
return diff(Acopy, P);
}
#endif /* LAPACK_LU_HH */
|
GitWeb
