[Trilinos-Users] [EXTERNAL] Tpetra::VbrMatrix interface

[Williams, Alan B]

Hello,

I wrote the VbrMatrix class, but I don’t think it was ever heavily used. I believe that a BlockCRS matrix capability is either planned or in development, and that’s probably a better bet. I haven’t been involved in tpetra development recently, and so it’s likely that the behavior and semantics of the VbrMatrix class have not kept up with best practices in other tpetra capabilities.

There are some unit-tests for VbrMatrix, in the directory tpetra/test/VbrMatrix. Any usage scenario that isn’t expressed in the unit-tests probably shouldn’t be counted on.

In particular, since the VbrMatrix header states that ‘getGlobalBlockRowView’ is only available prior to calling fillComplete, the pointers or views returned by that method are probably not still good after fillComplete has been called.

One thing you could try is to convert the global-row to a local-row (using the row-map) and then get a local row view after fillComplete has been called. I’m pretty sure that global-row will only convert to a *valid* local-row on 1 processor and *invalid* on the other procs. On the proc where it is valid, check the coefficient value. If that isn’t correct, then there is probably a bug.

Alan




From: "C.S. Natarajan" <csnataraj at gmail.com<mailto:csnataraj at gmail.com>>
Date: Sunday, August 3, 2014 at 8:49 PM
To: trilinos-users <trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>>
Subject: [EXTERNAL] [Trilinos-Users] Tpetra::VbrMatrix interface

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/20140804/fd66d741/attachment.html>