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