[Trilinos-Users] Fwd: extracting sub-matrix

[Hoang Giang Bui]

Yes, I try to add gids2 but it throws a segmentation fault error. Attach is
the full script. Maybe you take a look again.

BR
Bui



On Sat, Mar 2, 2013 at 9:24 AM, Bart Janssens <bart.janssens at lid.kviv.be>wrote:

> On Sat, Mar 2, 2013 at 6:42 AM, Hoang Giang Bui <hgbk2008 at gmail.com>
> wrote:
> > Thank you very much for the tip. Below is a snippet that I used to test
> the
> > method above, but it does not work as expected. Can you help to show me
> the
> > error in this. Basically, I try to extract the first 4x4 block of the
> > matrix.
> >
> >     Epetra_CrsMatrix * Mat = Gal.GetMatrix();
> >
> >     std::vector<int> gids1;
> >     gids1.push_back(0);
> >     gids1.push_back(1);
> >     gids1.push_back(2);
> >     gids1.push_back(3);
>
> I think the problem here is that you have to cover all indices of the
> matrix. So if your complete matrix is 12x12, you need to add a gids2
> with GIDs 4...11, and then GetBlock(0,0) should give you the first
> 4x4, GetBlock(0,1) the 4x8 and so on.
>
> Cheers,
>
> --
> Bart
>



-- 
With Best Regards !
Giang Bui
To learn and to excel



-- 
With Best Regards !
Giang Bui
To learn and to excel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://software.sandia.gov/pipermail/trilinos-users/attachments/20130302/dfc853f3/attachment.html 
-------------- next part --------------
#define HAVE_MPI

#define WATCH(a)    std::cout << #a << ": " << a << std::endl;

#ifdef HAVE_MPI
#include "mpi.h"
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_CrsMatrix.h"
#include "Epetra_Import.h"
//#include "Trilinos_Util.h"
#include "Trilinos_Util_CommandLineParser.h"
#include "Trilinos_Util_CrsMatrixGallery.h"
#include "Teko_BlockedEpetraOperator.hpp"
#include "Teuchos_RCP.hpp"

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

#ifdef HAVE_MPI
	MPI_Init(&argc,&argv);
	Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
	Epetra_SerialComm Comm;
#endif

	Trilinos_Util::CommandLineParser CLP(argc,argv);
	//CLP.ShowAll();
	CLP.ShowReallyAll();

	Trilinos_Util::CrsMatrixGallery Gal("",Comm,false);
	Gal.Set("problem_type","hb");
	Gal.Set("matrix_name","/media/data1/research/matrix-market/bcsstk01.rsa");
	Gal.Set("map_type","greedy");

	Gal.Set(CLP);

	Epetra_CrsMatrix * Mat = Gal.GetMatrix();
	
    WATCH(*Mat);

    std::vector<int> gids1;
    gids1.push_back(0);
    gids1.push_back(1);
    gids1.push_back(2);
    gids1.push_back(3);
    std::vector<int> gids2;
    for(int i = 4; i < 48; i++)
        gids2.push_back(i);

    std::vector<std::vector<int> > vars;
    vars.push_back(gids1);
    vars.push_back(gids2);
    
    Teuchos::RCP<const Epetra_Operator> matOp = Teuchos::rcp(Mat);
    Teko::Epetra::BlockedEpetraOperator blkExtractor(vars,matOp,"sub-matrix");

	Teuchos::RCP<const Epetra_Operator> blk1 = blkExtractor.GetBlock(0,0);
	
	WATCH(blk1);
	
	Comm.Barrier();
	
#ifdef HAVE_MPI
    MPI_Finalize();
#endif

	return 0;
}