2
// ************************************************************************
4
// Kokkos: A Fast Kernel Package
5
// Copyright (2004) Sandia Corporation
7
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8
// license for use of this work by or on behalf of the U.S. Government.
10
// This library is free software; you can redistribute it and/or modify
11
// it under the terms of the GNU Lesser General Public License as
12
// published by the Free Software Foundation; either version 2.1 of the
13
// License, or (at your option) any later version.
15
// This library is distributed in the hope that it will be useful, but
16
// WITHOUT ANY WARRANTY; without even the implied warranty of
17
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18
// Lesser General Public License for more details.
20
// You should have received a copy of the GNU Lesser General Public
21
// License along with this library; if not, write to the Free Software
22
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
24
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
26
// ************************************************************************
29
#ifndef KOKKOS_TEST_GENERATEOSKIPROBLEM_H
30
#define KOKKOS_TEST_GENERATEOSKIPROBLEM_H
31
#include "Kokkos_OskiVector.hpp"
32
#include "Kokkos_OskiMultiVector.hpp"
33
#include "Kokkos_OskiMatrix.hpp"
35
namespace KokkosTest {
38
//! KokkosTest::GenerateOskiProblem: Constructs a 2D PDE finite difference matrix using the list of x and y offsets.
41
nx (In) - number of grid points in x direction
43
ny (In) - number of grid points in y direction
45
The total number of equations will be nx*ny ordered such that the x direction changes
47
First equation is at point (0,0)
50
nx equation at (nx-1,0)
51
nx+1st equation at (0,1)
53
npoints (In) - number of points in finite difference stencil
54
xoff (In) - stencil offsets in x direction (of length npoints)
55
yoff (In) - stencil offsets in y direction (of length npoints)
56
A standard 5-point finite difference stencil would be described as:
58
xoff = [-1, 1, 0, 0, 0]
59
yoff = [ 0, 0, 0, -1, 1]
61
nrhs - Number of rhs to generate. (First interface produces vectors, so nrhs is not needed
63
A (Out) - Kokkos::OskiMatrix constructed for nx by ny grid using prescribed stencil
64
Off-diagonal values are random between 0 and 1. If diagonal is part of stencil,
65
diagonal will be slightly diag dominant.
66
x (Out) - Initial guess vector set to zero.
67
b (Out) - Generated RHS. Values satisfy b = A*xexact
68
xexact (Out) - Generated exact solution to Ax = b.
70
Note: Caller of this function is responsible for deleting all output objects.
73
template<typename OrdinalType, typename ScalarType>
74
class GenerateOskiProblem {
77
//@{ \name Constructors/Destructor.
79
//! Single RHS constuctor.
80
/* don't need this because Kokkos::OskiVector inherits from Kokkos::OskiMultivector
81
GenerateOskiProblem(bool isRowOriented,
82
OrdinalType nx, OrdinalType ny, OrdinalType npoints,
83
OrdinalType * xoff, OrdinalType * yoff,
84
Kokkos::CisMatrix<OrdinalType, ScalarType> *& A,
85
Kokkos::Vector<OrdinalType, ScalarType> *& x,
86
Kokkos::Vector<OrdinalType, ScalarType> *& b,
87
Kokkos::Vector<OrdinalType, ScalarType> *&xexact,
88
OrdinalType & numEntries);
90
//! Multi RHS constuctor.
91
GenerateOskiProblem(bool isRowOriented,
92
OrdinalType nx, OrdinalType ny, OrdinalType npoints,
93
OrdinalType * xoff, OrdinalType * yoff, OrdinalType nrhs,
94
Kokkos::CisMatrix<OrdinalType, ScalarType> *& A,
95
Kokkos::MultiVector<OrdinalType, ScalarType> *& x,
96
Kokkos::MultiVector<OrdinalType, ScalarType> *& b,
97
Kokkos::MultiVector<OrdinalType, ScalarType> *&xexact,
98
OrdinalType & numEntries);
102
//! OskiMatrix Destructor
103
virtual ~GenerateOskiProblem();
108
//! Copy constructor (Not implemented).
109
GenerateOskiProblem(const GenerateOskiProblem<OrdinalType, ScalarType>& source){};
111
/* void GenerateProblem(bool isRowOriented,
112
OrdinalType nx, OrdinalType ny, OrdinalType npoints,
113
OrdinalType * xoff, OrdinalType * yoff,
114
Kokkos::CisMatrix<OrdinalType, ScalarType> *& A,
115
Kokkos::Vector<OrdinalType, ScalarType> *& x,
116
Kokkos::Vector<OrdinalType, ScalarType> *& b,
117
Kokkos::Vector<OrdinalType, ScalarType> *&xexact,
118
OrdinalType & numEntries);
120
void GenerateProblem(bool isRowOriented,
121
OrdinalType nx, OrdinalType ny, OrdinalType npoints,
122
OrdinalType * xoff, OrdinalType * yoff, OrdinalType nrhs,
123
Kokkos::CisMatrix<OrdinalType, ScalarType> *& A,
124
Kokkos::MultiVector<OrdinalType, ScalarType> *& x,
125
Kokkos::MultiVector<OrdinalType, ScalarType> *& b,
126
Kokkos::MultiVector<OrdinalType, ScalarType> *&xexact,
127
OrdinalType & numEntries);
131
Kokkos::OskiVector<OrdinalType, ScalarType> * xd1;
132
Kokkos::OskiVector<OrdinalType, ScalarType> *bd1;
133
Kokkos::OskiVector<OrdinalType, ScalarType> *xexactd1;
134
Kokkos::OskiMultiVector<OrdinalType, ScalarType> * xd;
135
Kokkos::OskiMultiVector<OrdinalType, ScalarType> * bd;
136
Kokkos::OskiMultiVector<OrdinalType, ScalarType> * xexactd;
137
Kokkos::OskiMatrix<OrdinalType, ScalarType> * Ad; // Construct matrix
140
ScalarType ** xexactv;
141
OrdinalType numEquations;
144
OrdinalType ** indices;
145
//ScalarType ** values;
146
OrdinalType * pointers;
147
OrdinalType * allIndices;
148
ScalarType * allValues;
149
//OrdinalType * profiles;
151
} // namespace KokkosTest
153
using namespace KokkosTest;
154
//==============================================================================
155
/*template<typename OrdinalType, typename ScalarType>
156
GenerateOskiProblem<OrdinalType, ScalarType>::
157
GenerateOskiProblem(bool isRowOriented,
158
OrdinalType nx, OrdinalType ny, OrdinalType npoints,
159
OrdinalType * xoff, OrdinalType * yoff,
160
Kokkos::CisMatrix<OrdinalType, ScalarType> *& A,
161
Kokkos::Vector<OrdinalType, ScalarType> *& x,
162
Kokkos::Vector<OrdinalType, ScalarType> *& b,
163
Kokkos::Vector<OrdinalType, ScalarType> *&xexact,
164
OrdinalType & numEntries):
183
GenerateProblem(isRowOriented, nx, ny, npoints, xoff, yoff, A, x, b, xexact, numEntries);
186
//==============================================================================
187
template<typename OrdinalType, typename ScalarType>
188
GenerateOskiProblem<OrdinalType, ScalarType>::
189
GenerateOskiProblem(bool isRowOriented,
190
OrdinalType nx, OrdinalType ny, OrdinalType npoints,
191
OrdinalType * xoff, OrdinalType * yoff, OrdinalType nrhs,
192
Kokkos::CisMatrix<OrdinalType, ScalarType> *& A,
193
Kokkos::MultiVector<OrdinalType, ScalarType> *& x,
194
Kokkos::MultiVector<OrdinalType, ScalarType> *& b,
195
Kokkos::MultiVector<OrdinalType, ScalarType> *&xexact,
196
OrdinalType & numEntries):
215
GenerateProblem(isRowOriented, nx, ny, npoints, xoff, yoff, nrhs, A, x, b, xexact, numEntries);
218
//==============================================================================
219
template<typename OrdinalType, typename ScalarType>
220
void GenerateOskiProblem<OrdinalType, ScalarType>::
221
GenerateProblem(bool isRowOriented,
222
OrdinalType nx, OrdinalType ny, OrdinalType npoints,
223
OrdinalType * xoff, OrdinalType * yoff,
224
Kokkos::CisMatrix<OrdinalType, ScalarType> *& A,
225
Kokkos::Vector<OrdinalType, ScalarType> *& x,
226
Kokkos::Vector<OrdinalType, ScalarType> *& b,
227
Kokkos::Vector<OrdinalType, ScalarType> *&xexact,
228
OrdinalType & numEntries) {
230
Kokkos::MultiVector<OrdinalType, ScalarType> * x1, * b1, * xexact1;
232
GenerateProblem(isRowOriented, nx, ny, npoints, xoff, yoff, 1, A, x1, b1, xexact1, numEntries);
234
xd1 = new Kokkos::OskiVector<OrdinalType, ScalarType>();
235
bd1 = new Kokkos::OskiVector<OrdinalType, ScalarType>();
236
xexactd1 = new Kokkos::OskiVector<OrdinalType, ScalarType>();
238
xd1->initializeValues(x1->getNumRows(), x1->getValues(0), x1->getColInc());
239
bd1->initializeValues(b1->getNumRows(), b1->getValues(0), b1->getColInc());
240
xexactd1->initializeValues(xexact1->getNumRows(), xexact1->getValues(0), xexact1->getColInc());
242
x = dynamic_cast<Kokkos::Vector<OrdinalType, ScalarType> *>(xd1);
243
b = dynamic_cast<Kokkos::Vector<OrdinalType, ScalarType> *>(bd1);
244
xexact = dynamic_cast<Kokkos::Vector<OrdinalType, ScalarType> *>(xexactd1);
249
template<typename OrdinalType, typename ScalarType>
250
void GenerateOskiProblem<OrdinalType, ScalarType>::
251
GenerateProblem(bool isRowOriented,
252
OrdinalType nx, OrdinalType ny, OrdinalType npoints,
253
OrdinalType * xoff, OrdinalType * yoff, OrdinalType nrhs,
254
Kokkos::CisMatrix<OrdinalType, ScalarType> *& A,
255
Kokkos::MultiVector<OrdinalType, ScalarType> *& x,
256
Kokkos::MultiVector<OrdinalType, ScalarType> *& b,
257
Kokkos::MultiVector<OrdinalType, ScalarType> *&xexact,
258
OrdinalType & numEntries) {
260
// Number of equations is nx*ny.
261
numEquations = nx*ny;
265
xv = new ScalarType*[nrhs];
266
bv = new ScalarType*[nrhs];
267
xexactv = new ScalarType*[nrhs];
269
for (i=0; i<nrhs; i++) {
270
xv[i] = new ScalarType[numEquations];
271
bv[i] = new ScalarType[numEquations];
272
xexactv[i] = new ScalarType[numEquations];
275
for (i=0; i<nrhs; i++)
276
for (j=0; j<numEquations; j++) {
277
xexactv[i][j] = ((ScalarType) rand())/ ((ScalarType) RAND_MAX);
288
allIndices = new OrdinalType[numEquations*npoints];
289
allValues = new ScalarType[numEquations*npoints];
290
pointers = new OrdinalType[numEquations+1];
292
OrdinalType * curIndices;
293
ScalarType * curValues;
295
ScalarType dnpoints = (ScalarType) npoints;
298
for (i=0; i<numEquations; i++) {
300
OrdinalType rowID = i;
301
OrdinalType numIndices = 0;
303
curIndices = allIndices+numEntries;
304
curValues = allValues+numEntries;
306
for (j=0; j<npoints; j++) {
307
OrdinalType colID = rowID + xoff[j] + nx*yoff[j]; // Compute column ID based on stencil offsets
308
if (colID>-1 && colID<numEquations) {
309
curIndices[numIndices] = colID;
310
ScalarType value = - ((ScalarType) rand())/ ((ScalarType) RAND_MAX);
312
curValues[numIndices] = dnpoints - value; // Make diagonal dominant
314
curValues[numIndices] = -value;
316
for (k=0; k<nrhs; k++)
317
bv[k][i] += curValues[numIndices]*xexactv[k][curIndices[numIndices]];
319
for (k=0; k<nrhs; k++)
320
bv[k][curIndices[numIndices]] += curValues[numIndices]*xexactv[k][i];
325
pointers[i+1] = pointers[i]+numIndices;
327
// Now construct Kokkos objects and register pointers
328
xd = new Kokkos::OskiMultiVector<OrdinalType, ScalarType>();
329
bd = new Kokkos::OskiMultiVector<OrdinalType, ScalarType>();
330
xexactd = new Kokkos::OskiMultiVector<OrdinalType, ScalarType>();
331
Ad = new Kokkos::OskiMatrix<OrdinalType, ScalarType>(); // Construct matrix
333
xd->initializeValues(numEquations, nrhs, xv[0], numEquations, 1);
334
bd->initializeValues(numEquations, nrhs, bv[0], numEquations, 1);
335
xexactd->initializeValues(numEquations, nrhs, xexactv[0], numEquations, 1);
337
Ad->initializeStructure(numEquations, numEquations, isRowOriented, pointers, allIndices);
338
Ad->initializeValues(allValues);
340
x = dynamic_cast<Kokkos::MultiVector<OrdinalType, ScalarType> *>(xd);
341
b = dynamic_cast<Kokkos::MultiVector<OrdinalType, ScalarType> *>(bd);
342
xexact = dynamic_cast<Kokkos::MultiVector<OrdinalType, ScalarType> *>(xexactd);
343
A = dynamic_cast<Kokkos::CisMatrix<OrdinalType, ScalarType> *>(Ad);
348
template<typename OrdinalType, typename ScalarType>
349
GenerateOskiProblem<OrdinalType, ScalarType>::
350
~GenerateOskiProblem() {
353
if (xd1!=0) delete xd1;
354
if (bd1!=0) delete bd1;
355
if (xexactd1!=0) delete xexactd1;
356
if (xd!=0) delete xd;
357
if (bd!=0) delete bd;
358
if (xexactd!=0) delete xexactd;
359
if (Ad!=0) delete Ad;
362
for (i=0; i<nrhsv; i++) if (xv[i]!=0) delete [] xv[i];
366
for (i=0; i<nrhsv; i++) if (bv[i]!=0) delete [] bv[i];
370
for (i=0; i<nrhsv; i++) if (xexactv[i]!=0) delete [] xexactv[i];
374
for (i=0; i<numEquations; i++) if (indices[i]!=0) delete [] indices[i];
378
for (i=0; i<numEquations; i++) if (values[i]!=0) delete [] values[i];
381
if (pointers!=0) delete [] pointers;
382
if (allIndices!=0) delete [] allIndices;
383
if (allValues!=0) delete [] allValues;
385
#endif /* KOKKOS_TEST_GENERATEOSKIPROBLEM_H */