[Trilinos-Users] extracting sub-matrix

[Bart Janssens]

OK, attached is a slightly improved version of your program. It's
better to get the GIDs from the RangeMap. Other than that, your code
was OK, the crash was actually in the matrix gallery destructor. Also,
WATCH doesn't print the matrix in the case of a block, apparently, so
I'm using Thyra to print out the values.

Cheers,

Bart

On Sat, Mar 2, 2013 at 11:43 AM, Hoang Giang Bui <hgbk2008 at gmail.com> wrote:
>
>
> 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



-- 
Bart
-------------- 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"

#include <Thyra_EpetraLinearOp.hpp>
#include <Teuchos_VerboseObject.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);
	Gal.Set("problem_type","hb");
	Gal.Set("matrix_name","bcsstk01.rsa");
	Gal.Set("map_type","greedy");

	Gal.Set(CLP);

  Teuchos::RCP<Epetra_CrsMatrix> Mat = Teuchos::rcp(Gal.GetMatrix(),false);
	
  WATCH(*Mat);

  const Epetra_Map & range_map = Mat->OperatorRangeMap();
  const int num_elems = range_map.NumMyElements();
  int * matrix_global_ids = range_map.MyGlobalElements();
  
    std::vector<int> gids1;
    gids1.push_back(matrix_global_ids[0]);
    gids1.push_back(matrix_global_ids[1]);
    gids1.push_back(matrix_global_ids[2]);
    gids1.push_back(matrix_global_ids[3]);
    std::vector<int> gids2;
    for(int i = 4; i != num_elems; i++)
        gids2.push_back(matrix_global_ids[i]);

    std::vector<std::vector<int> > vars;
    vars.push_back(gids1);
    vars.push_back(gids2);
    
    Teko::Epetra::BlockedEpetraOperator blkExtractor(vars,Mat,"sub-matrix");

	Teuchos::RCP<const Epetra_Operator> blk1 = blkExtractor.GetBlock(0,0);

  
  Teuchos::RCP<Teuchos::FancyOStream> fancy_out = Teuchos::VerboseObjectBase::getDefaultOStream();
	Thyra::describeLinearOp(*Thyra::epetraLinearOp(blk1), *fancy_out, Teuchos::VERB_EXTREME);
	
	Comm.Barrier();
	
#ifdef HAVE_MPI
    MPI_Finalize();
#endif

	return 0;
}