[Trilinos-Users] Tpetra::VbrMatrix interface

C.S. Natarajan csnataraj at gmail.com
Sun Aug 3 20:49:44 MDT 2014


Hello,
        I am wondering if the Tpetra::Vbrmatrix interface is in regular
use. I would like to use that interface, however, seem to be running into
issues. I had brought this issue up last year but didn't get around to
filing a bug report.

 In the example below, each proc inserts a value of 1 into the global row
and column given by the variable *grow*. So, if the number of processors is
np, the value at the global row and column index grow should be equal to
np. While the crsMatrix interface produces the correct value, the Vbr
interface doesn't.

I would appreciate any input into what I might be doing wrong or whether
there really is a bug in the interface.

TIA,
C.S.N



#include <sstream>

//Teuchos Includes
#include "Teuchos_Array.hpp"
#include "Teuchos_ArrayView.hpp"
#include "Teuchos_ArrayViewDecl.hpp"
#include "Teuchos_BLAS.hpp"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_DefaultMpiComm.hpp"
#include "Teuchos_FancyOStream.hpp"
#include "Teuchos_GlobalMPISession.hpp"
#include "Teuchos_oblackholestream.hpp"
#include "Teuchos_OrdinalTraits.hpp"
#include "Teuchos_RCP.hpp"
#include "Teuchos_ScalarTraits.hpp"
#include "Teuchos_SerialDenseMatrix.hpp"
#include "Teuchos_ArrayRCPDecl.hpp"

//Tpetra Includea
#include "Tpetra_DefaultPlatform.hpp"
#include "Tpetra_Version.hpp"
#include "Tpetra_BlockMap_decl.hpp"
#include "Tpetra_BlockMap.hpp"
#include "Tpetra_Vector.hpp"
#include "Tpetra_VbrMatrix.hpp"
#include "Tpetra_BlockMultiVector.hpp"
#include "Tpetra_BlockMultiVector_decl.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Tpetra_CrsMatrix.hpp"

typedef Tpetra::DefaultPlatform::DefaultPlatformType            Platform;
typedef Tpetra::DefaultPlatform::DefaultPlatformType::NodeType  Node;


int main(int argc, char *argv[])
{

        using Teuchos::Array;
        typedef int global_size_t;
        typedef double Scalar;
        typedef int Ordinal;

        Teuchos::oblackholestream nullStream;
        Teuchos::GlobalMPISession mpiSession(&argc,&argv, NULL);

        Teuchos::RCP<const Teuchos::Comm<int> > comm =
Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
        Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
        Teuchos::RCP<Node> node = platform.getNode();

        size_t myRank = comm->getRank();
        size_t nProc = comm->getSize();
        std::ostream &out = (myRank == 0)?std::cout:nullStream;
        Teuchos::CommandLineProcessor clp;
        int nx = 64;
        clp.setOption("nx", &nx, "Number of nodes in X-direction.");
        int ny = 64;
        clp.setOption("ny", &ny, "Number of nodes in Y-direction.");
        int grow=12;
        clp.setOption("grow", &grow, "Global row to output");

        clp.recogniseAllOptions(true);
        clp.throwExceptions(false);

        Teuchos::CommandLineProcessor::EParseCommandLineReturn parseReturn
= clp.parse(argc, argv);
        if(parseReturn ==
Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED) return 0;
        if(parseReturn != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
return 1;


        const int blockSize = 2;

        Ordinal indexBase = 0;
        Teuchos::RCP<const Tpetra::BlockMap<Ordinal,Ordinal,Node> > tbmap =
Teuchos::rcp(new Tpetra::BlockMap<Ordinal,Ordinal,Node>(nx*ny, blockSize,
indexBase, comm));

        Teuchos::RCP<const Tpetra::Map<Ordinal> > tmap = Teuchos::rcp(new
Tpetra::Map<Ordinal>(nx*ny,indexBase, comm));

        Teuchos::RCP<Tpetra::VbrMatrix<Scalar, Ordinal> >A =
Teuchos::rcp(new Tpetra::VbrMatrix<Scalar, Ordinal>(tbmap, 9));

        Teuchos::RCP<Tpetra::CrsMatrix<Scalar, Ordinal> >Acrs =
Teuchos::rcp(new Tpetra::CrsMatrix<Scalar, Ordinal>(tmap, 9));

        Teuchos::SerialDenseMatrix <Ordinal, Scalar> S(2,2);

        Teuchos::Array<Scalar> value(1,0);

        Teuchos::Array<Ordinal> tcol(1,0);

        S(0,0) = value[0] = 1;
        S(1,1) = 1;
        long row, col;
        row = col = tcol[0] = grow;

        A->sumIntoGlobalBlockEntry(row, col, S);
        Acrs->insertGlobalValues(row, tcol (), value());
        A->globalAssemble();
        Acrs->globalAssemble();
        Acrs->fillComplete();

        Teuchos::Array< Teuchos::ArrayRCP< const Scalar > > blockEntries;
        Ordinal numptrows;
        Teuchos::ArrayView<const Ordinal> bcolv;
        Teuchos::Array<Ordinal> indices, pcolind;
        Teuchos::Array<Scalar> crsvals;
        size_t NumEntries = 0;

        A->getGlobalBlockRowView (grow, numptrows, bcolv, pcolind,
blockEntries);
        A->fillComplete();

        out << "\nOutputting global block row view " << std::endl;

        if (blockEntries.size() != 0)
        {
                std::cout << "block entries size is " <<
blockEntries.size() << std::endl;
                for(size_t i=0; i < blockEntries.size(); ++i)
                        for(int j = 0; j < blockEntries[i].size(); j++)
         std::cout << blockEntries[i][j] << std::endl;
        }

        NumEntries = Acrs->getNumEntriesInGlobalRow(grow);
        if (NumEntries != Teuchos::OrdinalTraits<size_t>::invalid())
        {
                std::cout << "outputting crs values " << std::endl;
                indices.resize(NumEntries); crsvals.resize(NumEntries);
                Acrs->getGlobalRowCopy (grow, indices(), crsvals(),
NumEntries);
                for(int i = 0; i < indices.size(); i++)
                        std::cout << indices[i] << " " << crsvals[i] <<
"\n";
        }
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://software.sandia.gov/pipermail/trilinos-users/attachments/20140803/f65c514c/attachment.html>


More information about the Trilinos-Users mailing list