FEI  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
test_LinearSystem.cpp
1 /*--------------------------------------------------------------------*/
2 /* Copyright 2005 Sandia Corporation. */
3 /* Under the terms of Contract DE-AC04-94AL85000, there is a */
4 /* non-exclusive license for use of this work by or on behalf */
5 /* of the U.S. Government. Export of this program may require */
6 /* a license from the United States Government. */
7 /*--------------------------------------------------------------------*/
8 
9 
10 #include <fei_macros.hpp>
11 
12 #include <test_utils/test_LinearSystem.hpp>
13 #include <test_utils/test_VectorSpace.hpp>
14 #include <test_utils/test_MatrixGraph.hpp>
15 #include <snl_fei_Factory.hpp>
16 #include <fei_Vector.hpp>
17 #include <fei_Matrix.hpp>
18 #include <snl_fei_LinearSystem_General.hpp>
19 
20 #include <test_utils/LibraryFactory.hpp>
21 
22 #ifdef HAVE_FEI_AZTECOO
23 #include <fei_Aztec_LinSysCore.hpp>
24 #endif
25 
26 #undef fei_file
27 #define fei_file "test_LinearSystem.cpp"
28 #include <fei_ErrMacros.hpp>
29 
30 test_LinearSystem::test_LinearSystem(MPI_Comm comm)
31  : tester(comm)
32 {
33 }
34 
35 test_LinearSystem::~test_LinearSystem()
36 {
37 }
38 
39 int test_LinearSystem::runtests()
40 {
41  CHK_ERR( test1() );
42  CHK_ERR( test2() );
43  CHK_ERR( test3() );
44  CHK_ERR( test4() );
45  CHK_ERR( test5() );
46 
47  return(0);
48 }
49 
50 int test_LinearSystem::test1()
51 {
52  return(0);
53 }
54 
55 int test_LinearSystem::test2()
56 {
57 #ifdef HAVE_FEI_AZTECOO
58  fei::SharedPtr<testData> testdata(new testData(localProc_, numProcs_));
59  std::vector<int>& fieldIDs = testdata->fieldIDs;
60  std::vector<int>& idTypes = testdata->idTypes;
61  std::vector<int>& ids = testdata->ids;
62 
63  fei::SharedPtr<LinearSystemCore> az_lsc(new fei_trilinos::Aztec_LinSysCore(comm_));
64 
65  char* param = new char[64];
66  sprintf(param,"debugOutput .");
67 
68  CHK_ERR( az_lsc->parameters(1, &param) );
69  delete [] param;
70 
71  fei::SharedPtr<fei::Factory> factory(new snl_fei::Factory(comm_, az_lsc));
72 
73  fei::SharedPtr<fei::VectorSpace> vectorSpacePtr =
74  test_VectorSpace::create_VectorSpace(comm_,
75  testdata.get(), localProc_, numProcs_,
76  false, false, "U_LS2", factory);
77 
78  fei::SharedPtr<fei::MatrixGraph> matrixGraphPtr =
79  test_MatrixGraph::create_MatrixGraph(testdata.get(), localProc_, numProcs_,
80  false, false, "U_LS2", vectorSpacePtr,
81  factory, path_);
82 
83  std::vector<int> crIDTypes(2);
84  std::vector<int> crFieldIDs(2);
85  crIDTypes[0] = idTypes[0]; crIDTypes[1] = idTypes[0];
86  crFieldIDs[0] = fieldIDs[0]; crFieldIDs[1] = fieldIDs[0];
87 
88  CHK_ERR( matrixGraphPtr->initLagrangeConstraint(0, idTypes[1],
89  2, //numIDs
90  &crIDTypes[0],
91  &(ids[1]),
92  &crFieldIDs[0]) );
93 
94  CHK_ERR( matrixGraphPtr->initComplete() );
95 
96  fei::SharedPtr<fei::Vector> vec_lsc = factory->createVector(vectorSpacePtr);
97 
98  fei::SharedPtr<fei::Vector> vec_lsc2 = factory->createVector(vectorSpacePtr, true);
99 
100  fei::SharedPtr<fei::Matrix> mat_lsc = factory->createMatrix(matrixGraphPtr);
101 
102  fei::SharedPtr<fei::LinearSystem> linsys = factory->createLinearSystem(matrixGraphPtr);
103  linsys->setMatrix(mat_lsc);
104  linsys->setSolutionVector(vec_lsc2);
105  linsys->setRHS(vec_lsc);
106 
107  int blockID=0;
108  int numIndices = matrixGraphPtr->getConnectivityNumIndices(blockID);
109 
110  std::vector<int> indicesArray(numIndices);
111  int* indicesPtr = &indicesArray[0];
112 
113  int checkNumIndices = 0;
114  CHK_ERR( matrixGraphPtr->getConnectivityIndices(blockID, 0,
115  numIndices, indicesPtr,
116  checkNumIndices) );
117 
118  std::vector<double> data(ids.size(), 1.0);
119  double* dptr = &data[0];
120  std::vector<double*> coefPtrs(ids.size());
121  std::vector<double> crdata(2);
122  crdata[0] = 1.0;
123  crdata[1] = -1.0;
124 
125  for(unsigned ii=0; ii<ids.size(); ++ii) coefPtrs[ii] = dptr;
126 
127  CHK_ERR( mat_lsc->sumIn(numIndices, indicesPtr, numIndices, indicesPtr,
128  &coefPtrs[0]) );
129 
130  CHK_ERR( vec_lsc->sumInFieldData(fieldIDs[0], idTypes[0],
131  ids.size(), &ids[0],
132  &data[0]) );
133 
134  CHK_ERR( linsys->loadLagrangeConstraint(0, &crdata[0], 0.0) );
135 
136  CHK_ERR( mat_lsc->gatherFromOverlap() );
137 
138  CHK_ERR( az_lsc->matrixLoadComplete() );
139 
140  CHK_ERR( linsys->loadComplete() );
141 
142  std::vector<int> crindices;
143  linsys->getConstrainedEqns(crindices);
144  if (crindices.size() != 2) {
145  ERReturn(-7);
146  }
147 
148  CHK_ERR( az_lsc->writeSystem("U_LS2") );
149 
150  MPI_Barrier(comm_);
151 #endif
152 
153  return(0);
154 }
155 
156 int test_LinearSystem::test3()
157 {
158 #ifdef HAVE_FEI_AZTECOO
159  fei::SharedPtr<testData> testdata(new testData(localProc_, numProcs_));
160  std::vector<int>& fieldIDs = testdata->fieldIDs;
161  std::vector<int>& idTypes = testdata->idTypes;
162  std::vector<int>& ids = testdata->ids;
163 
164  fei::SharedPtr<LinearSystemCore> az_lsc(new fei_trilinos::Aztec_LinSysCore(comm_));
165 
166  char* param = new char[64];
167  sprintf(param,"debugOutput .");
168 
169  CHK_ERR( az_lsc->parameters(1, &param) );
170 
171  fei::SharedPtr<fei::Factory> factory(new snl_fei::Factory(comm_, az_lsc));
172 
173  fei::SharedPtr<fei::VectorSpace> vectorSpacePtr =
174  test_VectorSpace::create_VectorSpace(comm_,
175  testdata.get(), localProc_, numProcs_,
176  false, false, "U_LS3", factory);
177 
178  fei::SharedPtr<fei::MatrixGraph> matrixGraphPtr =
179  test_MatrixGraph::create_MatrixGraph(testdata.get(), localProc_, numProcs_,
180  false, false, "U_LS3", vectorSpacePtr,
181  factory, path_);
182 
183  std::vector<int> crIDTypes(2);
184  std::vector<int> crFieldIDs(2);
185  crIDTypes[0] = idTypes[0]; crIDTypes[1] = idTypes[0];
186  crFieldIDs[0] = fieldIDs[0]; crFieldIDs[1] = fieldIDs[0];
187 
188  CHK_ERR( matrixGraphPtr->initPenaltyConstraint(0, idTypes[1],
189  2, //numIDs
190  &crIDTypes[0],
191  &(ids[1]),
192  &crFieldIDs[0]) );
193 
194  CHK_ERR( matrixGraphPtr->initComplete() );
195 
196  fei::SharedPtr<fei::Vector> vec_lsc = factory->createVector(vectorSpacePtr);
197 
198  fei::SharedPtr<fei::Vector> vec_lsc2 = factory->createVector(vectorSpacePtr, true);
199 
200  fei::SharedPtr<fei::Matrix> mat_lsc = factory->createMatrix(matrixGraphPtr);
201 
202  fei::SharedPtr<fei::LinearSystem> linsys = factory->createLinearSystem(matrixGraphPtr);
203  CHK_ERR( linsys->parameters(1, &param) );
204  delete [] param;
205 
206  linsys->setMatrix(mat_lsc);
207  linsys->setSolutionVector(vec_lsc2);
208  linsys->setRHS(vec_lsc);
209 
210  int blockID=0;
211  int numIndices = matrixGraphPtr->getConnectivityNumIndices(blockID);
212 
213  std::vector<int> indicesArray(numIndices);
214  int* indicesPtr = &indicesArray[0];
215 
216  int checkNumIndices = 0;
217  CHK_ERR( matrixGraphPtr->getConnectivityIndices(blockID, 0,
218  numIndices, indicesPtr,
219  checkNumIndices) );
220 
221  std::vector<double> data(ids.size(), 1.0);
222  double* dptr = &data[0];
223  std::vector<double*> coefPtrs(ids.size());
224  std::vector<double> crdata(2);
225  crdata[0] = 1.0;
226  crdata[1] = -1.0;
227 
228  for(unsigned ii=0; ii<ids.size(); ++ii) coefPtrs[ii] = dptr;
229 
230  CHK_ERR( mat_lsc->sumIn(numIndices, indicesPtr, numIndices, indicesPtr,
231  &coefPtrs[0]) );
232 
233  CHK_ERR( vec_lsc->sumInFieldData(fieldIDs[0], idTypes[0],
234  ids.size(), &ids[0],
235  &data[0]) );
236 
237  CHK_ERR( linsys->loadPenaltyConstraint(0, &crdata[0], 100.0, 0.0) );
238 
239  CHK_ERR( mat_lsc->gatherFromOverlap() );
240 
241  CHK_ERR( az_lsc->matrixLoadComplete() );
242 
243  CHK_ERR( linsys->loadComplete() );
244 
245  CHK_ERR( az_lsc->writeSystem("U_LS3") );
246 
247  MPI_Barrier(comm_);
248 #endif
249 
250  return(0);
251 }
252 
253 int test_LinearSystem::test4()
254 {
255 #ifdef HAVE_FEI_AZTECOO
256  fei::SharedPtr<testData> testdata(new testData(localProc_, numProcs_));
257  std::vector<int>& fieldIDs = testdata->fieldIDs;
258  std::vector<int>& idTypes = testdata->idTypes;
259  std::vector<int>& ids = testdata->ids;
260 
261  fei::SharedPtr<LinearSystemCore> az_lsc(new fei_trilinos::Aztec_LinSysCore(comm_));
262 
263  char* param = new char[64];
264  sprintf(param,"debugOutput .");
265 
266  CHK_ERR( az_lsc->parameters(1, &param) );
267  delete [] param;
268 
269  fei::SharedPtr<fei::Factory> factory(new snl_fei::Factory(comm_, az_lsc));
270 
271  fei::SharedPtr<fei::VectorSpace> vectorSpacePtr =
272  test_VectorSpace::create_VectorSpace(comm_,
273  testdata.get(), localProc_, numProcs_,
274  false, false, "U_LS4", factory);
275 
276  fei::SharedPtr<fei::MatrixGraph> matrixGraphPtr =
277  test_MatrixGraph::create_MatrixGraph(testdata.get(), localProc_, numProcs_,
278  false, false, "U_LS4", vectorSpacePtr,
279  factory, path_);
280 
281  std::vector<int> crIDTypes(2);
282  std::vector<int> crFieldIDs(2);
283  crIDTypes[0] = idTypes[0]; crIDTypes[1] = idTypes[0];
284  crFieldIDs[0] = fieldIDs[0]; crFieldIDs[1] = fieldIDs[0];
285 
286  CHK_ERR( matrixGraphPtr->initLagrangeConstraint(0, idTypes[1],
287  2, //numIDs
288  &crIDTypes[0],
289  &(ids[1]),
290  &crFieldIDs[0]) );
291 
292  CHK_ERR( matrixGraphPtr->initComplete() );
293 
294  fei::SharedPtr<fei::Vector> vec_lsc = factory->createVector(vectorSpacePtr);
295 
296  fei::SharedPtr<fei::Vector> vec_lsc2 = factory->createVector(vectorSpacePtr, true);
297 
298  fei::SharedPtr<fei::Matrix> mat_lsc = factory->createMatrix(matrixGraphPtr);
299 
300  fei::SharedPtr<fei::LinearSystem> linsys = factory->createLinearSystem(matrixGraphPtr);
301  linsys->setMatrix(mat_lsc);
302  linsys->setSolutionVector(vec_lsc2);
303  linsys->setRHS(vec_lsc);
304 
305  int blockID=0;
306  int numIndices = matrixGraphPtr->getConnectivityNumIndices(blockID);
307 
308  std::vector<int> indicesArray(numIndices);
309  int* indicesPtr = &indicesArray[0];
310 
311  int checkNumIndices = 0;
312  CHK_ERR( matrixGraphPtr->getConnectivityIndices(blockID, 0,
313  numIndices, indicesPtr,
314  checkNumIndices) );
315 
316  std::vector<double> data(ids.size(), 1.0);
317  double* dptr = &data[0];
318  std::vector<double*> coefPtrs(ids.size());
319  std::vector<double> crdata(2);
320  crdata[0] = 1.0;
321  crdata[1] = -1.0;
322 
323  for(unsigned ii=0; ii<ids.size(); ++ii) coefPtrs[ii] = dptr;
324 
325  CHK_ERR( mat_lsc->sumIn(numIndices, indicesPtr, numIndices, indicesPtr,
326  &coefPtrs[0]) );
327 
328  CHK_ERR( vec_lsc->sumInFieldData(fieldIDs[0], idTypes[0],
329  ids.size(), &ids[0], &data[0]) );
330 
331  CHK_ERR( linsys->loadLagrangeConstraint(0, &crdata[0], 0.0) );
332 
333  CHK_ERR( mat_lsc->gatherFromOverlap() );
334 
335  CHK_ERR( az_lsc->matrixLoadComplete() );
336 
337  CHK_ERR( linsys->loadComplete() );
338 
339  CHK_ERR( az_lsc->writeSystem("U_LS4") );
340 
341  MPI_Barrier(comm_);
342 #endif
343 
344  return(0);
345 }
346 
347 int test_LinearSystem::test5()
348 {
349  return(0);
350 }
351 
virtual void setMatrix(fei::SharedPtr< fei::Matrix > &matrix)
virtual int loadLagrangeConstraint(int constraintID, const double *weights, double rhsValue)=0
virtual int parameters(int numParams, const char *const *paramStrings)=0
virtual void setSolutionVector(fei::SharedPtr< fei::Vector > &soln)
virtual int loadPenaltyConstraint(int constraintID, const double *weights, double penaltyValue, double rhsValue)=0
virtual void setRHS(fei::SharedPtr< fei::Vector > &rhs)
virtual int loadComplete(bool applyBCs=true, bool globalAssemble=true)=0
virtual int sumIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)=0
virtual void getConstrainedEqns(std::vector< int > &crEqns) const =0
virtual int gatherFromOverlap(bool accumulate=true)=0