-
Notifications
You must be signed in to change notification settings - Fork 555
Expand file tree
/
Copy pathcg.cpp
More file actions
131 lines (109 loc) · 3.92 KB
/
Copy pathcg.cpp
File metadata and controls
131 lines (109 loc) · 3.92 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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
/*******************************************************
* Copyright (c) 2014, ArrayFire
* All rights reserved.
*
* This file is distributed under 3-clause BSD license.
* The complete license agreement can be obtained at:
* http://arrayfire.com/licenses/BSD-3-Clause
********************************************************/
#include <arrayfire.h>
#include <iostream>
using namespace af;
static size_t dimension = 4 * 1024;
static const int maxIter = 10;
static const int sparsityFactor = 7;
static array A;
static array spA; // Sparse A
static array x0;
static array b;
void setupInputs() {
// Generate a random input: A
array T = randu(dimension, dimension, f32);
// Create 0s in input.
// Anything that is no divisible by sparsityFactor will become 0.
A = floor(T * 1000);
A = A * ((A % sparsityFactor) == 0) / 1000;
// Make it positive definite
A = transpose(A) + A + A.dims(0) * identity(A.dims(0), A.dims(0), f32);
// Make A sparse as spA
spA = sparse(A);
// Generate x0: Random guess
x0 = randu(A.dims(0), f32);
// Generate b
b = matmul(A, x0);
std::cout << "Sparsity of A = "
<< 100.f * (float)sparseGetNNZ(spA) / (float)spA.elements() << "%"
<< std::endl;
std::cout << "Memory Usage of A = " << A.bytes() / (1024.f * 1024.f)
<< " MB" << std::endl;
std::cout << "Memory Usage of spA = "
<< (sparseGetValues(spA).bytes() + sparseGetRowIdx(spA).bytes() +
sparseGetColIdx(spA).bytes()) /
(1024.f * 1024.f)
<< " MB" << std::endl;
}
void sparseConjugateGradient(void) {
array x = constant(0, b.dims(), f32);
array r = b - matmul(spA, x);
array p = r;
for (int i = 0; i < maxIter; ++i) {
array Ap = matmul(spA, p);
array alpha_num = dot(r, r);
array alpha_den = dot(p, Ap);
array alpha = alpha_num / alpha_den;
r -= tile(alpha, Ap.dims()) * Ap;
x += tile(alpha, Ap.dims()) * p;
array beta_num = dot(r, r);
array beta = beta_num / alpha_num;
p = r + tile(beta, p.dims()) * p;
}
}
void denseConjugateGradient(void) {
array x = constant(0, b.dims(), f32);
array r = b - matmul(A, x);
array p = r;
for (int i = 0; i < maxIter; ++i) {
array Ap = matmul(A, p);
array alpha_num = dot(r, r);
array alpha_den = dot(p, Ap);
array alpha = alpha_num / alpha_den;
r -= tile(alpha, Ap.dims()) * Ap;
x += tile(alpha, Ap.dims()) * p;
array beta_num = dot(r, r);
array beta = beta_num / alpha_num;
p = r + tile(beta, p.dims()) * p;
}
}
void checkConjugateGradient(const af::array in) {
array x = constant(0, b.dims(), f32);
array r = b - matmul(in, x);
array p = r;
for (int i = 0; i < maxIter; ++i) {
array Ap = matmul(in, p);
array alpha_num = dot(r, r);
array alpha_den = dot(p, Ap);
array alpha = alpha_num / alpha_den;
r -= tile(alpha, Ap.dims()) * Ap;
x += tile(alpha, Ap.dims()) * p;
array beta_num = dot(r, r);
array beta = beta_num / alpha_num;
p = r + tile(beta, p.dims()) * p;
}
array res = x0 - x;
std::cout << "Final difference in solutions:\n";
af_print(dot(res, res));
}
int main(int, char **) {
af::info();
setupInputs();
std::cout << "Verifying Dense Conjugate Gradient:" << std::endl;
checkConjugateGradient(A);
std::cout << "Verifying Sparse Conjugate Gradient:" << std::endl;
checkConjugateGradient(spA);
af::sync();
std::cout << "Dense Conjugate Gradient Time: "
<< timeit(denseConjugateGradient) * 1000 << "ms" << std::endl;
std::cout << "Sparse Conjugate Gradient Time: "
<< timeit(sparseConjugateGradient) * 1000 << "ms" << std::endl;
return 0;
}