Re: [petsc-dev] Fortran binding for VecGetValuesSection

2022-06-25 Thread Jed Brown
Sorry about our late reply. Isn't VecGetValues similar? Note that 
VecGetValuesSection currently leaks a mutable pointer and really shouldn't 
exist like this without a matching VecRestoreValuesSection.

  PetscCall(VecGetArray(v, ));
  *values = [s->atlasOff[p]];
  PetscCall(VecRestoreArray(v, ));


Can you make a draft MR with what you've tried thus far?

Blaise Bourdin  writes:

> Hi,
>
> I’ve been trying to add a fortran90 binding for VecGetValuesSection, without 
> success, and could not find any function with a similar prototype to learn 
> from. Can anybody help?
>
> Regards,
> Blaise
>
> -- 
> Professor, Department of Mathematics & Statistics
> Hamilton Hall room 409A, McMaster University
> 1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
> https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243


Re: [petsc-dev] Problem with unnamed side sets in DMPlexCreateExodus

2022-06-25 Thread Jed Brown
Do you think this is more correct than just using "Face Sets", which gives you 
a numbered index?

There is curently some inconsistency between file formats in whether various 
sets show up as stand-alone labels or values under Face Sets. And there's this 
lingering issue to have a way to name Face Sets (and Cell Sets).

https://gitlab.com/petsc/petsc/-/issues/689

The issue is that if we just make named labels, there's no good way to do 
structured traversal (handle all Face Sets) and one also has to be careful to 
avoid namespace collisions between application sets and sets created by the 
mesh.

David Andrs  writes:

> Hello!
>
> The current behavior in DMPlexCreateExodus is assuming that side sets are 
> always named. However, that is not always true. If there are unnamed side 
> sets in an ExodusII file, the created DM will contain a label with an empty 
> name (“”) with as many stratas as there are side sets. See the sample output 
> from DMView below:
>
> Labels:
>   celltype: 3 strata with value/size (0 (138900), 4 (125514), 1 (265620))
>   depth: 3 strata with value/size (0 (138900), 1 (265620), 2 (125514))
>   Cell Sets: 4 strata with value/size (1 (57888), 2 (57888), 4 (1824), 5 
> (7914))
>   Face Sets: 6 strata with value/size (1000 (28944), 1001 (28944), 1003 
> (28944), 1005 (15144), 1100 (120), 1200 (120))
>   : 6 strata with value/size (1000 (28944), 1001 (28944), 1003 (28944), 1005 
> (15144), 1100 (120), 1200 (120))
>
> The attached patch fixes the behavior and names the side sets according to 
> their IDs if there is no name associated with it in the ExodusII file. See 
> the modified behavior below:
>
> Labels:
>   celltype: 3 strata with value/size (0 (138900), 4 (125514), 1 (265620))
>   depth: 3 strata with value/size (0 (138900), 1 (265620), 2 (125514))
>   Cell Sets: 4 strata with value/size (1 (57888), 2 (57888), 4 (1824), 5 
> (7914))
>   Face Sets: 6 strata with value/size (1000 (28944), 1001 (28944), 1003 
> (28944), 1005 (15144), 1100 (120), 1200 (120))
>   1000: 1 strata with value/size (1000 (28944))
>   1001: 1 strata with value/size (1001 (28944))
>   1003: 1 strata with value/size (1003 (28944))
>   1005: 1 strata with value/size (1005 (15144))
>   1100: 1 strata with value/size (1100 (120))
>   1200: 1 strata with value/size (1200 (120))
>
> The patch is against 3.17.0, but is seems it should apply cleanly on 3.17.2.
>
> With regards,
>
> David Andrs


Re: [petsc-dev] MatMPIAIJGetLocalMat problem with GPUs

2022-06-25 Thread Barry Smith


> On Jun 25, 2022, at 9:44 AM, Matthew Martineau  wrote:
> 
> Thanks - AmgX will accept unordered column indices.

  Any performance hit? 

  We can provide them sorted efficiently in the future but currently can only 
