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
|
#include "BLAS.h"
#include "Matrix.h"
#include "Random48.h"
#include <cstdlib> // for atoi
#define USE_PTHREADS_OR_NOT_NG
#include "Parallelizer.h"
typedef double MyRealType;
typedef PsimagLite::Matrix<MyRealType> MatrixType;
class MyBlasWrapper
{
public:
typedef PsimagLite::Vector<MatrixType*>::Type VectorMatrixType;
MyBlasWrapper(SizeType m, SizeType n, SizeType k, SizeType total)
: a_(total)
, b_(total)
, c_(total)
, myrandom48_(1234)
{
for (SizeType i = 0; i < total; ++i) {
SizeType factor = (i & 1) ? 2 : 1;
a_[i] = new MatrixType(m * factor, k * factor);
b_[i] = new MatrixType(k * factor, n * factor);
c_[i] = new MatrixType(m * factor, n * factor);
fillMatrix(*(a_[i]));
fillMatrix(*(b_[i]));
}
}
~MyBlasWrapper()
{
const SizeType total = a_.size();
for (SizeType i = 0; i < total; ++i) {
delete a_[i];
delete b_[i];
delete c_[i];
a_[i] = b_[i] = c_[i] = 0;
}
}
SizeType tasks() const { return a_.size(); }
void doTask(SizeType ind, SizeType)
{
assert(a_[ind] && b_[ind] && c_[ind]);
SizeType mm = a_[ind]->rows();
SizeType nn = b_[ind]->cols();
SizeType kk = a_[ind]->cols();
assert(kk == b_[ind]->rows());
assert(c_[ind]->rows() == mm && c_[ind]->cols() == nn);
SizeType lda = a_[ind]->rows();
SizeType ldb = b_[ind]->rows();
SizeType ldc = c_[ind]->rows();
MyRealType* aptr = &(a_[ind]->operator()(0, 0));
MyRealType* bptr = &(b_[ind]->operator()(0, 0));
MyRealType* cptr = &(c_[ind]->operator()(0, 0));
psimag::BLAS::GEMM('N', 'N',
mm, // rows of op(A)
nn, // columns of op(B)
kk, // columns of op(A)
1.0,
aptr,
lda, // first dimension of A
bptr,
ldb, // first dimension of B
0.0,
cptr,
ldc); // first dimension of C
*(a_[ind]) = *(b_[ind]);
*(b_[ind]) = *(c_[ind]);
}
private:
void fillMatrix(MatrixType& m)
{
for (SizeType j = 0; j < m.cols(); ++j)
for (SizeType i = 0; i < m.rows(); ++i)
m(i, j) = myrandom48_();
}
VectorMatrixType a_;
VectorMatrixType b_;
VectorMatrixType c_;
PsimagLite::Random48<MyRealType> myrandom48_;
};
int main(int argc, char** argv)
{
if (argc < 6) {
std::cerr << "USAGE: " << argv[0] << " m k n total threads\n";
return 1;
}
SizeType m = atoi(argv[1]);
SizeType k = atoi(argv[2]);
SizeType n = atoi(argv[3]);
SizeType total = atoi(argv[4]);
SizeType threads = atoi(argv[5]);
PsimagLite::CodeSectionParams codeSections(threads);
PsimagLite::Parallelizer<MyBlasWrapper> parallel(codeSections);
MyBlasWrapper myblasWrapper(m, n, k, total);
parallel.loopCreate(myblasWrapper);
}
|
