-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathIntegra.cpp
More file actions
65 lines (56 loc) · 2.45 KB
/
Integra.cpp
File metadata and controls
65 lines (56 loc) · 2.45 KB
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
#include <iostream>
#include <experimental/simd>
int main() {
int n = 24;
double* A = new double(n*n);
using Type = std::experimental::fixed_size_simd<double, 2>;
constexpr std::size_t K = 2;
for (int i = 1; i < n;i++) {
for (int j = 0; j < i;j++) {
Type norm(0), pro(0);// norm1(0), pro1(0);
for (int k = 0; k < n;k+=K) {
// store previous Gram-Schmidt in A
Type va(&A[j * n + k], std::experimental::element_aligned);
Type vb(&A[i * n + k], std::experimental::element_aligned);
//Type vc(&A[j * n + k + 2], std::experimental::element_aligned);
//Type vd(&A[i * n + k + 2], std::experimental::element_aligned);
// Type ve(&A[j * n + k + 4], std::experimental::element_aligned);
// Type vf(&A[i * n + k + 4], std::experimental::element_aligned);
// Type vg(&A[j * n + k + 6], std::experimental::element_aligned);
// Type vh(&A[i * n + k + 6], std::experimental::element_aligned);
norm += va * va;
pro += va * vb;
//norm += vc * vc;
//pro += vc * vd;
// norm += ve * ve;
// pro += ve * vf;
// norm += vg * vg;
// pro += vg * vh;
//norm += A[i * n + k] * A[i * n + k];
//in_pro += A[j * n + k] * A[i * n + k];
}
//double coef = (std::experimental::reduce(pro) + std::experimental::reduce(pro1)) /
// (std::experimental::reduce(norm) + std::experimental::reduce(norm1));
double coef = (std::experimental::reduce(pro)) /
(std::experimental::reduce(norm));
for (int k = 0;k < n;k+=K) {
//A[i * n + k] -= coef * A[j * n + k];
Type vi(&A[i * n + k], std::experimental::element_aligned);
Type vj(&A[j * n + k], std::experimental::element_aligned);
vi -= coef * vj;
vi.copy_to(&A[i * n + k], std::experimental::element_aligned);
}
}
}
for (int i = 0; i <n;i++) {
for (int j = i + 1; j < n;j++) {
double pro= 0.0;
for (int k = 0; k < n;k++) {
pro += A[i * n + k] * A[j * n + k];
}
if (abs(pro) > 1e-12) {
std::cout << pro << std::endl;
}
}
}
}