42 #ifndef KOKKOS_SEQUENTIAL_ADDSPARSEMATRICES_HPP
43 #define KOKKOS_SEQUENTIAL_ADDSPARSEMATRICES_HPP
45 #include <Kokkos_ConfigDefs.hpp>
49 namespace Sequential {
70 template<
class OffsetType,
class OrdinalType>
72 countMergedRowEntries (
const OffsetType ptr1[],
73 const OrdinalType ind1[],
74 const OffsetType ptr2[],
75 const OrdinalType ind2[],
104 template<
class OffsetType,
107 class BinaryFunctionType>
109 mergeRows (
const OffsetType ptrOut[],
110 OrdinalType indOut[],
112 const OffsetType ptr1[],
113 const OrdinalType ind1[],
114 const ScalarType val1[],
115 const OffsetType ptr2[],
116 const OrdinalType ind2[],
117 const ScalarType val2[],
119 BinaryFunctionType f);
126 template<
class ScalarType>
130 alpha_ (alpha), beta_ (beta) {}
133 operator () (
const ScalarType& A_ij,
const ScalarType& B_ij)
const {
134 return alpha_ * A_ij + beta_ * B_ij;
137 const ScalarType alpha_;
138 const ScalarType beta_;
160 template<
class OffsetType,
class OrdinalType,
class ScalarType>
162 addSparseMatrices (OffsetType*& ptrResult,
163 OrdinalType*& indResult,
164 ScalarType*& valResult,
165 const ScalarType alpha,
166 const OffsetType ptr1[],
167 const OrdinalType ind1[],
168 const ScalarType val1[],
169 const ScalarType beta,
170 const OffsetType ptr2[],
171 const OrdinalType ind2[],
172 const ScalarType val2[],
173 const OrdinalType numRows);
175 template<
class OffsetType,
class OrdinalType>
177 countMergedRowEntries (
const OffsetType ptr1[],
178 const OrdinalType ind1[],
179 const OffsetType ptr2[],
180 const OrdinalType ind2[],
183 const OffsetType start1 = ptr1[i];
184 const OffsetType end1 = ptr1[i+1];
185 const OffsetType start2 = ptr2[i];
186 const OffsetType end2 = ptr2[i+1];
188 OffsetType mark1 = start1, mark2 = start2;
189 OrdinalType count = 0;
190 while (mark1 < end1 && mark2 < end2) {
191 if (ind1[mark1] == ind2[mark2]) {
194 }
else if (ind1[mark1] < ind2[mark2]) {
202 count += end1 - mark1;
203 count += end2 - mark2;
207 template<
class OffsetType,
210 class BinaryFunctionType>
212 mergeRows (
const OffsetType ptrOut[],
213 OrdinalType indOut[],
215 const OffsetType ptr1[],
216 const OrdinalType ind1[],
217 const ScalarType val1[],
218 const OffsetType ptr2[],
219 const OrdinalType ind2[],
220 const ScalarType val2[],
222 BinaryFunctionType f)
224 const OffsetType startOut = ptrOut[i];
225 const OffsetType endOut = ptrOut[i+1];
226 const OffsetType start1 = ptr1[i];
227 const OffsetType end1 = ptr1[i+1];
228 const OffsetType start2 = ptr2[i];
229 const OffsetType end2 = ptr2[i+1];
232 OffsetType mark1 = start1, mark2 = start2, markOut = startOut;
233 while (mark1 < end1 && mark2 < end2 && markOut < endOut) {
234 if (ind1[mark1] == ind2[mark2]) {
235 indOut[markOut] = ind1[mark1];
236 valOut[markOut] = f (val1[mark1], val2[mark2]);
239 }
else if (ind1[mark1] < ind2[mark2]) {
240 indOut[markOut] = ind1[mark1];
242 valOut[markOut] = f (val1[mark1], zero);
245 indOut[markOut] = ind2[mark2];
247 valOut[markOut] = f (zero, val2[mark2]);
253 while (mark1 < end1 && markOut < endOut) {
254 indOut[markOut] = ind1[mark1];
256 valOut[markOut] = f (val1[mark1], zero);
260 while (mark2 < end2 && markOut < endOut) {
261 indOut[markOut] = ind2[mark2];
263 valOut[markOut] = f (zero, val2[mark2]);
271 markOut >= endOut && (mark1 < end1 || mark2 < end2),
273 "Kokkos::Sequential::mergeRows: Row " << i <<
" of the output array has "
274 << (end1 - mark1) <<
" + " << (end2 - mark2) <<
" too few entries.");
278 template<
class OffsetType,
class OrdinalType,
class ScalarType>
280 addSparseMatrices (OffsetType*& ptrResult,
281 OrdinalType*& indResult,
282 ScalarType*& valResult,
283 const ScalarType alpha,
284 const OffsetType ptr1[],
285 const OrdinalType ind1[],
286 const ScalarType val1[],
287 const ScalarType beta,
288 const OffsetType ptr2[],
289 const OrdinalType ind2[],
290 const ScalarType val2[],
291 const OrdinalType numRows)
299 OffsetType* ptrOut = NULL;
300 OrdinalType* indOut = NULL;
301 ScalarType* valOut = NULL;
308 ptrOut =
new OffsetType [numRows + 1];
315 if (alpha == STS::zero ()) {
316 if (beta == STS::zero ()) {
319 std::fill (ptrOut, ptrOut + numRows + 1, Teuchos::as<OffsetType> (0));
326 memcpy (ptrOut, ptr2, (numRows+1) *
sizeof (OffsetType));
327 const OffsetType numEntries = ptr2[numRows];
328 indOut =
new OrdinalType [numEntries];
329 memcpy (indOut, ind2, numEntries *
sizeof (OrdinalType));
330 valOut =
new ScalarType [numEntries];
331 if (beta == STS::one ()) {
332 memcpy (valOut, val2, numEntries *
sizeof (ScalarType));
334 for (OrdinalType k = 0; k < numEntries; ++k) {
335 valOut[k] = beta * val2[k];
344 else if (beta == STS::zero ()) {
346 memcpy (ptrOut, ptr1, (numRows+1) *
sizeof (OffsetType));
347 const OffsetType numEntries = ptr1[numRows];
348 indOut =
new OrdinalType [numEntries];
349 memcpy (indOut, ind1, numEntries *
sizeof (OrdinalType));
350 valOut =
new ScalarType [numEntries];
351 if (alpha == STS::one ()) {
352 memcpy (valOut, val1, numEntries *
sizeof (ScalarType));
354 for (OrdinalType k = 0; k < numEntries; ++k) {
355 valOut[k] = alpha * val1[k];
365 for (OrdinalType i = 0; i < numRows; ++i) {
366 ptrOut[i+1] = countMergedRowEntries (ptr1, ind1, ptr2, ind2, i);
370 for (OrdinalType i = 1; i < numRows; ++i) {
371 ptrOut[i+1] += ptrOut[i];
375 if (ptrOut != NULL) {
378 if (indOut != NULL) {
381 if (valOut != NULL) {
389 const OffsetType numEntries = ptrOut[numRows];
393 indOut =
new OrdinalType [numEntries];
394 valOut =
new ScalarType [numEntries];
400 if (alpha == STS::one () && beta == STS::one ()) {
401 for (OrdinalType i = 0; i < numRows; ++i) {
402 (void) mergeRows (ptrOut, indOut, valOut,
405 std::plus<ScalarType> ());
408 AddMatrixEntries<ScalarType> f (alpha, beta);
409 for (OrdinalType i = 0; i < numRows; ++i) {
410 (void) mergeRows (ptrOut, indOut, valOut,
417 if (indOut != NULL) {
420 if (valOut != NULL) {
423 if (ptrOut != NULL) {
438 #endif // KOKKOS_SEQUENTIAL_ADDSPARSEMATRICES_HPP
Compute C_ij = alpha * A_ij + beta * B_ij.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)