[Dune] Problem in implementing parallel linear solver

[Kumar, Paras]

Hi,

On 2017-12-07 17:00, Markus Blatt wrote:
> On Thu, Dec 07, 2017 at 02:49:13PM +0100, Kumar, Paras wrote:
>> >>
>> >>Here, are a few observations which might give some clues about the error
>> >>I
>> >>make, but I'm not sure.
>> >>1. As can be seen in the output file, A is tridiagonal where as the
>> >>distributed are not always symmetric.
>> >
>> >Please check where the unsymmetry happens. This should only happen
>> >rows whose index is not marked as owner.
>> >
>> >using OwnerSet = DuneCommunication::OwnerSet;
>> >using GlobalLookup = DuneCommunication::GlobalLookupIndexSet;
>> >GlobalLookup global(DuneComm.indexSet());
>> >auto pair = global.pair(index);
>> >if(!pair || OwnerSet::contains(pair.local().attribute())
>> >  cerr<<"Houston we have a problem at index "<<index<<" which is
>> >  owned.
>> >
>> Could you please explain how do we get the variable index. It seems to 
>> be
>> the local index (probably wrt *comm_redist in my example code). Also, 
>> do I
>> need to do this check for all dofs in an iterative manner.
>> 
> 
> Well if you have an iterator r over matrix rows, then you should use
> r.index() as index. You can also use the index of a column interator.
> 
> Only check the indices of rows, cols that exhibit unsymmetries.
> 
> Markus

When tried with the complete A matrix, I got the message for every row. 
When tried with the parallel_A matrices, I got a segmentation fault.
**************************************************************************
using OwnerSet = DuneCommunication::OwnerSet;
using GlobalLookup = DuneCommunication::GlobalLookupIndexSet;
GlobalLookup global(DuneComm.indexSet());
for(auto r = app->A.begin(); r !=app->A.end(); ++r){
//for(auto r = parallel_A.begin(); r !=parallel_A.end(); ++r){
	auto pair = global.pair(r.index());
		if(!pair || OwnerSet::contains(pair->local().attribute()))
   			std::cerr<<"Houston we have a problem at index "<<r.index()<<" 
which is owned. from rank "<<app->world_rank <<std::endl;
}
**************************************************************************
here app->A is the full A matrix which is only available at the master 
process while parallel_A is the redistributed part of the matrix on each 
of the processes.

Another query, should we have DuneComm.indexSet() or 
comm_redist->indexSet() [please refer the attached code file for 
details].

> 
> _______________________________________________
> Dune mailing list
> Dune at lists.dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune

With best regards,
Paras Kumar