FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
apply_mkl.cpp
1// FEAT3: Finite Element Analysis Toolbox, Version 3
2// Copyright (C) 2010 by Stefan Turek & the FEAT group
3// FEAT3 is released under the GNU General Public License version 3,
4// see the file 'copyright.txt' in the top level directory for details.
5
6// includes, FEAT
7#include <kernel/base_header.hpp>
8#include <kernel/lafem/arch/apply.hpp>
9#include <kernel/util/exception.hpp>
10
11#include <cstring>
12
13FEAT_DISABLE_WARNINGS
14#include <mkl.h>
15#include <mkl_spblas.h>
16FEAT_RESTORE_WARNINGS
17
19// disable 'deprecated' warnings for MKL functions, because this file won't compile otherwise
20#ifdef FEAT_COMPILER_MICROSOFT
21#pragma warning(disable: 4996)
22#endif
23
24using namespace FEAT;
25using namespace FEAT::LAFEM;
26using namespace FEAT::LAFEM::Arch;
27
28void Apply::csr_mkl(float * r, const float a, const float * const x, const float b, const float * const y, const float * const val, const Index * const col_ind, const Index * const row_ptr, const Index rows, const Index columns, const Index, const bool transposed)
29{
30 MKL_INT mrows = (MKL_INT)rows;
31 MKL_INT mcolumns = (MKL_INT)columns;
32
33 if (r != y)
34 {
35 MKL_INT one = 1;
36 if (transposed)
37 scopy(&mcolumns, (const float*)y, &one, r, &one);
38 else
39 scopy(&mrows, (const float*)y, &one, r, &one);
40 }
41
42#ifdef FEAT_USE_MKL_SPARSE_EXECUTOR
43
44 sparse_operation_t opt;
45 if (transposed)
46 opt = SPARSE_OPERATION_TRANSPOSE;
47 else
48 opt = SPARSE_OPERATION_NON_TRANSPOSE;
49
50 sparse_matrix_t A;
51 FEAT_DISABLE_WARNINGS
52 sparse_status_t status = mkl_sparse_s_create_csr(&A, SPARSE_INDEX_BASE_ZERO, mrows, mcolumns, (MKL_INT*)row_ptr, (MKL_INT*)(row_ptr + 1), (MKL_INT*)col_ind, (float*) val);
53 FEAT_RESTORE_WARNINGS
54 if (status != SPARSE_STATUS_SUCCESS)
55 XABORTM("MKL Sparse Error occurred in execution!\n");
56 matrix_descr md;
57 md.type = SPARSE_MATRIX_TYPE_GENERAL;
58 status = mkl_sparse_s_mv(opt, a, A, md, x, b, r);
59 if (status != SPARSE_STATUS_SUCCESS)
60 XABORTM("MKL Sparse Error occurred in execution!\n");
61 mkl_sparse_destroy(A);
62
63#else
64
65 char matdescra[6];
66 matdescra[0] = 'G';
67 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
68 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
69 matdescra[3] = 'C';
70
71 char trans;
72 if (transposed)
73 trans = 'T';
74 else
75 trans = 'N';
76
77 FEAT_DISABLE_WARNINGS
78 mkl_scsrmv(&trans, &mrows, &mcolumns, (const float*)&a, matdescra, (const float*) val, (const MKL_INT*)col_ind, (const MKL_INT*)row_ptr, (const MKL_INT*)(row_ptr) + 1, (const float*)x, (const float*)&b, r);
79 FEAT_RESTORE_WARNINGS
80
81#endif
82}
83
84void Apply::csr_mkl(double * r, const double a, const double * const x, const double b, const double * const y, const double * const val, const Index * const col_ind, const Index * const row_ptr, const Index rows, const Index columns, const Index, const bool transposed)
85{
86 MKL_INT mrows = (MKL_INT)rows;
87 MKL_INT mcolumns = (MKL_INT)columns;
88
89 if (r != y)
90 {
91 MKL_INT one = 1;
92 if (transposed)
93 dcopy(&mcolumns, (const double*)y, &one, r, &one);
94 else
95 dcopy(&mrows, (const double*)y, &one, r, &one);
96 }
97
98#ifdef FEAT_USE_MKL_SPARSE_EXECUTOR
99
100 sparse_operation_t opt;
101 if (transposed)
102 opt = SPARSE_OPERATION_TRANSPOSE;
103 else
104 opt = SPARSE_OPERATION_NON_TRANSPOSE;
105
106 sparse_matrix_t A;
107 FEAT_DISABLE_WARNINGS
108 sparse_status_t status = mkl_sparse_d_create_csr(&A, SPARSE_INDEX_BASE_ZERO, mrows, mcolumns, (MKL_INT*)row_ptr, (MKL_INT*)(row_ptr + 1), (MKL_INT*)col_ind, (double*) val);
109 FEAT_RESTORE_WARNINGS
110 if (status != SPARSE_STATUS_SUCCESS)
111 XABORTM("MKL Sparse Error occurred in execution!\n");
112 matrix_descr md;
113 md.type = SPARSE_MATRIX_TYPE_GENERAL;
114 status = mkl_sparse_d_mv(opt, a, A, md, x, b, r);
115 if (status != SPARSE_STATUS_SUCCESS)
116 XABORTM("MKL Sparse Error occurred in execution!\n");
117 mkl_sparse_destroy(A);
118
119#else
120
121 char matdescra[6];
122 matdescra[0] = 'G';
123 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
124 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
125 matdescra[3] = 'C';
126
127 char trans;
128 if (transposed)
129 trans = 'T';
130 else
131 trans = 'N';
132
133 FEAT_DISABLE_WARNINGS
134 mkl_dcsrmv(&trans, &mrows, &mcolumns, (const double*)&a, matdescra, (const double*) val, (const MKL_INT*)col_ind, (const MKL_INT*)row_ptr, (const MKL_INT*)(row_ptr) + 1, (const double*)x, (const double*)&b, r);
135 FEAT_RESTORE_WARNINGS
136
137#endif
138}
139
140void Apply::bcsr_mkl(float * r, const float a, const float * const x, const float b, const float * const y, const float * const val, const Index * const col_ind, const Index * const row_ptr, const Index rows, const Index columns, const Index, const int blocksize)
141{
142 MKL_INT mrows = (MKL_INT)rows;
143 MKL_INT mcolumns = (MKL_INT)columns;
144 MKL_INT mblocksize = (MKL_INT)blocksize;
145 MKL_INT mcopysize = mrows * mblocksize;
146
147 if (r != y)
148 {
149 MKL_INT one = 1;
150 scopy(&mcopysize, (const float*)y, &one, r, &one);
151 }
152
153
154#ifdef FEAT_USE_MKL_SPARSE_EXECUTOR
155
156 sparse_operation_t opt = SPARSE_OPERATION_NON_TRANSPOSE;
157
158 sparse_matrix_t A;
159 FEAT_DISABLE_WARNINGS
160 mkl_sparse_s_create_bsr(&A, SPARSE_INDEX_BASE_ZERO, SPARSE_LAYOUT_ROW_MAJOR, mrows, mcolumns, mblocksize, (MKL_INT*)row_ptr, (MKL_INT*)(row_ptr + 1), (MKL_INT*)col_ind, (float*) val);
161 FEAT_RESTORE_WARNINGS
162 matrix_descr md;
163 md.type = SPARSE_MATRIX_TYPE_GENERAL;
164 mkl_sparse_s_mv(opt, a, A, md, x, b, r);
165 mkl_sparse_destroy(A);
166
167#else
168
169 char trans = 'N';
170 char matdescra[6];
171 matdescra[0] = 'G';
172 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
173 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
174 matdescra[3] = 'C';
175
176 FEAT_DISABLE_WARNINGS
177 mkl_sbsrmv(&trans, &mrows, &mcolumns, &mblocksize, (const float*)&a, matdescra, (const float*) val, (const MKL_INT*)col_ind, (const MKL_INT*)row_ptr, (const MKL_INT*)(row_ptr) + 1, (const float*)x, (const float*)&b, r);
178 FEAT_RESTORE_WARNINGS
179
180#endif
181}
182
183void Apply::bcsr_mkl(double * r, const double a, const double * const x, const double b, const double * const y, const double * const val, const Index * const col_ind, const Index * const row_ptr, const Index rows, const Index columns, const Index, const int blocksize)
184{
185 MKL_INT mrows = (MKL_INT)rows;
186 MKL_INT mcolumns = (MKL_INT)columns;
187 MKL_INT mblocksize = (MKL_INT)blocksize;
188 MKL_INT mcopysize = mrows * mblocksize;
189
190 if (r != y)
191 {
192 MKL_INT one = 1;
193 dcopy(&mcopysize, (const double*)y, &one, r, &one);
194 }
195
196
197#ifdef FEAT_USE_MKL_SPARSE_EXECUTOR
198
199 sparse_operation_t opt = SPARSE_OPERATION_NON_TRANSPOSE;
200
201 sparse_matrix_t A;
202 FEAT_DISABLE_WARNINGS
203 mkl_sparse_d_create_bsr(&A, SPARSE_INDEX_BASE_ZERO, SPARSE_LAYOUT_ROW_MAJOR, mrows, mcolumns, mblocksize, (MKL_INT*)row_ptr, (MKL_INT*)(row_ptr + 1), (MKL_INT*)col_ind, (double*) val);
204 FEAT_RESTORE_WARNINGS
205 matrix_descr md;
206 md.type = SPARSE_MATRIX_TYPE_GENERAL;
207 mkl_sparse_d_mv(opt, a, A, md, x, b, r);
208 mkl_sparse_destroy(A);
209
210#else
211
212 char trans = 'N';
213 char matdescra[6];
214 matdescra[0] = 'G';
215 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
216 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
217 matdescra[3] = 'C';
218
219 FEAT_DISABLE_WARNINGS
220 mkl_dbsrmv(&trans, &mrows, &mcolumns, &mblocksize, (const double*)&a, matdescra, (const double*) val, (const MKL_INT*)col_ind, (const MKL_INT*)row_ptr, (const MKL_INT*)(row_ptr) + 1, (const double*)x, (const double*)&b, r);
221 FEAT_RESTORE_WARNINGS
222
223#endif
224}
225
226void Apply::bcsr_transposed_mkl(float * r, const float a, const float * const x, const float b, const float * const y, const float * const val, const Index * const col_ind, const Index * const row_ptr, const Index rows, const Index columns, const Index, const int blocksize)
227{
228 MKL_INT mrows = (MKL_INT)rows;
229 MKL_INT mcolumns = (MKL_INT)columns;
230 MKL_INT mblocksize = (MKL_INT)blocksize;
231 MKL_INT mcopysize = mcolumns * mblocksize;
232
233 if (r != y)
234 {
235 MKL_INT one = 1;
236 scopy(&mcopysize, (const float*)y, &one, r, &one);
237 }
238
239
240#ifdef FEAT_USE_MKL_SPARSE_EXECUTOR
241
242 sparse_operation_t opt = SPARSE_OPERATION_NON_TRANSPOSE;
243
244 sparse_matrix_t A;
245 FEAT_DISABLE_WARNINGS
246 mkl_sparse_s_create_bsr(&A, SPARSE_INDEX_BASE_ZERO, SPARSE_LAYOUT_ROW_MAJOR, mrows, mcolumns, mblocksize, (MKL_INT*)row_ptr, (MKL_INT*)(row_ptr + 1), (MKL_INT*)col_ind, (float*) val);
247 FEAT_RESTORE_WARNINGS
248 matrix_descr md;
249 md.type = SPARSE_MATRIX_TYPE_GENERAL;
250 mkl_sparse_s_mv(opt, a, A, md, x, b, r);
251 mkl_sparse_destroy(A);
252
253#else
254
255 char trans = 'T';
256 char matdescra[6];
257 matdescra[0] = 'G';
258 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
259 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
260 matdescra[3] = 'C';
261
262 FEAT_DISABLE_WARNINGS
263 mkl_sbsrmv(&trans, &mrows, &mcolumns, &mblocksize, (const float*)&a, matdescra, (const float*) val, (const MKL_INT*)col_ind, (const MKL_INT*)row_ptr, (const MKL_INT*)(row_ptr) + 1, (const float*)x, (const float*)&b, r);
264 FEAT_RESTORE_WARNINGS
265
266#endif
267}
268
269void Apply::bcsr_transposed_mkl(double * r, const double a, const double * const x, const double b, const double * const y, const double * const val, const Index * const col_ind, const Index * const row_ptr, const Index rows, const Index columns, const Index, const int blocksize)
270{
271 MKL_INT mrows = (MKL_INT)rows;
272 MKL_INT mcolumns = (MKL_INT)columns;
273 MKL_INT mblocksize = (MKL_INT)blocksize;
274 MKL_INT mcopysize = mcolumns * mblocksize;
275
276 if (r != y)
277 {
278 MKL_INT one = 1;
279 dcopy(&mcopysize, (const double*)y, &one, r, &one);
280 }
281
282
283#ifdef FEAT_USE_MKL_SPARSE_EXECUTOR
284
285 sparse_operation_t opt = SPARSE_OPERATION_NON_TRANSPOSE;
286
287 sparse_matrix_t A;
288 FEAT_DISABLE_WARNINGS
289 mkl_sparse_d_create_bsr(&A, SPARSE_INDEX_BASE_ZERO, SPARSE_LAYOUT_ROW_MAJOR, mrows, mcolumns, mblocksize, (MKL_INT*)row_ptr, (MKL_INT*)(row_ptr + 1), (MKL_INT*)col_ind, (double*) val);
290 FEAT_RESTORE_WARNINGS
291 matrix_descr md;
292 md.type = SPARSE_MATRIX_TYPE_GENERAL;
293 mkl_sparse_d_mv(opt, a, A, md, x, b, r);
294 mkl_sparse_destroy(A);
295
296#else
297
298 char trans = 'T';
299 char matdescra[6];
300 matdescra[0] = 'G';
301 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
302 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
303 matdescra[3] = 'C';
304
305 FEAT_DISABLE_WARNINGS
306 mkl_dbsrmv(&trans, &mrows, &mcolumns, &mblocksize, (const double*)&a, matdescra, (const double*) val, (const MKL_INT*)col_ind, (const MKL_INT*)row_ptr, (const MKL_INT*)(row_ptr) + 1, (const double*)x, (const double*)&b, r);
307 FEAT_RESTORE_WARNINGS
308
309#endif
310}
311
312void Apply::dense_mkl(float * r, const float alpha, const float beta, const float * const y, const float * const val, const float * const x, const Index rows, const Index columns)
313{
314 MKL_INT mrows = (MKL_INT)rows;
315 MKL_INT mcolumns = (MKL_INT)columns;
316
317 if (r != y)
318 {
319 MKL_INT one = 1;
320 scopy(&mrows, (const float*)y, &one, r, &one);
321 }
322 cblas_sgemv(CblasRowMajor, CblasNoTrans, mrows, mcolumns, alpha, val, mcolumns, x, 1, beta, r, 1);
323}
324
325void Apply::dense_mkl(double * r, const double alpha, const double beta, const double * const y, const double * const val, const double * const x, const Index rows, const Index columns)
326{
327 MKL_INT mrows = (MKL_INT)rows;
328 MKL_INT mcolumns = (MKL_INT)columns;
329
330 if (r != y)
331 {
332 MKL_INT one = 1;
333 dcopy(&mrows, (const double*)y, &one, r, &one);
334 }
335 cblas_dgemv(CblasRowMajor, CblasNoTrans, mrows, mcolumns, alpha, val, mcolumns, x, 1, beta, r, 1);
336}
XABORTM
#define XABORTM(msg)
Abortion macro definition with custom message.
Definition: assertion.hpp:192
base_header.hpp
FEAT Kernel base header.
FEAT::LAFEM
LAFEM namespace.
Definition: apply.hpp:22
FEAT
FEAT namespace.
Definition: adjactor.hpp:12
FEAT::Index
std::uint64_t Index
Index data type.
Definition: base_header.hpp:122