[Trilinos-Users] Tpetra::VbrMatrix interface
C.S. Natarajan
csnataraj at gmail.com
Mon Aug 4 22:31:11 MDT 2014
Alan, Mark,
Thank you very much for your responses. I do have constant
block sizes, however, for now I think I will just use the Crs interface.
Just curious, if you are doing in-house development (FEM), what is your
preferred methodology for problems with multiple DOF per node (Problems
span the range from electrostatics to fluid flow)?
I am going with the Tpetra, Intrepid (+shards), Muelu, Belos stack; is this
a reasonable choice? Also, does the Crs interface lend itself towards
block-preconditioning?
Thanks,
C.S.N
On Mon, Aug 4, 2014 at 4:02 PM, Hoemmen, Mark <mhoemme at sandia.gov> wrote:
> Greetings, and thanks for the VbrMatrix test code! There are some known
> issues with VbrMatrix. It should be considered broken until we have a bit
> more time to work on it. We do have a preliminary constant-block-size
> interface, Tpetra::Experimental::BlockCrsMatrix (NOT the same as
> Tpetra::BlockCrsMatrix), which optimizes for the common case where all
> blocks have the same dimensions. If you would like to experiment with
> this, please let me know and I'll show you how to use it. Thanks!
>
> mfh
>
> On 8/4/14, 12:00 PM, "trilinos-users-request at software.sandia.gov"
> <trilinos-users-request at software.sandia.gov> wrote:
>
> >Send Trilinos-Users mailing list submissions to
> > trilinos-users at software.sandia.gov
> >
> >To subscribe or unsubscribe via the World Wide Web, visit
> > https://software.sandia.gov/mailman/listinfo/trilinos-users
> >or, via email, send a message with subject or body 'help' to
> > trilinos-users-request at software.sandia.gov
> >
> >You can reach the person managing the list at
> > trilinos-users-owner at software.sandia.gov
> >
> >When replying, please edit your Subject line so it is more specific
> >than "Re: Contents of Trilinos-Users digest..."
> >
> >
> >Today's Topics:
> >
> > 1. Tpetra::VbrMatrix interface (C.S. Natarajan)
> > 2. Re: [EXTERNAL] Tpetra::VbrMatrix interface (Williams, Alan B)
> >
> >
> >----------------------------------------------------------------------
> >
> >Message: 1
> >Date: Sun, 3 Aug 2014 21:49:44 -0500
> >From: "C.S. Natarajan" <csnataraj at gmail.com>
> >To: "trilinos-users at software.sandia.gov"
> > <trilinos-users at software.sandia.gov>
> >Subject: [Trilinos-Users] Tpetra::VbrMatrix interface
> >Message-ID:
> > <
> CAO_E+w6YFVdsLoiQamKcv+mYBJRpsCMyeYnNUjdCL3QKdhkrRA at mail.gmail.com>
> >Content-Type: text/plain; charset="utf-8"
> >
> >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-0001.html>
> >
> >------------------------------
> >
> >Message: 2
> >Date: Mon, 4 Aug 2014 03:45:44 +0000
> >From: "Williams, Alan B" <william at sandia.gov>
> >To: "C.S. Natarajan" <csnataraj at gmail.com>,
> > "trilinos-users at software.sandia.gov"
> > <trilinos-users at software.sandia.gov>
> >Subject: Re: [Trilinos-Users] [EXTERNAL] Tpetra::VbrMatrix interface
> >Message-ID: <D0045827.1D524%william at sandia.gov>
> >Content-Type: text/plain; charset="windows-1252"
> >
> >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-0001.html>
> >
> >------------------------------
> >
> >_______________________________________________
> >Trilinos-Users mailing list
> >Trilinos-Users at software.sandia.gov
> >https://software.sandia.gov/mailman/listinfo/trilinos-users
> >
> >
> >End of Trilinos-Users Digest, Vol 108, Issue 2
> >**********************************************
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://software.sandia.gov/pipermail/trilinos-users/attachments/20140804/0b2f0458/attachment.html>
More information about the Trilinos-Users
mailing list