Question: Externally owned data in hypre_CSRMatrix

[ictmay]

I am working with the LANL AMP team. We have a parallel CSR matrix format that should be compatible with the Hypre format. Ideally we would wrap our format into a Hypre IJMatrix with a minimal amount of duplicate data.

Question #609 provides some initial information on how to approach this. I have looked at which (de)allocations are guarded by owns_data inside hypre_ParCSRMatrixCreate, hypre_CSRMatrixCreate, and the matching destroy functions to decide what we should own versus what Hypre should own. The heavy data all lives in the diag and offd blocks.

At present I do the following:

1. Create IJMatrix and set type to HYPRE_PARCSR, do not call initialize.
2. Call hypre_IJMatrixCreateParCSR. This sets owns_data in ParCSR to 1 and creates diag and offd structs.
3. Set need_aux to zero inside the IJMatrix->translator.
4. Grab the diag and offd blocks and call hypre_CSRMatrixInitialize manually. This allocates the ->i field which is not guarded by owns_data.
5. Set ownership to zero on diag and offd via hypre_CSRMatrixSetDataOwner.
6. Fill the ->i fields.
7. Set the ->big_j, ->j, and ->data fields of diag and offd to reference our data. These are all guarded by owns_data and are the heaviest memory requirements.

The issue comes from the ParCSR->col_map_offd field. I'd prefer to leave ParCSR->owns_data==1 so that it will manage the setup/teardown of the CSRMatrix structs themselves and then I can just set owns_data to zero in there. This means that Hypre needs to also own the col_map_offd field. The main benefit is that HYPRE_IJMatrixDestroy does (seem to) do all of the correct deallocations.

That is a long-winded set up to premise a couple of questions:

• Is this approach valid, or should I set ParCSR->owns_data=0 and manage more of the setup/teardown on my end?
• I can use HYPRE_IJMatrixAssemble to get col_map_offd created and filled. The internal function hypre_IJMatrixAssembleParCSR does not check the value of offd->owns_data and frees offd->big_j yielding an eventual double-free (line 2922 of IJMatrix_parcsr.c).
• I'm a little confused by the code that fills col_map_offd and rearranges offd->j. What are the exact requirements on col_map_offd? It must be more than just a map from local->global indices since offd->j also gets shuffled.

I did add a check for the offd->owns_data flag into hypre_IJMatrixAssembleParCSR, and that does avoid the double-free issue. Surprisingly this allows BoomerAMG to work as a solver, but I think that is a fluke. This doesn't really solve the overall issue. That function still modifies values in the big_j and j arrays that it doesn't own. As a minimal solution it should check the owns_data flag and throw an error.

That's quite the wall of text. I appreciate any help, and please let me know if any clarifications are needed.

[ictmay]

I am re-opening this to update it to what seems to be working. I have only medium confidence here, so grain of salt. Any input on whether this is valid would still be appreciated.

The setup is as follows. This is not a direct copy of the code but should get the point across.

1. Call HYPRE_IJMatrixCreate(..., &d_matrix) and HYPRE_IJMatrixSetObjectType to make a new parcsr matrix
2. Call hypre_IJMatrixCreateParCSR( d_matrix )
3. Extract the parcsr matrix via hypre_ParCSRMatrix *par_matrix = (hypre_ParCSRMatrix *) d_matrix->object
4. Extract on/off diagonal blocks hypre_CSRMatrix *diag = par_matrix->diag and hypre_CSRMatrix *offd = par_matrix->offd
5. Extract auxiliary matrix d_matrix->translator and set need_aux to zero inside it
6. Call hypre_CSRMatrixInitialize on both diag and offd
7. Fill diag->i according to my externally owned matrix
8. Fill offd->i if I have one, otherwise just set all zeros
9. Set ownership to zero with hypre_CSRMatrixSetDataOwner( diag, 0 ) and hypre_CSRMatrixSetDataOwner( offd, 0 )
10. Set diag->big_j, diag->j, and diag->data, to point at corresponding fields of external matrix. ->big_j presumably optional
11. Do the same for offd
12. Set diag->num_nonzeros and offd->num_nonzeros
13. Set par_matrix->col_map_offd to map converting from local to global indices in off diagonal block
14. Set offd->num_cols to match length of col_map_offd
15. Call hypre_CSRMatrixSetRownnz on diag and offd
16. Finally, set d_matrix->assemble_flag = 1

A few of these steps are non-trivial and potentially wrong.

Step 6 allocates the ->i field in diag and offd. These fields are not protected by ->owns_data so Hypre needs to allocate them (they are small so it doesn't really matter). Step 15 does the same thing for the same reason.

Step 13 and the layouts of ->j in steps 10 and 11 are also touchy. In diag the j field should be sorted to have the diagonal element occur first in each row, and the remaining indices should be sorted in ascending order. The j field in offd should have ascending indices in each row. Finally, col_map_offd should be in ascending order (changing the local indices to match).

[victorapm]

@ictmay thank you for the update and apologize for the delay.

On the surface, the logic seems correct including your comments about column ordering (except for remaining column indices being sorted in ascending order: that is not strictly necessary but advisable). Is this implementation openly available somewhere in AMP so we can take a look?

We can also discuss about adding owns field to some of the pointers in our matrix structure if you find it useful. Other suggestions to make the implementation of your interface easier are also welcomed!

Thank you!

[ictmay]

@victorapm No worries on the delay, and thank you for getting back to us.

For access to the source I'll have to talk to Bobby. I only work on the internal (restricted enclave) version, and I'm not sure what all has made it into the externally viewable version.

Regarding the owns flag(s), the main item is the col_map_offd member of hypre_ParCSRMatrix. I neglected to mention it above, but that one is also not guarded. Right now I am careful to null that field out before calling HYPRE_IJMatrixDestroy.

The column map isn't so large (except when it is...), so it could easily fall into the same category as the ->i and ->rownnz fields inside the on/off diagonal blocks. It also seems that Hypre will break in mysterious ways if this map is not ordered in exactly the right way. Having Hypre own and create this field would be fine and I could just make my own copy/reference of it after the fact for our communicators.

A nice solution that would avoid most of the touchy steps above, though one that would take some non-zero dev time and testing time, would be to make HYPRE_IJMatrixAssemble (or some similar function) callable in the non-owned data case. This would let Hypre set up all of the meta-data as it sees fit, and could remove many of the steps above.

I don't know how often people want to do shallow-copies into Hypre, and without broader demand it makes more sense to have this as an interface we maintain rather than something that gets added to Hypre.

[ictmay]

@victorapm All of the shallow copy code is present in our external code here:

https://github.com/AdvancedMultiPhysics/AMP/blob/master/src/matrices/data/hypre/HypreMatrixAdaptor.cpp

Lines 41 -- 43 create the IJMatrix, and the process described above starts on line 136 and continues to the end of the file.