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
18using namespace FEAT;
19using namespace FEAT::LAFEM;
20using namespace FEAT::LAFEM::Arch;
21
22void 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)
23{
24 MKL_INT mrows = (MKL_INT)rows;
25 MKL_INT mcolumns = (MKL_INT)columns;
26
27 if (r != y)
28 {
29 MKL_INT one = 1;
30 if (transposed)
31 scopy(&mcolumns, (const float*)y, &one, r, &one);
32 else
33 scopy(&mrows, (const float*)y, &one, r, &one);
34 }
35
36#ifndef FEAT_USE_MKL_LEGACY_SPMV
37
38 sparse_operation_t opt;
39 if (transposed)
40 opt = SPARSE_OPERATION_TRANSPOSE;
41 else
42 opt = SPARSE_OPERATION_NON_TRANSPOSE;
43
44 sparse_matrix_t A;
45 FEAT_DISABLE_WARNINGS
46 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);
47 FEAT_RESTORE_WARNINGS
48 if (status != SPARSE_STATUS_SUCCESS)
49 XABORTM("MKL Sparse Error occurred in execution!\n");
50 matrix_descr md;
51 md.type = SPARSE_MATRIX_TYPE_GENERAL;
52 status = mkl_sparse_s_mv(opt, a, A, md, x, b, r);
53 if (status != SPARSE_STATUS_SUCCESS)
54 XABORTM("MKL Sparse Error occurred in execution!\n");
55 mkl_sparse_destroy(A);
56
57#else
58
59 char matdescra[6];
60 matdescra[0] = 'G';
61 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
62 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
63 matdescra[3] = 'C';
64
65 char trans;
66 if (transposed)
67 trans = 'T';
68 else
69 trans = 'N';
70
71 FEAT_DISABLE_WARNINGS
72 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);
73 FEAT_RESTORE_WARNINGS
74
75#endif
76}
77
78void 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)
79{
80 MKL_INT mrows = (MKL_INT)rows;
81 MKL_INT mcolumns = (MKL_INT)columns;
82
83 if (r != y)
84 {
85 MKL_INT one = 1;
86 if (transposed)
87 dcopy(&mcolumns, (const double*)y, &one, r, &one);
88 else
89 dcopy(&mrows, (const double*)y, &one, r, &one);
90 }
91
92#ifndef FEAT_USE_MKL_LEGACY_SPMV
93
94 sparse_operation_t opt;
95 if (transposed)
96 opt = SPARSE_OPERATION_TRANSPOSE;
97 else
98 opt = SPARSE_OPERATION_NON_TRANSPOSE;
99
100 sparse_matrix_t A;
101 FEAT_DISABLE_WARNINGS
102 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);
103 FEAT_RESTORE_WARNINGS
104 if (status != SPARSE_STATUS_SUCCESS)
105 XABORTM("MKL Sparse Error occurred in execution!\n");
106 matrix_descr md;
107 md.type = SPARSE_MATRIX_TYPE_GENERAL;
108 status = mkl_sparse_d_mv(opt, a, A, md, x, b, r);
109 if (status != SPARSE_STATUS_SUCCESS)
110 XABORTM("MKL Sparse Error occurred in execution!\n");
111 mkl_sparse_destroy(A);
112
113#else
114
115 char matdescra[6];
116 matdescra[0] = 'G';
117 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
118 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
119 matdescra[3] = 'C';
120
121 char trans;
122 if (transposed)
123 trans = 'T';
124 else
125 trans = 'N';
126
127 FEAT_DISABLE_WARNINGS
128 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);
129 FEAT_RESTORE_WARNINGS
130
131#endif
132}
133
134void 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)
135{
136 MKL_INT mrows = (MKL_INT)rows;
137 MKL_INT mcolumns = (MKL_INT)columns;
138 MKL_INT mblocksize = (MKL_INT)blocksize;
139 MKL_INT mcopysize = mrows * mblocksize;
140
141 if (r != y)
142 {
143 MKL_INT one = 1;
144 scopy(&mcopysize, (const float*)y, &one, r, &one);
145 }
146
147
148#ifndef FEAT_USE_MKL_LEGACY_SPMV
149
150 sparse_operation_t opt = SPARSE_OPERATION_NON_TRANSPOSE;
151
152 sparse_matrix_t A;
153 FEAT_DISABLE_WARNINGS
154 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);
155 FEAT_RESTORE_WARNINGS
156 matrix_descr md;
157 md.type = SPARSE_MATRIX_TYPE_GENERAL;
158 mkl_sparse_s_mv(opt, a, A, md, x, b, r);
159 mkl_sparse_destroy(A);
160
161#else
162
163 char trans = 'N';
164 char matdescra[6];
165 matdescra[0] = 'G';
166 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
167 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
168 matdescra[3] = 'C';
169
170 FEAT_DISABLE_WARNINGS
171 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);
172 FEAT_RESTORE_WARNINGS
173
174#endif
175}
176
177void 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)
178{
179 MKL_INT mrows = (MKL_INT)rows;
180 MKL_INT mcolumns = (MKL_INT)columns;
181 MKL_INT mblocksize = (MKL_INT)blocksize;
182 MKL_INT mcopysize = mrows * mblocksize;
183
184 if (r != y)
185 {
186 MKL_INT one = 1;
187 dcopy(&mcopysize, (const double*)y, &one, r, &one);
188 }
189
190
191#ifndef FEAT_USE_MKL_LEGACY_SPMV
192
193 sparse_operation_t opt = SPARSE_OPERATION_NON_TRANSPOSE;
194
195 sparse_matrix_t A;
196 FEAT_DISABLE_WARNINGS
197 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);
198 FEAT_RESTORE_WARNINGS
199 matrix_descr md;
200 md.type = SPARSE_MATRIX_TYPE_GENERAL;
201 mkl_sparse_d_mv(opt, a, A, md, x, b, r);
202 mkl_sparse_destroy(A);
203
204#else
205
206 char trans = 'N';
207 char matdescra[6];
208 matdescra[0] = 'G';
209 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
210 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
211 matdescra[3] = 'C';
212
213 FEAT_DISABLE_WARNINGS
214 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);
215 FEAT_RESTORE_WARNINGS
216
217#endif
218}
219
220void 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)
221{
222 MKL_INT mrows = (MKL_INT)rows;
223 MKL_INT mcolumns = (MKL_INT)columns;
224 MKL_INT mblocksize = (MKL_INT)blocksize;
225 MKL_INT mcopysize = mcolumns * mblocksize;
226
227 if (r != y)
228 {
229 MKL_INT one = 1;
230 scopy(&mcopysize, (const float*)y, &one, r, &one);
231 }
232
233#ifndef FEAT_USE_MKL_LEGACY_SPMV
234
235 sparse_operation_t opt = SPARSE_OPERATION_NON_TRANSPOSE;
236
237 sparse_matrix_t A;
238 FEAT_DISABLE_WARNINGS
239 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);
240 FEAT_RESTORE_WARNINGS
241 matrix_descr md;
242 md.type = SPARSE_MATRIX_TYPE_GENERAL;
243 mkl_sparse_s_mv(opt, a, A, md, x, b, r);
244 mkl_sparse_destroy(A);
245
246#else
247
248 char trans = 'T';
249 char matdescra[6];
250 matdescra[0] = 'G';
251 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
252 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
253 matdescra[3] = 'C';
254
255 FEAT_DISABLE_WARNINGS
256 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);
257 FEAT_RESTORE_WARNINGS
258
259#endif
260}
261
262void 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)
263{
264 MKL_INT mrows = (MKL_INT)rows;
265 MKL_INT mcolumns = (MKL_INT)columns;
266 MKL_INT mblocksize = (MKL_INT)blocksize;
267 MKL_INT mcopysize = mcolumns * mblocksize;
268
269 if (r != y)
270 {
271 MKL_INT one = 1;
272 dcopy(&mcopysize, (const double*)y, &one, r, &one);
273 }
274
275#ifndef FEAT_USE_MKL_LEGACY_SPMV
276
277 sparse_operation_t opt = SPARSE_OPERATION_NON_TRANSPOSE;
278
279 sparse_matrix_t A;
280 FEAT_DISABLE_WARNINGS
281 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);
282 FEAT_RESTORE_WARNINGS
283 matrix_descr md;
284 md.type = SPARSE_MATRIX_TYPE_GENERAL;
285 mkl_sparse_d_mv(opt, a, A, md, x, b, r);
286 mkl_sparse_destroy(A);
287
288#else
289
290 char trans = 'T';
291 char matdescra[6];
292 matdescra[0] = 'G';
293 matdescra[1] = 0; //ingored by mkl, but valgrind complains otherwise
294 matdescra[2] = 0; //ingored by mkl, but valgrind complains otherwise
295 matdescra[3] = 'C';
296
297 FEAT_DISABLE_WARNINGS
298 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);
299 FEAT_RESTORE_WARNINGS
300
301#endif
302}
303
304void 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)
305{
306 MKL_INT mrows = (MKL_INT)rows;
307 MKL_INT mcolumns = (MKL_INT)columns;
308
309 if (r != y)
310 {
311 MKL_INT one = 1;
312 scopy(&mrows, (const float*)y, &one, r, &one);
313 }
314 cblas_sgemv(CblasRowMajor, CblasNoTrans, mrows, mcolumns, alpha, val, mcolumns, x, 1, beta, r, 1);
315}
316
317void 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)
318{
319 MKL_INT mrows = (MKL_INT)rows;
320 MKL_INT mcolumns = (MKL_INT)columns;
321
322 if (r != y)
323 {
324 MKL_INT one = 1;
325 dcopy(&mrows, (const double*)y, &one, r, &one);
326 }
327 cblas_dgemv(CblasRowMajor, CblasNoTrans, mrows, mcolumns, alpha, val, mcolumns, x, 1, beta, r, 1);
328}
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