[Trilinos-Users] FECrs Vector/Matrix

[David Neckels]

Hi Alan,

Okay. Using your previous explanation as a guide, I reworked things and 
figured out what I was doing wrong.

After creating the matrix and vectors with the unique map, and calling 
InsertValues, I was calling a version of
SumIntoGlobalValues
that exists only in the base class, Epetra_CrsMatrix:

int SumIntoGlobalValues(int GlobalRow, int NumEntries, double* Values, 
int* Indices);

Obviously, since the function is in the base class it didn't know about 
the extra (shared rows), and so it would return an error.
I now call:
int SumIntoGlobalValues(int numRows, const int* rows,
int numCols, const int* cols,
const double* values,
int format=Epetra_FECrsMatrix::COLUMN_MAJOR);

from Epetra_FECrsMatrix, which works correctly.

Also, after assembly I apply dirichlet boundary conditions by replacing 
the row for a dirichlet dof. I call

int ret = matrix->ExtractGlobalRowView(row_id, nentries, old_row, 
old_indices);

followed by

matrix->ReplaceGlobalValues(1, &row_id, nentries, old_indices, &new_row[0]);

At first this wasn't working correctly for shared dofs, since the 
ExtractView didn't know about the shared dofs on non-owning side (it is 
a function in the base class CrsMatrix).
What seems to work is to call:
matrix->GlobalAssemble(false /* don't fillcomplete */ );
and then extract the row and replace, on the owning processor.....

I don't know how I was getting away with my previous unique/shared map 
approach. It was giving the right answer, but I notice that my matrices 
are much better conditioned now.... A fluke, I guess.

Thanks a lot for the help.

-David N.

