[Dune] Problem in implementing parallel linear solver

[Kumar, Paras]

Hi again,

A small update to the previous email.

If, I use comm_redist instead of DuneComm (which kind of makes sense, 
once the redistribution has happened), it is observed that all indexes 
for parallel_A are owner for each process. Thus, the unsymmmetry also 
happens for rows whose index is Onwer according to the code snippet 
below. C

**************************************************************************
using OwnerSet = DuneCommunication::OwnerSet;
using GlobalLookup = DuneCommunication::GlobalLookupIndexSet;
GlobalLookup global(comm_redist->indexSet());
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;
}
**************************************************************************



On 2017-12-07 17:33, Kumar, Paras wrote:
> 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
> 
> _______________________________________________
> Dune mailing list
> Dune at lists.dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune