Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
UserInputForTests.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Zoltan2: A package of combinatorial algorithms for scientific computing
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Karen Devine (kddevin@sandia.gov)
39 // Erik Boman (egboman@sandia.gov)
40 // Siva Rajamanickam (srajama@sandia.gov)
41 //
42 // ***********************************************************************
43 //
44 // @HEADER
45 
50 #ifndef USERINPUTFORTESTS
51 #define USERINPUTFORTESTS
52 
53 #include "Zoltan2_TestHelpers.hpp"
54 #include <Zoltan2_XpetraTraits.hpp>
55 #include <Zoltan2_Typedefs.hpp>
56 
57 #include <Tpetra_MultiVector.hpp>
58 #include <Tpetra_CrsMatrix.hpp>
59 #include <Tpetra_Map.hpp>
60 #include <Xpetra_Vector.hpp>
61 #include <Xpetra_CrsMatrix.hpp>
62 #include <Xpetra_CrsGraph.hpp>
63 
64 #include <MatrixMarket_Tpetra.hpp>
65 
66 #ifdef HAVE_ZOLTAN2_GALERI
67 #include <Galeri_XpetraProblemFactory.hpp>
68 #include <Galeri_XpetraParameters.hpp>
69 #endif
70 
71 #include <Kokkos_DefaultNode.hpp>
72 
73 #include "GeometricGenerator.hpp"
74 #include <fstream>
75 #include <string>
76 
77 #include <TpetraExt_MatrixMatrix_def.hpp>
78 
79 // pamgen required includes
81 
82 
83 using Teuchos::RCP;
84 using Teuchos::ArrayRCP;
85 using Teuchos::ArrayView;
86 using Teuchos::Array;
87 using Teuchos::Comm;
88 using Teuchos::rcp;
89 using Teuchos::arcp;
90 using Teuchos::rcp_const_cast;
91 using Teuchos::ParameterList;
92 using namespace Zoltan2_TestingFramework;
93 
125 
127 {
128 public:
129 
130  typedef Tpetra::Map<zlno_t, zgno_t, znode_t> map_t;
131  typedef Tpetra::Export<zlno_t, zgno_t, znode_t> export_t;
132  typedef Tpetra::Import<zlno_t, zgno_t, znode_t> import_t;
133  typedef map_t::node_type default_znode_t;
134 
135 
155  UserInputForTests(string path, string testData,
156  const RCP<const Comm<int> > &c, bool debugInfo=false,
157  bool distributeInput=true);
158 
175  UserInputForTests(int x, int y, int z, string matrixType,
176  const RCP<const Comm<int> > &c, bool debugInfo=false,
177  bool distributeInput=true);
178 
194  UserInputForTests(const ParameterList &pList,
195  const RCP<const Comm<int> > &c);
196 
199  static void getUIRandomData(unsigned int seed, zlno_t length,
200  zscalar_t min, zscalar_t max, ArrayView<ArrayRCP<zscalar_t > > data);
201 
202  RCP<tMVector_t> getUICoordinates();
203 
204  RCP<tMVector_t> getUIWeights();
205 
206  RCP<tMVector_t> getUIEdgeWeights();
207 
208  RCP<tcrsMatrix_t> getUITpetraCrsMatrix();
209 
210  RCP<tcrsGraph_t> getUITpetraCrsGraph();
211 
212  RCP<tVector_t> getUITpetraVector();
213 
214  RCP<tMVector_t> getUITpetraMultiVector(int nvec);
215 
216  RCP<xcrsMatrix_t> getUIXpetraCrsMatrix();
217 
218  RCP<xcrsGraph_t> getUIXpetraCrsGraph();
219 
220  RCP<xVector_t> getUIXpetraVector();
221 
222  RCP<xMVector_t> getUIXpetraMultiVector(int nvec);
223 
224 #ifdef HAVE_ZOLTAN2_PAMGEN
225  PamgenMesh * getPamGenMesh(){return this->pamgen_mesh.operator->();}
226 #endif
227 
228 #ifdef HAVE_EPETRA_DATA_TYPES
229  RCP<Epetra_CrsGraph> getUIEpetraCrsGraph();
230 
231  RCP<Epetra_CrsMatrix> getUIEpetraCrsMatrix();
232 
233  RCP<Epetra_Vector> getUIEpetraVector();
234 
235  RCP<Epetra_MultiVector> getUIEpetraMultiVector(int nvec);
236 #endif
237  bool hasInput();
238 
239  bool hasInputDataType(const string &input_type);
240 
241  bool hasUICoordinates();
242 
243  bool hasUIWeights();
244 
245  bool hasUIEdgeWeights();
246 
247  bool hasUITpetraCrsMatrix();
248 
249  bool hasUITpetraCrsGraph();
250 
251  bool hasUITpetraVector();
252 
253  bool hasUITpetraMultiVector();
254 
255  bool hasUIXpetraCrsMatrix();
256 
257  bool hasUIXpetraCrsGraph();
258 
259  bool hasUIXpetraVector();
260 
261  bool hasUIXpetraMultiVector();
262 
263  bool hasPamgenMesh();
264 #ifdef HAVE_EPETRA_DATA_TYPES
265  bool hasUIEpetraCrsGraph();
266 
267  bool hasUIEpetraCrsMatrix();
268 
269  bool hasUIEpetraVector();
270 
271  bool hasUIEpetraMultiVector();
272 
273 #endif
274 
275 private:
276 
277  bool verbose_;
278 
279  const RCP<const Comm<int> > tcomm_;
280 
281  bool havePamgenMesh;
282 #ifdef HAVE_ZOLTAN2_PAMGEN
283  RCP<PamgenMesh> pamgen_mesh;
284 #endif
285 
286  RCP<tcrsMatrix_t> M_;
287  RCP<xcrsMatrix_t> xM_;
288 
289  RCP<tMVector_t> xyz_;
290  RCP<tMVector_t> vtxWeights_;
291  RCP<tMVector_t> edgWeights_;
292 
293 #ifdef HAVE_EPETRA_DATA_TYPES
294  RCP<const Epetra_Comm> ecomm_;
295  RCP<Epetra_CrsMatrix> eM_;
296  RCP<Epetra_CrsGraph> eG_;
297 #endif
298 
299  // Read a Matrix Market file into M_
300  // using Tpetra::MatrixMarket::Reader.
301  // If there are "Tim Davis" style coordinates
302  // that go with the file, read those into xyz_.
303 
304  void readMatrixMarketFile(string path, string testData,bool distributeInput = true);
305 
306  // Build matrix M_ from a mesh and a problem type
307  // with Galeri::Xpetra.
308 
309  void buildCrsMatrix(int xdim, int ydim, int zdim, string type,
310  bool distributeInput);
311 
312  // Read a Zoltan1 Chaco or Matrix Market file
313  // into M_. If it has geometric coordinates,
314  // read them into xyz_. If it has weights,
315  // read those into vtxWeights_ and edgWeights_.
316  void readZoltanTestData(string path, string testData,
317  bool distributeInput);
318 
319  // Modify the Maps of an input matrix to make them non-contiguous
320  RCP<tcrsMatrix_t> modifyMatrixGIDs(RCP<tcrsMatrix_t> &in);
321  inline zgno_t newID(const zgno_t id) { return id * 2 + 10001; }
322 
323  // Read Zoltan data that is in a .graph file.
324  void getUIChacoGraph(FILE *fptr, bool haveAssign, FILE *assignFile,
325  string name, bool distributeInput);
326 
327  // Read Zoltan data that is in a .coords file.
328  void getUIChacoCoords(FILE *fptr, string name);
329 
330  // Chaco reader code: This code is copied from zoltan/ch.
331  // It might benefit from a rewrite and simplification.
332 
333  // Chaco reader helper functions: copied from zoltan/ch
334  static const int CHACO_LINE_LENGTH=200;
335  char chaco_line[CHACO_LINE_LENGTH]; /* space to hold values */
336  int chaco_offset; /* offset into line for next data */
337  int chaco_break_pnt; /* place in sequence to pause */
338  int chaco_save_pnt; /* place in sequence to save */
339 
340  double chaco_read_val(FILE* infile, int *end_flag);
341  int chaco_read_int(FILE* infile, int *end_flag);
342  void chaco_flush_line(FILE*);
343 
344  // Chaco graph reader: copied from zoltan/ch
345  int chaco_input_graph(FILE *fin, const char *inname, int **start,
346  int **adjacency, int *nvtxs, int *nVwgts,
347  float **vweights, int *nEwgts, float **eweights);
348 
349  // Chaco coordinate reader: copied from zoltan/ch
350  int chaco_input_geom(FILE *fingeom, const char *geomname, int nvtxs,
351  int *igeom, double **x, double **y, double **z);
352 
353  // Chaco coordinate reader: copied from zoltan/ch
354  int chaco_input_assign(FILE *finassign, const char *assignname, int nvtxs,
355  short *assignments);
356 
357 
358  // Read a GeomGen.txt file into M_
359  // Read coordinates into xyz_.
360  // If iti has weights read those to vtxWeights_
361  // and edgeWeights_
362  void readGeometricGenTestData(string path, string testData);
363 
364  // Geometry Gnearatory helper function
365  void readGeoGenParams(string paramFileName,
366  ParameterList &geoparams);
367 
368  // utility methods used when reading geom gen files
369 
370  static string trim_right_copy(const string& s,
371  const string& delimiters = " \f\n\r\t\v" );
372 
373  static string trim_left_copy(const string& s,
374  const string& delimiters = " \f\n\r\t\v" );
375 
376  static string trim_copy(const string& s,
377  const string& delimiters = " \f\n\r\t\v" );
378 
379 
380  // Read a pamgen mesh
381  void readPamgenMeshFile(string path, string testData);
382 #ifdef HAVE_ZOLTAN2_PAMGEN
383  void setPamgenAdjacencyGraph();
384  void setPamgenCoordinateMV();
385 #endif
386 };
387 
388 UserInputForTests::UserInputForTests(string path, string testData,
389  const RCP<const Comm<int> > &c,
390  bool debugInfo, bool distributeInput):
391 verbose_(debugInfo), tcomm_(c), havePamgenMesh(false),
392 M_(), xM_(), xyz_(), vtxWeights_(), edgWeights_(),
394 ecomm_(), eM_(), eG_(),
395 #endif
396 chaco_offset(0), chaco_break_pnt(CHACO_LINE_LENGTH)
397 {
398  bool zoltan1 = false;
399  string::size_type loc = path.find("/zoltan/test/"); // Zoltan1 data
400  if (loc != string::npos)
401  zoltan1 = true;
402 
403  if (zoltan1)
404  readZoltanTestData(path, testData, distributeInput);
405  else
406  readMatrixMarketFile(path, testData);
407 
408 #ifdef HAVE_EPETRA_DATA_TYPES
409  ecomm_ = Xpetra::toEpetra(c);
410 #endif
411 }
412 
414  string matrixType,
415  const RCP<const Comm<int> > &c,
416  bool debugInfo,
417  bool distributeInput):
418 verbose_(debugInfo), tcomm_(c), havePamgenMesh(false),
419 M_(), xM_(), xyz_(), vtxWeights_(), edgWeights_(),
421 ecomm_(), eM_(), eG_(),
422 #endif
423 chaco_offset(0), chaco_break_pnt(CHACO_LINE_LENGTH)
424 {
425  if (matrixType.size() == 0){
426  int dim = 0;
427  if (x > 0) dim++;
428  if (y > 0) dim++;
429  if (z > 0) dim++;
430  if (dim == 1)
431  matrixType = string("Laplace1D");
432  else if (dim == 2)
433  matrixType = string("Laplace2D");
434  else if (dim == 3)
435  matrixType = string("Laplace3D");
436  else
437  throw std::runtime_error("input");
438 
439  if (verbose_ && tcomm_->getRank() == 0)
440  std::cout << "UserInputForTests, Matrix type : " << matrixType << std::endl;
441  }
442 
443  buildCrsMatrix(x, y, z, matrixType, distributeInput);
444 
445 #ifdef HAVE_EPETRA_DATA_TYPES
446  ecomm_ = Xpetra::toEpetra(c);
447 #endif
448 }
449 
450 UserInputForTests::UserInputForTests(const ParameterList &pList,
451  const RCP<const Comm<int> > &c):
452 tcomm_(c), havePamgenMesh(false),
453 M_(), xM_(), xyz_(), vtxWeights_(), edgWeights_(),
455 ecomm_(), eM_(), eG_(),
456 #endif
457 chaco_offset(0), chaco_break_pnt(CHACO_LINE_LENGTH)
458 {
459 
460  // get options
461  bool distributeInput = true, debugInfo = true;
462 
463  if(pList.isParameter("distribute input"))
464  distributeInput = pList.get<bool>("distribute input");
465 
466  if(pList.isParameter("debug"))
467  debugInfo = pList.get<bool>("debug");
468  this->verbose_ = debugInfo;
469 
470  if(pList.isParameter("input file"))
471  {
472 
473  // get input path
474  string path(".");
475  if(pList.isParameter("input path"))
476  path = pList.get<string>("input path");
477 
478  string testData = pList.get<string>("input file");
479 
480  // find out if we are working from the zoltan1 test diretory
482 
483  // find out if we are using the geometric generator
484  if(pList.isParameter("file type") && pList.get<string>("file type") == "Geometric Generator")
485  file_format = GEOMGEN;
486  else if(pList.isParameter("file type") && pList.get<string>("file type") == "Pamgen")
487  {
488  file_format = PAMGEN;
489  }
490  else if(pList.isParameter("file type") && pList.get<string>("file type") == "Chaco")
491  file_format = CHACO; // this flag calls read ZoltanTestData, which calls the chaco readers...
492 
493  // read the input file
494  switch (file_format) {
495  case GEOMGEN: readGeometricGenTestData(path,testData); break;
496  case PAMGEN: readPamgenMeshFile(path,testData); break;
497  case CHACO: readZoltanTestData(path, testData, distributeInput); break;
498  default: readMatrixMarketFile(path, testData, distributeInput); break;
499  }
500 
501  }else if(pList.isParameter("x") || pList.isParameter("y") || pList.isParameter("z")){
502 
503  int x,y,z;
504  x = y = z = 0;
505  if(pList.isParameter("x")) x = pList.get<int>("x");
506  if(pList.isParameter("y")) y = pList.get<int>("y");
507  if(pList.isParameter("z")) z = pList.get<int>("z");
508 
509  string problemType = "";
510  if(pList.isParameter("equation type")) problemType = pList.get<string>("equation type");
511 
512  if (problemType.size() == 0){
513  int dim = 0;
514  if (x > 0) dim++;
515  if (y > 0) dim++;
516  if (z > 0) dim++;
517  if (dim == 1)
518  problemType = string("Laplace1D");
519  else if (dim == 2)
520  problemType = string("Laplace2D");
521  else if (dim == 3)
522  problemType = string("Laplace3D");
523  else
524  throw std::runtime_error("input");
525 
526  if (verbose_ && tcomm_->getRank() == 0)
527  std::cout << "UserInputForTests, Matrix type : " << problemType << std::endl;
528  }
529 
530 
531  buildCrsMatrix(x, y, z, problemType, distributeInput);
532 
533  }else{
534  std::cerr << "Input file block undefined!" << std::endl;
535  }
536 
537 #ifdef HAVE_EPETRA_DATA_TYPES
538  ecomm_ = Xpetra::toEpetra(c);
539 #endif
540 
541 }
542 
543 
544 RCP<Zoltan2_TestingFramework::tMVector_t> UserInputForTests::getUICoordinates()
545 {
546  if (xyz_.is_null())
547  throw std::runtime_error("could not read coord file");
548  return xyz_;
549 }
550 
551 RCP<Zoltan2_TestingFramework::tMVector_t> UserInputForTests::getUIWeights()
552 {
553  return vtxWeights_;
554 }
555 
556 RCP<Zoltan2_TestingFramework::tMVector_t> UserInputForTests::getUIEdgeWeights()
557 {
558  return edgWeights_;
559 }
560 
561 RCP<Zoltan2_TestingFramework::tcrsMatrix_t> UserInputForTests::getUITpetraCrsMatrix()
562 {
563  if (M_.is_null())
564  throw std::runtime_error("could not read mtx file");
565  return M_;
566 }
567 
568 RCP<Zoltan2_TestingFramework::tcrsGraph_t> UserInputForTests::getUITpetraCrsGraph()
569 {
570  if (M_.is_null())
571  throw std::runtime_error("could not read mtx file");
572  return rcp_const_cast<tcrsGraph_t>(M_->getCrsGraph());
573 }
574 
575 RCP<Zoltan2_TestingFramework::tVector_t> UserInputForTests::getUITpetraVector()
576 {
577  RCP<tVector_t> V = rcp(new tVector_t(M_->getRowMap(), 1));
578  V->randomize();
579 
580  return V;
581 }
582 
583 RCP<Zoltan2_TestingFramework::tMVector_t> UserInputForTests::getUITpetraMultiVector(int nvec)
584 {
585  RCP<tMVector_t> mV = rcp(new tMVector_t(M_->getRowMap(), nvec));
586  mV->randomize();
587 
588  return mV;
589 }
590 
591 RCP<Zoltan2_TestingFramework::xcrsMatrix_t> UserInputForTests::getUIXpetraCrsMatrix()
592 {
593  if (M_.is_null())
594  throw std::runtime_error("could not read mtx file");
595  return xM_;
596 }
597 
598 RCP<Zoltan2_TestingFramework::xcrsGraph_t> UserInputForTests::getUIXpetraCrsGraph()
599 {
600  if (M_.is_null())
601  throw std::runtime_error("could not read mtx file");
602  return rcp_const_cast<xcrsGraph_t>(xM_->getCrsGraph());
603 }
604 
605 RCP<Zoltan2_TestingFramework::xVector_t> UserInputForTests::getUIXpetraVector()
606 {
608 }
609 
610 RCP<Zoltan2_TestingFramework::xMVector_t> UserInputForTests::getUIXpetraMultiVector(int nvec)
611 {
612  RCP<tMVector_t> tMV = getUITpetraMultiVector(nvec);
614 }
615 
616 #ifdef HAVE_EPETRA_DATA_TYPES
617 RCP<Epetra_CrsGraph> UserInputForTests::getUIEpetraCrsGraph()
618 {
619  if (M_.is_null())
620  throw std::runtime_error("could not read mtx file");
621  RCP<const tcrsGraph_t> tgraph = M_->getCrsGraph();
622  RCP<const Tpetra::Map<zlno_t, zgno_t> > trowMap = tgraph->getRowMap();
623  RCP<const Tpetra::Map<zlno_t, zgno_t> > tcolMap = tgraph->getColMap();
624 
625  int nElts = static_cast<int>(trowMap->getGlobalNumElements());
626  int nMyElts = static_cast<int>(trowMap->getNodeNumElements());
627  int base = 0;
628  ArrayView<const int> gids = trowMap->getNodeElementList();
629 
630  Epetra_BlockMap erowMap(nElts, nMyElts,
631  gids.getRawPtr(), 1, base, *ecomm_);
632 
633  Array<int> rowSize(nMyElts);
634  for (int i=0; i < nMyElts; i++){
635  rowSize[i] = static_cast<int>(M_->getNumEntriesInLocalRow(i));
636  }
637 
638  size_t maxRow = M_->getNodeMaxNumRowEntries();
639  Array<int> colGids(maxRow);
640  ArrayView<const int> colLid;
641 
642  eG_ = rcp(new Epetra_CrsGraph(Copy, erowMap,
643  rowSize.getRawPtr(), true));
644 
645  for (int i=0; i < nMyElts; i++){
646  tgraph->getLocalRowView(i, colLid);
647  for (int j=0; j < colLid.size(); j++)
648  colGids[j] = tcolMap->getGlobalElement(colLid[j]);
649  eG_->InsertGlobalIndices(gids[i], rowSize[i], colGids.getRawPtr());
650  }
651  eG_->FillComplete();
652  return eG_;
653 }
654 
655 RCP<Epetra_CrsMatrix> UserInputForTests::getUIEpetraCrsMatrix()
656 {
657  if (M_.is_null())
658  throw std::runtime_error("could not read mtx file");
659  RCP<Epetra_CrsGraph> egraph = getUIEpetraCrsGraph();
660  eM_ = rcp(new Epetra_CrsMatrix(Copy, *egraph));
661 
662  size_t maxRow = M_->getNodeMaxNumRowEntries();
663  int nrows = egraph->NumMyRows();
664  const Epetra_BlockMap &rowMap = egraph->RowMap();
665  const Epetra_BlockMap &colMap = egraph->ColMap();
666  Array<int> colGid(maxRow);
667 
668  for (int i=0; i < nrows; i++){
669  ArrayView<const int> colLid;
670  ArrayView<const zscalar_t> nz;
671  M_->getLocalRowView(i, colLid, nz);
672  size_t rowSize = colLid.size();
673  int rowGid = rowMap.GID(i);
674  for (size_t j=0; j < rowSize; j++){
675  colGid[j] = colMap.GID(colLid[j]);
676  }
677  eM_->InsertGlobalValues(rowGid, (int)rowSize, nz.getRawPtr(), colGid.getRawPtr());
678  }
679  eM_->FillComplete();
680  return eM_;
681 }
682 
683 RCP<Epetra_Vector> UserInputForTests::getUIEpetraVector()
684 {
685  RCP<Epetra_CrsGraph> egraph = getUIEpetraCrsGraph();
686  RCP<Epetra_Vector> V = rcp(new Epetra_Vector(egraph->RowMap()));
687  V->Random();
688  return V;
689 }
690 
691 RCP<Epetra_MultiVector> UserInputForTests::getUIEpetraMultiVector(int nvec)
692 {
693  RCP<Epetra_CrsGraph> egraph = getUIEpetraCrsGraph();
694  RCP<Epetra_MultiVector> mV =
695  rcp(new Epetra_MultiVector(egraph->RowMap(), nvec));
696  mV->Random();
697  return mV;
698 }
699 #endif
700 
702 {
703  // find out if an input source has been loaded
704  return this->hasUICoordinates() || \
705  this->hasUITpetraCrsMatrix() || \
706  this->hasUITpetraCrsGraph() || \
707  this->hasPamgenMesh();
708 }
709 
710 bool UserInputForTests::hasInputDataType(const string &input_type)
711 {
712  if(input_type == "coordinates")
713  return this->hasUICoordinates();
714  else if(input_type == "tpetra_vector")
715  return this->hasUITpetraVector();
716  else if(input_type == "tpetra_multivector")
717  return this->hasUITpetraMultiVector();
718  else if(input_type == "tpetra_crs_graph")
719  return this->hasUITpetraCrsGraph();
720  else if(input_type == "tpetra_crs_matrix")
721  return this->hasUITpetraCrsMatrix();
722  else if(input_type == "xpetra_vector")
723  return this->hasUIXpetraVector();
724  else if(input_type == "xpetra_multivector")
725  return this->hasUIXpetraMultiVector();
726  else if(input_type == "xpetra_crs_graph")
727  return this->hasUIXpetraCrsGraph();
728  else if(input_type == "xpetra_crs_matrix")
729  return this->hasUIXpetraCrsMatrix();
730 #ifdef HAVE_EPETRA_DATA_TYPES
731  else if(input_type == "epetra_vector")
732  return this->hasUIEpetraVector();
733  else if(input_type == "epetra_multivector")
734  return this->hasUIEpetraMultiVector();
735  else if(input_type == "epetra_crs_graph")
736  return this->hasUIEpetraCrsGraph();
737  else if(input_type == "epetra_crs_matrix")
738  return this->hasUIEpetraCrsMatrix();
739 #endif
740 
741  return false;
742 }
743 
745 {
746  return xyz_.is_null() ? false : true;
747 }
748 
750 {
751  return vtxWeights_.is_null() ? false : true;
752 }
753 
755 {
756  return edgWeights_.is_null() ? false : true;
757 }
758 
760 {
761  return M_.is_null() ? false : true;
762 }
763 
765 {
766  return M_.is_null() ? false : true;
767 }
768 
770 {
771  return true;
772 }
773 
775 {
776  return true;
777 }
778 
780 {
781  return M_.is_null() ? false : true;
782 }
783 
785 {
786  return M_.is_null() ? false : true;
787 }
788 
790 {
791  return true;
792 }
793 
795 {
796  return true;
797 }
798 
800 {
801  return this->havePamgenMesh;
802 }
803 
804 #ifdef HAVE_EPETRA_DATA_TYPES
805 bool UserInputForTests::hasUIEpetraCrsGraph()
806 {
807  return M_.is_null() ? false : true;
808 }
809 
810 bool UserInputForTests::hasUIEpetraCrsMatrix()
811 {
812  return hasUIEpetraCrsGraph();
813 }
814 
815 bool UserInputForTests::hasUIEpetraVector()
816 {
817  return hasUIEpetraCrsGraph();
818 }
819 
820 bool UserInputForTests::hasUIEpetraMultiVector()
821 {
822  return hasUIEpetraCrsGraph();
823 }
824 #endif
825 
826 void UserInputForTests::getUIRandomData(unsigned int seed, zlno_t length,
827  zscalar_t min, zscalar_t max,
828  ArrayView<ArrayRCP<zscalar_t > > data)
829 {
830  if (length < 1)
831  return;
832 
833  size_t dim = data.size();
834  for (size_t i=0; i < dim; i++){
835  zscalar_t *tmp = new zscalar_t [length];
836  if (!tmp)
837  throw (std::bad_alloc());
838  data[i] = Teuchos::arcp(tmp, 0, length, true);
839  }
840 
841  zscalar_t scalingFactor = (max-min) / RAND_MAX;
842  srand(seed);
843  for (size_t i=0; i < dim; i++){
844  zscalar_t *x = data[i].getRawPtr();
845  for (zlno_t j=0; j < length; j++)
846  *x++ = min + (zscalar_t(rand()) * scalingFactor);
847  }
848 }
849 
850 // utility methods used when reading geom gen files
851 
852 string UserInputForTests::trim_right_copy(
853  const string& s,
854  const string& delimiters)
855 {
856  return s.substr( 0, s.find_last_not_of( delimiters ) + 1 );
857 }
858 
859 string UserInputForTests::trim_left_copy(
860  const string& s,
861  const string& delimiters)
862 {
863  return s.substr( s.find_first_not_of( delimiters ) );
864 }
865 
866 string UserInputForTests::trim_copy(
867  const string& s,
868  const string& delimiters)
869 {
870  return trim_left_copy( trim_right_copy( s, delimiters ), delimiters );
871 }
872 
873 void UserInputForTests::readGeometricGenTestData(string path,
874  string testData)
875 {
876 
877  std::ostringstream fname;
878  fname << path << "/" << testData << ".txt";
879 
880  if (verbose_ && tcomm_->getRank() == 0)
881  std::cout << "UserInputForTests, Read: " << fname.str() << std::endl;
882 
883  Teuchos::ParameterList geoparams("geo params");
884  readGeoGenParams(fname.str(),geoparams);
885 
886  geometricgen_t * gg = new geometricgen_t(geoparams, this->tcomm_);
887 
888  // get coordinate and point info
889  int coord_dim = gg->getCoordinateDimension();
890  int numWeightsPerCoord = gg->getNumWeights();
891  zlno_t numLocalPoints = gg->getNumLocalCoords();
892  zgno_t numGlobalPoints = gg->getNumGlobalCoords();
893 
894  // allocate an array of coordinate arrays
895  zscalar_t **coords = new zscalar_t * [coord_dim];
896  for(int i = 0; i < coord_dim; ++i){
897  coords[i] = new zscalar_t[numLocalPoints];
898  }
899 
900  // get a copy of the data
901  gg->getLocalCoordinatesCopy(coords);
902 
903  // get an array of arrays of weight data (if any)
904  zscalar_t **weight = NULL;
905  if (numWeightsPerCoord) {
906  // memory allocation
907  weight = new zscalar_t * [numWeightsPerCoord];
908  for(int i = 0; i < numWeightsPerCoord; ++i){
909  weight[i] = new zscalar_t[numLocalPoints];
910  }
911 
912  // get a copy of the weight data
913  gg->getLocalWeightsCopy(weight);
914  }
915 
916  delete gg; // free up memory from geom gen
917 
918 
919  // make a Tpetra map
920  RCP<Tpetra::Map<zlno_t, zgno_t, znode_t> > mp =
921  rcp(new Tpetra::Map<zlno_t, zgno_t, znode_t>(numGlobalPoints, numLocalPoints, 0, this->tcomm_));
922 
923  // make an array of array views containing the coordinate data
924  Teuchos::Array<Teuchos::ArrayView<const zscalar_t> > coordView(coord_dim);
925  for (int i=0; i < coord_dim; i++){
926  if(numLocalPoints > 0){
927  Teuchos::ArrayView<const zscalar_t> a(coords[i], numLocalPoints);
928  coordView[i] = a;
929  }
930  else {
931  Teuchos::ArrayView<const zscalar_t> a;
932  coordView[i] = a;
933  }
934  }
935 
936  // set the xyz_ multivector
937  xyz_ = RCP<tMVector_t>(new
938  tMVector_t(mp, coordView.view(0, coord_dim),
939  coord_dim));
940 
941  // set the vtx weights
942  if (numWeightsPerCoord) {
943  // make an array of array views containing the weight data
944  Teuchos::Array<Teuchos::ArrayView<const zscalar_t> > weightView(numWeightsPerCoord);
945  for (int i=0; i < numWeightsPerCoord; i++){
946  if(numLocalPoints > 0){
947  Teuchos::ArrayView<const zscalar_t> a(weight[i], numLocalPoints);
948  weightView[i] = a;
949  }
950  else {
951  Teuchos::ArrayView<const zscalar_t> a;
952  weightView[i] = a;
953  }
954  }
955 
956  vtxWeights_ = RCP<tMVector_t>(new tMVector_t(mp, weightView.view(0, numWeightsPerCoord),
957  numWeightsPerCoord));
958  }
959 }
960 
961 void UserInputForTests::readGeoGenParams(string paramFileName,
962  ParameterList &geoparams){
963 
964  const char param_comment = '#';
965 
966  std::string input = "";
967  char inp[25000];
968  for(int i = 0; i < 25000; ++i){
969  inp[i] = 0;
970  }
971 
972  bool fail = false;
973  if(this->tcomm_->getRank() == 0){
974 
975  std::fstream inParam(paramFileName.c_str());
976  if (inParam.fail())
977  {
978  fail = true;
979  }
980  if(!fail)
981  {
982  std::string tmp = "";
983  getline (inParam,tmp);
984  while (!inParam.eof()){
985  if(tmp != ""){
986  tmp = trim_copy(tmp);
987  if(tmp != ""){
988  input += tmp + "\n";
989  }
990  }
991  getline (inParam,tmp);
992  }
993  inParam.close();
994  for (size_t i = 0; i < input.size(); ++i){
995  inp[i] = input[i];
996  }
997  }
998  }
999 
1000 
1001 
1002  int size = (int)input.size();
1003  if(fail){
1004  size = -1;
1005  }
1006  this->tcomm_->broadcast(0, sizeof(int), (char*) &size);
1007  if(size == -1){
1008  throw "File " + paramFileName + " cannot be opened.";
1009  }
1010  this->tcomm_->broadcast(0, size, inp);
1011  std::istringstream inParam(inp);
1012  string str;
1013  getline (inParam,str);
1014  while (!inParam.eof()){
1015  if(str[0] != param_comment){
1016  size_t pos = str.find('=');
1017  if(pos == string::npos){
1018  throw "Invalid Line:" + str + " in parameter file";
1019  }
1020  string paramname = trim_copy(str.substr(0,pos));
1021  string paramvalue = trim_copy(str.substr(pos + 1));
1022  geoparams.set(paramname, paramvalue);
1023  }
1024  getline (inParam,str);
1025  }
1026 }
1027 
1029 RCP<tcrsMatrix_t> UserInputForTests::modifyMatrixGIDs(
1030  RCP<tcrsMatrix_t> &inMatrix
1031 )
1032 {
1033  // Produce a new matrix with the same structure as inMatrix,
1034  // but whose row/column GIDs are non-contiguous values.
1035  // In this case GID g in inMatrix becomes g*2+1 in outMatrix.
1036 
1037  // Create the map for the new matrix: same structure as inMap but with
1038  // the GIDs modified.
1039  RCP<const map_t> inMap = inMatrix->getRowMap();
1040 
1041  size_t nRows = inMap->getNodeNumElements();
1042  auto inRows = inMap->getMyGlobalIndices();
1043  Teuchos::Array<zgno_t> outRows(nRows);
1044  for (size_t i = 0; i < nRows; i++) {
1045  outRows[i] = newID(inRows[i]);
1046  }
1047 
1048  Tpetra::global_size_t nGlobalRows = inMap->getGlobalNumElements();
1049  RCP<map_t> outMap = rcp(new map_t(nGlobalRows, outRows(), 0,
1050  inMap->getComm()));
1051 
1052 #ifdef INCLUDE_LENGTHY_OUTPUT
1053  // Sanity check output
1054  {
1055  std::cout << inMap->getComm()->getRank() << " KDDKDD "
1056  << "nGlobal " << inMap->getGlobalNumElements() << " "
1057  << outMap->getGlobalNumElements() << "; "
1058  << "nLocal " << inMap->getNodeNumElements() << " "
1059  << outMap->getNodeNumElements() << "; "
1060  << std::endl;
1061  std::cout << inMap->getComm()->getRank() << " KDDKDD ";
1062  for (size_t i = 0; i < nRows; i++)
1063  std::cout << "(" << inMap->getMyGlobalIndices()[i] << ", "
1064  << outMap->getMyGlobalIndices()[i] << ") ";
1065  std::cout << std::endl;
1066  }
1067 #endif // INCLUDE_LENGTHY_OUTPUT
1068 
1069  // Create a new matrix using the new map
1070  // Get the length of the longest row; allocate memory.
1071  size_t rowLen = inMatrix->getNodeMaxNumRowEntries();
1072  RCP<tcrsMatrix_t> outMatrix = rcp(new tcrsMatrix_t(outMap, rowLen));
1073 
1074  Teuchos::Array<zgno_t> indices(rowLen);
1075  Teuchos::Array<zscalar_t> values(rowLen);
1076 
1077  for (size_t i = 0; i < nRows; i++) {
1078  size_t nEntries;
1079  zgno_t inGid = inMap->getGlobalElement(i);
1080  inMatrix->getGlobalRowCopy(inGid, indices, values, nEntries);
1081  for (size_t j = 0; j < nEntries; j++)
1082  indices[j] = newID(indices[j]);
1083 
1084  zgno_t outGid = outMap->getGlobalElement(i);
1085  outMatrix->insertGlobalValues(outGid, indices(0, nEntries),
1086  values(0, nEntries));
1087  }
1088  outMatrix->fillComplete();
1089 
1090 #ifdef INCLUDE_LENGTHY_OUTPUT
1091  // Sanity check output
1092  {
1093  std::cout << inMap->getComm()->getRank() << " KDDKDD Rows "
1094  << "nGlobal " << inMatrix->getGlobalNumRows() << " "
1095  << outMatrix->getGlobalNumRows() << "; "
1096  << "nLocal " << inMatrix->getNodeNumRows() << " "
1097  << outMatrix->getNodeNumRows() << std::endl;
1098  std::cout << inMap->getComm()->getRank() << " KDDKDD NNZS "
1099  << "nGlobal " << inMatrix->getGlobalNumEntries() << " "
1100  << outMatrix->getGlobalNumEntries() << "; "
1101  << "nLocal " << inMatrix->getNodeNumEntries() << " "
1102  << outMatrix->getNodeNumEntries() << std::endl;
1103 
1104  size_t nIn, nOut;
1105  Teuchos::Array<zgno_t> in(rowLen), out(rowLen);
1106  Teuchos::Array<zscalar_t> inval(rowLen), outval(rowLen);
1107 
1108  for (size_t i = 0; i < nRows; i++) {
1109  std::cout << inMap->getComm()->getRank() << " KDDKDD " << i << " nnz(";
1110  inMatrix->getGlobalRowCopy(inMap->getGlobalElement(i), in, inval, nIn);
1111  outMatrix->getGlobalRowCopy(outMap->getGlobalElement(i), out, outval,
1112  nOut);
1113 
1114  std::cout << nIn << ", " << nOut << "): ";
1115  for (size_t j = 0; j < nIn; j++) {
1116  std::cout << "(" << in[j] << " " << inval[j] << ", "
1117  << out[j] << " " << outval[j] << ") ";
1118  }
1119  std::cout << std::endl;
1120  }
1121  }
1122 #endif // INCLUDE_LENGTHY_OUTPUT
1123 
1124  return outMatrix;
1125 }
1126 
1128 
1129 void UserInputForTests::readMatrixMarketFile(
1130  string path,
1131  string testData,
1132  bool distributeInput
1133 )
1134 {
1135  std::ostringstream fname;
1136  fname << path << "/" << testData << ".mtx";
1137 
1138  if (verbose_ && tcomm_->getRank() == 0)
1139  std::cout << "UserInputForTests, Read: " << fname.str() << std::endl;
1140 
1141  // FIXME (mfh 01 Aug 2016) Tpetra::MatrixMarket::Reader has a graph
1142  // ("pattern" matrix) reader. Call its readSparseGraphFile method.
1143 
1144  RCP<tcrsMatrix_t> toMatrix;
1145  RCP<tcrsMatrix_t> fromMatrix;
1146  bool aok = true;
1147  try{
1148  typedef Tpetra::MatrixMarket::Reader<tcrsMatrix_t> reader_type;
1149  fromMatrix = reader_type::readSparseFile(fname.str(), tcomm_,
1150  true, false, false);
1151 #ifdef KDD_NOT_READY_YET
1152  // See note below about modifying coordinate IDs as well.
1153  //if (makeNonContiguous)
1154  fromMatrix = modifyMatrixGIDs(fromMatrix);
1155 #endif
1156 
1157  if(!distributeInput)
1158  {
1159  if (verbose_ && tcomm_->getRank() == 0)
1160  std::cout << "Constructing serial distribution of matrix" << std::endl;
1161  // need to make a serial map and then import the data to redistribute it
1162  RCP<const map_t> fromMap = fromMatrix->getRowMap();
1163 
1164  size_t numGlobalCoords = fromMap->getGlobalNumElements();
1165  size_t numLocalCoords = this->tcomm_->getRank() == 0 ? numGlobalCoords : 0;
1166  RCP<const map_t> toMap = rcp(new map_t(numGlobalCoords,numLocalCoords, 0, tcomm_));
1167 
1168  RCP<import_t> importer = rcp(new import_t(fromMap, toMap));
1169  toMatrix = rcp(new tcrsMatrix_t(toMap,0));
1170  toMatrix->doImport(*fromMatrix, *importer, Tpetra::INSERT);
1171  toMatrix->fillComplete();
1172 
1173  }else{
1174  toMatrix = fromMatrix;
1175  }
1176  }catch (std::exception &e) {
1177  if (tcomm_->getRank() == 0) {
1178  std::cout << "UserInputForTests unable to read matrix market file:"
1179  << fname.str() << std::endl;
1180  std::cout << e.what() << std::endl;
1181  }
1182  aok = false;
1183  }
1184  TEST_FAIL_AND_THROW(*tcomm_, aok,
1185  "UserInputForTests unable to read matrix market file");
1186 
1187  M_ = toMatrix;
1188 #ifdef INCLUDE_LENGTHY_OUTPUT
1189  std::cout << tcomm_->getRank() << " KDDKDD " << M_->getNodeNumRows()
1190  << " " << M_->getGlobalNumRows()
1191  << " " << M_->getNodeNumEntries()
1192  << " " << M_->getGlobalNumEntries() << std::endl;
1193 #endif // INCLUDE_LENGTHY_OUTPUT
1194 
1196 
1197  // Open the coordinate file.
1198 
1199  fname.str("");
1200  fname << path << "/" << testData << "_coord.mtx";
1201 
1202  size_t coordDim = 0, numGlobalCoords = 0;
1203  size_t msg[2]={0,0};
1204  ArrayRCP<ArrayRCP<zscalar_t> > xyz;
1205  std::ifstream coordFile;
1206 
1207  if (tcomm_->getRank() == 0){
1208 
1209  if (verbose_)
1210  std::cout << "UserInputForTests, Read: " <<
1211  fname.str() << std::endl;
1212 
1213  int fail = 0;
1214  try{
1215  coordFile.open(fname.str().c_str());
1216  }
1217  catch (std::exception &e){ // there is no coordinate file
1218  fail = 1;
1219  }
1220 
1221  if (!fail){
1222 
1223  // Read past banner to number and dimension of coordinates.
1224 
1225  char c[256];
1226  bool done=false;
1227 
1228  while (!done && !fail && coordFile.good()){
1229  coordFile.getline(c, 256);
1230  if (!c[0])
1231  fail = 1;
1232  else if (c[0] == '%')
1233  continue;
1234  else {
1235  done=true;
1236  std::istringstream s(c);
1237  s >> numGlobalCoords >> coordDim;
1238  if (!s.eof() || numGlobalCoords < 1 || coordDim < 1)
1239  fail=1;
1240  }
1241  }
1242 
1243  if (done){
1244 
1245  // Read in the coordinates.
1246 
1247  xyz = Teuchos::arcp(new ArrayRCP<zscalar_t> [coordDim], 0, coordDim);
1248 
1249  for (size_t dim=0; !fail && dim < coordDim; dim++){
1250  size_t idx;
1251  zscalar_t *tmp = new zscalar_t [numGlobalCoords];
1252  if (!tmp)
1253  fail = 1;
1254  else{
1255  xyz[dim] = Teuchos::arcp(tmp, 0, numGlobalCoords);
1256 
1257  for (idx=0; !coordFile.eof() && idx < numGlobalCoords; idx++){
1258  coordFile.getline(c, 256);
1259  std::istringstream s(c);
1260  s >> tmp[idx];
1261  }
1262 
1263  if (idx < numGlobalCoords)
1264  fail = 1;
1265  }
1266  }
1267 
1268  if (fail){
1269  ArrayRCP<zscalar_t> emptyArray;
1270  for (size_t dim=0; dim < coordDim; dim++)
1271  xyz[dim] = emptyArray; // free the memory
1272 
1273  coordDim = 0;
1274  }
1275  }
1276  else{
1277  fail = 1;
1278  }
1279 
1280  coordFile.close();
1281  }
1282 
1283  msg[0] = coordDim;
1284  msg[1] = numGlobalCoords;
1285  }
1286 
1287  // Broadcast coordinate dimension
1288  Teuchos::broadcast<int, size_t>(*tcomm_, 0, 2, msg);
1289 
1290  coordDim = msg[0];
1291  numGlobalCoords = msg[1];
1292 
1293  if (coordDim == 0)
1294  return;
1295 
1296  zgno_t base = 0;
1297  RCP<const map_t> toMap;
1298 
1299  if (!M_.is_null()){
1300  const RCP<const map_t> &mapM = M_->getRowMap();
1301  toMap = mapM;
1302  }
1303  else{
1304  if (verbose_ && tcomm_->getRank() == 0)
1305  {
1306  std::cout << "Matrix was null. ";
1307  std::cout << "Constructing distribution map for coordinate vector."
1308  << std::endl;
1309  }
1310 
1311  if(!distributeInput)
1312  {
1313  if (verbose_ && tcomm_->getRank() == 0)
1314  std::cout << "Constructing serial distribution map for coordinates."
1315  << std::endl;
1316 
1317  size_t numLocalCoords = this->tcomm_->getRank()==0 ? numGlobalCoords : 0;
1318  toMap = rcp(new map_t(numGlobalCoords,numLocalCoords, base, tcomm_));
1319  }else{
1320  toMap = rcp(new map_t(numGlobalCoords, base, tcomm_));
1321  }
1322  }
1323 
1324  // Export coordinates to their owners
1325 
1326  xyz_ = rcp(new tMVector_t(toMap, coordDim));
1327 
1328  ArrayRCP<ArrayView<const zscalar_t> > coordLists(coordDim);
1329 
1330  if (tcomm_->getRank() == 0){
1331 
1332  for (size_t dim=0; dim < coordDim; dim++)
1333  coordLists[dim] = xyz[dim].view(0, numGlobalCoords);
1334 
1335  zgno_t *tmp = new zgno_t [numGlobalCoords];
1336  if (!tmp)
1337  throw std::bad_alloc();
1338 
1339  ArrayRCP<const zgno_t> rowIds = Teuchos::arcp(tmp, 0, numGlobalCoords);
1340 
1341 #ifdef KDD_NOT_READY_YET
1342  // TODO if modifyMatrixGIDs, we need to modify ids here as well
1343  for (zgno_t id=0; id < zgno_t(numGlobalCoords); id++)
1344  *tmp++ = newID(id);
1345 #else
1346  for (zgno_t id=0; id < zgno_t(numGlobalCoords); id++)
1347  *tmp++ = id;
1348 #endif
1349 
1350  RCP<const map_t> fromMap = rcp(new map_t(numGlobalCoords,
1351  rowIds.view(0, numGlobalCoords),
1352  base, tcomm_));
1353 
1354  tMVector_t allCoords(fromMap, coordLists.view(0, coordDim), coordDim);
1355 
1356  export_t exporter(fromMap, toMap);
1357 
1358  xyz_->doExport(allCoords, exporter, Tpetra::INSERT);
1359  }
1360  else{
1361 
1362  RCP<const map_t> fromMap = rcp(new map_t(numGlobalCoords,
1363  ArrayView<zgno_t>(), base, tcomm_));
1364 
1365  tMVector_t allCoords(fromMap, coordLists.view(0, coordDim), coordDim);
1366 
1367  export_t exporter(fromMap, toMap);
1368 
1369  xyz_->doExport(allCoords, exporter, Tpetra::INSERT);
1370  }
1371 }
1372 
1373 void UserInputForTests::buildCrsMatrix(int xdim, int ydim, int zdim,
1374  string problemType, bool distributeInput)
1375 {
1376 #ifdef HAVE_ZOLTAN2_GALERI
1377  Teuchos::CommandLineProcessor tclp;
1378  Galeri::Xpetra::Parameters<zgno_t> params(tclp,
1379  xdim, ydim, zdim, problemType);
1380 
1381  RCP<const Tpetra::Map<zlno_t, zgno_t> > map;
1382  if (distributeInput)
1383  map = rcp(new Tpetra::Map<zlno_t, zgno_t>(params.GetNumGlobalElements(),
1384  0, tcomm_));
1385  else {
1386  // All data initially on rank 0
1387  size_t nGlobalElements = params.GetNumGlobalElements();
1388  size_t nLocalElements = ((tcomm_->getRank() == 0) ? nGlobalElements : 0);
1389  map = rcp(new Tpetra::Map<zlno_t, zgno_t>(nGlobalElements, nLocalElements, 0,
1390  tcomm_));
1391  }
1392 
1393  if (verbose_ && tcomm_->getRank() == 0){
1394 
1395  std::cout << "Matrix is " << (distributeInput ? "" : "not");
1396  std::cout << "distributed." << std::endl;
1397 
1398  std::cout << "UserInputForTests, Create matrix with " << problemType;
1399  std::cout << " (and " << xdim;
1400  if (zdim > 0)
1401  std::cout << " x " << ydim << " x " << zdim;
1402  else if (ydim > 0)
1403  std::cout << " x" << ydim << " x 1";
1404  else
1405  std::cout << "x 1 x 1";
1406 
1407  std::cout << " mesh)" << std::endl;
1408 
1409  }
1410 
1411  bool aok = true;
1412  try{
1413  RCP<Galeri::Xpetra::Problem<Tpetra::Map<zlno_t, zgno_t>, Tpetra::CrsMatrix<zscalar_t, zlno_t, zgno_t>, Tpetra::MultiVector<zscalar_t, zlno_t, zgno_t> > > Pr =
1414  Galeri::Xpetra::BuildProblem<zscalar_t, zlno_t, zgno_t, Tpetra::Map<zlno_t, zgno_t>, Tpetra::CrsMatrix<zscalar_t, zlno_t, zgno_t>, Tpetra::MultiVector<zscalar_t, zlno_t, zgno_t> >
1415  (params.GetMatrixType(), map, params.GetParameterList());
1416  M_ = Pr->BuildMatrix();
1417  }
1418  catch (std::exception &e) { // Probably not enough memory
1419  aok = false;
1420  }
1421  TEST_FAIL_AND_THROW(*tcomm_, aok,
1422  "UserInputForTests Galeri::Xpetra::BuildProblem failed");
1423 
1425 
1426  // Compute the coordinates for the matrix rows.
1427 
1428  if (verbose_ && tcomm_->getRank() == 0)
1429  std::cout <<
1430  "UserInputForTests, Implied matrix row coordinates computed" <<
1431  std::endl;
1432 
1433  ArrayView<const zgno_t> gids = map->getNodeElementList();
1434  zlno_t count = static_cast<zlno_t>(gids.size());
1435  int dim = 3;
1436  size_t pos = problemType.find("2D");
1437  if (pos != string::npos)
1438  dim = 2;
1439  else if (problemType == string("Laplace1D") ||
1440  problemType == string("Identity"))
1441  dim = 1;
1442 
1443  Array<ArrayRCP<zscalar_t> > coordinates(dim);
1444 
1445  if (count > 0){
1446  for (int i=0; i < dim; i++){
1447  zscalar_t *c = new zscalar_t [count];
1448  if (!c)
1449  throw(std::bad_alloc());
1450  coordinates[i] = Teuchos::arcp(c, 0, count, true);
1451  }
1452 
1453  if (dim==3){
1454  zscalar_t *x = coordinates[0].getRawPtr();
1455  zscalar_t *y = coordinates[1].getRawPtr();
1456  zscalar_t *z = coordinates[2].getRawPtr();
1457  zgno_t xySize = xdim * ydim;
1458  for (zlno_t i=0; i < count; i++){
1459  zgno_t iz = gids[i] / xySize;
1460  zgno_t xy = gids[i] - iz*xySize;
1461  z[i] = zscalar_t(iz);
1462  y[i] = zscalar_t(xy / xdim);
1463  x[i] = zscalar_t(xy % xdim);
1464  }
1465  }
1466  else if (dim==2){
1467  zscalar_t *x = coordinates[0].getRawPtr();
1468  zscalar_t *y = coordinates[1].getRawPtr();
1469  for (zlno_t i=0; i < count; i++){
1470  y[i] = zscalar_t(gids[i] / xdim);
1471  x[i] = zscalar_t(gids[i] % xdim);
1472  }
1473  }
1474  else{
1475  zscalar_t *x = coordinates[0].getRawPtr();
1476  for (zlno_t i=0; i < count; i++)
1477  x[i] = zscalar_t(gids[i]);
1478  }
1479  }
1480 
1481  Array<ArrayView<const zscalar_t> > coordView(dim);
1482  if (count > 0)
1483  for (int i=0; i < dim; i++)
1484  coordView[i] = coordinates[i].view(0,count);
1485 
1486  xyz_ = rcp(new tMVector_t(map, coordView.view(0, dim), dim));
1487 #else
1488  throw std::runtime_error("Galeri input requested but Trilinos is "
1489  "not built with Galeri.");
1490 #endif
1491 }
1492 
1493 void UserInputForTests::readZoltanTestData(string path, string testData,
1494  bool distributeInput)
1495 {
1496  int rank = tcomm_->getRank();
1497  FILE *graphFile = NULL;
1498  FILE *coordFile = NULL;
1499  FILE *assignFile = NULL;
1500  int fileInfo[3];
1501 
1502  for (int i = 0; i < CHACO_LINE_LENGTH; i++) chaco_line[i] = '\0';
1503 
1504  if (rank == 0){
1505  // set chacho graph file name
1506  std::ostringstream chGraphFileName;
1507  chGraphFileName << path << "/" << testData << ".graph";
1508 
1509  // set chaco graph
1510  std::ostringstream chCoordFileName;
1511  chCoordFileName << path << "/" << testData << ".coords";
1512 
1513  // set chaco graph
1514  std::ostringstream chAssignFileName;
1515  chAssignFileName << path << "/" << testData << ".assign";
1516 
1517  // open file
1518  graphFile = fopen(chGraphFileName.str().c_str(), "r");
1519 
1520  if(!graphFile) // maybe the user is using the default zoltan1 path convention
1521  {
1522  chGraphFileName.str("");
1523  chCoordFileName.str("");
1524  // try constructing zoltan1 paths
1525  chGraphFileName << path << "/ch_" << testData << "/" << testData << ".graph";
1526  chCoordFileName << path << "/ch_" << testData << "/" << testData << ".coords";
1527  chAssignFileName << path << "/ch_" << testData << "/" << testData << ".assign";
1528  // try to open the graph file again, if this doesn't open
1529  // the user has not provided a valid path to the file
1530  graphFile = fopen(chGraphFileName.str().c_str(), "r");
1531  }
1532 
1533  memset(fileInfo, 0, sizeof(int) * 3); // set fileinfo to 0's
1534  if (graphFile){
1535  fileInfo[0] = 1;
1536  if (verbose_ && tcomm_->getRank() == 0)
1537  std::cout << "UserInputForTests, open " <<
1538  chGraphFileName.str () << std::endl;
1539 
1540  coordFile = fopen(chCoordFileName.str().c_str(), "r");
1541  if (coordFile){
1542  fileInfo[1] = 1;
1543  if (verbose_ && tcomm_->getRank() == 0)
1544  std::cout << "UserInputForTests, open " <<
1545  chCoordFileName.str () << std::endl;
1546  }
1547 
1548  assignFile = fopen(chAssignFileName.str().c_str(), "r");
1549  if (assignFile){
1550  fileInfo[2] = 1;
1551  if (verbose_ && tcomm_->getRank() == 0)
1552  std::cout << "UserInputForTests, open " <<
1553  chAssignFileName.str () << std::endl;
1554  }
1555  }else{
1556  if (verbose_ && tcomm_->getRank() == 0){
1557  std::cout << "UserInputForTests, unable to open file: ";
1558  std::cout << chGraphFileName.str() << std::endl;
1559  }
1560  }
1561  }
1562 
1563  // broadcast whether we have graphs and coords to all processes
1564  Teuchos::broadcast<int, int>(*tcomm_, 0, 3, fileInfo);
1565 
1566  bool haveGraph = (fileInfo[0] == 1);
1567  bool haveCoords = (fileInfo[1] == 1);
1568  bool haveAssign = (fileInfo[2] == 1);
1569 
1570  if (haveGraph){
1571  // builds M_, vtxWeights_, and edgWeights_ and closes file.
1572  try{
1573  getUIChacoGraph(graphFile, haveAssign, assignFile,
1574  testData, distributeInput);
1575  }
1577 
1578  if (haveCoords){
1579  // builds xyz_ and closes the file.
1580  try{
1581  getUIChacoCoords(coordFile, testData);
1582  }
1584  }
1585  }
1586 
1588 }
1589 
1590 void UserInputForTests::getUIChacoGraph(FILE *fptr, bool haveAssign,
1591  FILE *assignFile, string fname,
1592  bool distributeInput)
1593 {
1594  int rank = tcomm_->getRank();
1595  int graphCounts[5];
1596  int nvtxs=0, nedges=0;
1597  int nVwgts=0, nEwgts=0;
1598  int *start = NULL, *adj = NULL;
1599  float *ewgts = NULL, *vwgts = NULL;
1600  size_t *nzPerRow = NULL;
1601  size_t maxRowLen = 0;
1602  zgno_t base = 0;
1603  ArrayRCP<const size_t> rowSizes;
1604  int fail = 0;
1605  bool haveEdges = true;
1606 
1607  if (rank == 0){
1608 
1609  memset(graphCounts, 0, 5*sizeof(int));
1610 
1611  // Reads in the file and closes it when done.
1612  fail = chaco_input_graph(fptr, fname.c_str(), &start, &adj,
1613  &nvtxs, &nVwgts, &vwgts, &nEwgts, &ewgts);
1614 
1615  // There are Zoltan2 test graphs that have no edges.
1616 
1617  // nEwgts must be 1 or 0 - add error
1618 
1619  if (start == NULL)
1620  haveEdges = false;
1621 
1622  if (verbose_)
1623  {
1624  std::cout << "UserInputForTests, " << nvtxs << " vertices,";
1625  if (haveEdges)
1626  std::cout << start[nvtxs] << " edges,";
1627  else
1628  std::cout << "no edges,";
1629  std::cout << nVwgts << " vertex weights, ";
1630  std::cout << nEwgts << " edge weights" << std::endl;
1631  }
1632 
1633  if (nvtxs==0)
1634  fail = true;
1635 
1636  if (fail){
1637  Teuchos::broadcast<int, int>(*tcomm_, 0, 5, graphCounts);
1638  throw std::runtime_error("Unable to read chaco file");
1639  }
1640 
1641  if (haveEdges)
1642  nedges = start[nvtxs];
1643 
1644  nzPerRow = new size_t [nvtxs];
1645  if (!nzPerRow)
1646  throw std::bad_alloc();
1647  rowSizes = arcp(nzPerRow, 0, nvtxs, true);
1648 
1649  if (haveEdges){
1650  for (int i=0; i < nvtxs; i++){
1651  nzPerRow[i] = start[i+1] - start[i];
1652  if (nzPerRow[i] > maxRowLen)
1653  maxRowLen = nzPerRow[i];
1654  }
1655  }
1656  else{
1657  memset(nzPerRow, 0, sizeof(size_t) * nvtxs);
1658  }
1659 
1660  // Make sure base gid is zero.
1661 
1662  if (nedges){
1663  int chbase = adj[0];
1664  for (int i=1; i < nedges; i++)
1665  if (adj[i] < chbase)
1666  chbase = adj[i];
1667 
1668  if (chbase > 0){
1669  for (int i=0; i < nedges; i++)
1670  adj[i] -= chbase;
1671  }
1672  }
1673 
1674  graphCounts[0] = nvtxs;
1675  graphCounts[1] = nedges;
1676  graphCounts[2] = nVwgts;
1677  graphCounts[3] = nEwgts;
1678  graphCounts[4] = (int)maxRowLen; // size_t maxRowLen will fit; it is <= (int-int)
1679  }
1680 
1681  Teuchos::broadcast<int, int>(*tcomm_, 0, 5, graphCounts);
1682 
1683  if (graphCounts[0] == 0)
1684  throw std::runtime_error("Unable to read chaco file");
1685 
1686  haveEdges = (graphCounts[1] > 0);
1687 
1688  RCP<tcrsMatrix_t> fromMatrix;
1689  RCP<const map_t> fromMap;
1690 
1691  // Create a Tpetra::CrsMatrix where rank 0 has entire matrix.
1692  if (rank == 0){
1693  fromMap = rcp(new map_t(nvtxs, nvtxs, base, tcomm_));
1694 
1695  fromMatrix =
1696  rcp(new tcrsMatrix_t(fromMap, rowSizes(), Tpetra::StaticProfile));
1697 
1698  if (haveEdges){
1699 
1700  zgno_t *edgeIds = new zgno_t [nedges];
1701  if (nedges && !edgeIds)
1702  throw std::bad_alloc();
1703  for (int i=0; i < nedges; i++)
1704  edgeIds[i] = adj[i];
1705 
1706  free(adj);
1707  adj = NULL;
1708 
1709  zgno_t *nextId = edgeIds;
1710  Array<zscalar_t> values(nedges, 1.0);
1711  if (nedges > 0 && nEwgts > 0) {
1712  for (int i=0; i < nedges; i++)
1713  values[i] = ewgts[i];
1714  free(ewgts);
1715  ewgts = NULL;
1716  }
1717 
1718  for (int i=0; i < nvtxs; i++){
1719  if (nzPerRow[i] > 0){
1720  ArrayView<const zgno_t> rowNz(nextId, nzPerRow[i]);
1721  fromMatrix->insertGlobalValues(i, rowNz, values.view(start[i], start[i+1] - start[i]));
1722  nextId += nzPerRow[i];
1723  }
1724  }
1725 
1726  delete [] edgeIds;
1727  edgeIds = NULL;
1728  }
1729 
1730  fromMatrix->fillComplete();
1731  }
1732  else{
1733  nvtxs = graphCounts[0];
1734  nedges = graphCounts[1];
1735  nVwgts = graphCounts[2];
1736  nEwgts = graphCounts[3];
1737  maxRowLen = graphCounts[4];
1738 
1739  // Create a Tpetra::CrsMatrix where rank 0 has entire matrix.
1740 
1741  fromMap = rcp(new map_t(nvtxs, 0, base, tcomm_));
1742 
1743  fromMatrix =
1744  rcp(new tcrsMatrix_t(fromMap, rowSizes(), Tpetra::StaticProfile));
1745 
1746  fromMatrix->fillComplete();
1747  }
1748 
1749 #ifdef KDDKDDPRINT
1750  if (rank == 0) {
1751  size_t sz = fromMatrix->getNodeMaxNumRowEntries();
1752  Teuchos::Array<zgno_t> indices(sz);
1753  Teuchos::Array<zscalar_t> values(sz);
1754  for (size_t i = 0; i < fromMatrix->getNodeNumRows(); i++) {
1755  zgno_t gid = fromMatrix->getRowMap()->getGlobalElement(i);
1756  size_t num;
1757  fromMatrix->getGlobalRowCopy(gid, indices(), values(), num);
1758  std::cout << "ROW " << gid << ": ";
1759  for (size_t j = 0; j < num; j++)
1760  std::cout << indices[j] << " ";
1761  std::cout << std::endl;
1762  }
1763  }
1764 #endif
1765 
1766  RCP<const map_t> toMap;
1767  RCP<tcrsMatrix_t> toMatrix;
1768  RCP<import_t> importer;
1769 
1770  if (distributeInput) {
1771  if (haveAssign) {
1772  // Read assignments from Chaco assignment file
1773  short *assignments = new short[nvtxs];
1774  if (rank == 0) {
1775  fail = chaco_input_assign(assignFile, fname.c_str(), nvtxs, assignments);
1776  }
1777  // Broadcast coordinate dimension
1778  Teuchos::broadcast<int, short>(*tcomm_, 0, nvtxs, assignments);
1779 
1780  // Build map with my vertices
1781  Teuchos::Array<zgno_t> mine;
1782  for (int i = 0; i < nvtxs; i++) {
1783  if (assignments[i] == rank)
1784  mine.push_back(i);
1785  }
1786 
1787  Tpetra::global_size_t dummy =
1788  Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
1789  toMap = rcp(new map_t(dummy, mine(), base, tcomm_));
1790  delete [] assignments;
1791  }
1792  else {
1793  // Create a Tpetra::Map with default row distribution.
1794  toMap = rcp(new map_t(nvtxs, base, tcomm_));
1795  }
1796  toMatrix = rcp(new tcrsMatrix_t(toMap, maxRowLen));
1797 
1798  // Import the data.
1799  importer = rcp(new import_t(fromMap, toMap));
1800  toMatrix->doImport(*fromMatrix, *importer, Tpetra::INSERT);
1801  toMatrix->fillComplete();
1802  }
1803  else {
1804  toMap = fromMap;
1805  toMatrix = fromMatrix;
1806  }
1807 
1808  M_ = toMatrix;
1809 
1810  // Vertex weights, if any
1811 
1812  typedef ArrayRCP<const ArrayView<const zscalar_t> > arrayArray_t;
1813 
1814  if (nVwgts > 0){
1815 
1816  ArrayRCP<zscalar_t> weightBuf;
1817  ArrayView<const zscalar_t> *wgts = new ArrayView<const zscalar_t> [nVwgts];
1818 
1819  if (rank == 0){
1820  size_t len = nVwgts * nvtxs;
1821  zscalar_t *buf = new zscalar_t [len];
1822  if (!buf) throw std::bad_alloc();
1823  weightBuf = arcp(buf, 0, len, true);
1824 
1825  for (int widx=0; widx < nVwgts; widx++){
1826  wgts[widx] = ArrayView<const zscalar_t>(buf, nvtxs);
1827  float *vw = vwgts + widx;
1828  for (int i=0; i < nvtxs; i++, vw += nVwgts)
1829  buf[i] = *vw;
1830  buf += nvtxs;
1831  }
1832 
1833  free(vwgts);
1834  vwgts = NULL;
1835  }
1836 
1837  arrayArray_t vweights = arcp(wgts, 0, nVwgts, true);
1838 
1839  RCP<tMVector_t> fromVertexWeights =
1840  rcp(new tMVector_t(fromMap, vweights.view(0, nVwgts), nVwgts));
1841 
1842  RCP<tMVector_t> toVertexWeights;
1843  if (distributeInput) {
1844  toVertexWeights = rcp(new tMVector_t(toMap, nVwgts));
1845  toVertexWeights->doImport(*fromVertexWeights, *importer, Tpetra::INSERT);
1846  }
1847  else
1848  toVertexWeights = fromVertexWeights;
1849 
1850  vtxWeights_ = toVertexWeights;
1851  }
1852 
1853  // Edge weights, if any
1854 
1855  if (haveEdges && nEwgts > 0){
1856 
1857  // No longer distributing edge weights; they are the matrix values
1858  /*
1859  ArrayRCP<zscalar_t> weightBuf;
1860  ArrayView<const zscalar_t> *wgts = new ArrayView<const zscalar_t> [nEwgts];
1861 
1862  toMap = rcp(new map_t(nedges, M_->getNodeNumEntries(), base, tcomm_));
1863 
1864  if (rank == 0){
1865  size_t len = nEwgts * nedges;
1866  zscalar_t *buf = new zscalar_t [len];
1867  if (!buf) throw std::bad_alloc();
1868  weightBuf = arcp(buf, 0, len, true);
1869 
1870  for (int widx=0; widx < nEwgts; widx++){
1871  wgts[widx] = ArrayView<const zscalar_t>(buf, nedges);
1872  float *ew = ewgts + widx;
1873  for (int i=0; i < nedges; i++, ew += nEwgts)
1874  buf[i] = *ew;
1875  buf += nedges;
1876  }
1877 
1878  free(ewgts);
1879  ewgts = NULL;
1880  fromMap = rcp(new map_t(nedges, nedges, base, tcomm_));
1881  }
1882  else{
1883  fromMap = rcp(new map_t(nedges, 0, base, tcomm_));
1884  }
1885 
1886  arrayArray_t eweights = arcp(wgts, 0, nEwgts, true);
1887 
1888  RCP<tMVector_t> fromEdgeWeights;
1889  RCP<tMVector_t> toEdgeWeights;
1890  RCP<import_t> edgeImporter;
1891 
1892  if (distributeInput) {
1893  fromEdgeWeights =
1894  rcp(new tMVector_t(fromMap, eweights.view(0, nEwgts), nEwgts));
1895  toEdgeWeights = rcp(new tMVector_t(toMap, nEwgts));
1896  edgeImporter = rcp(new import_t(fromMap, toMap));
1897  toEdgeWeights->doImport(*fromEdgeWeights, *edgeImporter, Tpetra::INSERT);
1898  }
1899  else
1900  toEdgeWeights = fromEdgeWeights;
1901 
1902  edgWeights_ = toEdgeWeights;
1903  */
1904 
1905  toMap = rcp(new map_t(nedges, M_->getNodeNumEntries(), base, tcomm_));
1906  edgWeights_ = rcp(new tMVector_t(toMap, nEwgts));
1907 
1908  size_t maxSize = M_->getNodeMaxNumRowEntries();
1909  Array<zlno_t> colind(maxSize);
1910  Array<zscalar_t> vals(maxSize);
1911  size_t nEntries;
1912 
1913  for (size_t i = 0, idx = 0; i < M_->getNodeNumRows(); i++) {
1914  M_->getLocalRowCopy(i, colind, vals, nEntries);
1915  for (size_t j = 0; j < nEntries; j++) {
1916  edgWeights_->replaceLocalValue(idx, 0, vals[j]); // Assuming nEwgts==1
1917  idx++;
1918  }
1919  }
1920  }
1921 
1922  if (start) {
1923  free(start);
1924  start = NULL;
1925  }
1926 }
1927 
1928 void UserInputForTests::getUIChacoCoords(FILE *fptr, string fname)
1929 {
1930  int rank = tcomm_->getRank();
1931  int ndim=0;
1932  double *x=NULL, *y=NULL, *z=NULL;
1933  int fail = 0;
1934 
1935  size_t globalNumVtx = M_->getGlobalNumRows();
1936 
1937  if (rank == 0){
1938 
1939  // This function is from the Zoltan C-library.
1940 
1941  // Reads in the file and closes it when done.
1942  fail = chaco_input_geom(fptr, fname.c_str(), (int)globalNumVtx,
1943  &ndim, &x, &y, &z);
1944 
1945  if (fail)
1946  ndim = 0;
1947 
1948  if (verbose_){
1949  std::cout << "UserInputForTests, read " << globalNumVtx;
1950  std::cout << " " << ndim << "-dimensional coordinates." << std::endl;
1951  }
1952  }
1953 
1954  Teuchos::broadcast<int, int>(*tcomm_, 0, 1, &ndim);
1955 
1956  if (ndim == 0)
1957  throw std::runtime_error("Can't read coordinate file");
1958 
1959  ArrayRCP<ArrayRCP<const zscalar_t> > coords(ndim);
1960  zlno_t len = 0;
1961 
1962  if (rank == 0){
1963 
1964  for (int dim=0; dim < ndim; dim++){
1965  zscalar_t *v = new zscalar_t [globalNumVtx];
1966  if (!v)
1967  throw std::bad_alloc();
1968  coords[dim] = arcp<const zscalar_t>(v, 0, globalNumVtx, true);
1969  double *val = (dim==0 ? x : (dim==1 ? y : z));
1970  for (size_t i=0; i < globalNumVtx; i++)
1971  v[i] = zscalar_t(val[i]);
1972 
1973  free(val);
1974  }
1975 
1976  len = static_cast<zlno_t>(globalNumVtx);
1977  }
1978 
1979  RCP<const map_t> fromMap = rcp(new map_t(globalNumVtx, len, 0, tcomm_));
1980  RCP<const map_t> toMap = M_->getRowMap();
1981  RCP<import_t> importer = rcp(new import_t(fromMap, toMap));
1982 
1983  Array<ArrayView<const zscalar_t> > coordData;
1984  for (int dim=0; dim < ndim; dim++)
1985  coordData.push_back(coords[dim].view(0, len));
1986 
1987  RCP<tMVector_t> fromCoords =
1988  rcp(new tMVector_t(fromMap, coordData.view(0, ndim), ndim));
1989 
1990  RCP<tMVector_t> toCoords = rcp(new tMVector_t(toMap, ndim));
1991 
1992  toCoords->doImport(*fromCoords, *importer, Tpetra::INSERT);
1993 
1994  xyz_ = toCoords;
1995 
1996 }
1997 
2000 
2001 double UserInputForTests::chaco_read_val(
2002  FILE* infile, /* file to read value from */
2003  int *end_flag /* 0 => OK, 1 => EOL, -1 => EOF */
2004 )
2005 {
2006  double val; /* return value */
2007  char *ptr; /* ptr to next string to read */
2008  char *ptr2; /* ptr to next string to read */
2009  int length; /* length of line to read */
2010  int length_left;/* length of line still around */
2011  int white_seen; /* have I detected white space yet? */
2012  int done; /* checking for end of scan */
2013  int i; /* loop counter */
2014 
2015  *end_flag = 0;
2016 
2017  if (chaco_offset == 0 || chaco_offset >= chaco_break_pnt) {
2018  if (chaco_offset >= chaco_break_pnt) { /* Copy rest of line back to beginning. */
2019  length_left = CHACO_LINE_LENGTH - chaco_save_pnt - 1;
2020  ptr2 = chaco_line;
2021  ptr = &chaco_line[chaco_save_pnt];
2022  for (i=length_left; i; i--) *ptr2++ = *ptr++;
2023  length = chaco_save_pnt + 1;
2024  }
2025  else {
2026  length = CHACO_LINE_LENGTH;
2027  length_left = 0;
2028  }
2029 
2030  /* Now read next line, or next segment of current one. */
2031  ptr2 = fgets(&chaco_line[length_left], length, infile);
2032 
2033  if (ptr2 == (char *) NULL) { /* We've hit end of file. */
2034  *end_flag = -1;
2035  return((double) 0.0);
2036  }
2037 
2038  if ((chaco_line[CHACO_LINE_LENGTH - 2] != '\n') && (chaco_line[CHACO_LINE_LENGTH - 2] != '\f')
2039  && (strlen(chaco_line) == CHACO_LINE_LENGTH - 1)){
2040  /* Line too long. Find last safe place in chaco_line. */
2041  chaco_break_pnt = CHACO_LINE_LENGTH - 1;
2042  chaco_save_pnt = chaco_break_pnt;
2043  white_seen = 0;
2044  done = 0;
2045  while (!done) {
2046  --chaco_break_pnt;
2047  if (chaco_line[chaco_break_pnt] != '\0') {
2048  if (isspace((int)(chaco_line[chaco_break_pnt]))) {
2049  if (!white_seen) {
2050  chaco_save_pnt = chaco_break_pnt + 1;
2051  white_seen = 1;
2052  }
2053  }
2054  else if (white_seen) {
2055  done= 1;
2056  }
2057  }
2058  }
2059  }
2060  else {
2061  chaco_break_pnt = CHACO_LINE_LENGTH;
2062  }
2063 
2064  chaco_offset = 0;
2065  }
2066 
2067  while (isspace((int)(chaco_line[chaco_offset])) && chaco_offset < CHACO_LINE_LENGTH) chaco_offset++;
2068  if (chaco_line[chaco_offset] == '%' || chaco_line[chaco_offset] == '#') {
2069  *end_flag = 1;
2070  if (chaco_break_pnt < CHACO_LINE_LENGTH) {
2071  chaco_flush_line(infile);
2072  }
2073  return((double) 0.0);
2074  }
2075 
2076  ptr = &(chaco_line[chaco_offset]);
2077  val = strtod(ptr, &ptr2);
2078 
2079  if (ptr2 == ptr) { /* End of input line. */
2080  chaco_offset = 0;
2081  *end_flag = 1;
2082  return((double) 0.0);
2083  }
2084  else {
2085  chaco_offset = (int) (ptr2 - chaco_line) / sizeof(char);
2086  }
2087 
2088  return(val);
2089 }
2090 
2091 
2092 int UserInputForTests::chaco_read_int(
2093  FILE *infile, /* file to read value from */
2094  int *end_flag /* 0 => OK, 1 => EOL, -1 => EOF */
2095 )
2096 {
2097  int val; /* return value */
2098  char *ptr; /* ptr to next string to read */
2099  char *ptr2; /* ptr to next string to read */
2100  int length; /* length of line to read */
2101  int length_left; /* length of line still around */
2102  int white_seen; /* have I detected white space yet? */
2103  int done; /* checking for end of scan */
2104  int i; /* loop counter */
2105 
2106  *end_flag = 0;
2107 
2108  if (chaco_offset == 0 || chaco_offset >= chaco_break_pnt) {
2109  if (chaco_offset >= chaco_break_pnt) { /* Copy rest of line back to beginning. */
2110  length_left = CHACO_LINE_LENGTH - chaco_save_pnt - 1;
2111  ptr2 = chaco_line;
2112  ptr = &chaco_line[chaco_save_pnt];
2113  for (i=length_left; i; i--) *ptr2++ = *ptr++;
2114  length = chaco_save_pnt + 1;
2115  }
2116  else {
2117  length = CHACO_LINE_LENGTH;
2118  length_left = 0;
2119  }
2120 
2121  /* Now read next line, or next segment of current one. */
2122  ptr2 = fgets(&chaco_line[length_left], length, infile);
2123 
2124  if (ptr2 == (char *) NULL) { /* We've hit end of file. */
2125  *end_flag = -1;
2126  return(0);
2127  }
2128 
2129  if ((chaco_line[CHACO_LINE_LENGTH - 2] != '\n') && (chaco_line[CHACO_LINE_LENGTH - 2] != '\f')
2130  && (strlen(chaco_line) == CHACO_LINE_LENGTH - 1)){
2131  /* Line too long. Find last safe place in line. */
2132  chaco_break_pnt = CHACO_LINE_LENGTH - 1;
2133  chaco_save_pnt = chaco_break_pnt;
2134  white_seen = 0;
2135  done = 0;
2136  while (!done) {
2137  --chaco_break_pnt;
2138  if (chaco_line[chaco_break_pnt] != '\0') {
2139  if (isspace((int)(chaco_line[chaco_break_pnt]))) {
2140  if (!white_seen) {
2141  chaco_save_pnt = chaco_break_pnt + 1;
2142  white_seen = 1;
2143  }
2144  }
2145  else if (white_seen) {
2146  done= 1;
2147  }
2148  }
2149  }
2150  }
2151  else {
2152  chaco_break_pnt = CHACO_LINE_LENGTH;
2153  }
2154 
2155  chaco_offset = 0;
2156  }
2157 
2158  while (isspace((int)(chaco_line[chaco_offset])) && chaco_offset < CHACO_LINE_LENGTH) chaco_offset++;
2159  if (chaco_line[chaco_offset] == '%' || chaco_line[chaco_offset] == '#') {
2160  *end_flag = 1;
2161  if (chaco_break_pnt < CHACO_LINE_LENGTH) {
2162  chaco_flush_line(infile);
2163  }
2164  return(0);
2165  }
2166 
2167  ptr = &(chaco_line[chaco_offset]);
2168  val = (int) strtol(ptr, &ptr2, 10);
2169 
2170  if (ptr2 == ptr) { /* End of input chaco_line. */
2171  chaco_offset = 0;
2172  *end_flag = 1;
2173  return(0);
2174  }
2175  else {
2176  chaco_offset = (int) (ptr2 - chaco_line) / sizeof(char);
2177  }
2178 
2179  return(val);
2180 }
2181 
2182 void UserInputForTests::chaco_flush_line(
2183  FILE *infile /* file to read value from */
2184 )
2185 {
2186  char c; /* character being read */
2187 
2188  c = fgetc(infile);
2189  while (c != '\n' && c != '\f')
2190  c = fgetc(infile);
2191 }
2192 
2193 int UserInputForTests::chaco_input_graph(
2194  FILE *fin, /* input file */
2195  const char *inname, /* name of input file */
2196  int **start, /* start of edge list for each vertex */
2197  int **adjacency, /* edge list data */
2198  int *nvtxs, /* number of vertices in graph */
2199  int *nVwgts, /* # of vertex weights per node */
2200  float **vweights, /* vertex weight list data */
2201  int *nEwgts, /* # of edge weights per edge */
2202  float **eweights /* edge weight list data */
2203 )
2204 {
2205  int *adjptr; /* loops through adjacency data */
2206  float *ewptr; /* loops through edge weight data */
2207  int narcs; /* number of edges expected in graph */
2208  int nedges; /* twice number of edges really in graph */
2209  int nedge; /* loops through edges for each vertex */
2210  int flag; /* condition indicator */
2211  int skip_flag; /* should this edge be ignored? */
2212  int end_flag; /* indicates end of line or file */
2213  int vtx; /* vertex in graph */
2214  int line_num; /* line number in input file */
2215  int sum_edges; /* total number of edges read so far */
2216  int option = 0; /* input option */
2217  int using_ewgts; /* are edge weights in input file? */
2218  int using_vwgts; /* are vertex weights in input file? */
2219  int vtxnums; /* are vertex numbers in input file? */
2220  int vertex; /* current vertex being read */
2221  int new_vertex; /* new vertex being read */
2222  float weight; /* weight being read */
2223  float eweight; /* edge weight being read */
2224  int neighbor; /* neighbor of current vertex */
2225  int error_flag; /* error reading input? */
2226  int j; /* loop counters */
2227 
2228  /* Read first line of input (= nvtxs, narcs, option). */
2229  /* The (decimal) digits of the option variable mean: 1's digit not zero => input
2230  edge weights 10's digit not zero => input vertex weights 100's digit not zero
2231  => include vertex numbers */
2232 
2233  *start = NULL;
2234  *adjacency = NULL;
2235  *vweights = NULL;
2236  *eweights = NULL;
2237 
2238  error_flag = 0;
2239  line_num = 0;
2240 
2241  /* Read any leading comment lines */
2242  end_flag = 1;
2243  while (end_flag == 1) {
2244  *nvtxs = chaco_read_int(fin, &end_flag);
2245  ++line_num;
2246  }
2247  if (*nvtxs <= 0) {
2248  printf("ERROR in graph file `%s':", inname);
2249  printf(" Invalid number of vertices (%d).\n", *nvtxs);
2250  fclose(fin);
2251  return(1);
2252  }
2253 
2254  narcs = chaco_read_int(fin, &end_flag);
2255  if (narcs < 0) {
2256  printf("ERROR in graph file `%s':", inname);
2257  printf(" Invalid number of expected edges (%d).\n", narcs);
2258  fclose(fin);
2259  return(1);
2260  }
2261 
2262  /* Check if vertex or edge weights are used */
2263  if (!end_flag) {
2264  option = chaco_read_int(fin, &end_flag);
2265  }
2266  using_ewgts = option - 10 * (option / 10);
2267  option /= 10;
2268  using_vwgts = option - 10 * (option / 10);
2269  option /= 10;
2270  vtxnums = option - 10 * (option / 10);
2271 
2272  /* Get weight info from Chaco option */
2273  (*nVwgts) = using_vwgts;
2274  (*nEwgts) = using_ewgts;
2275 
2276  /* Read numbers of weights if they are specified separately */
2277  if (!end_flag && using_vwgts==1){
2278  j = chaco_read_int(fin, &end_flag);
2279  if (!end_flag) (*nVwgts) = j;
2280  }
2281  if (!end_flag && using_ewgts==1){
2282  j = chaco_read_int(fin, &end_flag);
2283  if (!end_flag) (*nEwgts) = j;
2284  }
2285 
2286  /* Discard rest of line */
2287  while (!end_flag)
2288  j = chaco_read_int(fin, &end_flag);
2289 
2290  /* Allocate space for rows and columns. */
2291  *start = (int *) malloc((unsigned) (*nvtxs + 1) * sizeof(int));
2292  if (narcs != 0)
2293  *adjacency = (int *) malloc((unsigned) (2 * narcs + 1) * sizeof(int));
2294  else
2295  *adjacency = NULL;
2296 
2297  if (using_vwgts)
2298  *vweights = (float *) malloc((unsigned) (*nvtxs) * (*nVwgts) * sizeof(float));
2299  else
2300  *vweights = NULL;
2301 
2302  if (using_ewgts)
2303  *eweights = (float *)
2304  malloc((unsigned) (2 * narcs + 1) * (*nEwgts) * sizeof(float));
2305  else
2306  *eweights = NULL;
2307 
2308  adjptr = *adjacency;
2309  ewptr = *eweights;
2310 
2311  sum_edges = 0;
2312  nedges = 0;
2313  (*start)[0] = 0;
2314  vertex = 0;
2315  vtx = 0;
2316  new_vertex = 1;
2317  while ((using_vwgts || vtxnums || narcs) && end_flag != -1) {
2318  ++line_num;
2319 
2320  /* If multiple input lines per vertex, read vertex number. */
2321  if (vtxnums) {
2322  j = chaco_read_int(fin, &end_flag);
2323  if (end_flag) {
2324  if (vertex == *nvtxs)
2325  break;
2326  printf("ERROR in graph file `%s':", inname);
2327  printf(" no vertex number in line %d.\n", line_num);
2328  fclose(fin);
2329  return (1);
2330  }
2331  if (j != vertex && j != vertex + 1) {
2332  printf("ERROR in graph file `%s':", inname);
2333  printf(" out-of-order vertex number in line %d.\n", line_num);
2334  fclose(fin);
2335  return (1);
2336  }
2337  if (j != vertex) {
2338  new_vertex = 1;
2339  vertex = j;
2340  }
2341  else
2342  new_vertex = 0;
2343  }
2344  else
2345  vertex = ++vtx;
2346 
2347  if (vertex > *nvtxs)
2348  break;
2349 
2350  /* If vertices are weighted, read vertex weight. */
2351  if (using_vwgts && new_vertex) {
2352  for (j=0; j<(*nVwgts); j++){
2353  weight = chaco_read_val(fin, &end_flag);
2354  if (end_flag) {
2355  printf("ERROR in graph file `%s':", inname);
2356  printf(" not enough weights for vertex %d.\n", vertex);
2357  fclose(fin);
2358  return (1);
2359  }
2360  (*vweights)[(vertex-1)*(*nVwgts)+j] = weight;
2361  }
2362  }
2363 
2364  nedge = 0;
2365 
2366  /* Read number of adjacent vertex. */
2367  neighbor = chaco_read_int(fin, &end_flag);
2368 
2369  while (!end_flag) {
2370  skip_flag = 0;
2371 
2372  if (using_ewgts) { /* Read edge weight if it's being input. */
2373  for (j=0; j<(*nEwgts); j++){
2374  eweight = chaco_read_val(fin, &end_flag);
2375 
2376  if (end_flag) {
2377  printf("ERROR in graph file `%s':", inname);
2378  printf(" not enough weights for edge (%d,%d).\n", vertex, neighbor);
2379  fclose(fin);
2380  return (1);
2381  }
2382 
2383  else {
2384  *ewptr++ = eweight;
2385  }
2386  }
2387  }
2388 
2389  /* Add edge to data structure. */
2390  if (!skip_flag) {
2391  if (++nedges > 2*narcs) {
2392  printf("ERROR in graph file `%s':", inname);
2393  printf(" at least %d adjacencies entered, but nedges = %d\n",
2394  nedges, narcs);
2395  fclose(fin);
2396  return (1);
2397  }
2398  *adjptr++ = neighbor;
2399  nedge++;
2400  }
2401 
2402  /* Read number of next adjacent vertex. */
2403  neighbor = chaco_read_int(fin, &end_flag);
2404  }
2405 
2406  sum_edges += nedge;
2407  (*start)[vertex] = sum_edges;
2408  }
2409 
2410  /* Make sure there's nothing else in file. */
2411  flag = 0;
2412  while (!flag && end_flag != -1) {
2413  chaco_read_int(fin, &end_flag);
2414  if (!end_flag)
2415  flag = 1;
2416  }
2417 
2418  (*start)[*nvtxs] = sum_edges;
2419 
2420  if (vertex != 0) { /* Normal file was read. */
2421  if (narcs) {
2422  }
2423  else { /* no edges, but did have vertex weights or vertex numbers */
2424  free(*start);
2425  *start = NULL;
2426  if (*adjacency != NULL)
2427  free(*adjacency);
2428  *adjacency = NULL;
2429  if (*eweights != NULL)
2430  free(*eweights);
2431  *eweights = NULL;
2432  }
2433  }
2434 
2435  else {
2436  /* Graph was empty */
2437  free(*start);
2438  if (*adjacency != NULL)
2439  free(*adjacency);
2440  if (*vweights != NULL)
2441  free(*vweights);
2442  if (*eweights != NULL)
2443  free(*eweights);
2444  *start = NULL;
2445  *adjacency = NULL;
2446  }
2447 
2448  fclose(fin);
2449 
2450  return (error_flag);
2451 }
2452 
2453 
2454 int UserInputForTests::chaco_input_geom(
2455  FILE *fingeom, /* geometry input file */
2456  const char *geomname, /* name of geometry file */
2457  int nvtxs, /* number of coordinates to read */
2458  int *igeom, /* dimensionality of geometry */
2459  double **x, /* coordinates of vertices */
2460  double **y,
2461  double **z
2462 )
2463 {
2464  double xc, yc, zc =0; /* first x, y, z coordinate */
2465  int nread; /* number of lines of coordinates read */
2466  int flag; /* any bad data at end of file? */
2467  int line_num; /* counts input lines in file */
2468  int end_flag; /* return conditional */
2469  int ndims; /* number of values in an input line */
2470  int i=0; /* loop counter */
2471 
2472  *x = *y = *z = NULL;
2473  line_num = 0;
2474  end_flag = 1;
2475  while (end_flag == 1) {
2476  xc = chaco_read_val(fingeom, &end_flag);
2477  ++line_num;
2478  }
2479 
2480  if (end_flag == -1) {
2481  printf("No values found in geometry file `%s'\n", geomname);
2482  fclose(fingeom);
2483  return (1);
2484  }
2485 
2486  ndims = 1;
2487  yc = chaco_read_val(fingeom, &end_flag);
2488  if (end_flag == 0) {
2489  ndims = 2;
2490  zc = chaco_read_val(fingeom, &end_flag);
2491  if (end_flag == 0) {
2492  ndims = 3;
2493  chaco_read_val(fingeom, &end_flag);
2494  if (!end_flag) {
2495  printf("Too many values on input line of geometry file `%s'\n",
2496  geomname);
2497 
2498  printf(" Maximum dimensionality is 3\n");
2499  fclose(fingeom);
2500  return (1);
2501  }
2502  }
2503  }
2504 
2505  *igeom = ndims;
2506 
2507  *x = (double *) malloc((unsigned) nvtxs * sizeof(double));
2508  (*x)[0] = xc;
2509  if (ndims > 1) {
2510  *y = (double *) malloc((unsigned) nvtxs * sizeof(double));
2511  (*y)[0] = yc;
2512  }
2513  if (ndims > 2) {
2514  *z = (double *) malloc((unsigned) nvtxs * sizeof(double));
2515  (*z)[0] = zc;
2516  }
2517 
2518  for (nread = 1; nread < nvtxs; nread++) {
2519  ++line_num;
2520  if (ndims == 1) {
2521  i = fscanf(fingeom, "%lf", &((*x)[nread]));
2522  }
2523  else if (ndims == 2) {
2524  i = fscanf(fingeom, "%lf%lf", &((*x)[nread]), &((*y)[nread]));
2525  }
2526  else if (ndims == 3) {
2527  i = fscanf(fingeom, "%lf%lf%lf", &((*x)[nread]), &((*y)[nread]),
2528  &((*z)[nread]));
2529  }
2530 
2531  if (i == EOF) {
2532  printf("Too few lines of values in geometry file; nvtxs=%d, but only %d read\n",
2533  nvtxs, nread);
2534  fclose(fingeom);
2535  return (1);
2536  }
2537  else if (i != ndims) {
2538  printf("Wrong number of values in line %d of geometry file `%s'\n",
2539  line_num, geomname);
2540  fclose(fingeom);
2541  return (1);
2542  }
2543  }
2544 
2545  /* Check for spurious extra stuff in file. */
2546  flag = 0;
2547  end_flag = 0;
2548  while (!flag && end_flag != -1) {
2549  chaco_read_val(fingeom, &end_flag);
2550  if (!end_flag)
2551  flag = 1;
2552  }
2553 
2554  fclose(fingeom);
2555 
2556  return (0);
2557 }
2558 
2559 // Chaco input assignments from filename.assign; copied from Zoltan
2560 
2561 int UserInputForTests::chaco_input_assign(
2562  FILE *finassign, /* input assignment file */
2563  const char *inassignname, /* name of input assignment file */
2564  int nvtxs, /* number of vertices to output */
2565  short *assignment) /* values to be printed */
2566 {
2567  int flag; /* logical conditional */
2568  int end_flag; /* return flag from read_int() */
2569  int i, j; /* loop counter */
2570 
2571  /* Get the assignment vector one line at a time, checking as you go. */
2572  /* First read past any comments at top. */
2573  end_flag = 1;
2574  while (end_flag == 1) {
2575  assignment[0] = chaco_read_int(finassign, &end_flag);
2576  }
2577 
2578  if (assignment[0] < 0) {
2579  printf("ERROR: Entry %d in assignment file `%s' less than zero (%d)\n",
2580  1, inassignname, assignment[0]);
2581  fclose(finassign);
2582  return (1);
2583  }
2584 
2585  if (end_flag == -1) {
2586  printf("ERROR: No values found in assignment file `%s'\n", inassignname);
2587  fclose(finassign);
2588  return (1);
2589  }
2590 
2591  flag = 0;
2592  int np = this->tcomm_->getSize();
2593  if (assignment[0] >= np) flag = assignment[0];
2594  srand(this->tcomm_->getRank());
2595  for (i = 1; i < nvtxs; i++) {
2596  j = fscanf(finassign, "%hd", &(assignment[i]));
2597  if (j != 1) {
2598  printf("ERROR: Too few values in assignment file `%s'.\n", inassignname);
2599  fclose(finassign);
2600  return (1);
2601  }
2602  if (assignment[i] < 0) {
2603  printf("ERROR: Entry %d in assignment file `%s' less than zero (%d)\n",
2604  i+1, inassignname, assignment[i]);
2605  fclose(finassign);
2606  return (1);
2607  }
2608  if (assignment[i] >= np) { // warn since perhaps an error -- initial part
2609  // assignment is greater than number of processors
2610  if (assignment[i] > flag)
2611  flag = assignment[i];
2612  assignment[i] = rand() % np; // randomly assign vtx to a proc in this case
2613  }
2614  }
2615  srand(rand());
2616 
2617  if (flag) {
2618  printf("WARNING: Possible error in assignment file `%s'\n",
2619  inassignname);
2620  printf(" Max assignment set (%d) greater than "
2621  "max processor rank (%d)\n", flag, np-1);
2622  printf(" Some vertices given random initial assignments\n");
2623  }
2624 
2625  /* Check for spurious extra stuff in file. */
2626  flag = 0;
2627  end_flag = 0;
2628  while (!flag && end_flag != -1) {
2629  chaco_read_int(finassign, &end_flag);
2630  if (!end_flag)
2631  flag = 1;
2632  }
2633  if (flag) {
2634  printf("WARNING: Possible error in assignment file `%s'\n", inassignname);
2635  printf(" Numerical data found after expected end of file\n");
2636  }
2637 
2638  fclose(finassign);
2639  return (0);
2640 }
2641 
2642 // Pamgen Reader
2643 void UserInputForTests::readPamgenMeshFile(string path, string testData)
2644 {
2645 #ifdef HAVE_ZOLTAN2_PAMGEN
2646  int rank = this->tcomm_->getRank();
2647  if (verbose_ && tcomm_->getRank() == 0)
2648  std::cout << "UserInputForTestsBD::readPamgenFile, Read: " << testData << std::endl;
2649 
2650  size_t len;
2651  std::fstream file;
2652  int dimension;
2653  if (rank == 0){
2654  // set file name
2655  std::ostringstream meshFileName;
2656  meshFileName << path << "/" << testData << ".pmgen";
2657  // open file
2658 
2659  file.open(meshFileName.str(), std::ios::in);
2660 
2661  if(!file.is_open()) // may be a problem with path or filename
2662  {
2663  if(verbose_ && tcomm_->getRank() == 0)
2664  {
2665  std::cout << "Unable to open pamgen mesh: ";
2666  std::cout << meshFileName.str();
2667  std::cout <<"\nPlease check file path and name." << std::endl;
2668  }
2669  len = 0; // broadcaset 0 length ->will cause exit
2670  }else{
2671  // write to character array
2672  // get size of file
2673  file.seekg (0,file.end);
2674  len = file.tellg();
2675  file.seekg (0);
2676 
2677  // get dimension
2678  dimension = 2;
2679  std::string line;
2680  while(std::getline(file,line))
2681  {
2682  if( line.find("nz") != std::string::npos ||
2683  line.find("nphi") != std::string::npos)
2684  {
2685  dimension = 3;
2686  break;
2687  }
2688  }
2689 
2690  file.clear();
2691  file.seekg(0, std::ios::beg);
2692  }
2693  }
2694 
2695  // broadcast the file size
2696  this->tcomm_->broadcast(0,sizeof(int), (char *)&dimension);
2697  this->tcomm_->broadcast(0,sizeof(size_t),(char *)&len);
2698  this->tcomm_->barrier();
2699 
2700  if(len == 0){
2701  if(verbose_ && tcomm_->getRank() == 0)
2702  std::cout << "Pamgen Mesh file size == 0, exiting UserInputForTests early." << std::endl;
2703  return;
2704  }
2705 
2706  char * file_data = new char[len+1];
2707  file_data[len] = '\0'; // critical to null terminate buffer
2708  if(rank == 0){
2709  file.read(file_data,len); // if proc 0 then read file
2710  }
2711 
2712  // broadcast the file to the world
2713  this->tcomm_->broadcast(0,(int)len+1,file_data);
2714  this->tcomm_->barrier();
2715 
2716  // Create the PamgenMesh
2717 
2718  this->pamgen_mesh = rcp(new PamgenMesh);
2719  this->havePamgenMesh = true;
2720  pamgen_mesh->createMesh(file_data,dimension,this->tcomm_);
2721 
2722  // save mesh info
2723  pamgen_mesh->storeMesh();
2724  this->tcomm_->barrier();
2725 
2726  // set coordinates
2727  this->setPamgenCoordinateMV();
2728 
2729  // set adjacency graph
2730  this->setPamgenAdjacencyGraph();
2731 
2732  this->tcomm_->barrier();
2733  if(rank == 0) file.close();
2734  delete [] file_data;
2735 #else
2736  throw std::runtime_error("Pamgen requested but Trilinos "
2737  "not built with Pamgen");
2738 #endif
2739 }
2740 
2741 #ifdef HAVE_ZOLTAN2_PAMGEN
2742 void UserInputForTests::setPamgenCoordinateMV()
2743 {
2744  int dimension = pamgen_mesh->num_dim;
2745  // get coordinate and point info;
2746 // zlno_t numLocalPoints = pamgen_mesh->num_nodes;
2747 // zgno_t numGlobalPoints = pamgen_mesh->num_nodes_global;
2748  zgno_t numelements = pamgen_mesh->num_elem;
2749  zgno_t numGlobalElements = pamgen_mesh->num_elems_global;
2750  // allocate and set an array of coordinate arrays
2751  zscalar_t **elem_coords = new zscalar_t * [dimension];
2752  for(int i = 0; i < dimension; ++i){
2753  elem_coords[i] = new zscalar_t[numelements];
2754  double *tmp = &pamgen_mesh->element_coord[i*numelements];
2755  for (int j = 0; j < numelements; j++) elem_coords[i][j] = tmp[j];
2756  }
2757 
2758  // make a Tpetra map
2759  typedef Tpetra::Map<zlno_t, zgno_t, znode_t> map_t;
2760  RCP<Tpetra::Map<zlno_t, zgno_t, znode_t> > mp;
2761  // mp = rcp(new map_t(numGlobalElements, numelements, 0, this->tcomm_)); // constructo 1
2762 
2763 // Array<zgno_t>::size_type numEltsPerProc = numelements;
2764  Array<zgno_t> elementList(numelements);
2765  for (Array<zgno_t>::size_type k = 0; k < numelements; ++k) {
2766  elementList[k] = pamgen_mesh->element_order_map[k];
2767  }
2768 
2769  mp = rcp (new map_t (numGlobalElements, elementList, 0, this->tcomm_)); // constructor 2
2770 
2771 
2772  // make an array of array views containing the coordinate data
2773  Teuchos::Array<Teuchos::ArrayView<const zscalar_t> > coordView(dimension);
2774  for (int i = 0; i < dimension; i++){
2775  if(numelements > 0){
2776  Teuchos::ArrayView<const zscalar_t> a(elem_coords[i], numelements);
2777  coordView[i] = a;
2778  }
2779  else {
2780  Teuchos::ArrayView<const zscalar_t> a;
2781  coordView[i] = a;
2782  }
2783  }
2784 
2785  // set the xyz_ multivector
2786  xyz_ = RCP<tMVector_t>(new
2787  tMVector_t(mp, coordView.view(0, dimension),
2788  dimension));
2789 }
2790 
2791 void UserInputForTests::setPamgenAdjacencyGraph()
2792 {
2793 // int rank = this->tcomm_->getRank();
2794 // if(rank == 0) std::cout << "Making a graph from our pamgen mesh...." << std::endl;
2795 
2796  // Define Types
2797 // typedef zlno_t lno_t;
2798 // typedef zgno_t gno_t;
2799  typedef Tpetra::Map<zlno_t, zgno_t, znode_t> map_t;
2800 
2801  // get info for setting up map
2802  size_t local_nodes = (size_t)this->pamgen_mesh->num_nodes;
2803  size_t local_els = (size_t)this->pamgen_mesh->num_elem;
2804 
2805  size_t global_els = (size_t)this->pamgen_mesh->num_elems_global; // global rows
2806  size_t global_nodes = (size_t)this->pamgen_mesh->num_nodes_global; //global columns
2807  // make map with global elements assigned to this mesh
2808  // make range map
2809 // if(rank == 0) std::cout << "Building Rowmap: " << std::endl;
2810  RCP<const map_t> rowMap = rcp(new map_t(global_els,0,this->tcomm_));
2811  RCP<const map_t> rangeMap = rowMap;
2812 
2813  // make domain map
2814  RCP<const map_t> domainMap = rcp(new map_t(global_nodes,0,this->tcomm_));
2815 
2816  // Get max number of nodes per element
2817  int blks = this->pamgen_mesh->num_elem_blk;
2818  int max_nodes_per_el = 0;
2819  for(int i = 0; i < blks; i++)
2820  if (this->pamgen_mesh->nodes_per_element[i] > max_nodes_per_el)
2821  max_nodes_per_el = this->pamgen_mesh->nodes_per_element[i];
2822 
2823  // make the element-node adjacency matrix
2824  Teuchos::RCP<tcrsMatrix_t> C = rcp(new tcrsMatrix_t(rowMap,max_nodes_per_el));
2825 
2826  Array<zgno_t> g_el_ids(local_els);
2827  for (size_t k = 0; k < local_els; ++k) {
2828  g_el_ids[k] = pamgen_mesh->global_element_numbers[k]-1;
2829  }
2830 
2831  Array<zgno_t> g_node_ids(local_nodes);
2832  for (size_t k = 0; k < local_nodes; ++k) {
2833  g_node_ids[k] = pamgen_mesh->global_node_numbers[k]-1;
2834  }
2835 
2836  zlno_t el_no = 0;
2837  zscalar_t one = static_cast<zscalar_t>(1);
2838 
2839 // if(rank == 0) std::cout << "Writing C... " << std::endl;
2840  for(int i = 0; i < blks; i++)
2841  {
2842  int el_per_block = this->pamgen_mesh->elements[i];
2843  int nodes_per_el = this->pamgen_mesh->nodes_per_element[i];
2844  int * connect = this->pamgen_mesh->elmt_node_linkage[i];
2845 
2846  for(int j = 0; j < el_per_block; j++)
2847  {
2848  const zgno_t gid = static_cast<zgno_t>(g_el_ids[el_no]);
2849  for(int k = 0; k < nodes_per_el; k++)
2850  {
2851  int g_node_i = g_node_ids[connect[j*nodes_per_el+k]-1];
2852  C->insertGlobalValues(gid,
2853  Teuchos::tuple<zgno_t>(g_node_i),
2854  Teuchos::tuple<zscalar_t>(one));
2855  }
2856  el_no++;
2857  }
2858  }
2859 
2860  C->fillComplete(domainMap, rangeMap);
2861 
2862 
2863  // Compute product C*C'
2864 // if(rank == 0) std::cout << "Compute Multiplication C... " << std::endl;
2865  RCP<tcrsMatrix_t> A = rcp(new tcrsMatrix_t(rowMap,0));
2866  Tpetra::MatrixMatrix::Multiply(*C, false, *C, true, *A);
2867 
2868  // remove entries not adjacent
2869  // make graph
2870 // if(rank == 0) std::cout << "Writing M_... " << std::endl;
2871  this->M_ = rcp(new tcrsMatrix_t(rowMap, A->getGlobalMaxNumRowEntries()));
2872 
2873 // if(rank == 0) std::cout << "\nSetting graph of connectivity..." << std::endl;
2874  Teuchos::ArrayView<const zgno_t> rowMapElementList =
2875  rowMap->getNodeElementList();
2876  for (Teuchos_Ordinal ii = 0; ii < rowMapElementList.size(); ii++)
2877  {
2878  zgno_t gid = rowMapElementList[ii];
2879  size_t numEntriesInRow = A->getNumEntriesInGlobalRow (gid);
2880  Array<zscalar_t> rowvals (numEntriesInRow);
2881  Array<zgno_t> rowinds (numEntriesInRow);
2882 
2883  // modified
2884  Array<zscalar_t> mod_rowvals;
2885  Array<zgno_t> mod_rowinds;
2886  A->getGlobalRowCopy (gid, rowinds (), rowvals (), numEntriesInRow);
2887  for (size_t i = 0; i < numEntriesInRow; i++) {
2888 // if (rowvals[i] == 2*(this->pamgen_mesh->num_dim-1))
2889 // {
2890  if (rowvals[i] >= 1)
2891  {
2892  mod_rowvals.push_back(one);
2893  mod_rowinds.push_back(rowinds[i]);
2894  }
2895  }
2896  this->M_->insertGlobalValues(gid, mod_rowinds, mod_rowvals);
2897  }
2898 
2899  this->M_->fillComplete();
2901  // if(rank == 0) std::cout << "Completed M... " << std::endl;
2902 
2903 }
2904 
2905 #endif
2906 
2907 #endif
keep typedefs that commonly appear in many places localized
#define TEST_FAIL_AND_THROW(comm, ok, s)
RCP< tMVector_t > getUITpetraMultiVector(int nvec)
Tpetra::CrsMatrix< zscalar_t, zlno_t, zgno_t, znode_t > tcrsMatrix_t
Definition: GraphModel.cpp:82
map_t::node_type default_znode_t
Tpetra::Import< zlno_t, zgno_t, znode_t > import_t
#define Z2_FORWARD_EXCEPTIONS
Forward an exception back through call stack.
Tpetra::Map< zlno_t, zgno_t, znode_t > map_t
RCP< tcrsMatrix_t > getUITpetraCrsMatrix()
string trim_right_copy(const string &s, const string &delimiters=" \f\n\r\t\v")
void readGeoGenParams(string paramFileName, Teuchos::ParameterList &geoparams, const RCP< const Teuchos::Comm< int > > &comm)
RCP< tVector_t > getUITpetraVector()
double zscalar_t
USERINPUT_FILE_FORMATS
A class that generates typical user input for testing.
static void getUIRandomData(unsigned int seed, zlno_t length, zscalar_t min, zscalar_t max, ArrayView< ArrayRCP< zscalar_t > > data)
Generate lists of random scalars.
Traits of Xpetra classes, including migration method.
int zlno_t
common code used by tests
#define HAVE_EPETRA_DATA_TYPES
RCP< tMVector_t > getUIWeights()
static RCP< User > convertToXpetra(const RCP< User > &a)
Convert the object to its Xpetra wrapped version.
UserInputForTests(string path, string testData, const RCP< const Comm< int > > &c, bool debugInfo=false, bool distributeInput=true)
Constructor that reads in a matrix/graph from disk.
Tpetra::Export< zlno_t, zgno_t, znode_t > export_t
static RCP< tMVector_t > coordinates
string trim_copy(const string &s, const string &delimiters=" \f\n\r\t\v")
RCP< xcrsMatrix_t > getUIXpetraCrsMatrix()
const char param_comment
RCP< xVector_t > getUIXpetraVector()
dictionary vals
Definition: xml2dox.py:186
bool hasInputDataType(const string &input_type)
GeometricGen::GeometricGenerator< zscalar_t, zlno_t, zgno_t, znode_t > geometricgen_t
RCP< tMVector_t > getUIEdgeWeights()
Tpetra::CrsGraph< zlno_t, zgno_t, znode_t > tcrsGraph_t
Definition: GraphModel.cpp:83
Xpetra::CrsGraph< zlno_t, zgno_t, znode_t > xcrsGraph_t
Tpetra::Vector< zscalar_t, zlno_t, zgno_t, znode_t > tVector_t
static const std::string fail
int zgno_t
string trim_left_copy(const string &s, const string &delimiters=" \f\n\r\t\v")
RCP< xMVector_t > getUIXpetraMultiVector(int nvec)
Tpetra::MultiVector< double, int, int > tMVector_t
RCP< tMVector_t > getUICoordinates()
Tpetra::global_size_t global_size_t
RCP< xcrsGraph_t > getUIXpetraCrsGraph()
RCP< tcrsGraph_t > getUITpetraCrsGraph()