Williams, Alan B wrote:
> If you can give me some example code that reproduces the -1
> return-errors you're seeing, I'll figure out what's going on.
>
> Perhaps there is a bug in FECrsMatrix.
>
> Alan
>
>
>   
>> -----Original Message-----
>> From: David Neckels [mailto:dneckels at ucar.edu] 
>> Sent: Wednesday, November 28, 2007 6:07 PM
>> To: Williams, Alan B
>> Cc: trilinos-users at software.sandia.gov
>> Subject: Re: [Trilinos-Users] FECrs Vector/Matrix
>>
>> Hi Alan,
>>
>> Thanks for responding, good to hear from you, hope all is well!
>>
>> I appreciate the basic explanation, I think it helps 
>> understand the theory on using the classes better.
>> However, I still don't seem to have luck with things. I have 
>> now tried:
>>
>> matrix = new Epetra_FECrsMatrix(Copy, *unique_map, max_df);  
>>
>> x = new Epetra_FEVector(*unique_map); //( solving Ax=b) b = 
>> new Epetra_FEVector(*unique_map);
>>
>> Call InsertGlobalValues for the shared set of dofs to build 
>> the matrix stencil.
>>
>> Then,
>> during assembly call SumIntoGlobalValues for the same set of 
>> dofs (the shared and unique dofs).
>> However, for shared dofs this function now returns -1 !!
>>
>> I then close the assembly:
>> matrix->GlobalAssembly();
>> b->GlobalAssembly()
>>
>> This gives the wrong result...
>> I notice that calling the InsertGlobalValues alone does 
>> appear to correctly sum the entries, but it doesn't seem to 
>> play well with SumIntoGlobalValues for some reason.
>>
>> For this example I am not using ReplaceGlobalValues; I have 
>> no constraints...
>>
>> If you like I can try to boil this down to a simple example....??
>>
>> Thanks,
>> David N.
>>
>>
>> // insert rows with all (shared) dofs
>> // sum into global with all (shared) dofs --here is I used 
>> unique_map to create the matrix, above, I would get bad 
>> return values??
>> // Call ExtractGlobalView and replace constrained rows of the 
>> matrix for shared dofs (unique_map in the matrix create 
>> caused this to fail) // Also sum into shared dofs of b
>>
>> // After assembly:
>> matrix->GlobalAssemble(*umap, *umap);
>> b->GlobalAssemble();
>>
>>
>>
>>
>> Williams, Alan B wrote:
>>     
>>> Hi David,
>>>
>>> The "FE" epetra objects are intended to handle the situation you 
>>> describe, where shared nodes exist on multiple processors, 
>>>       
>> so you're 
>>     
>>> on the right track.
>>>
>>> As you know, the non-FE epetra objects tend to be used for the 
>>> uniquely-owned case, where the Epetra_Map attribute has 
>>> non-overlapping entries.
>>>
>>> The general idea is that the FE objects allow a 
>>>       
>> finite-element caller 
>>     
>>> to assemble data on the local processor (i.e. overlapping data for 
>>> shared nodes), which may end up on another processor in the 
>>> non-overlapping distribution.
>>>
>>> The FE objects inherit their non-FE counterparts 
>>>       
>> (FECrsMatrix inherits 
>>     
>>> CrsMatrix, etc), and the Epetra_Map objects that are passed as 
>>> constructor arguments should describe the uniquely-owned, 
>>> non-overlapping distribution (which is what you want the underlying 
>>> CrsMatrix to end up with).
>>>
>>> When you assemble data into an FE object, it checks whether 
>>>       
>> the data 
>>     
>>> is local according to the Epetra_Map object it was 
>>>       
>> constructed with. 
>>     
>>> If so, that data is passed straight through to the non-FE 
>>>       
>> base-class 
>>     
>>> object. If the data belongs on another processor, then it is held 
>>> locally in "temporary" data structures until the 
>>>       
>> GlobalAssemble method is called.
>>     
>>> GlobalAssemble uses Import/Export operations to send that 
>>>       
>> data to the 
>>     
>>> appropriate owning processors, where it is placed into the non-FE 
>>> base-class object.
>>>
>>> So you should create your FECrsMatrix with the 
>>>       
>> non-overlapping row-map.
>>     
>>> FECrsMatrix shouldn't need your shared-map at all. Note 
>>>       
>> that the map 
>>     
>>> arguments which are passed to the GlobalAssemble method are the 
>>> domain-map and range-map for the underlying CrsMatrix, which is the 
>>> target of the global assembly. In general, a CrsMatrix has 4 maps:
>>> row-map, column-map, domain-map and range-map. In some cases the 
>>> domain-map is different than the column-map, and the range-map is 
>>> different than the row-map. See the doxygen documentation for 
>>> Epetra_CrsGraph for a description of these 4 map attributes. (At 
>>> trilinos.sandia.gov, go to packages => epetra => documentation.)
>>>
>>> The GlobalAssemble method should work without any map arguments in 
>>> most cases. (It internally obtains the domain-map and 
>>>       
>> range-map from 
>>     
>>> the underlying CrsMatrix.)
>>>
>>> Let me know if this doesn't clear things up for you, or if you are 
>>> still seeing errors.
>>>
>>> Note that non-zero return-values for SumIntoGlobalValues don't 
>>> necessarily indicate an error. If the return-value is a positive 
>>> number, it simply means that the matrix did an internal 
>>>       
>> allocation in 
>>     
>>> the process of storing the input values, because it wasn't 
>>>       
>> constructed 
>>     
>>> with enough information to pre-allocate all needed nonzero 
>>>       
>> positions.
>>     
>>> If you see a negative return-value, then that's an error.
>>>
>>> Alan
>>>
>>>
>>>   
>>>       
>>>> -----Original Message-----
>>>> From: trilinos-users-bounces at software.sandia.gov
>>>> [mailto:trilinos-users-bounces at software.sandia.gov] On Behalf Of 
>>>> David Neckels
>>>> Sent: Wednesday, November 28, 2007 3:12 PM
>>>> To: trilinos-users at software.sandia.gov
>>>> Subject: [Trilinos-Users] FECrs Vector/Matrix
>>>>
>>>> Hi,
>>>>
>>>> I installed the Epetra_FECrsMatrix and Epetra_FEVector 
>>>>         
>> classes in a 
>>     
>>>> parallel finite element code, but had some trouble figured 
>>>>         
>> out how to 
>>     
>>>> use the classes and their Epetra_Map's.  I settled on 
>>>>         
>> something that 
>>     
>>>> seems to work, but wasn't sure if it was correct?
>>>>
>>>> My code partitions on element boundary's, so faces, edges, 
>>>>         
>> and nodes 
>>     
>>>> on a shared element side exist on more than one processor.
>>>> Hence although each degree of freedom is owned by only one 
>>>>         
>> processor, 
>>     
>>>> it may reside on multiple processors.
>>>>
>>>> I have two corresponding Epetra_Map's:
>>>> unique_map: contains only the owned dofs on a processor 
>>>>         
>> (partitions 
>>     
>>>> the problem by processors).
>>>> shared_map: has overlap, i.e. a dof will appear on each processor 
>>>> where the dof resides.
>>>>
>>>> This was the sequence that seemed to work:
>>>>
>>>> matrix = new Epetra_FECrsMatrix(Copy, *shared_map, 
>>>>         
>> max_df); x = new 
>>     
>>>> Epetra_FEVector(*shared_map); //( solving Ax=b)  b = new 
>>>> Epetra_FEVector(*unique_map);
>>>>
>>>> // insert rows with all (shared) dofs // sum into global with all 
>>>> (shared) dofs --here is I used unique_map to create the matrix, 
>>>> above, I would get bad return values??
>>>> // Call ExtractGlobalView and replace constrained rows of 
>>>>         
>> the matrix 
>>     
>>>> for shared dofs (unique_map in the matrix create caused 
>>>>         
>> this to fail) 
>>     
>>>> // Also sum into shared dofs of b
>>>>
>>>> // After assembly:
>>>> matrix->GlobalAssemble(*umap, *umap);
>>>> b->GlobalAssemble();
>>>>
>>>> // Call Aztec solver....
>>>>
>>>> This dance between the unique and shared maps seemed strange. 
>>>>  I just wondered if a more thorough example for using the 
>>>>         
>> FE classes 
>>     
>>>> existed, or if someone knew what the correct approach was? 
>>>>         
>>  For now 
>>     
>>>> things seem to work....
>>>>
>>>> I have attached the routine where I set these maps up.
>>>>
>>>> Thanks for any help.
>>>>
>>>> -David N.
>>>>
>>>>
>>>>
>>>>     
>>>>         
>>>   
>>>       
>>
>>     
>
>