Zoltan2
Loading...
Searching...
No Matches
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
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 <KokkosCompat_DefaultNode.hpp>
72
74#include <fstream>
75#include <string>
76
77#include <TpetraExt_MatrixMatrix_def.hpp>
78
79// pamgen required includes
81
82
83using Teuchos::RCP;
84using Teuchos::ArrayRCP;
85using Teuchos::ArrayView;
86using Teuchos::Array;
87using Teuchos::Comm;
88using Teuchos::rcp;
89using Teuchos::arcp;
90using Teuchos::rcp_const_cast;
91using Teuchos::ParameterList;
92using namespace Zoltan2_TestingFramework;
93
125
127{
128public:
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
248
249 bool hasUITpetraCrsGraph();
250
251 bool hasUITpetraVector();
252
254
256
257 bool hasUIXpetraCrsGraph();
258
259 bool hasUIXpetraVector();
260
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
275private:
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
388UserInputForTests::UserInputForTests(string path, string testData,
389 const RCP<const Comm<int> > &c,
390 bool debugInfo, bool distributeInput):
391verbose_(debugInfo), tcomm_(c), havePamgenMesh(false),
392M_(), xM_(), xyz_(), vtxWeights_(), edgWeights_(),
393#ifdef HAVE_EPETRA_DATA_TYPES
394ecomm_(), eM_(), eG_(),
395#endif
396chaco_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):
418verbose_(debugInfo), tcomm_(c), havePamgenMesh(false),
419M_(), xM_(), xyz_(), vtxWeights_(), edgWeights_(),
420#ifdef HAVE_EPETRA_DATA_TYPES
421ecomm_(), eM_(), eG_(),
422#endif
423chaco_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
450UserInputForTests::UserInputForTests(const ParameterList &pList,
451 const RCP<const Comm<int> > &c):
452tcomm_(c), havePamgenMesh(false),
453M_(), xM_(), xyz_(), vtxWeights_(), edgWeights_(),
454#ifdef HAVE_EPETRA_DATA_TYPES
455ecomm_(), eM_(), eG_(),
456#endif
457chaco_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
544RCP<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
551RCP<Zoltan2_TestingFramework::tMVector_t> UserInputForTests::getUIWeights()
552{
553 return vtxWeights_;
554}
555
556RCP<Zoltan2_TestingFramework::tMVector_t> UserInputForTests::getUIEdgeWeights()
557{
558 return edgWeights_;
559}
560
561RCP<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
568RCP<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
575RCP<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
583RCP<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
591RCP<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
598RCP<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
605RCP<Zoltan2_TestingFramework::xVector_t> UserInputForTests::getUIXpetraVector()
606{
608}
609
610RCP<Zoltan2_TestingFramework::xMVector_t> UserInputForTests::getUIXpetraMultiVector(int nvec)
611{
612 RCP<tMVector_t> tMV = getUITpetraMultiVector(nvec);
614}
615
616#ifdef HAVE_EPETRA_DATA_TYPES
617RCP<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->getLocalNumElements());
627 int base = 0;
628 ArrayView<const int> gids = trowMap->getLocalElementList();
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_->getLocalMaxNumRowEntries();
639 Array<int> colGids(maxRow);
640
641 typename tcrsGraph_t::local_inds_host_view_type colLid;
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 (size_t j=0; j < colLid.extent(0); 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
655RCP<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_->getLocalMaxNumRowEntries();
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 for (int i=0; i < nrows; i++){
668 typename tcrsMatrix_t::local_inds_host_view_type colLid;
669 typename tcrsMatrix_t::values_host_view_type nz;
670 M_->getLocalRowView(i, colLid, nz);
671 size_t rowSize = colLid.size();
672 int rowGid = rowMap.GID(i);
673 for (size_t j=0; j < rowSize; j++){
674 colGid[j] = colMap.GID(colLid[j]);
675 }
676 eM_->InsertGlobalValues(rowGid, (int)rowSize, nz.data(), colGid.getRawPtr());
677 }
678 eM_->FillComplete();
679 return eM_;
680}
681
682RCP<Epetra_Vector> UserInputForTests::getUIEpetraVector()
683{
684 RCP<Epetra_CrsGraph> egraph = getUIEpetraCrsGraph();
685 RCP<Epetra_Vector> V = rcp(new Epetra_Vector(egraph->RowMap()));
686 V->Random();
687 return V;
688}
689
690RCP<Epetra_MultiVector> UserInputForTests::getUIEpetraMultiVector(int nvec)
691{
692 RCP<Epetra_CrsGraph> egraph = getUIEpetraCrsGraph();
693 RCP<Epetra_MultiVector> mV =
694 rcp(new Epetra_MultiVector(egraph->RowMap(), nvec));
695 mV->Random();
696 return mV;
697}
698#endif
699
701{
702 // find out if an input source has been loaded
703 return this->hasUICoordinates() || \
704 this->hasUITpetraCrsMatrix() || \
705 this->hasUITpetraCrsGraph() || \
706 this->hasPamgenMesh();
707}
708
709bool UserInputForTests::hasInputDataType(const string &input_type)
710{
711 if(input_type == "coordinates")
712 return this->hasUICoordinates();
713 else if(input_type == "tpetra_vector")
714 return this->hasUITpetraVector();
715 else if(input_type == "tpetra_multivector")
716 return this->hasUITpetraMultiVector();
717 else if(input_type == "tpetra_crs_graph")
718 return this->hasUITpetraCrsGraph();
719 else if(input_type == "tpetra_crs_matrix")
720 return this->hasUITpetraCrsMatrix();
721 else if(input_type == "xpetra_vector")
722 return this->hasUIXpetraVector();
723 else if(input_type == "xpetra_multivector")
724 return this->hasUIXpetraMultiVector();
725 else if(input_type == "xpetra_crs_graph")
726 return this->hasUIXpetraCrsGraph();
727 else if(input_type == "xpetra_crs_matrix")
728 return this->hasUIXpetraCrsMatrix();
729#ifdef HAVE_EPETRA_DATA_TYPES
730 else if(input_type == "epetra_vector")
731 return this->hasUIEpetraVector();
732 else if(input_type == "epetra_multivector")
733 return this->hasUIEpetraMultiVector();
734 else if(input_type == "epetra_crs_graph")
735 return this->hasUIEpetraCrsGraph();
736 else if(input_type == "epetra_crs_matrix")
737 return this->hasUIEpetraCrsMatrix();
738#endif
739
740 return false;
741}
742
744{
745 return xyz_.is_null() ? false : true;
746}
747
749{
750 return vtxWeights_.is_null() ? false : true;
751}
752
754{
755 return edgWeights_.is_null() ? false : true;
756}
757
759{
760 return M_.is_null() ? false : true;
761}
762
764{
765 return M_.is_null() ? false : true;
766}
767
769{
770 return true;
771}
772
774{
775 return true;
776}
777
779{
780 return M_.is_null() ? false : true;
781}
782
784{
785 return M_.is_null() ? false : true;
786}
787
789{
790 return true;
791}
792
794{
795 return true;
796}
797
799{
800 return this->havePamgenMesh;
801}
802
803#ifdef HAVE_EPETRA_DATA_TYPES
804bool UserInputForTests::hasUIEpetraCrsGraph()
805{
806 return M_.is_null() ? false : true;
807}
808
809bool UserInputForTests::hasUIEpetraCrsMatrix()
810{
811 return hasUIEpetraCrsGraph();
812}
813
814bool UserInputForTests::hasUIEpetraVector()
815{
816 return hasUIEpetraCrsGraph();
817}
818
819bool UserInputForTests::hasUIEpetraMultiVector()
820{
821 return hasUIEpetraCrsGraph();
822}
823#endif
824
825void UserInputForTests::getUIRandomData(unsigned int seed, zlno_t length,
826 zscalar_t min, zscalar_t max,
827 ArrayView<ArrayRCP<zscalar_t > > data)
828{
829 if (length < 1)
830 return;
831
832 size_t dim = data.size();
833 for (size_t i=0; i < dim; i++){
834 zscalar_t *tmp = new zscalar_t [length];
835 if (!tmp)
836 throw (std::bad_alloc());
837 data[i] = Teuchos::arcp(tmp, 0, length, true);
838 }
839
840 zscalar_t scalingFactor = (max-min) / RAND_MAX;
841 srand(seed);
842 for (size_t i=0; i < dim; i++){
843 zscalar_t *x = data[i].getRawPtr();
844 for (zlno_t j=0; j < length; j++)
845 *x++ = min + (zscalar_t(rand()) * scalingFactor);
846 }
847}
848
849// utility methods used when reading geom gen files
850
851string UserInputForTests::trim_right_copy(
852 const string& s,
853 const string& delimiters)
854{
855 return s.substr( 0, s.find_last_not_of( delimiters ) + 1 );
856}
857
858string UserInputForTests::trim_left_copy(
859 const string& s,
860 const string& delimiters)
861{
862 return s.substr( s.find_first_not_of( delimiters ) );
863}
864
865string UserInputForTests::trim_copy(
866 const string& s,
867 const string& delimiters)
868{
869 return trim_left_copy( trim_right_copy( s, delimiters ), delimiters );
870}
871
872void UserInputForTests::readGeometricGenTestData(string path,
873 string testData)
874{
875
876 std::ostringstream fname;
877 fname << path << "/" << testData << ".txt";
878
879 if (verbose_ && tcomm_->getRank() == 0)
880 std::cout << "UserInputForTests, Read: " << fname.str() << std::endl;
881
882 Teuchos::ParameterList geoparams("geo params");
883 readGeoGenParams(fname.str(),geoparams);
884
885 geometricgen_t * gg = new geometricgen_t(geoparams, this->tcomm_);
886
887 // get coordinate and point info
888 int coord_dim = gg->getCoordinateDimension();
889 int numWeightsPerCoord = gg->getNumWeights();
890 zlno_t numLocalPoints = gg->getNumLocalCoords();
891 zgno_t numGlobalPoints = gg->getNumGlobalCoords();
892
893 // allocate an array of coordinate arrays
894 zscalar_t **coords = new zscalar_t * [coord_dim];
895 for(int i = 0; i < coord_dim; ++i){
896 coords[i] = new zscalar_t[numLocalPoints];
897 }
898
899 // get a copy of the data
900 gg->getLocalCoordinatesCopy(coords);
901
902 // get an array of arrays of weight data (if any)
903 zscalar_t **weight = NULL;
904 if (numWeightsPerCoord) {
905 // memory allocation
906 weight = new zscalar_t * [numWeightsPerCoord];
907 for(int i = 0; i < numWeightsPerCoord; ++i){
908 weight[i] = new zscalar_t[numLocalPoints];
909 }
910
911 // get a copy of the weight data
912 gg->getLocalWeightsCopy(weight);
913 }
914
915 delete gg; // free up memory from geom gen
916
917
918 // make a Tpetra map
919 RCP<Tpetra::Map<zlno_t, zgno_t, znode_t> > mp =
920 rcp(new Tpetra::Map<zlno_t, zgno_t, znode_t>(numGlobalPoints, numLocalPoints, 0, this->tcomm_));
921
922 // make an array of array views containing the coordinate data
923 Teuchos::Array<Teuchos::ArrayView<const zscalar_t> > coordView(coord_dim);
924 for (int i=0; i < coord_dim; i++){
925 if(numLocalPoints > 0){
926 Teuchos::ArrayView<const zscalar_t> a(coords[i], numLocalPoints);
927 coordView[i] = a;
928 }
929 else {
930 Teuchos::ArrayView<const zscalar_t> a;
931 coordView[i] = a;
932 }
933 }
934
935 // set the xyz_ multivector
936 xyz_ = RCP<tMVector_t>(new
937 tMVector_t(mp, coordView.view(0, coord_dim),
938 coord_dim));
939
940 // set the vtx weights
941 if (numWeightsPerCoord) {
942 // make an array of array views containing the weight data
943 Teuchos::Array<Teuchos::ArrayView<const zscalar_t> > weightView(numWeightsPerCoord);
944 for (int i=0; i < numWeightsPerCoord; i++){
945 if(numLocalPoints > 0){
946 Teuchos::ArrayView<const zscalar_t> a(weight[i], numLocalPoints);
947 weightView[i] = a;
948 }
949 else {
950 Teuchos::ArrayView<const zscalar_t> a;
951 weightView[i] = a;
952 }
953 }
954
955 vtxWeights_ = RCP<tMVector_t>(new tMVector_t(mp, weightView.view(0, numWeightsPerCoord),
956 numWeightsPerCoord));
957 }
958}
959
960void UserInputForTests::readGeoGenParams(string paramFileName,
961 ParameterList &geoparams){
962
963 const char param_comment = '#';
964
965 std::string input = "";
966 char inp[25000];
967 for(int i = 0; i < 25000; ++i){
968 inp[i] = 0;
969 }
970
971 bool fail = false;
972 if(this->tcomm_->getRank() == 0){
973
974 std::fstream inParam(paramFileName.c_str());
975 if (inParam.fail())
976 {
977 fail = true;
978 }
979 if(!fail)
980 {
981 std::string tmp = "";
982 getline (inParam,tmp);
983 while (!inParam.eof()){
984 if(tmp != ""){
985 tmp = trim_copy(tmp);
986 if(tmp != ""){
987 input += tmp + "\n";
988 }
989 }
990 getline (inParam,tmp);
991 }
992 inParam.close();
993 for (size_t i = 0; i < input.size(); ++i){
994 inp[i] = input[i];
995 }
996 }
997 }
998
999
1000
1001 int size = (int)input.size();
1002 if(fail){
1003 size = -1;
1004 }
1005 this->tcomm_->broadcast(0, sizeof(int), (char*) &size);
1006 if(size == -1){
1007 throw "File " + paramFileName + " cannot be opened.";
1008 }
1009 this->tcomm_->broadcast(0, size, inp);
1010 std::istringstream inParam(inp);
1011 string str;
1012 getline (inParam,str);
1013 while (!inParam.eof()){
1014 if(str[0] != param_comment){
1015 size_t pos = str.find('=');
1016 if(pos == string::npos){
1017 throw "Invalid Line:" + str + " in parameter file";
1018 }
1019 string paramname = trim_copy(str.substr(0,pos));
1020 string paramvalue = trim_copy(str.substr(pos + 1));
1021 geoparams.set(paramname, paramvalue);
1022 }
1023 getline (inParam,str);
1024 }
1025}
1026
1028RCP<tcrsMatrix_t> UserInputForTests::modifyMatrixGIDs(
1029 RCP<tcrsMatrix_t> &inMatrix
1030)
1031{
1032 // Produce a new matrix with the same structure as inMatrix,
1033 // but whose row/column GIDs are non-contiguous values.
1034 // In this case GID g in inMatrix becomes g*2+1 in outMatrix.
1035
1036 // Create the map for the new matrix: same structure as inMap but with
1037 // the GIDs modified.
1038 RCP<const map_t> inMap = inMatrix->getRowMap();
1039
1040 size_t nRows = inMap->getLocalNumElements();
1041 auto inRows = inMap->getMyGlobalIndices();
1042 Teuchos::Array<zgno_t> outRows(nRows);
1043 for (size_t i = 0; i < nRows; i++) {
1044 outRows[i] = newID(inRows[i]);
1045 }
1046
1047 Tpetra::global_size_t nGlobalRows = inMap->getGlobalNumElements();
1048 RCP<map_t> outMap = rcp(new map_t(nGlobalRows, outRows(), 0,
1049 inMap->getComm()));
1050
1051#ifdef INCLUDE_LENGTHY_OUTPUT
1052 // Sanity check output
1053 {
1054 std::cout << inMap->getComm()->getRank() << " KDDKDD "
1055 << "nGlobal " << inMap->getGlobalNumElements() << " "
1056 << outMap->getGlobalNumElements() << "; "
1057 << "nLocal " << inMap->getLocalNumElements() << " "
1058 << outMap->getLocalNumElements() << "; "
1059 << std::endl;
1060 std::cout << inMap->getComm()->getRank() << " KDDKDD ";
1061 for (size_t i = 0; i < nRows; i++)
1062 std::cout << "(" << inMap->getMyGlobalIndices()[i] << ", "
1063 << outMap->getMyGlobalIndices()[i] << ") ";
1064 std::cout << std::endl;
1065 }
1066#endif // INCLUDE_LENGTHY_OUTPUT
1067
1068 // Create a new matrix using the new map
1069 // Get the length of the longest row; allocate memory.
1070 size_t rowLen = inMatrix->getLocalMaxNumRowEntries();
1071 RCP<tcrsMatrix_t> outMatrix = rcp(new tcrsMatrix_t(outMap, rowLen));
1072
1073 typename tcrsMatrix_t::nonconst_global_inds_host_view_type indices("Indices", rowLen);
1074 typename tcrsMatrix_t::nonconst_values_host_view_type values("Values", rowLen);
1075
1076 for (size_t i = 0; i < nRows; i++) {
1077 size_t nEntries;
1078 zgno_t inGid = inMap->getGlobalElement(i);
1079 inMatrix->getGlobalRowCopy(inGid, indices, values, nEntries);
1080 for (size_t j = 0; j < nEntries; j++)
1081 indices[j] = newID(indices[j]);
1082
1083 zgno_t outGid = outMap->getGlobalElement(i);
1084 ArrayView<zgno_t> Indices(indices.data(), nEntries);
1085 ArrayView<zscalar_t> Values(values.data(), nEntries);
1086 outMatrix->insertGlobalValues(outGid, Indices(0, nEntries),
1087 Values(0, nEntries));
1088 }
1089 outMatrix->fillComplete();
1090
1091#ifdef INCLUDE_LENGTHY_OUTPUT
1092 // Sanity check output
1093 {
1094 std::cout << inMap->getComm()->getRank() << " KDDKDD Rows "
1095 << "nGlobal " << inMatrix->getGlobalNumRows() << " "
1096 << outMatrix->getGlobalNumRows() << "; "
1097 << "nLocal " << inMatrix->getLocalNumRows() << " "
1098 << outMatrix->getLocalNumRows() << std::endl;
1099 std::cout << inMap->getComm()->getRank() << " KDDKDD NNZS "
1100 << "nGlobal " << inMatrix->getGlobalNumEntries() << " "
1101 << outMatrix->getGlobalNumEntries() << "; "
1102 << "nLocal " << inMatrix->getLocalNumEntries() << " "
1103 << outMatrix->getLocalNumEntries() << std::endl;
1104
1105 size_t nIn, nOut;
1106 Teuchos::Array<zgno_t> in(rowLen), out(rowLen);
1107 Teuchos::Array<zscalar_t> inval(rowLen), outval(rowLen);
1108
1109 for (size_t i = 0; i < nRows; i++) {
1110 std::cout << inMap->getComm()->getRank() << " KDDKDD " << i << " nnz(";
1111 inMatrix->getGlobalRowCopy(inMap->getGlobalElement(i), in, inval, nIn);
1112 outMatrix->getGlobalRowCopy(outMap->getGlobalElement(i), out, outval,
1113 nOut);
1114
1115 std::cout << nIn << ", " << nOut << "): ";
1116 for (size_t j = 0; j < nIn; j++) {
1117 std::cout << "(" << in[j] << " " << inval[j] << ", "
1118 << out[j] << " " << outval[j] << ") ";
1119 }
1120 std::cout << std::endl;
1121 }
1122 }
1123#endif // INCLUDE_LENGTHY_OUTPUT
1124
1125 return outMatrix;
1126}
1127
1129
1130void UserInputForTests::readMatrixMarketFile(
1131 string path,
1132 string testData,
1133 bool distributeInput
1134)
1135{
1136 std::ostringstream fname;
1137 fname << path << "/" << testData << ".mtx";
1138
1139 if (verbose_ && tcomm_->getRank() == 0)
1140 std::cout << "UserInputForTests, Read: " << fname.str() << std::endl;
1141
1142 // FIXME (mfh 01 Aug 2016) Tpetra::MatrixMarket::Reader has a graph
1143 // ("pattern" matrix) reader. Call its readSparseGraphFile method.
1144
1145 RCP<tcrsMatrix_t> toMatrix;
1146 RCP<tcrsMatrix_t> fromMatrix;
1147 bool aok = true;
1148 try{
1149 typedef Tpetra::MatrixMarket::Reader<tcrsMatrix_t> reader_type;
1150 fromMatrix = reader_type::readSparseFile(fname.str(), tcomm_,
1151 true, false, false);
1152#ifdef KDD_NOT_READY_YET
1153 // See note below about modifying coordinate IDs as well.
1154 //if (makeNonContiguous)
1155 fromMatrix = modifyMatrixGIDs(fromMatrix);
1156#endif
1157
1158 if(!distributeInput)
1159 {
1160 if (verbose_ && tcomm_->getRank() == 0)
1161 std::cout << "Constructing serial distribution of matrix" << std::endl;
1162 // need to make a serial map and then import the data to redistribute it
1163 RCP<const map_t> fromMap = fromMatrix->getRowMap();
1164
1165 size_t numGlobalCoords = fromMap->getGlobalNumElements();
1166 size_t numLocalCoords = this->tcomm_->getRank() == 0 ? numGlobalCoords : 0;
1167 RCP<const map_t> toMap = rcp(new map_t(numGlobalCoords,numLocalCoords, 0, tcomm_));
1168
1169 RCP<import_t> importer = rcp(new import_t(fromMap, toMap));
1170 toMatrix = rcp(new tcrsMatrix_t(toMap,0));
1171 toMatrix->doImport(*fromMatrix, *importer, Tpetra::INSERT);
1172 toMatrix->fillComplete();
1173
1174 }else{
1175 toMatrix = fromMatrix;
1176 }
1177 }catch (std::exception &e) {
1178 if (tcomm_->getRank() == 0) {
1179 std::cout << "UserInputForTests unable to read matrix market file:"
1180 << fname.str() << std::endl;
1181 std::cout << e.what() << std::endl;
1182 }
1183 aok = false;
1184 }
1185 TEST_FAIL_AND_THROW(*tcomm_, aok,
1186 "UserInputForTests unable to read matrix market file");
1187
1188 M_ = toMatrix;
1189#ifdef INCLUDE_LENGTHY_OUTPUT
1190 std::cout << tcomm_->getRank() << " KDDKDD " << M_->getLocalNumRows()
1191 << " " << M_->getGlobalNumRows()
1192 << " " << M_->getLocalNumEntries()
1193 << " " << M_->getGlobalNumEntries() << std::endl;
1194#endif // INCLUDE_LENGTHY_OUTPUT
1195
1197
1198 // Open the coordinate file.
1199
1200 fname.str("");
1201 fname << path << "/" << testData << "_coord.mtx";
1202
1203 size_t coordDim = 0, numGlobalCoords = 0;
1204 size_t msg[2]={0,0};
1205 ArrayRCP<ArrayRCP<zscalar_t> > xyz;
1206 std::ifstream coordFile;
1207
1208 if (tcomm_->getRank() == 0){
1209
1210 if (verbose_)
1211 std::cout << "UserInputForTests, Read: " <<
1212 fname.str() << std::endl;
1213
1214 int fail = 0;
1215 try{
1216 coordFile.open(fname.str().c_str());
1217 }
1218 catch (std::exception &){ // there is no coordinate file
1219 fail = 1;
1220 }
1221
1222 if (!fail){
1223
1224 // Read past banner to number and dimension of coordinates.
1225
1226 char c[256];
1227 bool done=false;
1228
1229 while (!done && !fail && coordFile.good()){
1230 coordFile.getline(c, 256);
1231 if (!c[0])
1232 fail = 1;
1233 else if (c[0] == '%')
1234 continue;
1235 else {
1236 done=true;
1237 std::istringstream s(c);
1238 s >> numGlobalCoords >> coordDim;
1239 if (!s.eof() || numGlobalCoords < 1 || coordDim < 1)
1240 fail=1;
1241 }
1242 }
1243
1244 if (done){
1245
1246 // Read in the coordinates.
1247
1248 xyz = Teuchos::arcp(new ArrayRCP<zscalar_t> [coordDim], 0, coordDim);
1249
1250 for (size_t dim=0; !fail && dim < coordDim; dim++){
1251 size_t idx;
1252 zscalar_t *tmp = new zscalar_t [numGlobalCoords];
1253 if (!tmp)
1254 fail = 1;
1255 else{
1256 xyz[dim] = Teuchos::arcp(tmp, 0, numGlobalCoords);
1257
1258 for (idx=0; !coordFile.eof() && idx < numGlobalCoords; idx++){
1259 coordFile.getline(c, 256);
1260 std::istringstream s(c);
1261 s >> tmp[idx];
1262 }
1263
1264 if (idx < numGlobalCoords)
1265 fail = 1;
1266 }
1267 }
1268
1269 if (fail){
1270 ArrayRCP<zscalar_t> emptyArray;
1271 for (size_t dim=0; dim < coordDim; dim++)
1272 xyz[dim] = emptyArray; // free the memory
1273
1274 coordDim = 0;
1275 }
1276 }
1277 else{
1278 fail = 1;
1279 }
1280
1281 coordFile.close();
1282 }
1283
1284 msg[0] = coordDim;
1285 msg[1] = numGlobalCoords;
1286 }
1287
1288 // Broadcast coordinate dimension
1289 Teuchos::broadcast<int, size_t>(*tcomm_, 0, 2, msg);
1290
1291 coordDim = msg[0];
1292 numGlobalCoords = msg[1];
1293
1294 if (coordDim == 0)
1295 return;
1296
1297 zgno_t base = 0;
1298 RCP<const map_t> toMap;
1299
1300 if (!M_.is_null()){
1301 const RCP<const map_t> &mapM = M_->getRowMap();
1302 toMap = mapM;
1303 }
1304 else{
1305 if (verbose_ && tcomm_->getRank() == 0)
1306 {
1307 std::cout << "Matrix was null. ";
1308 std::cout << "Constructing distribution map for coordinate vector."
1309 << std::endl;
1310 }
1311
1312 if(!distributeInput)
1313 {
1314 if (verbose_ && tcomm_->getRank() == 0)
1315 std::cout << "Constructing serial distribution map for coordinates."
1316 << std::endl;
1317
1318 size_t numLocalCoords = this->tcomm_->getRank()==0 ? numGlobalCoords : 0;
1319 toMap = rcp(new map_t(numGlobalCoords,numLocalCoords, base, tcomm_));
1320 }else{
1321 toMap = rcp(new map_t(numGlobalCoords, base, tcomm_));
1322 }
1323 }
1324
1325 // Export coordinates to their owners
1326
1327 xyz_ = rcp(new tMVector_t(toMap, coordDim));
1328
1329 ArrayRCP<ArrayView<const zscalar_t> > coordLists(coordDim);
1330
1331 if (tcomm_->getRank() == 0){
1332
1333 for (size_t dim=0; dim < coordDim; dim++)
1334 coordLists[dim] = xyz[dim].view(0, numGlobalCoords);
1335
1336 zgno_t *tmp = new zgno_t [numGlobalCoords];
1337 if (!tmp)
1338 throw std::bad_alloc();
1339
1340 ArrayRCP<const zgno_t> rowIds = Teuchos::arcp(tmp, 0, numGlobalCoords);
1341
1342#ifdef KDD_NOT_READY_YET
1343 // TODO if modifyMatrixGIDs, we need to modify ids here as well
1344 for (zgno_t id=0; id < zgno_t(numGlobalCoords); id++)
1345 *tmp++ = newID(id);
1346#else
1347 for (zgno_t id=0; id < zgno_t(numGlobalCoords); id++)
1348 *tmp++ = id;
1349#endif
1350
1351 RCP<const map_t> fromMap = rcp(new map_t(numGlobalCoords,
1352 rowIds.view(0, numGlobalCoords),
1353 base, tcomm_));
1354
1355 tMVector_t allCoords(fromMap, coordLists.view(0, coordDim), coordDim);
1356
1357 export_t exporter(fromMap, toMap);
1358
1359 xyz_->doExport(allCoords, exporter, Tpetra::INSERT);
1360 }
1361 else{
1362
1363 RCP<const map_t> fromMap = rcp(new map_t(numGlobalCoords,
1364 ArrayView<zgno_t>(), base, tcomm_));
1365
1366 tMVector_t allCoords(fromMap, coordLists.view(0, coordDim), coordDim);
1367
1368 export_t exporter(fromMap, toMap);
1369
1370 xyz_->doExport(allCoords, exporter, Tpetra::INSERT);
1371 }
1372}
1373
1374void UserInputForTests::buildCrsMatrix(int xdim, int ydim, int zdim,
1375 string problemType, bool distributeInput)
1376{
1377#ifdef HAVE_ZOLTAN2_GALERI
1378 Teuchos::CommandLineProcessor tclp;
1379 Galeri::Xpetra::Parameters<zgno_t> params(tclp,
1380 xdim, ydim, zdim, problemType);
1381
1382 RCP<const Tpetra::Map<zlno_t, zgno_t> > map;
1383 if (distributeInput)
1384 map = rcp(new Tpetra::Map<zlno_t, zgno_t>(params.GetNumGlobalElements(),
1385 0, tcomm_));
1386 else {
1387 // All data initially on rank 0
1388 size_t nGlobalElements = params.GetNumGlobalElements();
1389 size_t nLocalElements = ((tcomm_->getRank() == 0) ? nGlobalElements : 0);
1390 map = rcp(new Tpetra::Map<zlno_t, zgno_t>(nGlobalElements, nLocalElements, 0,
1391 tcomm_));
1392 }
1393
1394 if (verbose_ && tcomm_->getRank() == 0){
1395
1396 std::cout << "Matrix is " << (distributeInput ? "" : "not");
1397 std::cout << "distributed." << std::endl;
1398
1399 std::cout << "UserInputForTests, Create matrix with " << problemType;
1400 std::cout << " (and " << xdim;
1401 if (zdim > 0)
1402 std::cout << " x " << ydim << " x " << zdim;
1403 else if (ydim > 0)
1404 std::cout << " x" << ydim << " x 1";
1405 else
1406 std::cout << "x 1 x 1";
1407
1408 std::cout << " mesh)" << std::endl;
1409
1410 }
1411
1412 bool aok = true;
1413 try{
1414 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 =
1415 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> >
1416 (params.GetMatrixType(), map, params.GetParameterList());
1417 M_ = Pr->BuildMatrix();
1418 }
1419 catch (std::exception &) { // Probably not enough memory
1420 aok = false;
1421 }
1422 TEST_FAIL_AND_THROW(*tcomm_, aok,
1423 "UserInputForTests Galeri::Xpetra::BuildProblem failed");
1424
1426
1427 // Compute the coordinates for the matrix rows.
1428
1429 if (verbose_ && tcomm_->getRank() == 0)
1430 std::cout <<
1431 "UserInputForTests, Implied matrix row coordinates computed" <<
1432 std::endl;
1433
1434 ArrayView<const zgno_t> gids = map->getLocalElementList();
1435 zlno_t count = static_cast<zlno_t>(gids.size());
1436 int dim = 3;
1437 size_t pos = problemType.find("2D");
1438 if (pos != string::npos)
1439 dim = 2;
1440 else if (problemType == string("Laplace1D") ||
1441 problemType == string("Identity"))
1442 dim = 1;
1443
1444 Array<ArrayRCP<zscalar_t> > coordinates(dim);
1445
1446 if (count > 0){
1447 for (int i=0; i < dim; i++){
1448 zscalar_t *c = new zscalar_t [count];
1449 if (!c)
1450 throw(std::bad_alloc());
1451 coordinates[i] = Teuchos::arcp(c, 0, count, true);
1452 }
1453
1454 if (dim==3){
1455 zscalar_t *x = coordinates[0].getRawPtr();
1456 zscalar_t *y = coordinates[1].getRawPtr();
1457 zscalar_t *z = coordinates[2].getRawPtr();
1458 zgno_t xySize = xdim * ydim;
1459 for (zlno_t i=0; i < count; i++){
1460 zgno_t iz = gids[i] / xySize;
1461 zgno_t xy = gids[i] - iz*xySize;
1462 z[i] = zscalar_t(iz);
1463 y[i] = zscalar_t(xy / xdim);
1464 x[i] = zscalar_t(xy % xdim);
1465 }
1466 }
1467 else if (dim==2){
1468 zscalar_t *x = coordinates[0].getRawPtr();
1469 zscalar_t *y = coordinates[1].getRawPtr();
1470 for (zlno_t i=0; i < count; i++){
1471 y[i] = zscalar_t(gids[i] / xdim);
1472 x[i] = zscalar_t(gids[i] % xdim);
1473 }
1474 }
1475 else{
1476 zscalar_t *x = coordinates[0].getRawPtr();
1477 for (zlno_t i=0; i < count; i++)
1478 x[i] = zscalar_t(gids[i]);
1479 }
1480 }
1481
1482 Array<ArrayView<const zscalar_t> > coordView(dim);
1483 if (count > 0)
1484 for (int i=0; i < dim; i++)
1485 coordView[i] = coordinates[i].view(0,count);
1486
1487 xyz_ = rcp(new tMVector_t(map, coordView.view(0, dim), dim));
1488#else
1489 throw std::runtime_error("Galeri input requested but Trilinos is "
1490 "not built with Galeri.");
1491#endif
1492}
1493
1494void UserInputForTests::readZoltanTestData(string path, string testData,
1495 bool distributeInput)
1496{
1497 int rank = tcomm_->getRank();
1498 FILE *graphFile = NULL;
1499 FILE *coordFile = NULL;
1500 FILE *assignFile = NULL;
1501 int fileInfo[3];
1502
1503 for (int i = 0; i < CHACO_LINE_LENGTH; i++) chaco_line[i] = '\0';
1504
1505 if (rank == 0){
1506 // set chacho graph file name
1507 std::ostringstream chGraphFileName;
1508 chGraphFileName << path << "/" << testData << ".graph";
1509
1510 // set chaco graph
1511 std::ostringstream chCoordFileName;
1512 chCoordFileName << path << "/" << testData << ".coords";
1513
1514 // set chaco graph
1515 std::ostringstream chAssignFileName;
1516 chAssignFileName << path << "/" << testData << ".assign";
1517
1518 // open file
1519 graphFile = fopen(chGraphFileName.str().c_str(), "r");
1520
1521 if(!graphFile) // maybe the user is using the default zoltan1 path convention
1522 {
1523 chGraphFileName.str("");
1524 chCoordFileName.str("");
1525 // try constructing zoltan1 paths
1526 chGraphFileName << path << "/ch_" << testData << "/" << testData << ".graph";
1527 chCoordFileName << path << "/ch_" << testData << "/" << testData << ".coords";
1528 chAssignFileName << path << "/ch_" << testData << "/" << testData << ".assign";
1529 // try to open the graph file again, if this doesn't open
1530 // the user has not provided a valid path to the file
1531 graphFile = fopen(chGraphFileName.str().c_str(), "r");
1532 }
1533
1534 memset(fileInfo, 0, sizeof(int) * 3); // set fileinfo to 0's
1535 if (graphFile){
1536 fileInfo[0] = 1;
1537 if (verbose_ && tcomm_->getRank() == 0)
1538 std::cout << "UserInputForTests, open " <<
1539 chGraphFileName.str () << std::endl;
1540
1541 coordFile = fopen(chCoordFileName.str().c_str(), "r");
1542 if (coordFile){
1543 fileInfo[1] = 1;
1544 if (verbose_ && tcomm_->getRank() == 0)
1545 std::cout << "UserInputForTests, open " <<
1546 chCoordFileName.str () << std::endl;
1547 }
1548
1549 assignFile = fopen(chAssignFileName.str().c_str(), "r");
1550 if (assignFile){
1551 fileInfo[2] = 1;
1552 if (verbose_ && tcomm_->getRank() == 0)
1553 std::cout << "UserInputForTests, open " <<
1554 chAssignFileName.str () << std::endl;
1555 }
1556 }else{
1557 if (verbose_ && tcomm_->getRank() == 0){
1558 std::cout << "UserInputForTests, unable to open file: ";
1559 std::cout << chGraphFileName.str() << std::endl;
1560 }
1561 }
1562 }
1563
1564 // broadcast whether we have graphs and coords to all processes
1565 Teuchos::broadcast<int, int>(*tcomm_, 0, 3, fileInfo);
1566
1567 bool haveGraph = (fileInfo[0] == 1);
1568 bool haveCoords = (fileInfo[1] == 1);
1569 bool haveAssign = (fileInfo[2] == 1);
1570
1571 if (haveGraph){
1572 // builds M_, vtxWeights_, and edgWeights_ and closes file.
1573 try{
1574 getUIChacoGraph(graphFile, haveAssign, assignFile,
1575 testData, distributeInput);
1576 }
1578
1579 if (haveCoords){
1580 // builds xyz_ and closes the file.
1581 try{
1582 getUIChacoCoords(coordFile, testData);
1583 }
1585 }
1586 }
1587
1589}
1590
1591void UserInputForTests::getUIChacoGraph(FILE *fptr, bool haveAssign,
1592 FILE *assignFile, string fname,
1593 bool distributeInput)
1594{
1595 int rank = tcomm_->getRank();
1596 int graphCounts[5];
1597 int nvtxs=0, nedges=0;
1598 int nVwgts=0, nEwgts=0;
1599 int *start = NULL, *adj = NULL;
1600 float *ewgts = NULL, *vwgts = NULL;
1601 size_t *nzPerRow = NULL;
1602 size_t maxRowLen = 0;
1603 zgno_t base = 0;
1604 ArrayRCP<const size_t> rowSizes;
1605 int fail = 0;
1606 bool haveEdges = true;
1607
1608 if (rank == 0){
1609
1610 memset(graphCounts, 0, 5*sizeof(int));
1611
1612 // Reads in the file and closes it when done.
1613 fail = chaco_input_graph(fptr, fname.c_str(), &start, &adj,
1614 &nvtxs, &nVwgts, &vwgts, &nEwgts, &ewgts);
1615
1616 // There are Zoltan2 test graphs that have no edges.
1617
1618 // nEwgts must be 1 or 0 - add error
1619
1620 if (start == NULL)
1621 haveEdges = false;
1622
1623 if (verbose_)
1624 {
1625 std::cout << "UserInputForTests, " << nvtxs << " vertices,";
1626 if (haveEdges)
1627 std::cout << start[nvtxs] << " edges,";
1628 else
1629 std::cout << "no edges,";
1630 std::cout << nVwgts << " vertex weights, ";
1631 std::cout << nEwgts << " edge weights" << std::endl;
1632 }
1633
1634 if (nvtxs==0)
1635 fail = true;
1636
1637 if (fail){
1638 Teuchos::broadcast<int, int>(*tcomm_, 0, 5, graphCounts);
1639 throw std::runtime_error("Unable to read chaco file");
1640 }
1641
1642 if (haveEdges)
1643 nedges = start[nvtxs];
1644
1645 nzPerRow = new size_t [nvtxs];
1646 if (!nzPerRow)
1647 throw std::bad_alloc();
1648 rowSizes = arcp(nzPerRow, 0, nvtxs, true);
1649
1650 if (haveEdges){
1651 for (int i=0; i < nvtxs; i++){
1652 nzPerRow[i] = start[i+1] - start[i];
1653 if (nzPerRow[i] > maxRowLen)
1654 maxRowLen = nzPerRow[i];
1655 }
1656 }
1657 else{
1658 memset(nzPerRow, 0, sizeof(size_t) * nvtxs);
1659 }
1660
1661 // Make sure base gid is zero.
1662
1663 if (nedges){
1664 int chbase = 1; // chaco input files are one-based
1665 for (int i=0; i < nedges; i++)
1666 adj[i] -= chbase;
1667 }
1668
1669 graphCounts[0] = nvtxs;
1670 graphCounts[1] = nedges;
1671 graphCounts[2] = nVwgts;
1672 graphCounts[3] = nEwgts;
1673 graphCounts[4] = (int)maxRowLen; // size_t maxRowLen will fit; it is <= (int-int)
1674 }
1675
1676 Teuchos::broadcast<int, int>(*tcomm_, 0, 5, graphCounts);
1677
1678 if (graphCounts[0] == 0)
1679 throw std::runtime_error("Unable to read chaco file");
1680
1681 haveEdges = (graphCounts[1] > 0);
1682
1683 RCP<tcrsMatrix_t> fromMatrix;
1684 RCP<const map_t> fromMap;
1685
1686 // Create a Tpetra::CrsMatrix where rank 0 has entire matrix.
1687 if (rank == 0){
1688 fromMap = rcp(new map_t(nvtxs, nvtxs, base, tcomm_));
1689
1690 fromMatrix =
1691 rcp(new tcrsMatrix_t(fromMap, rowSizes()));
1692
1693 if (haveEdges){
1694
1695 zgno_t *edgeIds = new zgno_t [nedges];
1696 if (nedges && !edgeIds)
1697 throw std::bad_alloc();
1698 for (int i=0; i < nedges; i++)
1699 edgeIds[i] = adj[i];
1700
1701 free(adj);
1702 adj = NULL;
1703
1704 zgno_t *nextId = edgeIds;
1705 Array<zscalar_t> values(nedges, 1.0);
1706 if (nedges > 0 && nEwgts > 0) {
1707 for (int i=0; i < nedges; i++)
1708 values[i] = ewgts[i];
1709 free(ewgts);
1710 ewgts = NULL;
1711 }
1712
1713 for (int i=0; i < nvtxs; i++){
1714 if (nzPerRow[i] > 0){
1715 ArrayView<const zgno_t> rowNz(nextId, nzPerRow[i]);
1716 fromMatrix->insertGlobalValues(i, rowNz, values.view(start[i], start[i+1] - start[i]));
1717 nextId += nzPerRow[i];
1718 }
1719 }
1720
1721 delete [] edgeIds;
1722 edgeIds = NULL;
1723 }
1724
1725 fromMatrix->fillComplete();
1726 }
1727 else{
1728 nvtxs = graphCounts[0];
1729 nedges = graphCounts[1];
1730 nVwgts = graphCounts[2];
1731 nEwgts = graphCounts[3];
1732 maxRowLen = graphCounts[4];
1733
1734 // Create a Tpetra::CrsMatrix where rank 0 has entire matrix.
1735
1736 fromMap = rcp(new map_t(nvtxs, 0, base, tcomm_));
1737
1738 fromMatrix =
1739 rcp(new tcrsMatrix_t(fromMap, rowSizes()));
1740
1741 fromMatrix->fillComplete();
1742 }
1743
1744#ifdef KDDKDDPRINT
1745 if (rank == 0) {
1746 size_t sz = fromMatrix->getLocalMaxNumRowEntries();
1747 Teuchos::Array<zgno_t> indices(sz);
1748 Teuchos::Array<zscalar_t> values(sz);
1749 for (size_t i = 0; i < fromMatrix->getLocalNumRows(); i++) {
1750 zgno_t gid = fromMatrix->getRowMap()->getGlobalElement(i);
1751 size_t num;
1752 fromMatrix->getGlobalRowCopy(gid, indices(), values(), num);
1753 std::cout << "ROW " << gid << ": ";
1754 for (size_t j = 0; j < num; j++)
1755 std::cout << indices[j] << " ";
1756 std::cout << std::endl;
1757 }
1758 }
1759#endif
1760
1761 RCP<const map_t> toMap;
1762 RCP<tcrsMatrix_t> toMatrix;
1763 RCP<import_t> importer;
1764
1765 if (distributeInput) {
1766 if (haveAssign) {
1767 // Read assignments from Chaco assignment file
1768 short *assignments = new short[nvtxs];
1769 if (rank == 0) {
1770 fail = chaco_input_assign(assignFile, fname.c_str(), nvtxs, assignments);
1771 }
1772 // Broadcast coordinate dimension
1773 Teuchos::broadcast<int, short>(*tcomm_, 0, nvtxs, assignments);
1774
1775 // Build map with my vertices
1776 Teuchos::Array<zgno_t> mine;
1777 for (int i = 0; i < nvtxs; i++) {
1778 if (assignments[i] == rank)
1779 mine.push_back(i);
1780 }
1781
1782 Tpetra::global_size_t dummy =
1783 Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
1784 toMap = rcp(new map_t(dummy, mine(), base, tcomm_));
1785 delete [] assignments;
1786 }
1787 else {
1788 // Create a Tpetra::Map with default row distribution.
1789 toMap = rcp(new map_t(nvtxs, base, tcomm_));
1790 }
1791 toMatrix = rcp(new tcrsMatrix_t(toMap, maxRowLen));
1792
1793 // Import the data.
1794 importer = rcp(new import_t(fromMap, toMap));
1795 toMatrix->doImport(*fromMatrix, *importer, Tpetra::INSERT);
1796 toMatrix->fillComplete();
1797 }
1798 else {
1799 toMap = fromMap;
1800 toMatrix = fromMatrix;
1801 }
1802
1803 M_ = toMatrix;
1804
1805 // Vertex weights, if any
1806
1807 typedef ArrayRCP<const ArrayView<const zscalar_t> > arrayArray_t;
1808
1809 if (nVwgts > 0){
1810
1811 ArrayRCP<zscalar_t> weightBuf;
1812 ArrayView<const zscalar_t> *wgts = new ArrayView<const zscalar_t> [nVwgts];
1813
1814 if (rank == 0){
1815 size_t len = nVwgts * nvtxs;
1816 zscalar_t *buf = new zscalar_t [len];
1817 if (!buf) throw std::bad_alloc();
1818 weightBuf = arcp(buf, 0, len, true);
1819
1820 for (int widx=0; widx < nVwgts; widx++){
1821 wgts[widx] = ArrayView<const zscalar_t>(buf, nvtxs);
1822 float *vw = vwgts + widx;
1823 for (int i=0; i < nvtxs; i++, vw += nVwgts)
1824 buf[i] = *vw;
1825 buf += nvtxs;
1826 }
1827
1828 free(vwgts);
1829 vwgts = NULL;
1830 }
1831
1832 arrayArray_t vweights = arcp(wgts, 0, nVwgts, true);
1833
1834 RCP<tMVector_t> fromVertexWeights =
1835 rcp(new tMVector_t(fromMap, vweights.view(0, nVwgts), nVwgts));
1836
1837 RCP<tMVector_t> toVertexWeights;
1838 if (distributeInput) {
1839 toVertexWeights = rcp(new tMVector_t(toMap, nVwgts));
1840 toVertexWeights->doImport(*fromVertexWeights, *importer, Tpetra::INSERT);
1841 }
1842 else
1843 toVertexWeights = fromVertexWeights;
1844
1845 vtxWeights_ = toVertexWeights;
1846 }
1847
1848 // Edge weights, if any
1849
1850 if (haveEdges && nEwgts > 0){
1851
1852 // No longer distributing edge weights; they are the matrix values
1853 /*
1854 ArrayRCP<zscalar_t> weightBuf;
1855 ArrayView<const zscalar_t> *wgts = new ArrayView<const zscalar_t> [nEwgts];
1856
1857 toMap = rcp(new map_t(nedges, M_->getLocalNumEntries(), base, tcomm_));
1858
1859 if (rank == 0){
1860 size_t len = nEwgts * nedges;
1861 zscalar_t *buf = new zscalar_t [len];
1862 if (!buf) throw std::bad_alloc();
1863 weightBuf = arcp(buf, 0, len, true);
1864
1865 for (int widx=0; widx < nEwgts; widx++){
1866 wgts[widx] = ArrayView<const zscalar_t>(buf, nedges);
1867 float *ew = ewgts + widx;
1868 for (int i=0; i < nedges; i++, ew += nEwgts)
1869 buf[i] = *ew;
1870 buf += nedges;
1871 }
1872
1873 free(ewgts);
1874 ewgts = NULL;
1875 fromMap = rcp(new map_t(nedges, nedges, base, tcomm_));
1876 }
1877 else{
1878 fromMap = rcp(new map_t(nedges, 0, base, tcomm_));
1879 }
1880
1881 arrayArray_t eweights = arcp(wgts, 0, nEwgts, true);
1882
1883 RCP<tMVector_t> fromEdgeWeights;
1884 RCP<tMVector_t> toEdgeWeights;
1885 RCP<import_t> edgeImporter;
1886
1887 if (distributeInput) {
1888 fromEdgeWeights =
1889 rcp(new tMVector_t(fromMap, eweights.view(0, nEwgts), nEwgts));
1890 toEdgeWeights = rcp(new tMVector_t(toMap, nEwgts));
1891 edgeImporter = rcp(new import_t(fromMap, toMap));
1892 toEdgeWeights->doImport(*fromEdgeWeights, *edgeImporter, Tpetra::INSERT);
1893 }
1894 else
1895 toEdgeWeights = fromEdgeWeights;
1896
1897 edgWeights_ = toEdgeWeights;
1898 */
1899
1900 toMap = rcp(new map_t(nedges, M_->getLocalNumEntries(), base, tcomm_));
1901 edgWeights_ = rcp(new tMVector_t(toMap, nEwgts));
1902
1903 size_t maxSize = M_->getLocalMaxNumRowEntries();
1904 tcrsMatrix_t::nonconst_local_inds_host_view_type colind("colind", maxSize);
1905 tcrsMatrix_t::nonconst_values_host_view_type vals("values", maxSize);
1906 size_t nEntries;
1907
1908 for (size_t i = 0, idx = 0; i < M_->getLocalNumRows(); i++) {
1909 M_->getLocalRowCopy(i, colind, vals, nEntries);
1910 for (size_t j = 0; j < nEntries; j++) {
1911 edgWeights_->replaceLocalValue(idx, 0, vals[j]); // Assuming nEwgts==1
1912 idx++;
1913 }
1914 }
1915 }
1916
1917 if (start) {
1918 free(start);
1919 start = NULL;
1920 }
1921}
1922
1923void UserInputForTests::getUIChacoCoords(FILE *fptr, string fname)
1924{
1925 int rank = tcomm_->getRank();
1926 int ndim=0;
1927 double *x=NULL, *y=NULL, *z=NULL;
1928 int fail = 0;
1929
1930 size_t globalNumVtx = M_->getGlobalNumRows();
1931
1932 if (rank == 0){
1933
1934 // This function is from the Zoltan C-library.
1935
1936 // Reads in the file and closes it when done.
1937 fail = chaco_input_geom(fptr, fname.c_str(), (int)globalNumVtx,
1938 &ndim, &x, &y, &z);
1939
1940 if (fail)
1941 ndim = 0;
1942
1943 if (verbose_){
1944 std::cout << "UserInputForTests, read " << globalNumVtx;
1945 std::cout << " " << ndim << "-dimensional coordinates." << std::endl;
1946 }
1947 }
1948
1949 Teuchos::broadcast<int, int>(*tcomm_, 0, 1, &ndim);
1950
1951 if (ndim == 0)
1952 throw std::runtime_error("Can't read coordinate file");
1953
1954 ArrayRCP<ArrayRCP<const zscalar_t> > coords(ndim);
1955 zlno_t len = 0;
1956
1957 if (rank == 0){
1958
1959 for (int dim=0; dim < ndim; dim++){
1960 zscalar_t *v = new zscalar_t [globalNumVtx];
1961 if (!v)
1962 throw std::bad_alloc();
1963 coords[dim] = arcp<const zscalar_t>(v, 0, globalNumVtx, true);
1964 double *val = (dim==0 ? x : (dim==1 ? y : z));
1965 for (size_t i=0; i < globalNumVtx; i++)
1966 v[i] = zscalar_t(val[i]);
1967
1968 free(val);
1969 }
1970
1971 len = static_cast<zlno_t>(globalNumVtx);
1972 }
1973
1974 RCP<const map_t> fromMap = rcp(new map_t(globalNumVtx, len, 0, tcomm_));
1975 RCP<const map_t> toMap = M_->getRowMap();
1976 RCP<import_t> importer = rcp(new import_t(fromMap, toMap));
1977
1978 Array<ArrayView<const zscalar_t> > coordData;
1979 for (int dim=0; dim < ndim; dim++)
1980 coordData.push_back(coords[dim].view(0, len));
1981
1982 RCP<tMVector_t> fromCoords =
1983 rcp(new tMVector_t(fromMap, coordData.view(0, ndim), ndim));
1984
1985 RCP<tMVector_t> toCoords = rcp(new tMVector_t(toMap, ndim));
1986
1987 toCoords->doImport(*fromCoords, *importer, Tpetra::INSERT);
1988
1989 xyz_ = toCoords;
1990
1991}
1992
1995
1996double UserInputForTests::chaco_read_val(
1997 FILE* infile, /* file to read value from */
1998 int *end_flag /* 0 => OK, 1 => EOL, -1 => EOF */
1999)
2000{
2001 double val; /* return value */
2002 char *ptr; /* ptr to next string to read */
2003 char *ptr2; /* ptr to next string to read */
2004 int length; /* length of line to read */
2005 int length_left;/* length of line still around */
2006 int white_seen; /* have I detected white space yet? */
2007 int done; /* checking for end of scan */
2008 int i; /* loop counter */
2009
2010 *end_flag = 0;
2011
2012 if (chaco_offset == 0 || chaco_offset >= chaco_break_pnt) {
2013 if (chaco_offset >= chaco_break_pnt) { /* Copy rest of line back to beginning. */
2014 length_left = CHACO_LINE_LENGTH - chaco_save_pnt - 1;
2015 ptr2 = chaco_line;
2016 ptr = &chaco_line[chaco_save_pnt];
2017 for (i=length_left; i; i--) *ptr2++ = *ptr++;
2018 length = chaco_save_pnt + 1;
2019 }
2020 else {
2021 length = CHACO_LINE_LENGTH;
2022 length_left = 0;
2023 }
2024
2025 /* Now read next line, or next segment of current one. */
2026 ptr2 = fgets(&chaco_line[length_left], length, infile);
2027
2028 if (ptr2 == (char *) NULL) { /* We've hit end of file. */
2029 *end_flag = -1;
2030 return((double) 0.0);
2031 }
2032
2033 if ((chaco_line[CHACO_LINE_LENGTH - 2] != '\n') && (chaco_line[CHACO_LINE_LENGTH - 2] != '\f')
2034 && (strlen(chaco_line) == CHACO_LINE_LENGTH - 1)){
2035 /* Line too long. Find last safe place in chaco_line. */
2036 chaco_break_pnt = CHACO_LINE_LENGTH - 1;
2037 chaco_save_pnt = chaco_break_pnt;
2038 white_seen = 0;
2039 done = 0;
2040 while (!done) {
2041 --chaco_break_pnt;
2042 if (chaco_line[chaco_break_pnt] != '\0') {
2043 if (isspace((int)(chaco_line[chaco_break_pnt]))) {
2044 if (!white_seen) {
2045 chaco_save_pnt = chaco_break_pnt + 1;
2046 white_seen = 1;
2047 }
2048 }
2049 else if (white_seen) {
2050 done= 1;
2051 }
2052 }
2053 }
2054 }
2055 else {
2056 chaco_break_pnt = CHACO_LINE_LENGTH;
2057 }
2058
2059 chaco_offset = 0;
2060 }
2061
2062 while (isspace((int)(chaco_line[chaco_offset])) && chaco_offset < CHACO_LINE_LENGTH) chaco_offset++;
2063 if (chaco_line[chaco_offset] == '%' || chaco_line[chaco_offset] == '#') {
2064 *end_flag = 1;
2065 if (chaco_break_pnt < CHACO_LINE_LENGTH) {
2066 chaco_flush_line(infile);
2067 }
2068 return((double) 0.0);
2069 }
2070
2071 ptr = &(chaco_line[chaco_offset]);
2072 val = strtod(ptr, &ptr2);
2073
2074 if (ptr2 == ptr) { /* End of input line. */
2075 chaco_offset = 0;
2076 *end_flag = 1;
2077 return((double) 0.0);
2078 }
2079 else {
2080 chaco_offset = (int) (ptr2 - chaco_line) / sizeof(char);
2081 }
2082
2083 return(val);
2084}
2085
2086
2087int UserInputForTests::chaco_read_int(
2088 FILE *infile, /* file to read value from */
2089 int *end_flag /* 0 => OK, 1 => EOL, -1 => EOF */
2090)
2091{
2092 int val; /* return value */
2093 char *ptr; /* ptr to next string to read */
2094 char *ptr2; /* ptr to next string to read */
2095 int length; /* length of line to read */
2096 int length_left; /* length of line still around */
2097 int white_seen; /* have I detected white space yet? */
2098 int done; /* checking for end of scan */
2099 int i; /* loop counter */
2100
2101 *end_flag = 0;
2102
2103 if (chaco_offset == 0 || chaco_offset >= chaco_break_pnt) {
2104 if (chaco_offset >= chaco_break_pnt) { /* Copy rest of line back to beginning. */
2105 length_left = CHACO_LINE_LENGTH - chaco_save_pnt - 1;
2106 ptr2 = chaco_line;
2107 ptr = &chaco_line[chaco_save_pnt];
2108 for (i=length_left; i; i--) *ptr2++ = *ptr++;
2109 length = chaco_save_pnt + 1;
2110 }
2111 else {
2112 length = CHACO_LINE_LENGTH;
2113 length_left = 0;
2114 }
2115
2116 /* Now read next line, or next segment of current one. */
2117 ptr2 = fgets(&chaco_line[length_left], length, infile);
2118
2119 if (ptr2 == (char *) NULL) { /* We've hit end of file. */
2120 *end_flag = -1;
2121 return(0);
2122 }
2123
2124 if ((chaco_line[CHACO_LINE_LENGTH - 2] != '\n') && (chaco_line[CHACO_LINE_LENGTH - 2] != '\f')
2125 && (strlen(chaco_line) == CHACO_LINE_LENGTH - 1)){
2126 /* Line too long. Find last safe place in line. */
2127 chaco_break_pnt = CHACO_LINE_LENGTH - 1;
2128 chaco_save_pnt = chaco_break_pnt;
2129 white_seen = 0;
2130 done = 0;
2131 while (!done) {
2132 --chaco_break_pnt;
2133 if (chaco_line[chaco_break_pnt] != '\0') {
2134 if (isspace((int)(chaco_line[chaco_break_pnt]))) {
2135 if (!white_seen) {
2136 chaco_save_pnt = chaco_break_pnt + 1;
2137 white_seen = 1;
2138 }
2139 }
2140 else if (white_seen) {
2141 done= 1;
2142 }
2143 }
2144 }
2145 }
2146 else {
2147 chaco_break_pnt = CHACO_LINE_LENGTH;
2148 }
2149
2150 chaco_offset = 0;
2151 }
2152
2153 while (isspace((int)(chaco_line[chaco_offset])) && chaco_offset < CHACO_LINE_LENGTH) chaco_offset++;
2154 if (chaco_line[chaco_offset] == '%' || chaco_line[chaco_offset] == '#') {
2155 *end_flag = 1;
2156 if (chaco_break_pnt < CHACO_LINE_LENGTH) {
2157 chaco_flush_line(infile);
2158 }
2159 return(0);
2160 }
2161
2162 ptr = &(chaco_line[chaco_offset]);
2163 val = (int) strtol(ptr, &ptr2, 10);
2164
2165 if (ptr2 == ptr) { /* End of input chaco_line. */
2166 chaco_offset = 0;
2167 *end_flag = 1;
2168 return(0);
2169 }
2170 else {
2171 chaco_offset = (int) (ptr2 - chaco_line) / sizeof(char);
2172 }
2173
2174 return(val);
2175}
2176
2177void UserInputForTests::chaco_flush_line(
2178 FILE *infile /* file to read value from */
2179)
2180{
2181 char c; /* character being read */
2182
2183 c = fgetc(infile);
2184 while (c != '\n' && c != '\f')
2185 c = fgetc(infile);
2186}
2187
2188int UserInputForTests::chaco_input_graph(
2189 FILE *fin, /* input file */
2190 const char *inname, /* name of input file */
2191 int **start, /* start of edge list for each vertex */
2192 int **adjacency, /* edge list data */
2193 int *nvtxs, /* number of vertices in graph */
2194 int *nVwgts, /* # of vertex weights per node */
2195 float **vweights, /* vertex weight list data */
2196 int *nEwgts, /* # of edge weights per edge */
2197 float **eweights /* edge weight list data */
2198)
2199{
2200 int *adjptr; /* loops through adjacency data */
2201 float *ewptr; /* loops through edge weight data */
2202 int narcs; /* number of edges expected in graph */
2203 int nedges; /* twice number of edges really in graph */
2204 int nedge; /* loops through edges for each vertex */
2205 int flag; /* condition indicator */
2206 int skip_flag; /* should this edge be ignored? */
2207 int end_flag; /* indicates end of line or file */
2208 int vtx; /* vertex in graph */
2209 int line_num; /* line number in input file */
2210 int sum_edges; /* total number of edges read so far */
2211 int option = 0; /* input option */
2212 int using_ewgts; /* are edge weights in input file? */
2213 int using_vwgts; /* are vertex weights in input file? */
2214 int vtxnums; /* are vertex numbers in input file? */
2215 int vertex; /* current vertex being read */
2216 int new_vertex; /* new vertex being read */
2217 float weight; /* weight being read */
2218 float eweight; /* edge weight being read */
2219 int neighbor; /* neighbor of current vertex */
2220 int error_flag; /* error reading input? */
2221 int j; /* loop counters */
2222
2223 /* Read first line of input (= nvtxs, narcs, option). */
2224 /* The (decimal) digits of the option variable mean: 1's digit not zero => input
2225 edge weights 10's digit not zero => input vertex weights 100's digit not zero
2226 => include vertex numbers */
2227
2228 *start = NULL;
2229 *adjacency = NULL;
2230 *vweights = NULL;
2231 *eweights = NULL;
2232
2233 error_flag = 0;
2234 line_num = 0;
2235
2236 /* Read any leading comment lines */
2237 end_flag = 1;
2238 while (end_flag == 1) {
2239 *nvtxs = chaco_read_int(fin, &end_flag);
2240 ++line_num;
2241 }
2242 if (*nvtxs <= 0) {
2243 printf("ERROR in graph file `%s':", inname);
2244 printf(" Invalid number of vertices (%d).\n", *nvtxs);
2245 fclose(fin);
2246 return(1);
2247 }
2248
2249 narcs = chaco_read_int(fin, &end_flag);
2250 if (narcs < 0) {
2251 printf("ERROR in graph file `%s':", inname);
2252 printf(" Invalid number of expected edges (%d).\n", narcs);
2253 fclose(fin);
2254 return(1);
2255 }
2256
2257 /* Check if vertex or edge weights are used */
2258 if (!end_flag) {
2259 option = chaco_read_int(fin, &end_flag);
2260 }
2261 using_ewgts = option - 10 * (option / 10);
2262 option /= 10;
2263 using_vwgts = option - 10 * (option / 10);
2264 option /= 10;
2265 vtxnums = option - 10 * (option / 10);
2266
2267 /* Get weight info from Chaco option */
2268 (*nVwgts) = using_vwgts;
2269 (*nEwgts) = using_ewgts;
2270
2271 /* Read numbers of weights if they are specified separately */
2272 if (!end_flag && using_vwgts==1){
2273 j = chaco_read_int(fin, &end_flag);
2274 if (!end_flag) (*nVwgts) = j;
2275 }
2276 if (!end_flag && using_ewgts==1){
2277 j = chaco_read_int(fin, &end_flag);
2278 if (!end_flag) (*nEwgts) = j;
2279 }
2280
2281 /* Discard rest of line */
2282 while (!end_flag)
2283 j = chaco_read_int(fin, &end_flag);
2284
2285 /* Allocate space for rows and columns. */
2286 *start = (int *) malloc((unsigned) (*nvtxs + 1) * sizeof(int));
2287 if (narcs != 0)
2288 *adjacency = (int *) malloc((unsigned) (2 * narcs + 1) * sizeof(int));
2289 else
2290 *adjacency = NULL;
2291
2292 if (using_vwgts)
2293 *vweights = (float *) malloc((unsigned) (*nvtxs) * (*nVwgts) * sizeof(float));
2294 else
2295 *vweights = NULL;
2296
2297 if (using_ewgts)
2298 *eweights = (float *)
2299 malloc((unsigned) (2 * narcs + 1) * (*nEwgts) * sizeof(float));
2300 else
2301 *eweights = NULL;
2302
2303 adjptr = *adjacency;
2304 ewptr = *eweights;
2305
2306 sum_edges = 0;
2307 nedges = 0;
2308 (*start)[0] = 0;
2309 vertex = 0;
2310 vtx = 0;
2311 new_vertex = 1;
2312 while ((using_vwgts || vtxnums || narcs) && end_flag != -1) {
2313 ++line_num;
2314
2315 /* If multiple input lines per vertex, read vertex number. */
2316 if (vtxnums) {
2317 j = chaco_read_int(fin, &end_flag);
2318 if (end_flag) {
2319 if (vertex == *nvtxs)
2320 break;
2321 printf("ERROR in graph file `%s':", inname);
2322 printf(" no vertex number in line %d.\n", line_num);
2323 fclose(fin);
2324 return (1);
2325 }
2326 if (j != vertex && j != vertex + 1) {
2327 printf("ERROR in graph file `%s':", inname);
2328 printf(" out-of-order vertex number in line %d.\n", line_num);
2329 fclose(fin);
2330 return (1);
2331 }
2332 if (j != vertex) {
2333 new_vertex = 1;
2334 vertex = j;
2335 }
2336 else
2337 new_vertex = 0;
2338 }
2339 else
2340 vertex = ++vtx;
2341
2342 if (vertex > *nvtxs)
2343 break;
2344
2345 /* If vertices are weighted, read vertex weight. */
2346 if (using_vwgts && new_vertex) {
2347 for (j=0; j<(*nVwgts); j++){
2348 weight = chaco_read_val(fin, &end_flag);
2349 if (end_flag) {
2350 printf("ERROR in graph file `%s':", inname);
2351 printf(" not enough weights for vertex %d.\n", vertex);
2352 fclose(fin);
2353 return (1);
2354 }
2355 (*vweights)[(vertex-1)*(*nVwgts)+j] = weight;
2356 }
2357 }
2358
2359 nedge = 0;
2360
2361 /* Read number of adjacent vertex. */
2362 neighbor = chaco_read_int(fin, &end_flag);
2363
2364 while (!end_flag) {
2365 skip_flag = 0;
2366
2367 if (using_ewgts) { /* Read edge weight if it's being input. */
2368 for (j=0; j<(*nEwgts); j++){
2369 eweight = chaco_read_val(fin, &end_flag);
2370
2371 if (end_flag) {
2372 printf("ERROR in graph file `%s':", inname);
2373 printf(" not enough weights for edge (%d,%d).\n", vertex, neighbor);
2374 fclose(fin);
2375 return (1);
2376 }
2377
2378 else {
2379 *ewptr++ = eweight;
2380 }
2381 }
2382 }
2383
2384 /* Add edge to data structure. */
2385 if (!skip_flag) {
2386 if (++nedges > 2*narcs) {
2387 printf("ERROR in graph file `%s':", inname);
2388 printf(" at least %d adjacencies entered, but nedges = %d\n",
2389 nedges, narcs);
2390 fclose(fin);
2391 return (1);
2392 }
2393 *adjptr++ = neighbor;
2394 nedge++;
2395 }
2396
2397 /* Read number of next adjacent vertex. */
2398 neighbor = chaco_read_int(fin, &end_flag);
2399 }
2400
2401 sum_edges += nedge;
2402 (*start)[vertex] = sum_edges;
2403 }
2404
2405 /* Make sure there's nothing else in file. */
2406 flag = 0;
2407 while (!flag && end_flag != -1) {
2408 chaco_read_int(fin, &end_flag);
2409 if (!end_flag)
2410 flag = 1;
2411 }
2412
2413 (*start)[*nvtxs] = sum_edges;
2414
2415 if (vertex != 0) { /* Normal file was read. */
2416 if (narcs) {
2417 }
2418 else { /* no edges, but did have vertex weights or vertex numbers */
2419 free(*start);
2420 *start = NULL;
2421 if (*adjacency != NULL)
2422 free(*adjacency);
2423 *adjacency = NULL;
2424 if (*eweights != NULL)
2425 free(*eweights);
2426 *eweights = NULL;
2427 }
2428 }
2429
2430 else {
2431 /* Graph was empty */
2432 free(*start);
2433 if (*adjacency != NULL)
2434 free(*adjacency);
2435 if (*vweights != NULL)
2436 free(*vweights);
2437 if (*eweights != NULL)
2438 free(*eweights);
2439 *start = NULL;
2440 *adjacency = NULL;
2441 }
2442
2443 fclose(fin);
2444
2445 return (error_flag);
2446}
2447
2448
2449int UserInputForTests::chaco_input_geom(
2450 FILE *fingeom, /* geometry input file */
2451 const char *geomname, /* name of geometry file */
2452 int nvtxs, /* number of coordinates to read */
2453 int *igeom, /* dimensionality of geometry */
2454 double **x, /* coordinates of vertices */
2455 double **y,
2456 double **z
2457)
2458{
2459 double xc, yc, zc =0; /* first x, y, z coordinate */
2460 int nread; /* number of lines of coordinates read */
2461 int flag; /* any bad data at end of file? */
2462 int line_num; /* counts input lines in file */
2463 int end_flag; /* return conditional */
2464 int ndims; /* number of values in an input line */
2465 int i=0; /* loop counter */
2466
2467 *x = *y = *z = NULL;
2468 line_num = 0;
2469 end_flag = 1;
2470 while (end_flag == 1) {
2471 xc = chaco_read_val(fingeom, &end_flag);
2472 ++line_num;
2473 }
2474
2475 if (end_flag == -1) {
2476 printf("No values found in geometry file `%s'\n", geomname);
2477 fclose(fingeom);
2478 return (1);
2479 }
2480
2481 ndims = 1;
2482 yc = chaco_read_val(fingeom, &end_flag);
2483 if (end_flag == 0) {
2484 ndims = 2;
2485 zc = chaco_read_val(fingeom, &end_flag);
2486 if (end_flag == 0) {
2487 ndims = 3;
2488 chaco_read_val(fingeom, &end_flag);
2489 if (!end_flag) {
2490 printf("Too many values on input line of geometry file `%s'\n",
2491 geomname);
2492
2493 printf(" Maximum dimensionality is 3\n");
2494 fclose(fingeom);
2495 return (1);
2496 }
2497 }
2498 }
2499
2500 *igeom = ndims;
2501
2502 *x = (double *) malloc((unsigned) nvtxs * sizeof(double));
2503 (*x)[0] = xc;
2504 if (ndims > 1) {
2505 *y = (double *) malloc((unsigned) nvtxs * sizeof(double));
2506 (*y)[0] = yc;
2507 }
2508 if (ndims > 2) {
2509 *z = (double *) malloc((unsigned) nvtxs * sizeof(double));
2510 (*z)[0] = zc;
2511 }
2512
2513 for (nread = 1; nread < nvtxs; nread++) {
2514 ++line_num;
2515 if (ndims == 1) {
2516 i = fscanf(fingeom, "%lf", &((*x)[nread]));
2517 }
2518 else if (ndims == 2) {
2519 i = fscanf(fingeom, "%lf%lf", &((*x)[nread]), &((*y)[nread]));
2520 }
2521 else if (ndims == 3) {
2522 i = fscanf(fingeom, "%lf%lf%lf", &((*x)[nread]), &((*y)[nread]),
2523 &((*z)[nread]));
2524 }
2525
2526 if (i == EOF) {
2527 printf("Too few lines of values in geometry file; nvtxs=%d, but only %d read\n",
2528 nvtxs, nread);
2529 fclose(fingeom);
2530 return (1);
2531 }
2532 else if (i != ndims) {
2533 printf("Wrong number of values in line %d of geometry file `%s'\n",
2534 line_num, geomname);
2535 fclose(fingeom);
2536 return (1);
2537 }
2538 }
2539
2540 /* Check for spurious extra stuff in file. */
2541 flag = 0;
2542 end_flag = 0;
2543 while (!flag && end_flag != -1) {
2544 chaco_read_val(fingeom, &end_flag);
2545 if (!end_flag)
2546 flag = 1;
2547 }
2548
2549 fclose(fingeom);
2550
2551 return (0);
2552}
2553
2554// Chaco input assignments from filename.assign; copied from Zoltan
2555
2556int UserInputForTests::chaco_input_assign(
2557 FILE *finassign, /* input assignment file */
2558 const char *inassignname, /* name of input assignment file */
2559 int nvtxs, /* number of vertices to output */
2560 short *assignment) /* values to be printed */
2561{
2562 int flag; /* logical conditional */
2563 int end_flag; /* return flag from read_int() */
2564 int i, j; /* loop counter */
2565
2566 /* Get the assignment vector one line at a time, checking as you go. */
2567 /* First read past any comments at top. */
2568 end_flag = 1;
2569 while (end_flag == 1) {
2570 assignment[0] = chaco_read_int(finassign, &end_flag);
2571 }
2572
2573 if (assignment[0] < 0) {
2574 printf("ERROR: Entry %d in assignment file `%s' less than zero (%d)\n",
2575 1, inassignname, assignment[0]);
2576 fclose(finassign);
2577 return (1);
2578 }
2579
2580 if (end_flag == -1) {
2581 printf("ERROR: No values found in assignment file `%s'\n", inassignname);
2582 fclose(finassign);
2583 return (1);
2584 }
2585
2586 flag = 0;
2587 int np = this->tcomm_->getSize();
2588 if (assignment[0] >= np) flag = assignment[0];
2589 srand(this->tcomm_->getRank());
2590 for (i = 1; i < nvtxs; i++) {
2591 j = fscanf(finassign, "%hd", &(assignment[i]));
2592 if (j != 1) {
2593 printf("ERROR: Too few values in assignment file `%s'.\n", inassignname);
2594 fclose(finassign);
2595 return (1);
2596 }
2597 if (assignment[i] < 0) {
2598 printf("ERROR: Entry %d in assignment file `%s' less than zero (%d)\n",
2599 i+1, inassignname, assignment[i]);
2600 fclose(finassign);
2601 return (1);
2602 }
2603 if (assignment[i] >= np) { // warn since perhaps an error -- initial part
2604 // assignment is greater than number of processors
2605 if (assignment[i] > flag)
2606 flag = assignment[i];
2607 assignment[i] = rand() % np; // randomly assign vtx to a proc in this case
2608 }
2609 }
2610 srand(rand());
2611
2612 if (flag) {
2613 printf("WARNING: Possible error in assignment file `%s'\n",
2614 inassignname);
2615 printf(" Max assignment set (%d) greater than "
2616 "max processor rank (%d)\n", flag, np-1);
2617 printf(" Some vertices given random initial assignments\n");
2618 }
2619
2620 /* Check for spurious extra stuff in file. */
2621 flag = 0;
2622 end_flag = 0;
2623 while (!flag && end_flag != -1) {
2624 chaco_read_int(finassign, &end_flag);
2625 if (!end_flag)
2626 flag = 1;
2627 }
2628 if (flag) {
2629 printf("WARNING: Possible error in assignment file `%s'\n", inassignname);
2630 printf(" Numerical data found after expected end of file\n");
2631 }
2632
2633 fclose(finassign);
2634 return (0);
2635}
2636
2637// Pamgen Reader
2638void UserInputForTests::readPamgenMeshFile(string path, string testData)
2639{
2640#ifdef HAVE_ZOLTAN2_PAMGEN
2641 int rank = this->tcomm_->getRank();
2642 if (verbose_ && tcomm_->getRank() == 0)
2643 std::cout << "UserInputForTestsBD::readPamgenFile, Read: " << testData << std::endl;
2644
2645 size_t len;
2646 std::fstream file;
2647 int dimension;
2648 if (rank == 0){
2649 // set file name
2650 std::ostringstream meshFileName;
2651 meshFileName << path << "/" << testData << ".pmgen";
2652 // open file
2653
2654 file.open(meshFileName.str(), std::ios::in);
2655
2656 if(!file.is_open()) // may be a problem with path or filename
2657 {
2658 if(verbose_ && tcomm_->getRank() == 0)
2659 {
2660 std::cout << "Unable to open pamgen mesh: ";
2661 std::cout << meshFileName.str();
2662 std::cout <<"\nPlease check file path and name." << std::endl;
2663 }
2664 len = 0; // broadcaset 0 length ->will cause exit
2665 }else{
2666 // write to character array
2667 // get size of file
2668 file.seekg (0,file.end);
2669 len = file.tellg();
2670 file.seekg (0);
2671
2672 // get dimension
2673 dimension = 2;
2674 std::string line;
2675 while(std::getline(file,line))
2676 {
2677 if( line.find("nz") != std::string::npos ||
2678 line.find("nphi") != std::string::npos)
2679 {
2680 dimension = 3;
2681 break;
2682 }
2683 }
2684
2685 file.clear();
2686 file.seekg(0, std::ios::beg);
2687 }
2688 }
2689
2690 // broadcast the file size
2691 this->tcomm_->broadcast(0,sizeof(int), (char *)&dimension);
2692 this->tcomm_->broadcast(0,sizeof(size_t),(char *)&len);
2693 this->tcomm_->barrier();
2694
2695 if(len == 0){
2696 if(verbose_ && tcomm_->getRank() == 0)
2697 std::cout << "Pamgen Mesh file size == 0, exiting UserInputForTests early." << std::endl;
2698 return;
2699 }
2700
2701 char * file_data = new char[len+1];
2702 file_data[len] = '\0'; // critical to null terminate buffer
2703 if(rank == 0){
2704 file.read(file_data,len); // if proc 0 then read file
2705 }
2706
2707 // broadcast the file to the world
2708 this->tcomm_->broadcast(0,(int)len+1,file_data);
2709 this->tcomm_->barrier();
2710
2711 // Create the PamgenMesh
2712
2713 this->pamgen_mesh = rcp(new PamgenMesh);
2714 this->havePamgenMesh = true;
2715 pamgen_mesh->createMesh(file_data,dimension,this->tcomm_);
2716
2717 // save mesh info
2718 pamgen_mesh->storeMesh();
2719 this->tcomm_->barrier();
2720
2721 // set coordinates
2722 this->setPamgenCoordinateMV();
2723
2724 // set adjacency graph
2725 this->setPamgenAdjacencyGraph();
2726
2727 this->tcomm_->barrier();
2728 if(rank == 0) file.close();
2729 delete [] file_data;
2730#else
2731 throw std::runtime_error("Pamgen requested but Trilinos "
2732 "not built with Pamgen");
2733#endif
2734}
2735
2736#ifdef HAVE_ZOLTAN2_PAMGEN
2737void UserInputForTests::setPamgenCoordinateMV()
2738{
2739 int dimension = pamgen_mesh->num_dim;
2740 // get coordinate and point info;
2741// zlno_t numLocalPoints = pamgen_mesh->num_nodes;
2742// zgno_t numGlobalPoints = pamgen_mesh->num_nodes_global;
2743 zgno_t numelements = pamgen_mesh->num_elem;
2744 zgno_t numGlobalElements = pamgen_mesh->num_elems_global;
2745 // allocate and set an array of coordinate arrays
2746 zscalar_t **elem_coords = new zscalar_t * [dimension];
2747 for(int i = 0; i < dimension; ++i){
2748 elem_coords[i] = new zscalar_t[numelements];
2749 double *tmp = &pamgen_mesh->element_coord[i*numelements];
2750 for (int j = 0; j < numelements; j++) elem_coords[i][j] = tmp[j];
2751 }
2752
2753 // make a Tpetra map
2754 RCP<Tpetra::Map<zlno_t, zgno_t, znode_t> > mp;
2755 // mp = rcp(new map_t(numGlobalElements, numelements, 0, this->tcomm_)); // constructo 1
2756
2757// Array<zgno_t>::size_type numEltsPerProc = numelements;
2758 Array<zgno_t> elementList(numelements);
2759 for (Array<zgno_t>::size_type k = 0; k < numelements; ++k) {
2760 elementList[k] = pamgen_mesh->element_order_map[k];
2761 }
2762
2763 mp = rcp (new map_t (numGlobalElements, elementList, 0, this->tcomm_)); // constructor 2
2764
2765
2766 // make an array of array views containing the coordinate data
2767 Teuchos::Array<Teuchos::ArrayView<const zscalar_t> > coordView(dimension);
2768 for (int i = 0; i < dimension; i++){
2769 if(numelements > 0){
2770 Teuchos::ArrayView<const zscalar_t> a(elem_coords[i], numelements);
2771 coordView[i] = a;
2772 }
2773 else {
2774 Teuchos::ArrayView<const zscalar_t> a;
2775 coordView[i] = a;
2776 }
2777 }
2778
2779 // set the xyz_ multivector
2780 xyz_ = RCP<tMVector_t>(new
2781 tMVector_t(mp, coordView.view(0, dimension),
2782 dimension));
2783}
2784
2785void UserInputForTests::setPamgenAdjacencyGraph()
2786{
2787// int rank = this->tcomm_->getRank();
2788// if(rank == 0) std::cout << "Making a graph from our pamgen mesh...." << std::endl;
2789
2790 // Define Types
2791// typedef zlno_t lno_t;
2792// typedef zgno_t gno_t;
2793
2794 // get info for setting up map
2795 size_t local_nodes = (size_t)this->pamgen_mesh->num_nodes;
2796 size_t local_els = (size_t)this->pamgen_mesh->num_elem;
2797
2798 size_t global_els = (size_t)this->pamgen_mesh->num_elems_global; // global rows
2799 size_t global_nodes = (size_t)this->pamgen_mesh->num_nodes_global; //global columns
2800 // make map with global elements assigned to this mesh
2801 // make range map
2802// if(rank == 0) std::cout << "Building Rowmap: " << std::endl;
2803 RCP<const map_t> rowMap = rcp(new map_t(global_els,0,this->tcomm_));
2804 RCP<const map_t> rangeMap = rowMap;
2805
2806 // make domain map
2807 RCP<const map_t> domainMap = rcp(new map_t(global_nodes,0,this->tcomm_));
2808
2809 // Get max number of nodes per element
2810 int blks = this->pamgen_mesh->num_elem_blk;
2811 int max_nodes_per_el = 0;
2812 for(int i = 0; i < blks; i++)
2813 if (this->pamgen_mesh->nodes_per_element[i] > max_nodes_per_el)
2814 max_nodes_per_el = this->pamgen_mesh->nodes_per_element[i];
2815
2816 // make the element-node adjacency matrix
2817 Teuchos::RCP<tcrsMatrix_t> C = rcp(new tcrsMatrix_t(rowMap,max_nodes_per_el));
2818
2819 Array<zgno_t> g_el_ids(local_els);
2820 for (size_t k = 0; k < local_els; ++k) {
2821 g_el_ids[k] = pamgen_mesh->global_element_numbers[k]-1;
2822 }
2823
2824 Array<zgno_t> g_node_ids(local_nodes);
2825 for (size_t k = 0; k < local_nodes; ++k) {
2826 g_node_ids[k] = pamgen_mesh->global_node_numbers[k]-1;
2827 }
2828
2829 zlno_t el_no = 0;
2830 zscalar_t one = static_cast<zscalar_t>(1);
2831
2832// if(rank == 0) std::cout << "Writing C... " << std::endl;
2833 for(int i = 0; i < blks; i++)
2834 {
2835 int el_per_block = this->pamgen_mesh->elements[i];
2836 int nodes_per_el = this->pamgen_mesh->nodes_per_element[i];
2837 int * connect = this->pamgen_mesh->elmt_node_linkage[i];
2838
2839 for(int j = 0; j < el_per_block; j++)
2840 {
2841 const zgno_t gid = static_cast<zgno_t>(g_el_ids[el_no]);
2842 for(int k = 0; k < nodes_per_el; k++)
2843 {
2844 int g_node_i = g_node_ids[connect[j*nodes_per_el+k]-1];
2845 C->insertGlobalValues(gid,
2846 Teuchos::tuple<zgno_t>(g_node_i),
2847 Teuchos::tuple<zscalar_t>(one));
2848 }
2849 el_no++;
2850 }
2851 }
2852
2853 C->fillComplete(domainMap, rangeMap);
2854
2855
2856 // Compute product C*C'
2857// if(rank == 0) std::cout << "Compute Multiplication C... " << std::endl;
2858 RCP<tcrsMatrix_t> A = rcp(new tcrsMatrix_t(rowMap,0));
2859 Tpetra::MatrixMatrix::Multiply(*C, false, *C, true, *A);
2860
2861 // remove entries not adjacent
2862 // make graph
2863// if(rank == 0) std::cout << "Writing M_... " << std::endl;
2864 this->M_ = rcp(new tcrsMatrix_t(rowMap, A->getGlobalMaxNumRowEntries()));
2865
2866// if(rank == 0) std::cout << "\nSetting graph of connectivity..." << std::endl;
2867 Teuchos::ArrayView<const zgno_t> rowMapElementList =
2868 rowMap->getLocalElementList();
2869 for (Teuchos_Ordinal ii = 0; ii < rowMapElementList.size(); ii++)
2870 {
2871 zgno_t gid = rowMapElementList[ii];
2872 size_t numEntriesInRow = A->getNumEntriesInGlobalRow (gid);
2873 typename tcrsMatrix_t::nonconst_global_inds_host_view_type rowinds("Indices", numEntriesInRow);
2874 typename tcrsMatrix_t::nonconst_values_host_view_type rowvals("Values", numEntriesInRow);
2875
2876 // modified
2877 Array<zscalar_t> mod_rowvals;
2878 Array<zgno_t> mod_rowinds;
2879 A->getGlobalRowCopy (gid, rowinds, rowvals, numEntriesInRow);
2880 for (size_t i = 0; i < numEntriesInRow; i++) {
2881// if (rowvals[i] == 2*(this->pamgen_mesh->num_dim-1))
2882// {
2883 if (rowvals[i] >= 1)
2884 {
2885 mod_rowvals.push_back(one);
2886 mod_rowinds.push_back(rowinds[i]);
2887 }
2888 }
2889 this->M_->insertGlobalValues(gid, mod_rowinds, mod_rowvals);
2890 }
2891
2892 this->M_->fillComplete();
2894 // if(rank == 0) std::cout << "Completed M... " << std::endl;
2895
2896}
2897
2898#endif
2899
2900#endif
#define TEST_FAIL_AND_THROW(comm, ok, s)
const char param_comment
USERINPUT_FILE_FORMATS
A class that generates typical user input for testing.
@ MATRIX_MARKET
#define Z2_FORWARD_EXCEPTIONS
Forward an exception back through call stack.
common code used by tests
float zscalar_t
Tpetra::Map ::local_ordinal_type zlno_t
Tpetra::Map ::global_ordinal_type zgno_t
keep typedefs that commonly appear in many places localized
Traits of Xpetra classes, including migration method.
map_t::node_type default_znode_t
RCP< xcrsGraph_t > getUIXpetraCrsGraph()
Tpetra::Map< zlno_t, zgno_t, znode_t > map_t
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.
Tpetra::Export< zlno_t, zgno_t, znode_t > export_t
RCP< tVector_t > getUITpetraVector()
Tpetra::Import< zlno_t, zgno_t, znode_t > import_t
RCP< tMVector_t > getUIEdgeWeights()
RCP< xMVector_t > getUIXpetraMultiVector(int nvec)
bool hasInputDataType(const string &input_type)
RCP< tMVector_t > getUITpetraMultiVector(int nvec)
RCP< tMVector_t > getUICoordinates()
RCP< tcrsGraph_t > getUITpetraCrsGraph()
RCP< xVector_t > getUIXpetraVector()
RCP< xcrsMatrix_t > getUIXpetraCrsMatrix()
RCP< tcrsMatrix_t > getUITpetraCrsMatrix()
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.
RCP< tMVector_t > getUIWeights()
static const std::string fail
Tpetra::Map map_t
GeometricGen::GeometricGenerator< zscalar_t, zlno_t, zgno_t, znode_t > geometricgen_t
Tpetra::CrsMatrix< zscalar_t, zlno_t, zgno_t, znode_t > tcrsMatrix_t
Tpetra::MultiVector< zscalar_t, zlno_t, zgno_t, znode_t > tMVector_t
Tpetra::Vector< zscalar_t, zlno_t, zgno_t, znode_t > tVector_t
fname
Begin.
Definition validXML.py:19
static RCP< tMVector_t > coordinates
static RCP< User > convertToXpetra(const RCP< User > &a)
Convert the object to its Xpetra wrapped version.
Tpetra::MultiVector< zscalar_t, zlno_t, zgno_t, znode_t > tMVector_t