provide unordered column indices efficiently on the GPU.

  Barry
> 
> From: Barry Smith 
> Sent: Saturday, 25 June 2022, 14:39
> To: Mark Adams 
> Cc: Matthew Martineau ; For users of the development 
> version of PETSc 
> Subject: Re: [petsc-dev] MatMPIAIJGetLocalMat problem with GPUs
> 
> External email: Use caution opening links or attachments 
> 
> 
>   Does AMGX require sorted column indices? (Python indentation notation below)
> 
>   If not 
>  just use MatMPIAIJGetLocalMatMerge instead of MatMPIAIJGetLocalMat.
> 
>   If yes, 
>  on the first call
> 
>  Mat tmplocal;
>  PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_INITIAL_MATRIX, ));
>  PetscCall(MatConvert(tmplocal,MATSEQAIJCUSPARSE,>localA));
>  PetscCall(MatDestroy());
> 
>  leave the later calls as is with
> 
>  PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_REUSE_MATRIX, >localA);
> 
> Eventually, someone will need to buckle down and write 
> MatMPIAIJGetLocalMat_SeqAIJCUSPARSE, but that can be done later.
> 
> Barry
> 
> 
>> On Jun 25, 2022, at 9:13 AM, Mark Adams > > wrote:
>> 
>> 
>> 
>> On Fri, Jun 24, 2022 at 1:54 PM Barry Smith > > wrote:
>> 
>> 
>>> On Jun 24, 2022, at 1:38 PM, Mark Adams >> > wrote:
>>> 
>>> I am rearranging the code for clarity from the repo but I have:
>>> 
>>>   PetscBool is_dev_ptrs;
>>>   PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_INITIAL_MATRIX, >localA));
>>>   PetscCall(PetscObjectTypeCompareAny((PetscObject)amgx->localA, 
>>> _dev_ptrs, MATAIJCUSPARSE, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
>>>   PetscPrintf(PETSC_COMM_SELF,"checking against mataijcusparse amgx->localA 
>>> = %d\n",is_dev_ptrs);
>>>   PetscCall(PetscObjectTypeCompareAny((PetscObject)Pmat, _dev_ptrs, 
>>> MATAIJCUSPARSE, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
>>>   PetscPrintf(PETSC_COMM_SELF,"checking against mataijcusparse Pmat = 
>>> %d\n",is_dev_ptrs);
>>> 
>>> And it seems to show that MatMPIAIJGetLocalMat takes a mpiaijcusparse Mat 
>>> and returns an seqaij mat (see below):
>> 
>>  Yes, this is how it currently behaves as Stefano has indicated. Thus it is 
>> not currently directly suitable for use with GPUs. As Stefano has indicated 
>> it has be revised to handle mpiaijcusparse matrices correctly in the same 
>> way that MatMPIAIJGetLocalMatMerge has been revised for GPUs.
>> 
>> 
>> OK, sorry, I did not understand that this is not supported. We need a 
>> MatMPIAIJCusparseGetLocalMatMerge (I read this as supported with "hstack" 
>> format, unsorted?, by Stefano)
>> 
>> What is the best way to proceed? Should we just convert to amgx->localA to 
>> mpiaijcusparse if Pmat is a cusparse matrix?
>> If so, should this code go in amgx or 
>> MatMPIAIJGetLocalMat(MAT_INITIAL_MATRIX) ?
>> Or should I add a MatMPIAIJCusparseGetLocalMatMerge that simply wraps these 
>> two calls for now?
>> 
>> Thanks,
>> Mark
>> 
>>  
>> 
>>> 
>>> AMGX version 2.2.0.132-opensource
>>> Built on Jun 24 2022, 09:21:43
>>> Compiled with CUDA Runtime 11.5, using CUDA driver 11.5
>>> checking against mataijcusparse amgx->localA = 0
>>> checking against mataijcusparse Pmat = 1
>>> localA_name seqaij
>>> Pmat_name mpiaijcusparse
>>> 
>>> Matt's existing testing code (below) then shows the types that conform with 
>>> these tests and prints that I added.
>>> 
>>>   // XXX DEBUG REMOVE
>>>   const char* localA_name;
>>>   PetscObjectGetType((PetscObject)amgx->localA, _name);
>>>   PetscPrintf(PETSC_COMM_SELF,"localA_name %s\n", localA_name);
>>>   const char* Pmat_name;
>>>   PetscObjectGetType((PetscObject)Pmat, _name);
>>>   PetscPrintf(PETSC_COMM_SELF,"Pmat_name %s\n", Pmat_name);
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> On Fri, Jun 24, 2022 at 10:00 AM Barry Smith >> > wrote:
>>> 
>>> 
 On Jun 24, 2022, at 8:58 AM, Mark Adams >>> > wrote:
 
 And before we move to the MR, I think Matt found a clear problem:
 
 * PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_REUSE_MATRIX, >localA));
 returns "localA seqaij"
 
 * And, oddly, PetscCall(MatSeqAIJGetArrayRead(amgx->localA, 
 >values)); returns:
 "it seems to detect that the pointer is a device mapped pointer but that 
 it is invalid"
>>> 
>>>   It does not return a device mapped pointer, it returns a valid host 
>>> pointer only. MatSeqAIJGetArrayRead() is intended to only return a host 
>>> pointer, it cannot return a device pointer. MatSeqAIJCusparseGetArrayRead() 
>>> returns device pointers and should be used for this purpose.
>>> 
 
 Matt, lets just comment out the REUSE line and add another INITIAL line 
 (destroying the old Mat of course), and lets press on.
>>> 
>>>   Looking at the code there is no way 

Re: [petsc-dev] MatMPIAIJGetLocalMat problem with GPUs

2022-06-25 Thread Matthew Martineau via petsc-dev
Thanks - AmgX will accept unordered column indices.


From: Barry Smith 
Sent: Saturday, 25 June 2022, 14:39
To: Mark Adams 
Cc: Matthew Martineau ; For users of the development 
version of PETSc 
Subject: Re: [petsc-dev] MatMPIAIJGetLocalMat problem with GPUs

External email: Use caution opening links or attachments


  Does AMGX require sorted column indices? (Python indentation notation below)

  If not
 just use MatMPIAIJGetLocalMatMerge instead of MatMPIAIJGetLocalMat.

  If yes,
 on the first call

 Mat tmplocal;
 PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_INITIAL_MATRIX, ));
 PetscCall(MatConvert(tmplocal,MATSEQAIJCUSPARSE,>localA));
 PetscCall(MatDestroy());

 leave the later calls as is with

 PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_REUSE_MATRIX, >localA);

Eventually, someone will need to buckle down and write 
MatMPIAIJGetLocalMat_SeqAIJCUSPARSE, but that can be done later.

Barry


On Jun 25, 2022, at 9:13 AM, Mark Adams 
mailto:mfad...@lbl.gov>> wrote:



On Fri, Jun 24, 2022 at 1:54 PM Barry Smith 
mailto:bsm...@petsc.dev>> wrote:


On Jun 24, 2022, at 1:38 PM, Mark Adams 
mailto:mfad...@lbl.gov>> wrote:

I am rearranging the code for clarity from the repo but I have:

  PetscBool is_dev_ptrs;
  PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_INITIAL_MATRIX, >localA));
  PetscCall(PetscObjectTypeCompareAny((PetscObject)amgx->localA, _dev_ptrs, 
MATAIJCUSPARSE, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
  PetscPrintf(PETSC_COMM_SELF,"checking against mataijcusparse amgx->localA = 
%d\n",is_dev_ptrs);
  PetscCall(PetscObjectTypeCompareAny((PetscObject)Pmat, _dev_ptrs, 
MATAIJCUSPARSE, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
  PetscPrintf(PETSC_COMM_SELF,"checking against mataijcusparse Pmat = 
%d\n",is_dev_ptrs);

And it seems to show that MatMPIAIJGetLocalMat takes a mpiaijcusparse Mat and 
returns an seqaij mat (see below):

 Yes, this is how it currently behaves as Stefano has indicated. Thus it is not 
currently directly suitable for use with GPUs. As Stefano has indicated it has 
be revised to handle mpiaijcusparse matrices correctly in the same way that 
MatMPIAIJGetLocalMatMerge has been revised for GPUs.


OK, sorry, I did not understand that this is not supported. We need a 
MatMPIAIJCusparseGetLocalMatMerge (I read this as supported with "hstack" 
format, unsorted?, by Stefano)

What is the best way to proceed? Should we just convert to amgx->localA to 
mpiaijcusparse if Pmat is a cusparse matrix?
If so, should this code go in amgx or MatMPIAIJGetLocalMat(MAT_INITIAL_MATRIX) ?
Or should I add a MatMPIAIJCusparseGetLocalMatMerge that simply wraps these two 
calls for now?

Thanks,
Mark




AMGX version 2.2.0.132-opensource
Built on Jun 24 2022, 09:21:43
Compiled with CUDA Runtime 11.5, using CUDA driver 11.5
checking against mataijcusparse amgx->localA = 0
checking against mataijcusparse Pmat = 1
localA_name seqaij
Pmat_name mpiaijcusparse

Matt's existing testing code (below) then shows the types that conform with 
these tests and prints that I added.

  // XXX DEBUG REMOVE
  const char* localA_name;
  PetscObjectGetType((PetscObject)amgx->localA, _name);
  PetscPrintf(PETSC_COMM_SELF,"localA_name %s\n", localA_name);
  const char* Pmat_name;
  PetscObjectGetType((PetscObject)Pmat, _name);
  PetscPrintf(PETSC_COMM_SELF,"Pmat_name %s\n", Pmat_name);






On Fri, Jun 24, 2022 at 10:00 AM Barry Smith 
mailto:bsm...@petsc.dev>> wrote:


On Jun 24, 2022, at 8:58 AM, Mark Adams 
mailto:mfad...@lbl.gov>> wrote:

And before we move to the MR, I think Matt found a clear problem:

* PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_REUSE_MATRIX, >localA));
returns "localA seqaij"

* And, oddly, PetscCall(MatSeqAIJGetArrayRead(amgx->localA, >values)); 
returns:
"it seems to detect that the pointer is a device mapped pointer but that it is 
invalid"

  It does not return a device mapped pointer, it returns a valid host pointer 
only. MatSeqAIJGetArrayRead() is intended to only return a host pointer, it 
cannot return a device pointer. MatSeqAIJCusparseGetArrayRead() returns device 
pointers and should be used for this purpose.


Matt, lets just comment out the REUSE line and add another INITIAL line 
(destroying the old Mat of course), and lets press on.

  Looking at the code there is no way that simply using INITIAL instead of 
REUSE will make a code that does not work on the GPU run on the GPU. The 
MatMPIAIJGetLocalMat() returns only a MATSEQAIJ matrix regardless of the 
INITIAL versus REUSE and one can never get a device pointer from a non-GPU 
matrix.

   As noted by Stefano, the code either needs to use MatMPIAIJGetLocalMatMerge 
() which does return a CUSPARSE matrix (but the columns are not supported) or 
MatMPIAIJGetLocalMat()
needs to be updated to return a CUSPARSE matrix when the input MPI matrix is a 
CUSPARSE matrix.




We can keep the debugging code for now.

We (PETSc) can work on this independently,

Thanks,
Mark

On Fri, 

Re: [petsc-dev] MatMPIAIJGetLocalMat problem with GPUs

2022-06-25 Thread Barry Smith

  Does AMGX require sorted column indices? (Python indentation notation below)

  If not 
 just use MatMPIAIJGetLocalMatMerge instead of MatMPIAIJGetLocalMat.

  If yes, 
 on the first call

 Mat tmplocal;
 PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_INITIAL_MATRIX, ));
 PetscCall(MatConvert(tmplocal,MATSEQAIJCUSPARSE,>localA));
 PetscCall(MatDestroy());

 leave the later calls as is with

 PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_REUSE_MATRIX, >localA);

Eventually, someone will need to buckle down and write 
MatMPIAIJGetLocalMat_SeqAIJCUSPARSE, but that can be done later.

Barry


> On Jun 25, 2022, at 9:13 AM, Mark Adams  wrote:
> 
> 
> 
> On Fri, Jun 24, 2022 at 1:54 PM Barry Smith  > wrote:
> 
> 
>> On Jun 24, 2022, at 1:38 PM, Mark Adams > > wrote:
>> 
>> I am rearranging the code for clarity from the repo but I have:
>> 
>>   PetscBool is_dev_ptrs;
>>   PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_INITIAL_MATRIX, >localA));
>>   PetscCall(PetscObjectTypeCompareAny((PetscObject)amgx->localA, 
>> _dev_ptrs, MATAIJCUSPARSE, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
>>   PetscPrintf(PETSC_COMM_SELF,"checking against mataijcusparse amgx->localA 
>> = %d\n",is_dev_ptrs);
>>   PetscCall(PetscObjectTypeCompareAny((PetscObject)Pmat, _dev_ptrs, 
>> MATAIJCUSPARSE, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
>>   PetscPrintf(PETSC_COMM_SELF,"checking against mataijcusparse Pmat = 
>> %d\n",is_dev_ptrs);
>> 
>> And it seems to show that MatMPIAIJGetLocalMat takes a mpiaijcusparse Mat 
>> and returns an seqaij mat (see below):
> 
>  Yes, this is how it currently behaves as Stefano has indicated. Thus it is 
> not currently directly suitable for use with GPUs. As Stefano has indicated 
> it has be revised to handle mpiaijcusparse matrices correctly in the same way 
> that MatMPIAIJGetLocalMatMerge has been revised for GPUs.
> 
> 
> OK, sorry, I did not understand that this is not supported. We need a 
> MatMPIAIJCusparseGetLocalMatMerge (I read this as supported with "hstack" 
> format, unsorted?, by Stefano)
> 
> What is the best way to proceed? Should we just convert to amgx->localA to 
> mpiaijcusparse if Pmat is a cusparse matrix?
> If so, should this code go in amgx or 
> MatMPIAIJGetLocalMat(MAT_INITIAL_MATRIX) ?
> Or should I add a MatMPIAIJCusparseGetLocalMatMerge that simply wraps these 
> two calls for now?
> 
> Thanks,
> Mark
> 
>  
> 
>> 
>> AMGX version 2.2.0.132-opensource
>> Built on Jun 24 2022, 09:21:43
>> Compiled with CUDA Runtime 11.5, using CUDA driver 11.5
>> checking against mataijcusparse amgx->localA = 0
>> checking against mataijcusparse Pmat = 1
>> localA_name seqaij
>> Pmat_name mpiaijcusparse
>> 
>> Matt's existing testing code (below) then shows the types that conform with 
>> these tests and prints that I added.
>> 
>>   // XXX DEBUG REMOVE
>>   const char* localA_name;
>>   PetscObjectGetType((PetscObject)amgx->localA, _name);
>>   PetscPrintf(PETSC_COMM_SELF,"localA_name %s\n", localA_name);
>>   const char* Pmat_name;
>>   PetscObjectGetType((PetscObject)Pmat, _name);
>>   PetscPrintf(PETSC_COMM_SELF,"Pmat_name %s\n", Pmat_name);
>> 
>> 
>> 
>> 
>> 
>> 
>> On Fri, Jun 24, 2022 at 10:00 AM Barry Smith > > wrote:
>> 
>> 
>>> On Jun 24, 2022, at 8:58 AM, Mark Adams >> > wrote:
>>> 
>>> And before we move to the MR, I think Matt found a clear problem:
>>> 
>>> * PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_REUSE_MATRIX, >localA));
>>> returns "localA seqaij"
>>> 
>>> * And, oddly, PetscCall(MatSeqAIJGetArrayRead(amgx->localA, 
>>> >values)); returns:
>>> "it seems to detect that the pointer is a device mapped pointer but that it 
>>> is invalid"
>> 
>>   It does not return a device mapped pointer, it returns a valid host 
>> pointer only. MatSeqAIJGetArrayRead() is intended to only return a host 
>> pointer, it cannot return a device pointer. MatSeqAIJCusparseGetArrayRead() 
>> returns device pointers and should be used for this purpose.
>> 
>>> 
>>> Matt, lets just comment out the REUSE line and add another INITIAL line 
>>> (destroying the old Mat of course), and lets press on.
>> 
>>   Looking at the code there is no way that simply using INITIAL instead of 
>> REUSE will make a code that does not work on the GPU run on the GPU. The 
>> MatMPIAIJGetLocalMat() returns only a MATSEQAIJ matrix regardless of the 
>> INITIAL versus REUSE and one can never get a device pointer from a non-GPU 
>> matrix.
>> 
>>As noted by Stefano, the code either needs to use 
>> MatMPIAIJGetLocalMatMerge () which does return a CUSPARSE matrix (but the 
>> columns are not supported) or MatMPIAIJGetLocalMat()
>> needs to be updated to return a CUSPARSE matrix when the input MPI matrix is 
>> a CUSPARSE matrix.
>> 
>> 
>> 
>> 
>>> We can keep the debugging code for now.
>>> 
>>> We (PETSc) can work on this independently,
>>> 
>>> Thanks,
>>> Mark
>>> 
>>> On Fri, Jun 24, 2022 at 

Re: [petsc-dev] MatMPIAIJGetLocalMat problem with GPUs

2022-06-25 Thread Mark Adams
On Fri, Jun 24, 2022 at 1:54 PM Barry Smith  wrote:

>
>
> On Jun 24, 2022, at 1:38 PM, Mark Adams  wrote:
>
> I am rearranging the code for clarity from the repo but I have:
>
>   PetscBool is_dev_ptrs;
>   PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_INITIAL_MATRIX, >localA));
>   PetscCall(PetscObjectTypeCompareAny((PetscObject)amgx->localA,
> _dev_ptrs, MATAIJCUSPARSE, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
>   PetscPrintf(PETSC_COMM_SELF,"*checking against mataijcusparse
> amgx->localA = %d*\n",is_dev_ptrs);
>   PetscCall(PetscObjectTypeCompareAny((PetscObject)Pmat, _dev_ptrs,
> MATAIJCUSPARSE, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
>   PetscPrintf(PETSC_COMM_SELF,"*checking against mataijcusparse Pmat = %d*
> \n",is_dev_ptrs);
>
> And it seems to show that MatMPIAIJGetLocalMat takes a mpiaijcusparse Mat
> and returns an seqaij mat (see below):
>
>
>  Yes, this is how it currently behaves as Stefano has indicated. Thus it
> is not currently directly suitable for use with GPUs. As Stefano has
> indicated it has be revised to handle mpiaijcusparse matrices correctly
> in the same way that MatMPIAIJGetLocalMatMerge has been revised for GPUs.
>
>
OK, sorry, I did not understand that this is not supported. We need a
MatMPIAIJCusparseGetLocalMatMerge (I read this as supported with "hstack"
format, unsorted?, by Stefano)

What is the best way to proceed? Should we just convert to amgx->localA
to mpiaijcusparse if Pmat is a cusparse matrix?
If so, should this code go in amgx or
MatMPIAIJGetLocalMat(MAT_INITIAL_MATRIX) ?
Or should I add a MatMPIAIJCusparseGetLocalMatMerge that simply wraps these
two calls for now?

Thanks,
Mark



>
>
> AMGX version 2.2.0.132-opensource
> Built on Jun 24 2022, 09:21:43
> Compiled with CUDA Runtime 11.5, using CUDA driver 11.5
>
> *checking against mataijcusparse amgx->localA = 0checking against
> mataijcusparse Pmat = 1*
> localA_name seqaij
> Pmat_name mpiaijcusparse
>
> Matt's existing testing code (below) then shows the types that conform
> with these tests and prints that I added.
>
>   // XXX DEBUG REMOVE
>   const char* localA_name;
>   PetscObjectGetType((PetscObject)amgx->localA, _name);
>   PetscPrintf(PETSC_COMM_SELF,"localA_name %s\n", localA_name);
>   const char* Pmat_name;
>   PetscObjectGetType((PetscObject)Pmat, _name);
>   PetscPrintf(PETSC_COMM_SELF,"Pmat_name %s\n", Pmat_name);
>
>
>
>
>
>
> On Fri, Jun 24, 2022 at 10:00 AM Barry Smith  wrote:
>
>>
>>
>> On Jun 24, 2022, at 8:58 AM, Mark Adams  wrote:
>>
>> And before we move to the MR, I think Matt found a clear problem:
>>
>> * PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_REUSE_MATRIX, >localA));
>> returns "localA seqaij"
>>
>> * And, oddly, PetscCall(MatSeqAIJGetArrayRead(amgx->localA,
>> >values)); returns:
>> "it seems to detect that the pointer is a device mapped pointer but that
>> it is invalid"
>>
>>
>>   It does not return a device mapped pointer, it returns a valid host
>> pointer only. MatSeqAIJGetArrayRead() is intended to only return a host
>> pointer, it cannot return a device pointer. MatSeqAIJCusparseGetArrayRead()
>> returns device pointers and should be used for this purpose.
>>
>>
>> Matt, lets just comment out the REUSE line and add another INITIAL line
>> (destroying the old Mat of course), and lets press on.
>>
>>
>>   Looking at the code there is no way that simply using INITIAL instead
>> of REUSE will make a code that does not work on the GPU run on the GPU. The
>> MatMPIAIJGetLocalMat() returns only a MATSEQAIJ matrix regardless of
>> the INITIAL versus REUSE and one can never get a device pointer from a
>> non-GPU matrix.
>>
>>As noted by Stefano, the code either needs to use
>> MatMPIAIJGetLocalMatMerge () which does return a CUSPARSE matrix (but
>> the columns are not supported) or MatMPIAIJGetLocalMat()
>> needs to be updated to return a CUSPARSE matrix when the input MPI matrix
>> is a CUSPARSE matrix.
>>
>>
>>
>>
>> We can keep the debugging code for now.
>>
>> We (PETSc) can work on this independently,
>>
>> Thanks,
>> Mark
>>
>> On Fri, Jun 24, 2022 at 8:51 AM Mark Adams  wrote:
>>
>>> I am not seeing this response, I see my "hstack" comment last.
>>> https://gitlab.com/petsc/petsc/-/merge_requests/4323
>>>
>>> On Thu, Jun 23, 2022 at 4:37 PM Barry Smith  wrote:
>>>

  I have responded in the MR, which has all the context and the code.
 Please move this conversation from petsc-dev to the MR. Note you can use
 the little cartoon cloud symbol (upper write of the sub window with my
 text)  to reply to my post and keep everything in a thread for clarity.

   We are confused because it seems you are trying a variety of things
 and we don't know how the different things you tried resulted in the
 multiple errors you reported.



 On Jun 23, 2022, at 3:59 PM, Matthew Martineau 
 wrote:

 I checked in the changes and some debugging statements.

 PetscCall(MatMPIAIJGetLocalMat(Pmat, MAT_REUSE_MATRIX,