Re: [petsc-users] Question about petsc4py createWithArray function

[Samar Khatiwala]

Hi Stefano,

Thank you for this. Your example was invaluable in helping me figure out what 
is going on. (I did not know about about sys.getrefcount - I’m only just 
starting to learn python.) It was indeed the refcount. Until the Vec is 
destroyed it can’t go to zero and my array destructor method won’t be triggered 
(which is not really a problem because the memory is freed elsewhere). I take 
back what I said about there being a memory leak in petsc4py. Sorry.

Thanks again for your help!

Best,

Samar

On May 6, 2024, at 2:09 PM, Stefano Zampini  wrote:

Samar

After a second look, I believe the petsc4py code is correct. You can test it 
using the script below.
After destroy is called (or del), the reference count of the numpy array is 
back to its initial state.
Maybe you are not calling del or destroy? a MWE would help to understand your 
use case

kl-18448:~ szampini$ cat t.py
import sys
import numpy as np
from petsc4py import PETSc

x = np.zeros(4, dtype=PETSc.ScalarType)
print('initial ref count',sys.getrefcount(x))

v = PETSc.Vec().createWithArray(x)
print('after create',sys.getrefcount(x))

# check if they share memory
v.view()
x[1] = 2
v.view()

# free
v.destroy()
# you can also call del
# del v
print('after destroy',sys.getrefcount(x))

kl-18448:~ szampini$ python t.py
initial ref count 2
after create 3
Vec Object: 1 MPI process
  type: seq
0.
0.
0.
0.
Vec Object: 1 MPI process
  type: seq
0.
2.
0.
0.
after destroy 2



Il giorno lun 6 mag 2024 alle ore 08:49 Samar Khatiwala 
mailto:samar.khatiw...@earth.ox.ac.uk>> ha 
scritto:
Hi Stefano,

Thanks for looking into this. Since createWithArray calls VecCreateMPIWithArray 
which, as Matt noted and is documented 
(https://urldefense.us/v3/__https://petsc.org/main/manualpages/Vec/VecCreateMPIWithArray/__;!!G_uCfscf7eWS!fl-7I7xkuEmTFH8M9EaARIZL5vhGFWkIk5bVIGt0bkt_cYcuwPSwgOWaec1-DvYwiKWxA5q0HS0VTLBmD4NnSe1PfG057Y4SBd-8qs8$
 ) doesn’t free the memory, then there’s a memory leak (and, furthermore, 
calling del on the original array will have no effect).

Lisandro: would be great if you can provide some guidance.

Thanks,

Samar

On May 3, 2024, at 12:45 PM, Stefano Zampini 
mailto:stefano.zamp...@gmail.com>> wrote:

While waiting for our Python wizard to shed light on this, I note that, from 
the documentation of PyArray_FROM_OTF 
https://urldefense.us/v3/__https://numpy.org/devdocs/user/c-info.how-to-extend.html*converting-an-arbitrary-sequence-object__;Iw!!G_uCfscf7eWS!fl-7I7xkuEmTFH8M9EaARIZL5vhGFWkIk5bVIGt0bkt_cYcuwPSwgOWaec1-DvYwiKWxA5q0HS0VTLBmD4NnSe1PfG057Y4S1ro4qy0$
 , we have

The object can be any Python object convertible to an ndarray. If the object is 
already (a subclass of) the ndarray that satisfies the requirements then a new 
reference is returned.

I guess we should call "del" on the ndarray returned by iarray_s after having 
called  self.set_attr('__array__', array) in this case, but let's wait for 
Lisandro to confirm




Il giorno ven 3 mag 2024 alle ore 11:42 Samar Khatiwala 
mailto:samar.khatiw...@earth.ox.ac.uk>> ha 
scritto:
This Message Is From an External Sender
This message came from outside your organization.

Hi Matt,

Thanks so much for the quick reply!

Regarding #2, I put some debug statement in my code and what I find is that 
when I use createWithArray on my Cython-allocated numpy array, the destructor I 
set for it is no longer called when I delete the array. (If I don’t use 
createWithArray then the destructor is triggered.) I interpret that to suggest 
that the petsc4py Vec is somehow ’taking over’ management of the numpy array. 
But I don’t understand where that could be happening. (I don’t think it has to 
do with the actual freeing of memory by PETSc's VecDestroy.)

createWithArray calls iarray_s which in turn calls PyArray_FROM_OTF. Could it 
be there’s something going on there? The numpy documentation is unclear.

Lisandro: do you have any thoughts on this?

Thanks,

Samar

On May 2, 2024, at 11:56 PM, Matthew Knepley 
mailto:knep...@gmail.com>> wrote:

On Thu, May 2, 2024 at 12:53 PM Samar Khatiwala 
mailto:samar.khatiw...@earth.ox.ac.uk>> wrote:
This Message Is From an External Sender
This message came from outside your organization.


Hello,

I have a couple of questions about createWithArray in petsc4py:

1) What is the correct usage for creating a standard MPI Vec with it? Something 
like this seems to work but is it right?:

On each rank do:
a = np.zeros(localSize)
v = PETSc.Vec().createWithArray(a, comm=PETSc.COMM_WORLD)

Is that all it takes?

That looks right to me.

2) Who ‘owns’ the underlying memory for a Vec created with the createWithArray 
method, i.e., who is responsible for managing it and doing garbage collection? 
In my problem, the numpy array is created in a Cython module where memory is 
allocated, and a pointer to it is associated with a numpy ndarray via 
PyArray_SimpleNewFromData and PyArray_SetBaseObject. I have a deallocator 
method of my own that is called when the numpy array is 

Re: [petsc-users] Question about petsc4py createWithArray function

[Stefano Zampini]

Samar

After a second look, I believe the petsc4py code is correct. You can test
it using the script below.
After destroy is called (or del), the reference count of the numpy array is
back to its initial state.
Maybe you are not calling del or destroy? a MWE would help to understand
your use case


kl-18448:~ szampini$ cat t.py

import sys

import numpy as np

from petsc4py import PETSc


x = np.zeros(4, dtype=PETSc.ScalarType)

print('initial ref count',sys.getrefcount(x))


v = PETSc.Vec().createWithArray(x)

print('after create',sys.getrefcount(x))


# check if they share memory

v.view()

x[1] = 2

v.view()


# free

v.destroy()

# you can also call del

# del v

print('after destroy',sys.getrefcount(x))


kl-18448:~ szampini$ python t.py

initial ref count 2

after create 3

Vec Object: 1 MPI process

  type: seq

0.

0.

0.

0.

Vec Object: 1 MPI process

  type: seq

0.

2.

0.

0.

after destroy 2



Il giorno lun 6 mag 2024 alle ore 08:49 Samar Khatiwala <
samar.khatiw...@earth.ox.ac.uk> ha scritto:

> Hi Stefano,
>
> Thanks for looking into this. Since createWithArray calls
> VecCreateMPIWithArray which, as Matt noted and is documented (
> https://urldefense.us/v3/__https://petsc.org/main/manualpages/Vec/VecCreateMPIWithArray/__;!!G_uCfscf7eWS!Z9RkH5ffTuJDtBUf8_Gk0BHuG__BKv4jPYeg89Rp6_GcS9VTcFs2J8uyLd5_wiqdDmy9ABXjczA3PYB1raRliVrs4DeDixY$
>  ) doesn’t
> free the memory, then there’s a memory leak (and, furthermore, calling del
> on the original array will have no effect).
>
> Lisandro: would be great if you can provide some guidance.
>
> Thanks,
>
> Samar
>
> On May 3, 2024, at 12:45 PM, Stefano Zampini 
> wrote:
>
> While waiting for our Python wizard to shed light on this, I note that,
> from the documentation of PyArray_FROM_OTF
> https://urldefense.us/v3/__https://numpy.org/devdocs/user/c-info.how-to-extend.html*converting-an-arbitrary-sequence-object__;Iw!!G_uCfscf7eWS!Z9RkH5ffTuJDtBUf8_Gk0BHuG__BKv4jPYeg89Rp6_GcS9VTcFs2J8uyLd5_wiqdDmy9ABXjczA3PYB1raRliVrsskf0ioU$
>  ,
> we have
>
> The object can be any Python object convertible to an ndarray. If the
> object is already (a subclass of) the ndarray that satisfies the
> requirements then a new reference is returned.
>
> I guess we should call "del" on the ndarray returned by iarray_s after
> having called  self.set_attr('__array__', array) in this case, but let's
> wait for Lisandro to confirm
>
>
>
>
> Il giorno ven 3 mag 2024 alle ore 11:42 Samar Khatiwala <
> samar.khatiw...@earth.ox.ac.uk> ha scritto:
>
>> Hi Matt, Thanks so much for the quick reply! Regarding #2, I put some
>> debug statement in my code and what I find is that when I use
>> createWithArray on my Cython-allocated numpy array, the destructor I set
>> for it is no longer called when I delete
>> ZjQcmQRYFpfptBannerStart
>> This Message Is From an External Sender
>> This message came from outside your organization.
>>
>> ZjQcmQRYFpfptBannerEnd
>> Hi Matt,
>>
>> Thanks so much for the quick reply!
>>
>> Regarding #2, I put some debug statement in my code and what I find is
>> that when I use createWithArray on my Cython-allocated numpy array, the
>> destructor I set for it is no longer called when I delete the array. (If I
>> don’t use createWithArray then the destructor is triggered.) I interpret
>> that to suggest that the petsc4py Vec is somehow ’taking over’ management
>> of the numpy array. But I don’t understand where that could be
>> happening. (I don’t think it has to do with the actual freeing of memory by
>> PETSc's VecDestroy.)
>>
>> createWithArray calls iarray_s which in turn calls PyArray_FROM_OTF.
>> Could it be there’s something going on there? The numpy documentation is
>> unclear.
>>
>> Lisandro: do you have any thoughts on this?
>>
>> Thanks,
>>
>> Samar
>>
>> On May 2, 2024, at 11:56 PM, Matthew Knepley  wrote:
>>
>> On Thu, May 2, 2024 at 12:53 PM Samar Khatiwala <
>> samar.khatiw...@earth.ox.ac.uk> wrote:
>>
>>> This Message Is From an External Sender
>>> This message came from outside your organization.
>>>
>>>
>>> Hello,
>>>
>>> I have a couple of questions about createWithArray in petsc4py:
>>>
>>> 1) What is the correct usage for creating a standard MPI Vec with it? 
>>> Something like this seems to work but is it right?:
>>>
>>> On each rank do:
>>> a = np.zeros(localSize)
>>> v = PETSc.Vec().createWithArray(a, comm=PETSc.COMM_WORLD)
>>>
>>> Is that all it takes?
>>>
>>>
>> That looks right to me.
>>
>>> 2) Who ‘owns’ the underlying memory for a Vec created with the 
>>> createWithArray method, i.e., who is responsible for managing it and doing 
>>> garbage collection? In my problem, the numpy array is created in a Cython 
>>> module where memory is allocated, and a pointer to it is associated with a 
>>> numpy ndarray via PyArray_SimpleNewFromData and PyArray_SetBaseObject. I 
>>> have a deallocator method of my own that is called when the numpy array is 
>>> deleted/goes out of scope/whenever python does garbage collection. All of 
>>> that works 

Re: [petsc-users] Question about petsc4py createWithArray function

[Samar Khatiwala]

Hi Stefano,

Thanks for looking into this. Since createWithArray calls VecCreateMPIWithArray 
which, as Matt noted and is documented 
(https://urldefense.us/v3/__https://petsc.org/main/manualpages/Vec/VecCreateMPIWithArray/__;!!G_uCfscf7eWS!YpH_OGif6W8Ql9KZDYBiRMKyALSRhMZOHRnkCkBHqBcVB4ef_UBHFbT_EuTgx8ivGmDNQU8bRVaOc-B_tqQ2heGjdvoE9tp8ZbPi3QQ$
 ) doesn’t free the memory, then there’s a memory leak (and, furthermore, 
calling del on the original array will have no effect).

Lisandro: would be great if you can provide some guidance.

Thanks,

Samar

On May 3, 2024, at 12:45 PM, Stefano Zampini  wrote:

While waiting for our Python wizard to shed light on this, I note that, from 
the documentation of PyArray_FROM_OTF 
https://urldefense.us/v3/__https://numpy.org/devdocs/user/c-info.how-to-extend.html*converting-an-arbitrary-sequence-object__;Iw!!G_uCfscf7eWS!YpH_OGif6W8Ql9KZDYBiRMKyALSRhMZOHRnkCkBHqBcVB4ef_UBHFbT_EuTgx8ivGmDNQU8bRVaOc-B_tqQ2heGjdvoE9tp8_QyxtQk$
 , we have

The object can be any Python object convertible to an ndarray. If the object is 
already (a subclass of) the ndarray that satisfies the requirements then a new 
reference is returned.

I guess we should call "del" on the ndarray returned by iarray_s after having 
called  self.set_attr('__array__', array) in this case, but let's wait for 
Lisandro to confirm




Il giorno ven 3 mag 2024 alle ore 11:42 Samar Khatiwala 
mailto:samar.khatiw...@earth.ox.ac.uk>> ha 
scritto:
Hi Matt, Thanks so much for the quick reply! Regarding #2, I put some debug 
statement in my code and what I find is that when I use createWithArray on my 
Cython-allocated numpy array, the destructor I set for it is no longer called 
when I delete
ZjQcmQRYFpfptBannerStart
This Message Is From an External Sender
This message came from outside your organization.

ZjQcmQRYFpfptBannerEnd
Hi Matt,

Thanks so much for the quick reply!

Regarding #2, I put some debug statement in my code and what I find is that 
when I use createWithArray on my Cython-allocated numpy array, the destructor I 
set for it is no longer called when I delete the array. (If I don’t use 
createWithArray then the destructor is triggered.) I interpret that to suggest 
that the petsc4py Vec is somehow ’taking over’ management of the numpy array. 
But I don’t understand where that could be happening. (I don’t think it has to 
do with the actual freeing of memory by PETSc's VecDestroy.)

createWithArray calls iarray_s which in turn calls PyArray_FROM_OTF. Could it 
be there’s something going on there? The numpy documentation is unclear.

Lisandro: do you have any thoughts on this?

Thanks,

Samar

On May 2, 2024, at 11:56 PM, Matthew Knepley 
mailto:knep...@gmail.com>> wrote:

On Thu, May 2, 2024 at 12:53 PM Samar Khatiwala 
mailto:samar.khatiw...@earth.ox.ac.uk>> wrote:
This Message Is From an External Sender
This message came from outside your organization.


Hello,

I have a couple of questions about createWithArray in petsc4py:

1) What is the correct usage for creating a standard MPI Vec with it? Something 
like this seems to work but is it right?:

On each rank do:
a = np.zeros(localSize)
v = PETSc.Vec().createWithArray(a, comm=PETSc.COMM_WORLD)

Is that all it takes?

That looks right to me.

2) Who ‘owns’ the underlying memory for a Vec created with the createWithArray 
method, i.e., who is responsible for managing it and doing garbage collection? 
In my problem, the numpy array is created in a Cython module where memory is 
allocated, and a pointer to it is associated with a numpy ndarray via 
PyArray_SimpleNewFromData and PyArray_SetBaseObject. I have a deallocator 
method of my own that is called when the numpy array is deleted/goes out of 
scope/whenever python does garbage collection. All of that works fine. But if I 
use this array to create a Vec with createWithArray what happens when the Vec 
is, e.g., destroyed? Will my deallocator be called?

No. The PETSc struct will be deallocated, but the storage will not be touched.

  Thanks,

 Matt

Or does petsc4py know that it doesn’t own the memory and won’t attempt to free 
it? I can’t quite figure out from the petsc4py code what is going on. And help 
would be appreciated.

Thanks very much.

Samar




--
What most experimenters take for granted before they begin their experiments is 
infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YpH_OGif6W8Ql9KZDYBiRMKyALSRhMZOHRnkCkBHqBcVB4ef_UBHFbT_EuTgx8ivGmDNQU8bRVaOc-B_tqQ2heGjdvoE9tp8YS3M6pE$
 




--
Stefano

Re: [petsc-users] Question about petsc4py createWithArray function

[Stefano Zampini]

While waiting for our Python wizard to shed light on this, I note that,
from the documentation of PyArray_FROM_OTF
https://urldefense.us/v3/__https://numpy.org/devdocs/user/c-info.how-to-extend.html*converting-an-arbitrary-sequence-object__;Iw!!G_uCfscf7eWS!cPqvd3brsxbc9HOSop19DutpF2hSdn6iPY382KBUd45kJQBA2AyAU5Neifq0WTDf49xH3CybbVSfg7gg6NOaKMOEYdkmd8Y$
 ,
we have

The object can be any Python object convertible to an ndarray. If the
object is already (a subclass of) the ndarray that satisfies the
requirements then a new reference is returned.

I guess we should call "del" on the ndarray returned by iarray_s after
having called  self.set_attr('__array__', array) in this case, but let's
wait for Lisandro to confirm




Il giorno ven 3 mag 2024 alle ore 11:42 Samar Khatiwala <
samar.khatiw...@earth.ox.ac.uk> ha scritto:

> Hi Matt, Thanks so much for the quick reply! Regarding #2, I put some
> debug statement in my code and what I find is that when I use
> createWithArray on my Cython-allocated numpy array, the destructor I set
> for it is no longer called when I delete
> ZjQcmQRYFpfptBannerStart
> This Message Is From an External Sender
> This message came from outside your organization.
>
> ZjQcmQRYFpfptBannerEnd
> Hi Matt,
>
> Thanks so much for the quick reply!
>
> Regarding #2, I put some debug statement in my code and what I find is
> that when I use createWithArray on my Cython-allocated numpy array, the
> destructor I set for it is no longer called when I delete the array. (If I
> don’t use createWithArray then the destructor is triggered.) I interpret
> that to suggest that the petsc4py Vec is somehow ’taking over’ management
> of the numpy array. But I don’t understand where that could be
> happening. (I don’t think it has to do with the actual freeing of memory by
> PETSc's VecDestroy.)
>
> createWithArray calls iarray_s which in turn calls PyArray_FROM_OTF.
> Could it be there’s something going on there? The numpy documentation is
> unclear.
>
> Lisandro: do you have any thoughts on this?
>
> Thanks,
>
> Samar
>
> On May 2, 2024, at 11:56 PM, Matthew Knepley  wrote:
>
> On Thu, May 2, 2024 at 12:53 PM Samar Khatiwala <
> samar.khatiw...@earth.ox.ac.uk> wrote:
>
>> This Message Is From an External Sender
>> This message came from outside your organization.
>>
>>
>> Hello,
>>
>> I have a couple of questions about createWithArray in petsc4py:
>>
>> 1) What is the correct usage for creating a standard MPI Vec with it? 
>> Something like this seems to work but is it right?:
>>
>> On each rank do:
>> a = np.zeros(localSize)
>> v = PETSc.Vec().createWithArray(a, comm=PETSc.COMM_WORLD)
>>
>> Is that all it takes?
>>
>>
> That looks right to me.
>
>> 2) Who ‘owns’ the underlying memory for a Vec created with the 
>> createWithArray method, i.e., who is responsible for managing it and doing 
>> garbage collection? In my problem, the numpy array is created in a Cython 
>> module where memory is allocated, and a pointer to it is associated with a 
>> numpy ndarray via PyArray_SimpleNewFromData and PyArray_SetBaseObject. I 
>> have a deallocator method of my own that is called when the numpy array is 
>> deleted/goes out of scope/whenever python does garbage collection. All of 
>> that works fine. But if I use this array to create a Vec with 
>> createWithArray what happens when the Vec is, e.g., destroyed? Will my 
>> deallocator be called?
>>
>> No. The PETSc struct will be deallocated, but the storage will not be
> touched.
>
>   Thanks,
>
>  Matt
>
>> Or does petsc4py know that it doesn’t own the memory and won’t attempt to 
>> free it? I can’t quite figure out from the petsc4py code what is going on. 
>> And help would be appreciated.
>>
>> Thanks very much.
>>
>> Samar
>>
>>
>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!cPqvd3brsxbc9HOSop19DutpF2hSdn6iPY382KBUd45kJQBA2AyAU5Neifq0WTDf49xH3CybbVSfg7gg6NOaKMOEOy05fZU$
>  
> 
>
>
>

-- 
Stefano

Re: [petsc-users] Question about petsc4py createWithArray function

[Samar Khatiwala]

Hi Matt,

Thanks so much for the quick reply!

Regarding #2, I put some debug statement in my code and what I find is that 
when I use createWithArray on my Cython-allocated numpy array, the destructor I 
set for it is no longer called when I delete the array. (If I don’t use 
createWithArray then the destructor is triggered.) I interpret that to suggest 
that the petsc4py Vec is somehow ’taking over’ management of the numpy array. 
But I don’t understand where that could be happening. (I don’t think it has to 
do with the actual freeing of memory by PETSc's VecDestroy.)

createWithArray calls iarray_s which in turn calls PyArray_FROM_OTF. Could it 
be there’s something going on there? The numpy documentation is unclear.

Lisandro: do you have any thoughts on this?

Thanks,

Samar

On May 2, 2024, at 11:56 PM, Matthew Knepley  wrote:

On Thu, May 2, 2024 at 12:53 PM Samar Khatiwala 
mailto:samar.khatiw...@earth.ox.ac.uk>> wrote:
This Message Is From an External Sender
This message came from outside your organization.


Hello,

I have a couple of questions about createWithArray in petsc4py:

1) What is the correct usage for creating a standard MPI Vec with it? Something 
like this seems to work but is it right?:

On each rank do:
a = np.zeros(localSize)
v = PETSc.Vec().createWithArray(a, comm=PETSc.COMM_WORLD)

Is that all it takes?

That looks right to me.

2) Who ‘owns’ the underlying memory for a Vec created with the createWithArray 
method, i.e., who is responsible for managing it and doing garbage collection? 
In my problem, the numpy array is created in a Cython module where memory is 
allocated, and a pointer to it is associated with a numpy ndarray via 
PyArray_SimpleNewFromData and PyArray_SetBaseObject. I have a deallocator 
method of my own that is called when the numpy array is deleted/goes out of 
scope/whenever python does garbage collection. All of that works fine. But if I 
use this array to create a Vec with createWithArray what happens when the Vec 
is, e.g., destroyed? Will my deallocator be called?

No. The PETSc struct will be deallocated, but the storage will not be touched.

  Thanks,

 Matt

Or does petsc4py know that it doesn’t own the memory and won’t attempt to free 
it? I can’t quite figure out from the petsc4py code what is going on. And help 
would be appreciated.

Thanks very much.

Samar




--
What most experimenters take for granted before they begin their experiments is 
infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!dXZuKR18IO2hZ4uPEqHaG0Jb0ALPgw1T0dm4HAWfyBfsCvkYTmq7aG7yw3yc-5CVhcFMaUL3WzGM8YHlBTTA_JMs_dYUCTX8wgHdY_Y$
 

Re: [petsc-users] Question about petsc4py createWithArray function

[Matthew Knepley]

On Thu, May 2, 2024 at 12:53 PM Samar Khatiwala <
samar.khatiw...@earth.ox.ac.uk> wrote:

> Hello, I have a couple of questions about createWithArray in petsc4py: 1)
> What is the correct usage for creating a standard MPI Vec with it?
> Something like this seems to work but is it right?: On each rank do: a =
> np. zeros(localSize) v = PETSc. Vec(). createWithArray(a,
> ZjQcmQRYFpfptBannerStart
> This Message Is From an External Sender
> This message came from outside your organization.
>
> ZjQcmQRYFpfptBannerEnd
>
> Hello,
>
> I have a couple of questions about createWithArray in petsc4py:
>
> 1) What is the correct usage for creating a standard MPI Vec with it? 
> Something like this seems to work but is it right?:
>
> On each rank do:
> a = np.zeros(localSize)
> v = PETSc.Vec().createWithArray(a, comm=PETSc.COMM_WORLD)
>
> Is that all it takes?
>
>
That looks right to me.

> 2) Who ‘owns’ the underlying memory for a Vec created with the 
> createWithArray method, i.e., who is responsible for managing it and doing 
> garbage collection? In my problem, the numpy array is created in a Cython 
> module where memory is allocated, and a pointer to it is associated with a 
> numpy ndarray via PyArray_SimpleNewFromData and PyArray_SetBaseObject. I have 
> a deallocator method of my own that is called when the numpy array is 
> deleted/goes out of scope/whenever python does garbage collection. All of 
> that works fine. But if I use this array to create a Vec with createWithArray 
> what happens when the Vec is, e.g., destroyed? Will my deallocator be called?
>
> No. The PETSc struct will be deallocated, but the storage will not be
touched.

  Thanks,

 Matt

> Or does petsc4py know that it doesn’t own the memory and won’t attempt to 
> free it? I can’t quite figure out from the petsc4py code what is going on. 
> And help would be appreciated.
>
> Thanks very much.
>
> Samar
>
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!anRKQn_8v127_CSkzk22FfTFRodKT0G2BLwgi_kPAUQt_eqCiySYWSg3ctwewCceXLeJY4FU01ONUG3xfSeD$
  

Re: [petsc-users] Question about how to use DS to do the discretization

[Matthew Knepley]

On Sat, Mar 23, 2024 at 8:03 AM Gong Yujie 
wrote:

> Dear PETSc group, I'm reading the DS part for the discretization start
> from SNES ex17. c which is a demo for solving linear elasticity problem. I
> have two questions for the details. The first question is for the residual
> function. Is the residual
> ZjQcmQRYFpfptBannerStart
> This Message Is From an External Sender
> This message came from outside your organization.
>
> ZjQcmQRYFpfptBannerEnd
> Dear PETSc group,
>
> I'm reading the DS part for the discretization start from SNES ex17.c
> which is a demo for solving linear elasticity problem. I have two questions
> for the details.
>
> The first question is for the residual function. Is the residual
> calculated as this? The dot product is a little weird because of the
> dimension of the result.
> Here \sigma is the stress tensor, \phi_i is the test function for the i-th
> function (Linear elasticity in 3D contains three equations).
>

The stress term in the momentum equation is

  (-div sigma) . psi = d_i sigma_{ij} psi_j

which is then integrated by parts

  sigma_{ij} d_i psi_j

This is linear isotropic elasticity, so

  sigma = \mu (d_i u_j + d_j u_i) + \lambda \delta_{ij} sigma_{kk}

In PETSc, we phrase the term in the weak form as

  f^1_{ij} d_i psi_j

so f1[i * dim + j] below is sigma_{ij}, from line 298 of ex17.c

  for (PetscInt c = 0; c < Nc; ++c) {
for (PetscInt d = 0; d < dim; ++d) {
  f1[c * dim + d] += mu * (u_x[c * dim + d] + u_x[d * dim + c]);
  f1[c * dim + c] += lambda * u_x[d * dim + d];
}
  }


> The second question is how to derive the Jacobian of the system (line 330
> in ex17.c). As shown in the PetscDSSetJacobian, we need to provide function
> g3() which I think is a 4th-order tensor with size 3*3*3*3 in this linear
> elasticity case. I'm not sure how to get it. Are there any references on
> how to get this Jacobian?
>

The Jacobian indices are shown here:
https://urldefense.us/v3/__https://petsc.org/main/manualpages/FE/PetscFEIntegrateJacobian/__;!!G_uCfscf7eWS!eNdXZ9QpzXPTMCtLneeS7td5OYVeKmhdOrw_fswDze2u7v2JaO7kBmFTVakyHDmWJOCHjf-0GrqGc53QKQec$
 
where the g3 term is

  \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u)
\nabla\phi^{gc}_g(q)

To get the Jacobian, we use u = \sum_i u_i psi_i, where psi_i is a vector,
and then differentiate the expression with respect to the coefficient u_i.
Since the operator is already linear,this is just matching indices

  for (PetscInt c = 0; c < Nc; ++c) {
for (PetscInt d = 0; d < dim; ++d) {
  g3[((c * Nc + c) * dim + d) * dim + d] += mu;
  g3[((c * Nc + d) * dim + d) * dim + c] += mu;
  g3[((c * Nc + d) * dim + c) * dim + d] += lambda;
}
  }

Take the first mu term, mu (d_c u_d ) (d_c psi_d). We know that fc == gc
and df == dg, so we get

  g3[((c * Nc + c) * dim + d) * dim + d] += mu;

For the second term, we transpose grad u, so fc == dg and gc == df,

  g3[((c * Nc + d) * dim + d) * dim + c] += mu;

Finally, for the lambda term, fc == df and gc == dg because we are matching
terms in the sum

  g3[((c * Nc + d) * dim + c) * dim + d] += lambda;

  Thanks,

 Matt


> I've checked about the comment before this Jacobian function (line 330 in
> ex17.c) but don't know how to get this.
>
> Thanks in advance!
>
> Best Regards,
> Yujie
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!eNdXZ9QpzXPTMCtLneeS7td5OYVeKmhdOrw_fswDze2u7v2JaO7kBmFTVakyHDmWJOCHjf-0GrqGc4DkrV1c$
  

Re: [petsc-users] Question about how to use DS to do the discretization

[Stefano Zampini]

Take a look at
https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/snes/tutorials/ex11.c?ref_type=heads__;!!G_uCfscf7eWS!er4CI8GIe7OCWvCmRKQpZt6FOz1QYvbuZOdf2Fm7pvMGee3I9M5bhjNytv42F9C17NpBy0i6mTfgEmQfUR_QOqwC7gC6pYk$
 
and the discussion at the beginning (including the reference to the
original paper)

On Sat, Mar 23, 2024, 15:03 Gong Yujie  wrote:

> Dear PETSc group, I'm reading the DS part for the discretization start
> from SNES ex17. c which is a demo for solving linear elasticity problem. I
> have two questions for the details. The first question is for the residual
> function. Is the residual
> ZjQcmQRYFpfptBannerStart
> This Message Is From an External Sender
> This message came from outside your organization.
>
> ZjQcmQRYFpfptBannerEnd
> Dear PETSc group,
>
> I'm reading the DS part for the discretization start from SNES ex17.c
> which is a demo for solving linear elasticity problem. I have two questions
> for the details.
>
> The first question is for the residual function. Is the residual
> calculated as this? The dot product is a little weird because of the
> dimension of the result.
> Here \sigma is the stress tensor, \phi_i is the test function for the i-th
> function (Linear elasticity in 3D contains three equations).
>
> The second question is how to derive the Jacobian of the system (line 330
> in ex17.c). As shown in the PetscDSSetJacobian, we need to provide function
> g3() which I think is a 4th-order tensor with size 3*3*3*3 in this linear
> elasticity case. I'm not sure how to get it. Are there any references on
> how to get this Jacobian?
>
> I've checked about the comment before this Jacobian function (line 330 in
> ex17.c) but don't know how to get this.
>
> Thanks in advance!
>
> Best Regards,
> Yujie
>

Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST

[Zhang, Hong via petsc-users]

Pierre,
Please review 
https://gitlab.com/petsc/petsc/-/merge_requests/7292#5479b9c001ff78d18f5236b7e21d8b1cb9acf476
I moved your checking PetscCheck((*C)->ops->productsymbolic,...) from 
MatProduct_Private() to MatProductSymbolic() .
With this fix, you can run ex195 with
./ex195 -test_nest*nest or ./ex195 -test_matproduct
...
[0]PETSC ERROR: No support for this operation for this object type
[0]PETSC ERROR: Unspecified symbolic phase for product AB with A nest, B nest. 
Call MatProductSetFromOptions() first or the product is not supported

Both tests give the same error output.

In the design of MatProduct routines, we set API calling process as
MatProductCreate()
MatProductSetType()
MatProductSetFromOptions() //determines the implementations of symbolic and 
numeric products, must be called!
MatProductSymbolic()
MatProductNumeric()

Hong


From: Pierre Jolivet 
Sent: Tuesday, February 13, 2024 2:48 PM
To: Zhang, Hong 
Cc: petsc-users@mcs.anl.gov ; Hana Honnerová 

Subject: Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST



On 13 Feb 2024, at 9:21 PM, Zhang, Hong via petsc-users 
 wrote:

Pierre,
I can repeat your change in ex27.c on petsc-release. However, replacing
+Mat D;
+PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DECIDE, ));
+PetscCall(MatDestroy());
with
 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, array, ));
+Mat D;
+PetscCall(MatProductCreate(C, C, NULL, ));
+PetscCall(MatProductSetType(D, MATPRODUCT_AB));
+PetscCall(MatProductSetFromOptions(D));
+PetscCall(MatProductSymbolic(D));
+PetscCall(MatProductNumeric(D));
+PetscCall(MatDestroy());

Sure, I fixed a single MatProduct bug, my guts tell me there are other 
unhandled corner cases, and you indeed found another one.

./ex27 -f  farzad_B_rhs -truncate -solve_augmented
...
[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Unspecified symbolic phase for product AB with A nest, B nest. 
Call MatProductSetFromOptions() first or the product is not supported
...
[0]PETSC ERROR: #1 MatProductSymbolic() at 
/Users/hongzhang-sun/soft/petsc/src/mat/interface/matproduct.c:807
[0]PETSC ERROR: #2 main() at ex27.c:250

i.e., same confusing error message as reported by Hana, because this calling 
process does not call MatProduct_Private() with your fix. A fix to this is to 
modify the error message in MatProductSymbolic():

--- a/src/mat/interface/matproduct.c
+++ b/src/mat/interface/matproduct.c
@@ -804,7 +804,7 @@ PetscErrorCode MatProductSymbolic(Mat mat)
 ...
-PetscCheck(!missing, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, 
"Unspecified symbolic phase for product %s. Call MatProductSetFromOptions() 
first", errstr);
+PetscCheck(!missing, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, 
"Unspecified symbolic phase for product %s. Call MatProductSetFromOptions() 
first or the product is not supported", errstr);

I don’t see how this is less confusing.
In fact, to me, error message with conditionals in the sentence but without 
PETSc telling the value of each expression of the conditionals is infuriating.
How do I know if MatProductSetFromOptions() has not been called, or if the 
product is not supported?
This is very difficult to debug, to me, and if it would be possible to catch 
this with two different checks, it would be much better.
But the current design of MatProduct may not allow us to do it, so I will not 
be opposed to such a change.
Maybe add a Boolean à la pc->setupcalled or pc->setfromoptionscalled in the 
MatProduct structure to be able to distinguish better the cause of the failure?

Thanks,
Pierre

with this fix, I get
./ex27 -f farzad_B_rhs -truncate -solve_augmented
...
[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Unspecified symbolic phase for product AB with A nest, B nest. 
Call MatProductSetFromOptions() first or the product is not supported

If you agree with this fix, I'll create a MR for it.
Hong



From: Pierre Jolivet 
Sent: Tuesday, February 13, 2024 12:08 AM
To: Zhang, Hong 
Cc: Hana Honnerová ; petsc-users@mcs.anl.gov 

Subject: Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST



On 13 Feb 2024, at 12:33 AM, Zhang, Hong  wrote:

Pierre,
I just modified the error message of MatProductSymbolic() and added a testing 
segment in src/mat/tests/ex195.c. I have not pushed my change yet.

Your fix at 
https://gitlab.com/petsc/petsc/-/commit/9dcea022de3b0309e5c16b8c554ad9c85dea29cf?merge_request_iid=7283
 is more general. Has this fix merged to release and main? With latest main and 
release, I get same previous error message.

I don’t (anymore, but could prior to my fix).
The trigger is MatMatMult() with MAT_INITIAL_MATRIX in PCSetUp_LSC().
Reproducible with:
diff --git a/src/ksp/ksp/tutorials/ex27.c b/src/ksp/ksp/tutorials/ex27.c
index 116b7df8522..9bdf4d7334a 10

Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST

[Pierre Jolivet]

> On 13 Feb 2024, at 9:21 PM, Zhang, Hong via petsc-users 
>  wrote:
> 
> Pierre,
> I can repeat your change in ex27.c on petsc-release. However, replacing 
> +Mat D;
> +PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DECIDE, ));
> +PetscCall(MatDestroy());
> with
>  PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, array, ));
> +Mat D;
> +PetscCall(MatProductCreate(C, C, NULL, ));
> +PetscCall(MatProductSetType(D, MATPRODUCT_AB));
> +PetscCall(MatProductSetFromOptions(D));
> +PetscCall(MatProductSymbolic(D));
> +PetscCall(MatProductNumeric(D));
> +PetscCall(MatDestroy());

Sure, I fixed a single MatProduct bug, my guts tell me there are other 
unhandled corner cases, and you indeed found another one.

> ./ex27 -f  farzad_B_rhs -truncate -solve_augmented
> ...
> [0]PETSC ERROR: Petsc has generated inconsistent data
> [0]PETSC ERROR: Unspecified symbolic phase for product AB with A nest, B 
> nest. Call MatProductSetFromOptions() first or the product is not supported
> ...
> [0]PETSC ERROR: #1 MatProductSymbolic() at 
> /Users/hongzhang-sun/soft/petsc/src/mat/interface/matproduct.c:807
> [0]PETSC ERROR: #2 main() at ex27.c:250
> 
> i.e., same confusing error message as reported by Hana, because this calling 
> process does not call MatProduct_Private() with your fix. A fix to this is to 
> modify the error message in MatProductSymbolic():
> 
> --- a/src/mat/interface/matproduct.c
> +++ b/src/mat/interface/matproduct.c
> @@ -804,7 +804,7 @@ PetscErrorCode MatProductSymbolic(Mat mat)
>  ...
> -PetscCheck(!missing, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, 
> "Unspecified symbolic phase for product %s. Call MatProductSetFromOptions() 
> first", errstr);
> +PetscCheck(!missing, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, 
> "Unspecified symbolic phase for product %s. Call MatProductSetFromOptions() 
> first or the product is not supported", errstr);

I don’t see how this is less confusing.
In fact, to me, error message with conditionals in the sentence but without 
PETSc telling the value of each expression of the conditionals is infuriating.
How do I know if MatProductSetFromOptions() has not been called, or if the 
product is not supported?
This is very difficult to debug, to me, and if it would be possible to catch 
this with two different checks, it would be much better.
But the current design of MatProduct may not allow us to do it, so I will not 
be opposed to such a change.
Maybe add a Boolean à la pc->setupcalled or pc->setfromoptionscalled in the 
MatProduct structure to be able to distinguish better the cause of the failure?

Thanks,
Pierre

> with this fix, I get 
> ./ex27 -f farzad_B_rhs -truncate -solve_augmented
> ...
> [0]PETSC ERROR: Petsc has generated inconsistent data
> [0]PETSC ERROR: Unspecified symbolic phase for product AB with A nest, B 
> nest. Call MatProductSetFromOptions() first or the product is not supported
> 
> If you agree with this fix, I'll create a MR for it.
> Hong
> 
> 
> From: Pierre Jolivet 
> Sent: Tuesday, February 13, 2024 12:08 AM
> To: Zhang, Hong 
> Cc: Hana Honnerová ; petsc-users@mcs.anl.gov 
> 
> Subject: Re: [petsc-users] question on PCLSC with matrix blocks of type 
> MATNEST
>  
> 
> 
>> On 13 Feb 2024, at 12:33 AM, Zhang, Hong  wrote:
>> 
>> Pierre,
>> I just modified the error message of MatProductSymbolic() and added a 
>> testing segment in src/mat/tests/ex195.c. I have not pushed my change yet.
>> 
>> Your fix at 
>> https://gitlab.com/petsc/petsc/-/commit/9dcea022de3b0309e5c16b8c554ad9c85dea29cf?merge_request_iid=7283
>>  is more general. Has this fix merged to release and main? With latest main 
>> and release, I get same previous error message.
> 
> I don’t (anymore, but could prior to my fix).
> The trigger is MatMatMult() with MAT_INITIAL_MATRIX in PCSetUp_LSC().
> Reproducible with:
> diff --git a/src/ksp/ksp/tutorials/ex27.c b/src/ksp/ksp/tutorials/ex27.c
> index 116b7df8522..9bdf4d7334a 100644
> --- a/src/ksp/ksp/tutorials/ex27.c
> +++ b/src/ksp/ksp/tutorials/ex27.c
> @@ -245,2 +245,5 @@ int main(int argc, char **args)
>  PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, array, ));
> +Mat D;
> +PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DECIDE, ));
> +PetscCall(MatDestroy());
>  if (!sbaij) PetscCall(MatNestSetVecType(C, VECNEST));
> Which now generates:
> $ ./ex27 -f ${DATAFILESPATH}/matrices/farzad_B_rhs -truncate -solve_augmented 
> Failed to load RHS, so use a vector of all ones.
> Failed to load initial guess, so use a vector of all zeros.
> [0]PETSC ERROR: - Error Messa

Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST

[Zhang, Hong via petsc-users]

Pierre,
I can repeat your change in ex27.c on petsc-release. However, replacing
+Mat D;
+PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DECIDE, ));
+PetscCall(MatDestroy());
with
 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, array, ));
+Mat D;
+PetscCall(MatProductCreate(C, C, NULL, ));
+PetscCall(MatProductSetType(D, MATPRODUCT_AB));
+PetscCall(MatProductSetFromOptions(D));
+PetscCall(MatProductSymbolic(D));
+PetscCall(MatProductNumeric(D));
+PetscCall(MatDestroy());
./ex27 -f  farzad_B_rhs -truncate -solve_augmented
...
[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Unspecified symbolic phase for product AB with A nest, B nest. 
Call MatProductSetFromOptions() first or the product is not supported
...
[0]PETSC ERROR: #1 MatProductSymbolic() at 
/Users/hongzhang-sun/soft/petsc/src/mat/interface/matproduct.c:807
[0]PETSC ERROR: #2 main() at ex27.c:250

i.e., same confusing error message as reported by Hana, because this calling 
process does not call MatProduct_Private() with your fix. A fix to this is to 
modify the error message in MatProductSymbolic():

--- a/src/mat/interface/matproduct.c
+++ b/src/mat/interface/matproduct.c
@@ -804,7 +804,7 @@ PetscErrorCode MatProductSymbolic(Mat mat)
 ...
-PetscCheck(!missing, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, 
"Unspecified symbolic phase for product %s. Call MatProductSetFromOptions() 
first", errstr);
+PetscCheck(!missing, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, 
"Unspecified symbolic phase for product %s. Call MatProductSetFromOptions() 
first or the product is not supported", errstr);

with this fix, I get
./ex27 -f farzad_B_rhs -truncate -solve_augmented
...
[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Unspecified symbolic phase for product AB with A nest, B nest. 
Call MatProductSetFromOptions() first or the product is not supported

If you agree with this fix, I'll create a MR for it.
Hong



From: Pierre Jolivet 
Sent: Tuesday, February 13, 2024 12:08 AM
To: Zhang, Hong 
Cc: Hana Honnerová ; petsc-users@mcs.anl.gov 

Subject: Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST



On 13 Feb 2024, at 12:33 AM, Zhang, Hong  wrote:

Pierre,
I just modified the error message of MatProductSymbolic() and added a testing 
segment in src/mat/tests/ex195.c. I have not pushed my change yet.

Your fix at 
https://gitlab.com/petsc/petsc/-/commit/9dcea022de3b0309e5c16b8c554ad9c85dea29cf?merge_request_iid=7283
 is more general. Has this fix merged to release and main? With latest main and 
release, I get same previous error message.

I don’t (anymore, but could prior to my fix).
The trigger is MatMatMult() with MAT_INITIAL_MATRIX in PCSetUp_LSC().
Reproducible with:
diff --git a/src/ksp/ksp/tutorials/ex27.c b/src/ksp/ksp/tutorials/ex27.c
index 116b7df8522..9bdf4d7334a 100644
--- a/src/ksp/ksp/tutorials/ex27.c
+++ b/src/ksp/ksp/tutorials/ex27.c
@@ -245,2 +245,5 @@ int main(int argc, char **args)
 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, array, ));
+Mat D;
+PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DECIDE, ));
+PetscCall(MatDestroy());
 if (!sbaij) PetscCall(MatNestSetVecType(C, VECNEST));
Which now generates:
$ ./ex27 -f ${DATAFILESPATH}/matrices/farzad_B_rhs -truncate -solve_augmented
Failed to load RHS, so use a vector of all ones.
Failed to load initial guess, so use a vector of all zeros.
[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR: No support for this operation for this object type
[0]PETSC ERROR: MatProduct AB not supported for nest and nest
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.

Thanks,
Pierre

Hong


From: Pierre Jolivet 
Sent: Sunday, February 11, 2024 7:43 AM
To: Zhang, Hong 
Cc: Hana Honnerová ; petsc-users@mcs.anl.gov 

Subject: Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST


On 8 Feb 2024, at 5:37 PM, Zhang, Hong via petsc-users 
 wrote:

Hana,
"product AB with A nest, B nest" is not supported by PETSc. I do not know why 
PETSc does not display such an error message. I'll check it.

Did you?
A naive fix is to simply add the missing PetscCheck() in MatProduct_Private() 
right after 
MatProductSetFromOptions()https://petsc.org/release/src/mat/interface/matrix.c.html#line10026
 (notice that it is there line 10048 in the other code branch)
I have this at 
https://gitlab.com/petsc/petsc/-/commit/9dcea022de3b0309e5c16b8c554ad9c85dea29cf?merge_request_iid=7283
 (coupled with code refactoring to avoid missing any other operations), but 
maybe we could do things more elegantly.

Thanks,
Pierre

Hong

From: petsc-users  on behalf of Hana Honnerová 

Sent: Thursday, February 8

Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST

[Pierre Jolivet]

> On 13 Feb 2024, at 12:33 AM, Zhang, Hong  wrote:
> 
> Pierre,
> I just modified the error message of MatProductSymbolic() and added a testing 
> segment in src/mat/tests/ex195.c. I have not pushed my change yet.
> 
> Your fix at 
> https://gitlab.com/petsc/petsc/-/commit/9dcea022de3b0309e5c16b8c554ad9c85dea29cf?merge_request_iid=7283
>  is more general. Has this fix merged to release and main? With latest main 
> and release, I get same previous error message.

I don’t (anymore, but could prior to my fix).
The trigger is MatMatMult() with MAT_INITIAL_MATRIX in PCSetUp_LSC().
Reproducible with:
diff --git a/src/ksp/ksp/tutorials/ex27.c b/src/ksp/ksp/tutorials/ex27.c
index 116b7df8522..9bdf4d7334a 100644
--- a/src/ksp/ksp/tutorials/ex27.c
+++ b/src/ksp/ksp/tutorials/ex27.c
@@ -245,2 +245,5 @@ int main(int argc, char **args)
 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, array, ));
+Mat D;
+PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DECIDE, ));
+PetscCall(MatDestroy());
 if (!sbaij) PetscCall(MatNestSetVecType(C, VECNEST));
Which now generates:
$ ./ex27 -f ${DATAFILESPATH}/matrices/farzad_B_rhs -truncate -solve_augmented 
Failed to load RHS, so use a vector of all ones.
Failed to load initial guess, so use a vector of all zeros.
[0]PETSC ERROR: - Error Message 
--
[0]PETSC ERROR: No support for this operation for this object type
[0]PETSC ERROR: MatProduct AB not supported for nest and nest
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.

Thanks,
Pierre

> Hong
> 
> From: Pierre Jolivet 
> Sent: Sunday, February 11, 2024 7:43 AM
> To: Zhang, Hong 
> Cc: Hana Honnerová ; petsc-users@mcs.anl.gov 
> 
> Subject: Re: [petsc-users] question on PCLSC with matrix blocks of type 
> MATNEST
>  
> 
> On 8 Feb 2024, at 5:37 PM, Zhang, Hong via petsc-users 
>  wrote:
> 
> Hana,
> "product AB with A nest, B nest" is not supported by PETSc. I do not know why 
> PETSc does not display such an error message. I'll check it.
> 
> Did you?
> A naive fix is to simply add the missing PetscCheck() in MatProduct_Private() 
> right after 
> MatProductSetFromOptions()https://petsc.org/release/src/mat/interface/matrix.c.html#line10026
>  (notice that it is there line 10048 in the other code branch)
> I have this at 
> https://gitlab.com/petsc/petsc/-/commit/9dcea022de3b0309e5c16b8c554ad9c85dea29cf?merge_request_iid=7283
>  (coupled with code refactoring to avoid missing any other operations), but 
> maybe we could do things more elegantly.
> 
> Thanks,
> Pierre
> 
> Hong
> From: petsc-users  on behalf of Hana 
> Honnerová 
> Sent: Thursday, February 8, 2024 4:45 AM
> To: petsc-users@mcs.anl.gov 
> Subject: [petsc-users] question on PCLSC with matrix blocks of type MATNEST
>  
> Hi all,
> I am trying to solve linear systems arising from isogeometric discretization 
> (similar to FEM) of the Navier-Stokes equations in parallel using PETSc. The 
> linear systems are of saddle-point type, so I would like to use the 
> PCFIELDSPLIT preconditioner with the -pc_fieldsplit_detect_saddle_point 
> option, Schur complement factorization and the LSC Schur complement 
> preconditioner. I do not provide any user-defined operators for PCLSC in my 
> codes (at least for now).
> I store the matrix as a MATNEST consisting of 4 blocks (F for 
> velocity-velocity part, Bt for velocity-pressure part, B for 
> pressure-velocity part and NULL for pressure-pressure part). It is also 
> convenient for me to store the blocks F, Bt and B as another MATNEST 
> consisting of blocks corresponding to individual velocity components. 
> 
> However, in this setting, I get the following error message:
> [0]PETSC ERROR: Petsc has generated inconsistent data
> [0]PETSC ERROR: Unspecified symbolic phase for product AB with A nest, B 
> nest. Call MatProductSetFromOptions() first
> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.20.4, unknown 
> [0]PETSC ERROR: 
> /home/hhornik/gismo/build-petsc-mpi/RelWithDebInfo/bin/gsINSSolverPETScTest 
> on a arch-debug named ThinkPad-T14 by hhornik Thu Feb  8 11:04:17 2024
> [0]PETSC ERROR: Configure options PETSC_ARCH=arch-debug --with-debugging=1 
> --download-mumps --download-scalapack
> [0]PETSC ERROR: #1 MatProductSymbolic() at 
> /home/hhornik/Software/PETSc/src/mat/interface/matproduct.c:807
> [0]PETSC ERROR: #2 MatProduct_Private() at 
> /home/hhornik/Software/PETSc/src/mat/interface/matrix.c:10027
> [0]PETSC ERROR: #3 MatMatMult() at 
> /home/hhornik/Software/PETSc/src/mat/interface/matrix.c:10103
> [0]PETSC ERROR: #4 PCSetUp_LSC() at 
> /home/hhornik/Softw

Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST

[Zhang, Hong via petsc-users]

Pierre,
I just modified the error message of MatProductSymbolic() and added a testing 
segment in src/mat/tests/ex195.c. I have not pushed my change yet.

Your fix at 
https://gitlab.com/petsc/petsc/-/commit/9dcea022de3b0309e5c16b8c554ad9c85dea29cf?merge_request_iid=7283
 is more general. Has this fix merged to release and main? With latest main and 
release, I get same previous error message.
Hong


From: Pierre Jolivet 
Sent: Sunday, February 11, 2024 7:43 AM
To: Zhang, Hong 
Cc: Hana Honnerová ; petsc-users@mcs.anl.gov 

Subject: Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST


On 8 Feb 2024, at 5:37 PM, Zhang, Hong via petsc-users 
 wrote:

Hana,
"product AB with A nest, B nest" is not supported by PETSc. I do not know why 
PETSc does not display such an error message. I'll check it.

Did you?
A naive fix is to simply add the missing PetscCheck() in MatProduct_Private() 
right after MatProductSetFromOptions() 
https://petsc.org/release/src/mat/interface/matrix.c.html#line10026 (notice 
that it is there line 10048 in the other code branch)
I have this at 
https://gitlab.com/petsc/petsc/-/commit/9dcea022de3b0309e5c16b8c554ad9c85dea29cf?merge_request_iid=7283
 (coupled with code refactoring to avoid missing any other operations), but 
maybe we could do things more elegantly.

Thanks,
Pierre

Hong

From: petsc-users  on behalf of Hana Honnerová 

Sent: Thursday, February 8, 2024 4:45 AM
To: petsc-users@mcs.anl.gov 
Subject: [petsc-users] question on PCLSC with matrix blocks of type MATNEST

Hi all,
I am trying to solve linear systems arising from isogeometric discretization 
(similar to FEM) of the Navier-Stokes equations in parallel using PETSc. The 
linear systems are of saddle-point type, so I would like to use the 
PCFIELDSPLIT preconditioner with the -pc_fieldsplit_detect_saddle_point option, 
Schur complement factorization and the LSC Schur complement preconditioner. I 
do not provide any user-defined operators for PCLSC in my codes (at least for 
now).
I store the matrix as a MATNEST consisting of 4 blocks (F for velocity-velocity 
part, Bt for velocity-pressure part, B for pressure-velocity part and NULL for 
pressure-pressure part). It is also convenient for me to store the blocks F, Bt 
and B as another MATNEST consisting of blocks corresponding to individual 
velocity components.

However, in this setting, I get the following error message:
[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Unspecified symbolic phase for product AB with A nest, B nest. 
Call MatProductSetFromOptions() first
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.20.4, unknown
[0]PETSC ERROR: 
/home/hhornik/gismo/build-petsc-mpi/RelWithDebInfo/bin/gsINSSolverPETScTest on 
a arch-debug named ThinkPad-T14 by hhornik Thu Feb  8 11:04:17 2024
[0]PETSC ERROR: Configure options PETSC_ARCH=arch-debug --with-debugging=1 
--download-mumps --download-scalapack
[0]PETSC ERROR: #1 MatProductSymbolic() at 
/home/hhornik/Software/PETSc/src/mat/interface/matproduct.c:807
[0]PETSC ERROR: #2 MatProduct_Private() at 
/home/hhornik/Software/PETSc/src/mat/interface/matrix.c:10027
[0]PETSC ERROR: #3 MatMatMult() at 
/home/hhornik/Software/PETSc/src/mat/interface/matrix.c:10103
[0]PETSC ERROR: #4 PCSetUp_LSC() at 
/home/hhornik/Software/PETSc/src/ksp/pc/impls/lsc/lsc.c:79
[0]PETSC ERROR: #5 PCSetUp() at 
/home/hhornik/Software/PETSc/src/ksp/pc/interface/precon.c:1080
[0]PETSC ERROR: #6 KSPSetUp() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:415
[0]PETSC ERROR: #7 KSPSolve_Private() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:832
[0]PETSC ERROR: #8 KSPSolve() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:1079
[0]PETSC ERROR: #9 PCApply_FieldSplit_Schur() at 
/home/hhornik/Software/PETSc/src/ksp/pc/impls/fieldsplit/fieldsplit.c:1165
[0]PETSC ERROR: #10 PCApply() at 
/home/hhornik/Software/PETSc/src/ksp/pc/interface/precon.c:498
[0]PETSC ERROR: #11 KSP_PCApply() at 
/home/hhornik/Software/PETSc/include/petsc/private/kspimpl.h:383
[0]PETSC ERROR: #12 KSPFGMRESCycle() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/impls/gmres/fgmres/fgmres.c:123
[0]PETSC ERROR: #13 KSPSolve_FGMRES() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/impls/gmres/fgmres/fgmres.c:235
[0]PETSC ERROR: #14 KSPSolve_Private() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:906
[0]PETSC ERROR: #15 KSPSolve() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:1079
[0]PETSC ERROR: #16 applySolver() at 
/home/hhornik/gismo/optional/gsIncompressibleFlow/src/gsINSSolver.hpp:531

I could not find any solution for this so far. My question is: Is it possible 
to use the LSC preconditioner in such case, where the matrix blocks are of type 
MATNEST? And if so, how?
Thank you,
Hana Honnerova

Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST

[Pierre Jolivet]

> On 8 Feb 2024, at 5:37 PM, Zhang, Hong via petsc-users 
>  wrote:
> 
> Hana,
> "product AB with A nest, B nest" is not supported by PETSc. I do not know why 
> PETSc does not display such an error message. I'll check it.

Did you?
A naive fix is to simply add the missing PetscCheck() in MatProduct_Private() 
right after MatProductSetFromOptions() 
https://petsc.org/release/src/mat/interface/matrix.c.html#line10026 (notice 
that it is there line 10048 in the other code branch)
I have this at 
https://gitlab.com/petsc/petsc/-/commit/9dcea022de3b0309e5c16b8c554ad9c85dea29cf?merge_request_iid=7283
 (coupled with code refactoring to avoid missing any other operations), but 
maybe we could do things more elegantly.

Thanks,
Pierre

> Hong
> From: petsc-users  on behalf of Hana 
> Honnerová 
> Sent: Thursday, February 8, 2024 4:45 AM
> To: petsc-users@mcs.anl.gov 
> Subject: [petsc-users] question on PCLSC with matrix blocks of type MATNEST
>  
> Hi all,
> I am trying to solve linear systems arising from isogeometric discretization 
> (similar to FEM) of the Navier-Stokes equations in parallel using PETSc. The 
> linear systems are of saddle-point type, so I would like to use the 
> PCFIELDSPLIT preconditioner with the -pc_fieldsplit_detect_saddle_point 
> option, Schur complement factorization and the LSC Schur complement 
> preconditioner. I do not provide any user-defined operators for PCLSC in my 
> codes (at least for now).
> I store the matrix as a MATNEST consisting of 4 blocks (F for 
> velocity-velocity part, Bt for velocity-pressure part, B for 
> pressure-velocity part and NULL for pressure-pressure part). It is also 
> convenient for me to store the blocks F, Bt and B as another MATNEST 
> consisting of blocks corresponding to individual velocity components. 
> 
> However, in this setting, I get the following error message:
> [0]PETSC ERROR: Petsc has generated inconsistent data
> [0]PETSC ERROR: Unspecified symbolic phase for product AB with A nest, B 
> nest. Call MatProductSetFromOptions() first
> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.20.4, unknown 
> [0]PETSC ERROR: 
> /home/hhornik/gismo/build-petsc-mpi/RelWithDebInfo/bin/gsINSSolverPETScTest 
> on a arch-debug named ThinkPad-T14 by hhornik Thu Feb  8 11:04:17 2024
> [0]PETSC ERROR: Configure options PETSC_ARCH=arch-debug --with-debugging=1 
> --download-mumps --download-scalapack
> [0]PETSC ERROR: #1 MatProductSymbolic() at 
> /home/hhornik/Software/PETSc/src/mat/interface/matproduct.c:807
> [0]PETSC ERROR: #2 MatProduct_Private() at 
> /home/hhornik/Software/PETSc/src/mat/interface/matrix.c:10027
> [0]PETSC ERROR: #3 MatMatMult() at 
> /home/hhornik/Software/PETSc/src/mat/interface/matrix.c:10103
> [0]PETSC ERROR: #4 PCSetUp_LSC() at 
> /home/hhornik/Software/PETSc/src/ksp/pc/impls/lsc/lsc.c:79
> [0]PETSC ERROR: #5 PCSetUp() at 
> /home/hhornik/Software/PETSc/src/ksp/pc/interface/precon.c:1080
> [0]PETSC ERROR: #6 KSPSetUp() at 
> /home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:415
> [0]PETSC ERROR: #7 KSPSolve_Private() at 
> /home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:832
> [0]PETSC ERROR: #8 KSPSolve() at 
> /home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:1079
> [0]PETSC ERROR: #9 PCApply_FieldSplit_Schur() at 
> /home/hhornik/Software/PETSc/src/ksp/pc/impls/fieldsplit/fieldsplit.c:1165
> [0]PETSC ERROR: #10 PCApply() at 
> /home/hhornik/Software/PETSc/src/ksp/pc/interface/precon.c:498
> [0]PETSC ERROR: #11 KSP_PCApply() at 
> /home/hhornik/Software/PETSc/include/petsc/private/kspimpl.h:383
> [0]PETSC ERROR: #12 KSPFGMRESCycle() at 
> /home/hhornik/Software/PETSc/src/ksp/ksp/impls/gmres/fgmres/fgmres.c:123
> [0]PETSC ERROR: #13 KSPSolve_FGMRES() at 
> /home/hhornik/Software/PETSc/src/ksp/ksp/impls/gmres/fgmres/fgmres.c:235
> [0]PETSC ERROR: #14 KSPSolve_Private() at 
> /home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:906
> [0]PETSC ERROR: #15 KSPSolve() at 
> /home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:1079
> [0]PETSC ERROR: #16 applySolver() at 
> /home/hhornik/gismo/optional/gsIncompressibleFlow/src/gsINSSolver.hpp:531
> 
> I could not find any solution for this so far. My question is: Is it possible 
> to use the LSC preconditioner in such case, where the matrix blocks are of 
> type MATNEST? And if so, how?
> Thank you,
> Hana Honnerova

Re: [petsc-users] question on PCLSC with matrix blocks of type MATNEST

[Zhang, Hong via petsc-users]

Hana,
"product AB with A nest, B nest" is not supported by PETSc. I do not know why 
PETSc does not display such an error message. I'll check it.
Hong

From: petsc-users  on behalf of Hana Honnerová 

Sent: Thursday, February 8, 2024 4:45 AM
To: petsc-users@mcs.anl.gov 
Subject: [petsc-users] question on PCLSC with matrix blocks of type MATNEST


Hi all,

I am trying to solve linear systems arising from isogeometric discretization 
(similar to FEM) of the Navier-Stokes equations in parallel using PETSc. The 
linear systems are of saddle-point type, so I would like to use the 
PCFIELDSPLIT preconditioner with the -pc_fieldsplit_detect_saddle_point option, 
Schur complement factorization and the LSC Schur complement preconditioner. I 
do not provide any user-defined operators for PCLSC in my codes (at least for 
now).

I store the matrix as a MATNEST consisting of 4 blocks (F for velocity-velocity 
part, Bt for velocity-pressure part, B for pressure-velocity part and NULL for 
pressure-pressure part). It is also convenient for me to store the blocks F, Bt 
and B as another MATNEST consisting of blocks corresponding to individual 
velocity components.

However, in this setting, I get the following error message:

[0]PETSC ERROR: Petsc has generated inconsistent data
[0]PETSC ERROR: Unspecified symbolic phase for product AB with A nest, B nest. 
Call MatProductSetFromOptions() first
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.20.4, unknown
[0]PETSC ERROR: 
/home/hhornik/gismo/build-petsc-mpi/RelWithDebInfo/bin/gsINSSolverPETScTest on 
a arch-debug named ThinkPad-T14 by hhornik Thu Feb  8 11:04:17 2024
[0]PETSC ERROR: Configure options PETSC_ARCH=arch-debug --with-debugging=1 
--download-mumps --download-scalapack
[0]PETSC ERROR: #1 MatProductSymbolic() at 
/home/hhornik/Software/PETSc/src/mat/interface/matproduct.c:807
[0]PETSC ERROR: #2 MatProduct_Private() at 
/home/hhornik/Software/PETSc/src/mat/interface/matrix.c:10027
[0]PETSC ERROR: #3 MatMatMult() at 
/home/hhornik/Software/PETSc/src/mat/interface/matrix.c:10103
[0]PETSC ERROR: #4 PCSetUp_LSC() at 
/home/hhornik/Software/PETSc/src/ksp/pc/impls/lsc/lsc.c:79
[0]PETSC ERROR: #5 PCSetUp() at 
/home/hhornik/Software/PETSc/src/ksp/pc/interface/precon.c:1080
[0]PETSC ERROR: #6 KSPSetUp() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:415
[0]PETSC ERROR: #7 KSPSolve_Private() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:832
[0]PETSC ERROR: #8 KSPSolve() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:1079
[0]PETSC ERROR: #9 PCApply_FieldSplit_Schur() at 
/home/hhornik/Software/PETSc/src/ksp/pc/impls/fieldsplit/fieldsplit.c:1165
[0]PETSC ERROR: #10 PCApply() at 
/home/hhornik/Software/PETSc/src/ksp/pc/interface/precon.c:498
[0]PETSC ERROR: #11 KSP_PCApply() at 
/home/hhornik/Software/PETSc/include/petsc/private/kspimpl.h:383
[0]PETSC ERROR: #12 KSPFGMRESCycle() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/impls/gmres/fgmres/fgmres.c:123
[0]PETSC ERROR: #13 KSPSolve_FGMRES() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/impls/gmres/fgmres/fgmres.c:235
[0]PETSC ERROR: #14 KSPSolve_Private() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:906
[0]PETSC ERROR: #15 KSPSolve() at 
/home/hhornik/Software/PETSc/src/ksp/ksp/interface/itfunc.c:1079
[0]PETSC ERROR: #16 applySolver() at 
/home/hhornik/gismo/optional/gsIncompressibleFlow/src/gsINSSolver.hpp:531

I could not find any solution for this so far. My question is: Is it possible 
to use the LSC preconditioner in such case, where the matrix blocks are of type 
MATNEST? And if so, how?

Thank you,
Hana Honnerova

Re: [petsc-users] Question on using petsc matrix with hip

[Runfeng Jin]

Thank you!

Junchao Zhang  于2024年2月5日周一 12:30写道：

>
>
>
> On Sun, Feb 4, 2024 at 9:57 PM Runfeng Jin  wrote:
>
>> Hi,
>>
>> I see in document the mat status in hip is still in development, but
>> I see some hip interface in MAT. I want to use petsc with primme(eigenvalue
>> solver) to diagonalize the matrix in AMD GPU. The primme need user to
>> provide MV (matrix-vector operation in y=AX) in GPU, and I want to use
>> PETSc to realize it. Can PETSc support this in hip?  thank you!
>>
> petsc/hip backend supports sparse matrix vector multiplication (via AMD
> hipsparce).  The matrix has to be in the CSR format (i.e., MATAIJ in
> petsc).
>
>>
>> Runfeng Jin
>>
>

Re: [petsc-users] Question on using petsc matrix with hip

[Junchao Zhang]

On Sun, Feb 4, 2024 at 9:57 PM Runfeng Jin  wrote:

> Hi,
>
> I see in document the mat status in hip is still in development, but I
> see some hip interface in MAT. I want to use petsc with primme(eigenvalue
> solver) to diagonalize the matrix in AMD GPU. The primme need user to
> provide MV (matrix-vector operation in y=AX) in GPU, and I want to use
> PETSc to realize it. Can PETSc support this in hip?  thank you!
>
petsc/hip backend supports sparse matrix vector multiplication (via AMD
hipsparce).  The matrix has to be in the CSR format (i.e., MATAIJ in
petsc).

>
> Runfeng Jin
>

Re: [petsc-users] Question about a parallel implementation of PCFIELDSPLIT

[Matthew Knepley]

On Tue, Jan 23, 2024 at 11:06 AM Pantelis Moschopoulos <
pmoschopou...@outlook.com> wrote:

> Dear Matt,
>
> Thank you for your explanation. The new methodology is straightforward to
> implement.
> Still, I have one more question . When I use the option
> -pc_fieldsplit_schur_precondition full, PETSc computes internally the exact
> Schur complement matrix representation. Based on the example matrix that
> you send, the Schur complement is: S = -v^t (A^-1) u. How will PETSc will
> calculate the vector (A^-1) u ? Or it calculates the exact Schur complement
> matrix differently?
>

FULL calls MatSchurComplementComputeExplicitOperator(), which
calls KSPMatSolve() to compute A^{-1} u, which default to KSPSolve for each
column if no specialized code is available. So it should just use your
solver for the (0,0) block, and then take the dot product with v.

  Thanks,

 Matt


> Thanks,
> Pantelis
> --
> *Από:* Matthew Knepley 
> *Στάλθηκε:* Τρίτη, 23 Ιανουαρίου 2024 5:21 μμ
> *Προς:* Pantelis Moschopoulos 
> *Κοιν.:* Barry Smith ; petsc-users@mcs.anl.gov <
> petsc-users@mcs.anl.gov>
> *Θέμα:* Re: [petsc-users] Question about a parallel implementation of
> PCFIELDSPLIT
>
> On Tue, Jan 23, 2024 at 9:45 AM Pantelis Moschopoulos <
> pmoschopou...@outlook.com> wrote:
>
> Dear Matt,
>
> I read about the MATLRC. However, its correct usage is not clear to me so
> I have the following questions:
>
>1. The U and V input matrices should be created as dense using
>MatCreateDense?
>
>
> Yes. If you have one row, it looks like a vector, or a matrix with one
> column.
>
> If you have 1 row on the bottom, then
>
>   U = [0, 0, ..., 0, 1]
>   V = [the row]
>   C = [1]
>
> will give you that. However, you have an extra row and column?
>
>
>
>1. I use the command MatCreateLRC just to declare the matrix and then
>MatLRCSetMats to pass the values of the constituents?
>
> You can use
>
>   MatCreate(comm, )
>   MatSetSizes(M, ...)
>   MatSetType(M, MATLRC)
>   MatLRCSetMats(M, ...)
>
> However, you are right that it is a little more complicated, because A is
> not just the upper block here.
>
>
>1.
>
> Then, how do I proceed? How I apply the step of Sherman-Morrison-Woobury
> formula? I intend to use iterative solvers for A (main matrix) so I will
> not have its A^-1 at hand which I think is what the
> Sherman-Morrison-Woobury formula needs.
>
>
> I think I was wrong. MatLRC is not the best fit. We should use MatNest
> instead. Then you could have
>
>   A   u
>   v^t 0
>
> as your matrix. We could still get an explicit Schur complement
> automatically using nested FieldSplit. So
>
>   1) Make the MATNEST matrix as shown above
>
>   2) Use PCFIELDSPLIT  for it. This should be an easy IS, since there is
> only one row in the second one.
>
>   3) Select the full Schur complement
>
> -pc_type fieldsplit
> -pc_fieldsplit_type schur
> -pc_fieldsplit_schur_fact_type full
> -pc_fieldsplit_schur_precondition full
>
>   4) Use a recursive FieldSplit (might be able to use
> -fieldsplit_0_pc_fieldsplit_detect_saddle_point)
>
> -fieldsplit_0_pc_type fieldsplit
> -fieldsplit_0_pc_fieldsplit_0_fields 0,1,2
> -fieldsplit_0_pc_fieldsplit_1_fields 3
>
> I think this does what you want, and should be much easier than getting
> MatLRC to do it.
>
>   Thanks,
>
>  Matt
>
>
> Thanks,
> Pantelis
> --
> *Από:* Matthew Knepley 
> *Στάλθηκε:* Τρίτη, 23 Ιανουαρίου 2024 3:20 μμ
> *Προς:* Pantelis Moschopoulos 
> *Κοιν.:* Barry Smith ; petsc-users@mcs.anl.gov <
> petsc-users@mcs.anl.gov>
> *Θέμα:* Re: [petsc-users] Question about a parallel implementation of
> PCFIELDSPLIT
>
> On Tue, Jan 23, 2024 at 8:16 AM Pantelis Moschopoulos <
> pmoschopou...@outlook.com> wrote:
>
> Dear Matt,
>
> Thank you for your response. This is an idealized setup where I have only
> one row/column. Sometimes we might need two or even three constraints based
> on the application. Thus, I will pursue the user-defined IS.
>
>
> Anything < 50 I would use MatLRC. The bottleneck is the inversion of a
> dense matrix of size k x k, where k is the number of constraints. Using an
> IS is definitely fine, but dense rows can detract from iterative
> convergence.
>
>
> When I supply the IS using the command PCFieldSplitSetIS, I do not specify
> anything in the matrix set up right?
>
>
> You should just need to specify the rows for each field as an IS.
>
>   Thanks,
>
>  Matt
>
>
> Thanks,
> Pantelis
> --

Re: [petsc-users] Question about a parallel implementation of PCFIELDSPLIT

[Matthew Knepley]

On Tue, Jan 23, 2024 at 9:45 AM Pantelis Moschopoulos <
pmoschopou...@outlook.com> wrote:

> Dear Matt,
>
> I read about the MATLRC. However, its correct usage is not clear to me so
> I have the following questions:
>
>1. The U and V input matrices should be created as dense using
>MatCreateDense?
>
>
Yes. If you have one row, it looks like a vector, or a matrix with one
column.

If you have 1 row on the bottom, then

  U = [0, 0, ..., 0, 1]
  V = [the row]
  C = [1]

will give you that. However, you have an extra row and column?


>
>1. I use the command MatCreateLRC just to declare the matrix and then
>MatLRCSetMats to pass the values of the constituents?
>
> You can use

  MatCreate(comm, )
  MatSetSizes(M, ...)
  MatSetType(M, MATLRC)
  MatLRCSetMats(M, ...)

However, you are right that it is a little more complicated, because A is
not just the upper block here.


>1.
>
> Then, how do I proceed? How I apply the step of Sherman-Morrison-Woobury
> formula? I intend to use iterative solvers for A (main matrix) so I will
> not have its A^-1 at hand which I think is what the
> Sherman-Morrison-Woobury formula needs.
>

I think I was wrong. MatLRC is not the best fit. We should use MatNest
instead. Then you could have

  A   u
  v^t 0

as your matrix. We could still get an explicit Schur complement
automatically using nested FieldSplit. So

  1) Make the MATNEST matrix as shown above

  2) Use PCFIELDSPLIT  for it. This should be an easy IS, since there is
only one row in the second one.

  3) Select the full Schur complement

-pc_type fieldsplit
-pc_fieldsplit_type schur
-pc_fieldsplit_schur_fact_type full
-pc_fieldsplit_schur_precondition full

  4) Use a recursive FieldSplit (might be able to use
-fieldsplit_0_pc_fieldsplit_detect_saddle_point)

-fieldsplit_0_pc_type fieldsplit
-fieldsplit_0_pc_fieldsplit_0_fields 0,1,2
-fieldsplit_0_pc_fieldsplit_1_fields 3

I think this does what you want, and should be much easier than getting
MatLRC to do it.

  Thanks,

 Matt


> Thanks,
> Pantelis
> --
> *Από:* Matthew Knepley 
> *Στάλθηκε:* Τρίτη, 23 Ιανουαρίου 2024 3:20 μμ
> *Προς:* Pantelis Moschopoulos 
> *Κοιν.:* Barry Smith ; petsc-users@mcs.anl.gov <
> petsc-users@mcs.anl.gov>
> *Θέμα:* Re: [petsc-users] Question about a parallel implementation of
> PCFIELDSPLIT
>
> On Tue, Jan 23, 2024 at 8:16 AM Pantelis Moschopoulos <
> pmoschopou...@outlook.com> wrote:
>
> Dear Matt,
>
> Thank you for your response. This is an idealized setup where I have only
> one row/column. Sometimes we might need two or even three constraints based
> on the application. Thus, I will pursue the user-defined IS.
>
>
> Anything < 50 I would use MatLRC. The bottleneck is the inversion of a
> dense matrix of size k x k, where k is the number of constraints. Using an
> IS is definitely fine, but dense rows can detract from iterative
> convergence.
>
>
> When I supply the IS using the command PCFieldSplitSetIS, I do not specify
> anything in the matrix set up right?
>
>
> You should just need to specify the rows for each field as an IS.
>
>   Thanks,
>
>  Matt
>
>
> Thanks,
> Pantelis
> ----------
> *Από:* Matthew Knepley 
> *Στάλθηκε:* Τρίτη, 23 Ιανουαρίου 2024 2:51 μμ
> *Προς:* Pantelis Moschopoulos 
> *Κοιν.:* Barry Smith ; petsc-users@mcs.anl.gov <
> petsc-users@mcs.anl.gov>
> *Θέμα:* Re: [petsc-users] Question about a parallel implementation of
> PCFIELDSPLIT
>
> On Tue, Jan 23, 2024 at 4:23 AM Pantelis Moschopoulos <
> pmoschopou...@outlook.com> wrote:
>
> Dear Matt and Dear Barry,
>
> I have some follow up questions regarding FieldSplit.
> Let's assume that I solve again the 3D Stokes flow but now I have also a
> global constraint that controls the flow rate at the inlet. Now, the matrix
> has the same unknowns as before, i.e. ux0,uy0,uz0,p0//ux1,uy1,uz1,p1//...,
> but the last line (and the last column) corresponds to the contribution of
> the global constraint equation. I want to incorporate the last line (and
> last column)  into the local block of velocities (split 0) and the
> pressure. The problem is how I do that. I have two questions:
>
>1. Now, the block size should be 5 in the matrix and vector creation
>for this problem?
>
> No. Blocksize is only useful when the vector/matrix layout is completely
> regular, meaning _every_ block looks the same. Here you have a single row
> to be added in.
>
>
>1. I have to rely entirely on PCFieldSplitSetIS to create the two
>blocks? Can I augment simply the previously defined block 0 with the last
>line of the matrix?
>
> If you want to add in a s

Re: [petsc-users] Question about a parallel implementation of PCFIELDSPLIT

[Matthew Knepley]

On Tue, Jan 23, 2024 at 8:16 AM Pantelis Moschopoulos <
pmoschopou...@outlook.com> wrote:

> Dear Matt,
>
> Thank you for your response. This is an idealized setup where I have only
> one row/column. Sometimes we might need two or even three constraints based
> on the application. Thus, I will pursue the user-defined IS.
>

Anything < 50 I would use MatLRC. The bottleneck is the inversion of a
dense matrix of size k x k, where k is the number of constraints. Using an
IS is definitely fine, but dense rows can detract from iterative
convergence.


> When I supply the IS using the command PCFieldSplitSetIS, I do not specify
> anything in the matrix set up right?
>

You should just need to specify the rows for each field as an IS.

  Thanks,

 Matt


> Thanks,
> Pantelis
> --
> *Από:* Matthew Knepley 
> *Στάλθηκε:* Τρίτη, 23 Ιανουαρίου 2024 2:51 μμ
> *Προς:* Pantelis Moschopoulos 
> *Κοιν.:* Barry Smith ; petsc-users@mcs.anl.gov <
> petsc-users@mcs.anl.gov>
> *Θέμα:* Re: [petsc-users] Question about a parallel implementation of
> PCFIELDSPLIT
>
> On Tue, Jan 23, 2024 at 4:23 AM Pantelis Moschopoulos <
> pmoschopou...@outlook.com> wrote:
>
> Dear Matt and Dear Barry,
>
> I have some follow up questions regarding FieldSplit.
> Let's assume that I solve again the 3D Stokes flow but now I have also a
> global constraint that controls the flow rate at the inlet. Now, the matrix
> has the same unknowns as before, i.e. ux0,uy0,uz0,p0//ux1,uy1,uz1,p1//...,
> but the last line (and the last column) corresponds to the contribution of
> the global constraint equation. I want to incorporate the last line (and
> last column)  into the local block of velocities (split 0) and the
> pressure. The problem is how I do that. I have two questions:
>
>1. Now, the block size should be 5 in the matrix and vector creation
>for this problem?
>
> No. Blocksize is only useful when the vector/matrix layout is completely
> regular, meaning _every_ block looks the same. Here you have a single row
> to be added in.
>
>
>1. I have to rely entirely on PCFieldSplitSetIS to create the two
>blocks? Can I augment simply the previously defined block 0 with the last
>line of the matrix?
>
> If you want to add in a single row, then you have to specify the IS
> yourself since we cannot generate it from the regular pattern.
>
> However, if you know that you will only ever have a single constraint row
> (which I assume is fairly dense), then I would suggest instead using
> MatLRC, which Jose developed for SLEPc. This handles the last row/col as a
> low-rank correction. One step of Sherman-Morrison-Woobury solves this
> exactly. It requires a solve for A, for which you can use FieldSplit as
> normal.
>
>   Thanks,
>
>  Matt
>
>
> Up to this moment, I use the following commands to create the Field split:
> ufields(3) = [0, 1, 2]
> pfields(1) = [3]
>
> call PCSetType(pc, PCFIELDSPLIT, ierr)
> call PCFieldSplitSetBlockSize(pc, 4,ierr)
>  call PCFieldSplitSetFields(pc, "0", 3, ufields, ufields,ierr)
>  call PCFieldSplitSetFields(pc, "1", 1, pfields, pfields,ierr)
>
> Thanks,
> Pantelis
>
>
> --
> *Από:* Matthew Knepley 
> *Στάλθηκε:* Παρασκευή, 19 Ιανουαρίου 2024 11:31 μμ
> *Προς:* Barry Smith 
> *Κοιν.:* Pantelis Moschopoulos ;
> petsc-users@mcs.anl.gov 
> *Θέμα:* Re: [petsc-users] Question about a parallel implementation of
> PCFIELDSPLIT
>
> On Fri, Jan 19, 2024 at 4:25 PM Barry Smith  wrote:
>
>
>Generally fieldsplit is used on problems that have a natural "split" of
> the variables into two or more subsets. For example u0,v0,u1,v1,u2,v2,u3,v4
> This is often indicated in the vectors and matrices with the "blocksize"
> argument, 2 in this case. DM also often provides this information.
>
>When laying out a vector/matrix with a blocksize one must ensure that
> an equal number of of the subsets appears on each MPI process. So, for
> example, if the above vector is distributed over 3 MPI processes one could
> use   u0,v0,u1,v1   u2,v2  u3,v3  but one cannot use u0,v0,u1
>  v1,u2,v2   u3,v3.  Another way to think about it is that one must split up
> the vector as indexed by block among the processes. For most multicomponent
> problems this type of decomposition is very natural in the logic of the
> code.
>
>
> This blocking is only convenient, not necessary. You can specify your own
> field division using PCFieldSplitSetIS().
>
>   Thanks,
>
>  Matt
>
>
>   Barry
>
>
> On Jan 19, 2024, at 3:19 AM, Pantelis Moschopoulos <
> pmoschopou...@outlook.com> wrot

Re: [petsc-users] Question about a parallel implementation of PCFIELDSPLIT

[Matthew Knepley]

On Tue, Jan 23, 2024 at 4:23 AM Pantelis Moschopoulos <
pmoschopou...@outlook.com> wrote:

> Dear Matt and Dear Barry,
>
> I have some follow up questions regarding FieldSplit.
> Let's assume that I solve again the 3D Stokes flow but now I have also a
> global constraint that controls the flow rate at the inlet. Now, the matrix
> has the same unknowns as before, i.e. ux0,uy0,uz0,p0//ux1,uy1,uz1,p1//...,
> but the last line (and the last column) corresponds to the contribution of
> the global constraint equation. I want to incorporate the last line (and
> last column)  into the local block of velocities (split 0) and the
> pressure. The problem is how I do that. I have two questions:
>
>1. Now, the block size should be 5 in the matrix and vector creation
>for this problem?
>
> No. Blocksize is only useful when the vector/matrix layout is completely
regular, meaning _every_ block looks the same. Here you have a single row
to be added in.

>
>1. I have to rely entirely on PCFieldSplitSetIS to create the two
>blocks? Can I augment simply the previously defined block 0 with the last
>line of the matrix?
>
> If you want to add in a single row, then you have to specify the IS
yourself since we cannot generate it from the regular pattern.

However, if you know that you will only ever have a single constraint row
(which I assume is fairly dense), then I would suggest instead using
MatLRC, which Jose developed for SLEPc. This handles the last row/col as a
low-rank correction. One step of Sherman-Morrison-Woobury solves this
exactly. It requires a solve for A, for which you can use FieldSplit as
normal.

  Thanks,

 Matt


> Up to this moment, I use the following commands to create the Field split:
> ufields(3) = [0, 1, 2]
> pfields(1) = [3]
>
> call PCSetType(pc, PCFIELDSPLIT, ierr)
> call PCFieldSplitSetBlockSize(pc, 4,ierr)
>  call PCFieldSplitSetFields(pc, "0", 3, ufields, ufields,ierr)
>  call PCFieldSplitSetFields(pc, "1", 1, pfields, pfields,ierr)
>
> Thanks,
> Pantelis
>
>
> --
> *Από:* Matthew Knepley 
> *Στάλθηκε:* Παρασκευή, 19 Ιανουαρίου 2024 11:31 μμ
> *Προς:* Barry Smith 
> *Κοιν.:* Pantelis Moschopoulos ;
> petsc-users@mcs.anl.gov 
> *Θέμα:* Re: [petsc-users] Question about a parallel implementation of
> PCFIELDSPLIT
>
> On Fri, Jan 19, 2024 at 4:25 PM Barry Smith  wrote:
>
>
>Generally fieldsplit is used on problems that have a natural "split" of
> the variables into two or more subsets. For example u0,v0,u1,v1,u2,v2,u3,v4
> This is often indicated in the vectors and matrices with the "blocksize"
> argument, 2 in this case. DM also often provides this information.
>
>When laying out a vector/matrix with a blocksize one must ensure that
> an equal number of of the subsets appears on each MPI process. So, for
> example, if the above vector is distributed over 3 MPI processes one could
> use   u0,v0,u1,v1   u2,v2  u3,v3  but one cannot use u0,v0,u1
>  v1,u2,v2   u3,v3.  Another way to think about it is that one must split up
> the vector as indexed by block among the processes. For most multicomponent
> problems this type of decomposition is very natural in the logic of the
> code.
>
>
> This blocking is only convenient, not necessary. You can specify your own
> field division using PCFieldSplitSetIS().
>
>   Thanks,
>
>  Matt
>
>
>   Barry
>
>
> On Jan 19, 2024, at 3:19 AM, Pantelis Moschopoulos <
> pmoschopou...@outlook.com> wrote:
>
> Dear all,
>
> When I am using PCFIELDSPLIT and pc type "schur" in serial mode everything
> works fine. When I turn now to parallel, I observe that the number of ranks
> that I can use must divide the number of N without any remainder, where N
> is the number of unknowns. Otherwise, an error of the following form
> emerges: "Local columns of A10 3473 do not equal local rows of A00 3471".
>
> Can I do something to overcome this?
>
> Thanks,
> Pantelis
>
>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/~knepley/>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>

Re: [petsc-users] Question about a parallel implementation of PCFIELDSPLIT

[Matthew Knepley]

On Fri, Jan 19, 2024 at 4:25 PM Barry Smith  wrote:

>
>Generally fieldsplit is used on problems that have a natural "split" of
> the variables into two or more subsets. For example u0,v0,u1,v1,u2,v2,u3,v4
> This is often indicated in the vectors and matrices with the "blocksize"
> argument, 2 in this case. DM also often provides this information.
>
>When laying out a vector/matrix with a blocksize one must ensure that
> an equal number of of the subsets appears on each MPI process. So, for
> example, if the above vector is distributed over 3 MPI processes one could
> use   u0,v0,u1,v1   u2,v2  u3,v3  but one cannot use u0,v0,u1
>  v1,u2,v2   u3,v3.  Another way to think about it is that one must split up
> the vector as indexed by block among the processes. For most multicomponent
> problems this type of decomposition is very natural in the logic of the
> code.
>

This blocking is only convenient, not necessary. You can specify your own
field division using PCFieldSplitSetIS().

  Thanks,

 Matt


>   Barry
>
>
> On Jan 19, 2024, at 3:19 AM, Pantelis Moschopoulos <
> pmoschopou...@outlook.com> wrote:
>
> Dear all,
>
> When I am using PCFIELDSPLIT and pc type "schur" in serial mode everything
> works fine. When I turn now to parallel, I observe that the number of ranks
> that I can use must divide the number of N without any remainder, where N
> is the number of unknowns. Otherwise, an error of the following form
> emerges: "Local columns of A10 3473 do not equal local rows of A00 3471".
>
> Can I do something to overcome this?
>
> Thanks,
> Pantelis
>
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about a parallel implementation of PCFIELDSPLIT

[Barry Smith]

   Generally fieldsplit is used on problems that have a natural "split" of the 
variables into two or more subsets. For example u0,v0,u1,v1,u2,v2,u3,v4 This is 
often indicated in the vectors and matrices with the "blocksize" argument, 2 in 
this case. DM also often provides this information. 

   When laying out a vector/matrix with a blocksize one must ensure that an 
equal number of of the subsets appears on each MPI process. So, for example, if 
the above vector is distributed over 3 MPI processes one could use   
u0,v0,u1,v1   u2,v2  u3,v3  but one cannot use u0,v0,u1v1,u2,v2   
u3,v3.  Another way to think about it is that one must split up the vector as 
indexed by block among the processes. For most multicomponent problems this 
type of decomposition is very natural in the logic of the code.

  Barry
 

> On Jan 19, 2024, at 3:19 AM, Pantelis Moschopoulos 
>  wrote:
> 
> Dear all,
> 
> When I am using PCFIELDSPLIT and pc type "schur" in serial mode everything 
> works fine. When I turn now to parallel, I observe that the number of ranks 
> that I can use must divide the number of N without any remainder, where N is 
> the number of unknowns. Otherwise, an error of the following form emerges: 
> "Local columns of A10 3473 do not equal local rows of A00 3471".
> 
> Can I do something to overcome this?
> 
> Thanks,
> Pantelis

Re: [petsc-users] Question about petsc4py with cuda

[Matthew Knepley]

On Mon, Jan 15, 2024 at 11:57 AM MIA via petsc-users <
petsc-users@mcs.anl.gov> wrote:

> Hi,
>
> I am trying to solve a large linear equation, which needs a GPU solver as
> comparison. I install a CUDA-enabled PETSc and petsc4py from sources using
> the release tarball. According to the test results after installation, the
> PETSc can successfully work with cuda.
>

Here is a How-To for GPUs: https://petsc.org/main/faq/#doc-faq-gpuhowto

  Thanks,

 Matt


> All my programs are written in python, so I turn to petsc4py. But I do not
> find any commands that define variables on coda device or define where the
> direct solver is executed. I check `nvidia-smi` and find my cuda does not
> work at all when executing my python script:
>
> from petsc4py import PETSc
> import numpy as np
>
> n = 1000
>
> nnz = 3 * np.ones(n, dtype=np.int32)
> nnz[0] = nnz[-1] = 2
>
> A = PETSc.Mat()
> A.createAIJ([n, n], nnz=nnz)
>
> # First set the first row
> A.setValue(0, 0, 2)
> A.setValue(0, 1, -1)
> # Now we fill the last row
> A.setValue(999, 998, -1)
> A.setValue(999, 999, 2)
>
> # And now everything else
> for index in range(1, n - 1):
> A.setValue(index, index - 1, -1)
> A.setValue(index, index, 2)
> A.setValue(index, index + 1, -1)
>
> A.assemble()
>
> indexptr, indices, data = A.getValuesCSR()
> b = A.createVecLeft()
> b.array[:] = 1
> for i in range(10):
> ksp = PETSc.KSP().create()
> ksp.setOperators(A)
> ksp.setType('preonly')
> ksp.setConvergenceHistory()
> ksp.getPC().setType('lu')
> x = A.createVecRight()
> ksp.solve(2*b, x)
> residual = A * x - 2*b
> if i % 10 == 0:
> print(f"The relative residual is: {residual.norm() / b.norm()}.”)
>
> What should I do to utlize GPU to execute the KSP task? Are there some
> settings to be modified?
>
> Looking forward to your early reply. Thanks a lot.
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about MATMPIAIJ and column indexes ordering

[Barry Smith]

   Added clarification to the man pages in 
https://gitlab.com/petsc/petsc/-/merge_requests/7170


> On Jan 8, 2024, at 4:31 AM, Deuse, Mathieu via petsc-users 
>  wrote:
> 
> Hello,
>  
> I have a piece of code which generates a matrix in CSR format, but the 
> without sorting the column indexes in increasing order within each row. This 
> seems not to be 100% compatible with the MATMPIAIJ format: the documentation 
> of MatCreateMPIAIJWithArrays indeed mentions 'row-major ordering'.
>  
> For example, consider the 2x2 matrix (1 2; 3 4), which in my code could be 
> stored as i=[0, 2, 4], j=[1, 0, 0, 1], v=[2, 1, 3, 4]. I can generate the 
> matrix as follows (on 1 proc): MatCreateMPIAIJWithArrays(PETSC_COMM_SELF, 2, 
> 2, 2, 2, i, j, v, ). This appears to work fine, and I can then use the 
> matrix in a KSP for example. However, if I try to update the entry values 
> (same order and values v=[2, 1, 3, 4]) with MatUpdateMPIAIJWithArray(matrix, 
> v), it seems that PETSc does not memorize the order of the column indexes and 
> the matrix that I get now is (2 1; 3 4). I get the same result with 
> MatUpdateMPIAIJWithArrays(matrix, 2, 2, 2, 2, i, j, v). On the other hand, if 
> the column indexes are sorted within each row (i=[0, 2, 4], j=[0, 1, 0, 1], 
> v=[1, 2, 3, 4]), then it works fine. I have attached a minimal working 
> example (C++).
>  
> Can I safely rely on MatCreateMPIAIJWithArrays working fine with unsorted 
> column indexes as long as I do not use MatUpdateMPIAIJWithArray(s)? Or should 
> I do the sorting myself before calling MatCreateMPIAIJWithArrays? (or 
> alternatively use another matrix format).
>  
> Thanks in advance for the help.
>  
> Kind regards,
>  
> Mathieu Deuse
> 

Re: [petsc-users] Question about MATMPIAIJ and column indexes ordering

[Barry Smith]

> On Jan 8, 2024, at 4:31 AM, Deuse, Mathieu via petsc-users 
>  wrote:
> 
> Hello,
>  
> I have a piece of code which generates a matrix in CSR format, but the 
> without sorting the column indexes in increasing order within each row. This 
> seems not to be 100% compatible with the MATMPIAIJ format: the documentation 
> of MatCreateMPIAIJWithArrays indeed mentions 'row-major ordering'.
>  
> For example, consider the 2x2 matrix (1 2; 3 4), which in my code could be 
> stored as i=[0, 2, 4], j=[1, 0, 0, 1], v=[2, 1, 3, 4]. I can generate the 
> matrix as follows (on 1 proc): MatCreateMPIAIJWithArrays(PETSC_COMM_SELF, 2, 
> 2, 2, 2, i, j, v, ). This appears to work fine, and I can then use the 
> matrix in a KSP for example. However, if I try to update the entry values 
> (same order and values v=[2, 1, 3, 4]) with MatUpdateMPIAIJWithArray(matrix, 
> v), it seems that PETSc does not memorize the order of the column indexes and 
> the matrix that I get now is (2 1; 3 4). I get the same result with 
> MatUpdateMPIAIJWithArrays(matrix, 2, 2, 2, 2, i, j, v). On the other hand, if 
> the column indexes are sorted within each row (i=[0, 2, 4], j=[0, 1, 0, 1], 
> v=[1, 2, 3, 4]), then it works fine. I have attached a minimal working 
> example (C++).
>  
> Can I safely rely on MatCreateMPIAIJWithArrays working fine with unsorted 
> column indexes as long as I do not use MatUpdateMPIAIJWithArray(s)?

   Yes, this is correct. The column indices do not need to be sorted if you 
never call MatUpdateMPIAIJWithArray().



> Or should I do the sorting myself before calling MatCreateMPIAIJWithArrays? 
> (or alternatively use another matrix format).
>  
> Thanks in advance for the help.
>  
> Kind regards,
>  
> Mathieu Deuse
> 

Re: [petsc-users] Question on output vector in vtk file

[Matthew Knepley]

On Mon, Dec 11, 2023 at 4:32 AM Gong Yujie 
wrote:

> Dear PETSc developers,
>
> I have a DMPlex DM with 1 field 1dof. I'd like to output a vector with
> block size equals to 3. It seems that there is no response using command
> line option or using some code about PetscViewer.
>

I am not sure how we can do this. If you only have 1 dof per cell (I
assume), how can we have a blocksize of 3?

  Thanks,

 Matt


> The DM is generated with (note that I'm not using PetscFE for
> discretization, just for allocate dof.)
>
> *PetscCall(DMPlexCreateExodusFromFile(PETSC_COMM_WORLD,"tube.exo",interpolate,));*
>
> *PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF,dim,1,PETSC_TRUE,1,PETSC_DETERMINE,));*
> *PetscCall(PetscObjectSetName((PetscObject)fe,"potential_field"));*
> *PetscCall(DMSetField(dm,0,NULL,(PetscObject)fe));*
> *PetscCall(DMPlexDistribute(dm,0,,));*
>
> The Vector is created using
> *PetscCall(DMCreateGlobalVector(dm,_1));*
> *PetscCall(VecSetLocalToGlobalMapping(phi_1,Itog));*
> *PetscCall(VecGetLocalSize(phi_1,_local_size_test));*
> *PetscCall(VecCreateMPI(PETSC_COMM_WORLD, vec_local_size_test*3,
> PETSC_DETERMINE, _grad_psi));*
> *PetscCall(VecSetBlockSize(u_grad_psi, 3));*
> *PetscCall(VecSetLocalToGlobalMapping(u_grad_psi,Itog));*
>
> The output command line option is just --vec_view vtk:test.vtk. The PETSc
> version I'm using is 3.19.5.
>
> Could you please give me some advice?
>
> Best Regards,
> Yujie
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about using PETSc to generate matrix preconditioners

[Jed Brown]

I think random matrices will produce misleading results. The chance of randomly 
generating a matrix that resembles an application is effectively zero. I think 
you'd be better off with some model problems varying parameters that control 
the physical regime (e.g., shifts to a Laplacian, advection with/without 
upwinding at varying cell Peclet numbers, varying nondimensional parameters, 
time step, and boundary layers for fluids, varying Poisson ratio and material 
contrast in elasticity).

"Sun, Yixuan via petsc-users"  writes:

> Hello PETSc team,
>
> My name is Yixuan, and I am a postdoc at MCS. We are working on a project 
> where we want to predict the preconditioners from the corresponding matrices. 
> For preliminary exploration purposes, we generated random matrices and their 
> preconditioners using PETSc. The current deep learning models look promising 
> with the generated data. However, we are not sure about the correctness of 
> the data generation code and wanted to ask if you could help check the 
> correctness of our script (< 100 lines) and let us know if there were any 
> issues. Here is the 
> link
>  to the script.
>
> Please let me know if you have any questions or concerns. Thank you in 
> advance!
>
>
> Warm regards,
> Yixuan
> 
> Yixuan Sun
> Postdoctoral Researcher
> Mathematics and Computer Science
> Argonne National Laboratory

Re: [petsc-users] Question about residue norm in PETSc

[Barry Smith]

  Also look at 
https://petsc.org/release/manualpages/KSP/KSPSetPCSide/#kspsetpcside and 
https://petsc.org/release/manualpages/KSP/KSPSetNormType/#kspsetnormtype  in 
PETSc different Krylov solvers have different default values for this. 



> On Jul 2, 2023, at 1:47 AM, 王赫萌  wrote:
> 
> Dear PETSc Team,
> 
> Sorry to bother! My name is Hemeng Wang, and I am currently learning the use 
> of PETSc software package. I am confused while calculating the norm of 
> residue.
> 
> I calculated residue norm by myself with:
> ```
>   PetscCall(VecNorm(b, NORM_2, _b)); // (main.c, line 74)
> 
>   PetscCall(VecDuplicate(b, )); // (main.c, line 105)
>   PetscCall(MatMult(A, x, u));
>   PetscCall(VecAXPY(b, -1.0, u));
>   PetscCall(VecNorm(b, NORM_2, _delta));
> ```
> and check the (norm_delta) / (norm_b). It seems not the `atol` which set by 
> `KSPSetTolerances()`.
> (I set atol as 1e-12, but got 5e-7 finally)
> (options:  -ksp_type cg -pc_type gamg -ksp_converged_reason -ksp_norm_type 
> unpreconditioned -ksp_monitor_true_residual)
> 
> I also check the soure code of `KSPSolve_CG` in 
> `petsc/src/ksp/ksp/impls/cg/cg.c`. And could not figure out where is the 
> difference.
> 
> I will really really appreciated if someone can explain the calculation of 
> `resid norm` in petsc. And where is my mistake.
> 
> To provide you with more context, here are the source code about my 
> implementation. And the output of my test.
> 
> main.c
> Main code of my program
> 
> mmio.h  mmloader.h
> Headers for matrix read
> 
> Makefile
> For compiling, same as sample provided
> 
> task.sh
> A script for running program in `slurm`
> 
> slurm-4803840.out
> Output of my test
> 
> Thank you very much for your time and attention. I greatly appreciate your 
> support and look forward to hearing from you soon.
> Best regards,
> Hemeng Wang
> 
> 

Re: [petsc-users] Question about residue norm in PETSc

[Matthew Knepley]

On Sun, Jul 2, 2023 at 8:45 AM 王赫萌  wrote:

> Thank so much for your patience! I'm really grateful for that!
>
> Could you explain the calculation of "1.31278e+06"
> It appears that my calculation of this value is "18.0468"
>

100 KSP unpreconditioned resid norm 1.312782529439e+06 true resid norm
1.312782529439e+06 ||r(i)||/||b|| 1.804684612214e+01
norm_delta 1312782.529439, norm_b 72743.044439
Norm of error 18.0468 iterations 100
Residue of error 1.31278e+06 iterations 100

Here norm_delta matches true resid norm, because they are both A x - b.
What you are calling
"Norm of error" is actually the "Relative Residual" ||r_100|| / ||b||
= 18.0468, which is correct since
that is 1312782.529439 / 72743.044439. What you call "Residue of error" is
the "Residual" just
as reported in the KSP monitor output.

  Thanks,

  Matt


> Please forgive my ignorance. I just begin my learning on iterative method.
>
> Thank you so much for the reply!
>
> Welcome to China if you have time! ^_^
>
>
> <https://maas.mail.163.com/dashi-web-extend/html/proSignature.html?ftlId=1=%E7%8E%8B%E8%B5%AB%E8%90%8C=wanghemeng%40163.com=https%3A%2F%2Fmail-online.nosdn.127.net%2Fqiyelogo%2FdefaultAvatar.png=%5B%22wanghemeng%40163.com%22%5D>
> ---- Replied Message ----
> From Matthew Knepley 
> Date 7/2/2023 20:32
> To 王赫萌 
> Cc PETSc 
> Subject Re: [petsc-users] Question about residue norm in PETSc
> On Sun, Jul 2, 2023 at 8:19 AM 王赫萌  wrote:
>
>> Here is the mat and rhs used in code! (May need to change the data path)
>>
>> mat:
>>
>> https://studentcupeducn-my.sharepoint.com/:u:/g/personal/wanghemeng_student_cup_edu_cn/Ed76oGtC1ttDriZsObbPR74BCnDPUP8aicVXQEL4sO1AyQ?e=zeszik
>> rhs:
>>
>> https://studentcupeducn-my.sharepoint.com/:u:/g/personal/wanghemeng_student_cup_edu_cn/EdHRqWbzVmtIkAppOLL1UMIBM7tK7ws0gEASESGHuGC3yw?e=SMQSmY
>>
>> I tried using
>> PetscCall(VecAXPY(u, -1.0, b));
>> but is just as same as
>> PetscCall(VecAXPY(b, -1.0, u));
>>
>> Thank you so much for that!!!
>>
>
>  90 KSP unpreconditioned resid norm 9.749157899195e+05 true resid norm
> 9.749157899195e+05 ||r(i)||/||b|| 1.340218569960e+01
>  91 KSP unpreconditioned resid norm 1.073123446417e+06 true resid norm
> 1.073123446417e+06 ||r(i)||/||b|| 1.475224820050e+01
>  92 KSP unpreconditioned resid norm 1.170251286554e+06 true resid norm
> 1.170251286554e+06 ||r(i)||/||b|| 1.608746644557e+01
>  93 KSP unpreconditioned resid norm 1.264719067990e+06 true resid norm
> 1.264719067990e+06 ||r(i)||/||b|| 1.738611681365e+01
>  94 KSP unpreconditioned resid norm 1.329446257320e+06 true resid norm
> 1.329446257320e+06 ||r(i)||/||b|| 1.827592270272e+01
>  95 KSP unpreconditioned resid norm 1.365944956372e+06 true resid norm
> 1.365944956372e+06 ||r(i)||/||b|| 1.877767100504e+01
>  96 KSP unpreconditioned resid norm 1.369513563400e+06 true resid norm
> 1.369513563400e+06 ||r(i)||/||b|| 1.882672871297e+01
>  97 KSP unpreconditioned resid norm 1.364905651654e+06 true resid norm
> 1.364905651654e+06 ||r(i)||/||b|| 1.876338366353e+01
>  98 KSP unpreconditioned resid norm 1.352584030803e+06 true resid norm
> 1.352584030803e+06 ||r(i)||/||b|| 1.859399810996e+01
>  99 KSP unpreconditioned resid norm 1.330589478009e+06 true resid norm
> 1.330589478009e+06 ||r(i)||/||b|| 1.829163857903e+01
> 100 KSP unpreconditioned resid norm 1.312782529439e+06 true resid norm
> 1.312782529439e+06 ||r(i)||/||b|| 1.804684612214e+01
> Linear solve did not converge due to DIVERGED_ITS iterations 100
> KSPSolve Time: 7579.364000 ms
> norm_delta 1312782.529439, norm_b 72743.044439
> Norm of error 18.0468 iterations 100
> Residue of error 1.31278e+06 iterations 100
>
> I ran with
>
>   -pc_type jacobi -ksp_max_it 100
>
> because GAMG takes a long time to setup on my laptop. Those numbers match
> exactly.
>
>THanks,
>
>  Matt
>
> I'm such a beginner T_T
>>  Replied Message 
>> From Matthew Knepley 
>> Date 7/2/2023 20:10
>> To 王赫萌 
>> Cc PETSc 
>> Subject Re: [petsc-users] Question about residue norm in PETSc
>> On Sun, Jul 2, 2023 at 8:05 AM Matthew Knepley  wrote:
>>
>>> On Sun, Jul 2, 2023 at 7:53 AM 王赫萌  wrote:
>>>
>>>> Thanks for your reply!
>>>> So sorry that I made a mistake in the description.
>>>> I set the tolerances by:
>>>> PetscCall(KSPSetTolerances(ksp, 1e-12, DBL_MIN, PETSC_DEFAULT,
>>>> PETSC_DEFAULT));
>>>> and got (by passing `-ksp_norm_type unpreconditioned
>>>> -ksp_monitor_true_residual`)
>>>> 74 KSP unpreconditioned resid norm 7.256655641876e-08 true resid norm
>>

Re: [petsc-users] Question about residue norm in PETSc

[Matthew Knepley]

On Sun, Jul 2, 2023 at 8:19 AM 王赫萌  wrote:

> Here is the mat and rhs used in code! (May need to change the data path)
>
> mat:
>
> https://studentcupeducn-my.sharepoint.com/:u:/g/personal/wanghemeng_student_cup_edu_cn/Ed76oGtC1ttDriZsObbPR74BCnDPUP8aicVXQEL4sO1AyQ?e=zeszik
> rhs:
>
> https://studentcupeducn-my.sharepoint.com/:u:/g/personal/wanghemeng_student_cup_edu_cn/EdHRqWbzVmtIkAppOLL1UMIBM7tK7ws0gEASESGHuGC3yw?e=SMQSmY
>
> I tried using
> PetscCall(VecAXPY(u, -1.0, b));
> but is just as same as
> PetscCall(VecAXPY(b, -1.0, u));
>
> Thank you so much for that!!!
>

 90 KSP unpreconditioned resid norm 9.749157899195e+05 true resid norm
9.749157899195e+05 ||r(i)||/||b|| 1.340218569960e+01
 91 KSP unpreconditioned resid norm 1.073123446417e+06 true resid norm
1.073123446417e+06 ||r(i)||/||b|| 1.475224820050e+01
 92 KSP unpreconditioned resid norm 1.170251286554e+06 true resid norm
1.170251286554e+06 ||r(i)||/||b|| 1.608746644557e+01
 93 KSP unpreconditioned resid norm 1.264719067990e+06 true resid norm
1.264719067990e+06 ||r(i)||/||b|| 1.738611681365e+01
 94 KSP unpreconditioned resid norm 1.329446257320e+06 true resid norm
1.329446257320e+06 ||r(i)||/||b|| 1.827592270272e+01
 95 KSP unpreconditioned resid norm 1.365944956372e+06 true resid norm
1.365944956372e+06 ||r(i)||/||b|| 1.877767100504e+01
 96 KSP unpreconditioned resid norm 1.369513563400e+06 true resid norm
1.369513563400e+06 ||r(i)||/||b|| 1.882672871297e+01
 97 KSP unpreconditioned resid norm 1.364905651654e+06 true resid norm
1.364905651654e+06 ||r(i)||/||b|| 1.876338366353e+01
 98 KSP unpreconditioned resid norm 1.352584030803e+06 true resid norm
1.352584030803e+06 ||r(i)||/||b|| 1.859399810996e+01
 99 KSP unpreconditioned resid norm 1.330589478009e+06 true resid norm
1.330589478009e+06 ||r(i)||/||b|| 1.829163857903e+01
100 KSP unpreconditioned resid norm 1.312782529439e+06 true resid norm
1.312782529439e+06 ||r(i)||/||b|| 1.804684612214e+01
Linear solve did not converge due to DIVERGED_ITS iterations 100
KSPSolve Time: 7579.364000 ms
norm_delta 1312782.529439, norm_b 72743.044439
Norm of error 18.0468 iterations 100
Residue of error 1.31278e+06 iterations 100

I ran with

  -pc_type jacobi -ksp_max_it 100

because GAMG takes a long time to setup on my laptop. Those numbers match
exactly.

   THanks,

 Matt

I'm such a beginner T_T
>  Replied Message 
> From Matthew Knepley 
> Date 7/2/2023 20:10
> To 王赫萌 
> Cc PETSc 
> Subject Re: [petsc-users] Question about residue norm in PETSc
> On Sun, Jul 2, 2023 at 8:05 AM Matthew Knepley  wrote:
>
>> On Sun, Jul 2, 2023 at 7:53 AM 王赫萌  wrote:
>>
>>> Thanks for your reply!
>>> So sorry that I made a mistake in the description.
>>> I set the tolerances by:
>>> PetscCall(KSPSetTolerances(ksp, 1e-12, DBL_MIN, PETSC_DEFAULT,
>>> PETSC_DEFAULT));
>>> and got (by passing `-ksp_norm_type unpreconditioned
>>> -ksp_monitor_true_residual`)
>>> 74 KSP unpreconditioned resid norm 7.256655641876e-08 true resid norm
>>> 7.256655641876e-08 ||r(i)||/||b|| 9.975738158726e-13
>>> I'm wondering why the ` ||r(i)||/||b||` is different with mine which
>>> calculated by:
>>> ```
>>>   PetscCall(VecNorm(b, NORM_2, _b)); // (main.c, line 74) (which is
>>> 72743.044439)
>>>   PetscCall(VecDuplicate(b, )); // (main.c, line 105)
>>>   PetscCall(MatMult(A, x, u));
>>>   PetscCall(VecAXPY(b, -1.0, u));
>>>   PetscCall(VecNorm(b, NORM_2, _delta)); // (which is 0.039608)
>>> ```
>>> and (norm_delta) / (norm_b) = 5.44496e-07 which is higher and different
>>> with the rtol I set (1e-12).
>>> Sorry again for the waste of your time. I would really appreciated if
>>> you could help me again!
>>>
>>
>> 1) 7.256655641876e-08 / 72743.044439 = 9.975738158726e-13 so
>> ||r_i||/||b|| is correct in the output
>>
>> 2) You are asking why you calculate a different residual? I will have to
>> run your code.
>>
>
> I built your code, but you did not send the matrix and rhs.
>
> I suggest using
>
> PetscCall(VecAXPY(u, -1.0, b));
>
> instead so that you do not change b, and keep the residual in u.
>
>   Thanks,
>
>  Matt
>
>   Thanks,
>
>  Matt
>
>
>>   Thanks,
>>
>>  Matt
>>
>>
>>> Best regards!
>>> Hemeng Wang
>>>
>>>
>>>  Replied Message 
>>> From Matthew Knepley 
>>> Date 7/2/2023 18:51
>>> To 王赫萌 
>>> Cc petsc-users@mcs.anl.gov
>>> 
>>> Subject Re: [petsc-users] Question about residue norm in PETSc
>>> On Sun, Jul 2, 2023 at 2:24 AM 王赫萌  wrote:
>>>
>>>> Dear PETS

Re: [petsc-users] Question about residue norm in PETSc

[Matthew Knepley]

On Sun, Jul 2, 2023 at 8:05 AM Matthew Knepley  wrote:

> On Sun, Jul 2, 2023 at 7:53 AM 王赫萌  wrote:
>
>> Thanks for your reply!
>> So sorry that I made a mistake in the description.
>> I set the tolerances by:
>> PetscCall(KSPSetTolerances(ksp, 1e-12, DBL_MIN, PETSC_DEFAULT,
>> PETSC_DEFAULT));
>> and got (by passing `-ksp_norm_type unpreconditioned
>> -ksp_monitor_true_residual`)
>> 74 KSP unpreconditioned resid norm 7.256655641876e-08 true resid norm
>> 7.256655641876e-08 ||r(i)||/||b|| 9.975738158726e-13
>> I'm wondering why the ` ||r(i)||/||b||` is different with mine which
>> calculated by:
>> ```
>>   PetscCall(VecNorm(b, NORM_2, _b)); // (main.c, line 74) (which is
>> 72743.044439)
>>   PetscCall(VecDuplicate(b, )); // (main.c, line 105)
>>   PetscCall(MatMult(A, x, u));
>>   PetscCall(VecAXPY(b, -1.0, u));
>>   PetscCall(VecNorm(b, NORM_2, _delta)); // (which is 0.039608)
>> ```
>> and (norm_delta) / (norm_b) = 5.44496e-07 which is higher and different
>> with the rtol I set (1e-12).
>> Sorry again for the waste of your time. I would really appreciated if you
>> could help me again!
>>
>
> 1) 7.256655641876e-08 / 72743.044439 = 9.975738158726e-13 so ||r_i||/||b||
> is correct in the output
>
> 2) You are asking why you calculate a different residual? I will have to
> run your code.
>

I built your code, but you did not send the matrix and rhs.

I suggest using

PetscCall(VecAXPY(u, -1.0, b));

instead so that you do not change b, and keep the residual in u.

  Thanks,

 Matt

  Thanks,

 Matt


>   Thanks,
>
>  Matt
>
>
>> Best regards!
>> Hemeng Wang
>>
>>
>>  Replied Message 
>> From Matthew Knepley 
>> Date 7/2/2023 18:51
>> To 王赫萌 
>> Cc petsc-users@mcs.anl.gov
>> 
>> Subject Re: [petsc-users] Question about residue norm in PETSc
>> On Sun, Jul 2, 2023 at 2:24 AM 王赫萌  wrote:
>>
>>> Dear PETSc Team,
>>>
>>> Sorry to bother! My name is Hemeng Wang, and I am currently learning the
>>> use of PETSc software package. I am confused while calculating the norm of
>>> residue.
>>>
>>> I calculated residue norm by myself with:
>>> ```
>>>   PetscCall(VecNorm(b, NORM_2, _b)); // (main.c, line 74)
>>>
>>>   PetscCall(VecDuplicate(b, )); // (main.c, line 105)
>>>   PetscCall(MatMult(A, x, u));
>>>   PetscCall(VecAXPY(b, -1.0, u));
>>>   PetscCall(VecNorm(b, NORM_2, _delta));
>>> ```
>>> and check the (norm_delta) / (norm_b). It seems not the `atol` which set
>>> by `KSPSetTolerances()`.
>>> (I set atol as 1e-12, but got 5e-7 finally)
>>> (options:  -ksp_type cg -pc_type gamg -ksp_converged_reason
>>> -ksp_norm_type unpreconditioned -ksp_monitor_true_residual)
>>>
>>
>> If you are using the default convergence test, there is an absolute
>> tolerance (atol) _and_ a relative tolerance (rtol). It seems likely you hit
>> the relative tolerance. You can check this using
>>
>>   -ksp_converged_reason
>>
>> You could make rtol really small if you want to just see the atol
>>
>>   -ksp_rtol 1e-20
>>
>>   Thanks,
>>
>>  Matt
>>
>>
>>> I also check the soure code of `KSPSolve_CG` in
>>> `petsc/src/ksp/ksp/impls/cg/cg.c`. And could not figure out where is the
>>> difference.
>>>
>>> I will really really appreciated if someone can explain the calculation
>>> of `resid norm` in petsc. And where is my mistake.
>>>
>>> To provide you with more context, here are the source code about my
>>> implementation. And the output of my test.
>>>
>>> main.c
>>> Main code of my program
>>>
>>> mmio.h  mmloader.h
>>> Headers for matrix read
>>>
>>> Makefile
>>> For compiling, same as sample provided
>>>
>>> task.sh
>>> A script for running program in `slurm`
>>>
>>> slurm-4803840.out
>>> Output of my test
>>>
>>> Thank you very much for your time and attention. I greatly appreciate
>>> your support and look forward to hearing from you soon.
>>> Best regards,
>>> Hemeng Wang
>>>
>>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> <http://www.cse.buffalo.edu/~knepley/>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/~knepley/>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>

Re: [petsc-users] Question about residue norm in PETSc

[Matthew Knepley]

On Sun, Jul 2, 2023 at 7:53 AM 王赫萌  wrote:

> Thanks for your reply!
> So sorry that I made a mistake in the description.
> I set the tolerances by:
> PetscCall(KSPSetTolerances(ksp, 1e-12, DBL_MIN, PETSC_DEFAULT,
> PETSC_DEFAULT));
> and got (by passing `-ksp_norm_type unpreconditioned
> -ksp_monitor_true_residual`)
> 74 KSP unpreconditioned resid norm 7.256655641876e-08 true resid norm
> 7.256655641876e-08 ||r(i)||/||b|| 9.975738158726e-13
> I'm wondering why the ` ||r(i)||/||b||` is different with mine which
> calculated by:
> ```
>   PetscCall(VecNorm(b, NORM_2, _b)); // (main.c, line 74) (which is
> 72743.044439)
>   PetscCall(VecDuplicate(b, )); // (main.c, line 105)
>   PetscCall(MatMult(A, x, u));
>   PetscCall(VecAXPY(b, -1.0, u));
>   PetscCall(VecNorm(b, NORM_2, _delta)); // (which is 0.039608)
> ```
> and (norm_delta) / (norm_b) = 5.44496e-07 which is higher and different
> with the rtol I set (1e-12).
> Sorry again for the waste of your time. I would really appreciated if you
> could help me again!
>

1) 7.256655641876e-08 / 72743.044439 = 9.975738158726e-13 so ||r_i||/||b||
is correct in the output

2) You are asking why you calculate a different residual? I will have to
run your code.

  Thanks,

 Matt


> Best regards!
> Hemeng Wang
>
>
>  Replied Message 
> From Matthew Knepley 
> Date 7/2/2023 18:51
> To 王赫萌 
> Cc petsc-users@mcs.anl.gov
> 
> Subject Re: [petsc-users] Question about residue norm in PETSc
> On Sun, Jul 2, 2023 at 2:24 AM 王赫萌  wrote:
>
>> Dear PETSc Team,
>>
>> Sorry to bother! My name is Hemeng Wang, and I am currently learning the
>> use of PETSc software package. I am confused while calculating the norm of
>> residue.
>>
>> I calculated residue norm by myself with:
>> ```
>>   PetscCall(VecNorm(b, NORM_2, _b)); // (main.c, line 74)
>>
>>   PetscCall(VecDuplicate(b, )); // (main.c, line 105)
>>   PetscCall(MatMult(A, x, u));
>>   PetscCall(VecAXPY(b, -1.0, u));
>>   PetscCall(VecNorm(b, NORM_2, _delta));
>> ```
>> and check the (norm_delta) / (norm_b). It seems not the `atol` which set
>> by `KSPSetTolerances()`.
>> (I set atol as 1e-12, but got 5e-7 finally)
>> (options:  -ksp_type cg -pc_type gamg -ksp_converged_reason
>> -ksp_norm_type unpreconditioned -ksp_monitor_true_residual)
>>
>
> If you are using the default convergence test, there is an absolute
> tolerance (atol) _and_ a relative tolerance (rtol). It seems likely you hit
> the relative tolerance. You can check this using
>
>   -ksp_converged_reason
>
> You could make rtol really small if you want to just see the atol
>
>   -ksp_rtol 1e-20
>
>   Thanks,
>
>  Matt
>
>
>> I also check the soure code of `KSPSolve_CG` in
>> `petsc/src/ksp/ksp/impls/cg/cg.c`. And could not figure out where is the
>> difference.
>>
>> I will really really appreciated if someone can explain the calculation
>> of `resid norm` in petsc. And where is my mistake.
>>
>> To provide you with more context, here are the source code about my
>> implementation. And the output of my test.
>>
>> main.c
>> Main code of my program
>>
>> mmio.h  mmloader.h
>> Headers for matrix read
>>
>> Makefile
>> For compiling, same as sample provided
>>
>> task.sh
>> A script for running program in `slurm`
>>
>> slurm-4803840.out
>> Output of my test
>>
>> Thank you very much for your time and attention. I greatly appreciate
>> your support and look forward to hearing from you soon.
>> Best regards,
>> Hemeng Wang
>>
>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/~knepley/>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>

Re: [petsc-users] Question about residue norm in PETSc

[Matthew Knepley]

On Sun, Jul 2, 2023 at 2:24 AM 王赫萌  wrote:

> Dear PETSc Team,
>
> Sorry to bother! My name is Hemeng Wang, and I am currently learning the
> use of PETSc software package. I am confused while calculating the norm of
> residue.
>
> I calculated residue norm by myself with:
> ```
>   PetscCall(VecNorm(b, NORM_2, _b)); // (main.c, line 74)
>
>   PetscCall(VecDuplicate(b, )); // (main.c, line 105)
>   PetscCall(MatMult(A, x, u));
>   PetscCall(VecAXPY(b, -1.0, u));
>   PetscCall(VecNorm(b, NORM_2, _delta));
> ```
> and check the (norm_delta) / (norm_b). It seems not the `atol` which set
> by `KSPSetTolerances()`.
> (I set atol as 1e-12, but got 5e-7 finally)
> (options:  -ksp_type cg -pc_type gamg -ksp_converged_reason -ksp_norm_type
> unpreconditioned -ksp_monitor_true_residual)
>

If you are using the default convergence test, there is an absolute
tolerance (atol) _and_ a relative tolerance (rtol). It seems likely you hit
the relative tolerance. You can check this using

  -ksp_converged_reason

You could make rtol really small if you want to just see the atol

  -ksp_rtol 1e-20

  Thanks,

 Matt


> I also check the soure code of `KSPSolve_CG` in
> `petsc/src/ksp/ksp/impls/cg/cg.c`. And could not figure out where is the
> difference.
>
> I will really really appreciated if someone can explain the calculation of
> `resid norm` in petsc. And where is my mistake.
>
> To provide you with more context, here are the source code about my
> implementation. And the output of my test.
>
> main.c
> Main code of my program
>
> mmio.h  mmloader.h
> Headers for matrix read
>
> Makefile
> For compiling, same as sample provided
>
> task.sh
> A script for running program in `slurm`
>
> slurm-4803840.out
> Output of my test
>
> Thank you very much for your time and attention. I greatly appreciate your
> support and look forward to hearing from you soon.
> Best regards,
> Hemeng Wang
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about using MatCreateAIJ

[Barry Smith]

  I am concerned this is not good advice being provided.

  Let's back up and look more closely at your use case.

  * What is the ratio of new nonzero locations added compared to the initial 
number of nonzeros for your code, for each of your iterations?

  * Is it possible for many iterations, no or very few new nonzeros are being 
added?

  * Are many previous nonzero values becoming zero (or unneeded) later? Again 
as a ratio compared to the initial number of nonzeros?

  * Can you quantify the difference in time between initially filling the 
matrix and then refilling it using the reset preallocation and not using the 
reset preallocation?

The effect you report that resetting the preallocation results in slower code 
is possible if relatively few additional nonzero locations are being created.


  After a matrix is assembled with a given nonzero structure (regardless of how 
it was filled it, using preallocation or not), setting nonzero values into new 
locations will be slow due to needing to do possibly many mallocs() (as much as 
one for each new nonzero introduced). Resetting the initially provided 
preallocation removes the need for all the new mallocs(), but at the expense of 
needing additional bookkeeping while setting the values. If you did not 
preallocate originally, then there is no way to prevent the additional 
mallocs() the second time through, so if you never preallocate but need to add 
many new nonzero locations adding the new nonzeros will be time consuming; 
hence in that situation providing the initial preallocation (taking into 
account all future nonzeros appearing) will pay off in possibly a big way.

I will look at the bug you report for when MatResetPreallocation()is called 
before the first matrix assembly and see if it can be fixed.

Barry



> On Jun 18, 2023, at 9:01 AM, Junchao Zhang  wrote:
> 
> yeah,  see a C example at 
> https://gitlab.com/petsc/petsc/-/blob/main/src/mat/tests/ex259.c
> 
> I guess you can code in this outline with petsc-3.19
> 
> MatCreate
> MatSetSizes
> MatSetFromOptions
> iteration-loop start
> 
> MatResetPreallocation(A,...)
> 
> Fill Matrix A with MatSetValues
> 
> MatAssemblyXX(A_WAVE,MAT_FINAL_ASSEMBLY,IERR)
> 
>  iteration-loop end
> 
> 
> --Junchao Zhang
> 
> 
> On Sun, Jun 18, 2023 at 7:37 AM Matthew Knepley  > wrote:
>> On Sun, Jun 18, 2023 at 4:26 AM Karsten Lettmann 
>> > > wrote:
>>> Dear Matthew,
>>> 
>>> 
>>> 
>>> thanks for you help.
>>> 
>>> 
>>> 
>>> 1) I tested your suggestion to pass NULL as the arguments for the 
>>> MatXAIJSetPreallocation.
>>> 
>>> So the old version was:
>>> 
>>> CALL 
>>> MatCreateMPIAIJ(MPI_GROUP,N_local,N_local,N_global,N_global,0,DNNZ,0,ONNZ,A,IERR)
>>> 
>>> And after you suggestion it is now:
>>> 
>>>   CALL MatCreate(MPI_GROUP,A,IERR)
>>>   CALL MatSetType(A,MATAIJ,IERR)
>>>   CALL MatSetSizes(A,N_local,N_local,N_global,N_global,IERR)
>>>   CALL 
>>> MatXAIJSetPreallocation(A,1,DNNZ,ONNZ,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,IERR)
>>> 
>>> 
>>> 
>>>   Setting block-size = 1.
>>> 
>>> 
>>> 
>>> 2) Concerning the error with MatResetPreallocation:
>>> 
>>> We have an iterative loop, in which the matrix A is filled very often with 
>>> different non-zero structure.
>>> Further, I read in the manual pages that due to performance issues, one 
>>> should preallocate enough space, as operations as matsetvalues might be 
>>> time consuming due to additional 
>>> on-demand allocations.
>>> 
>>> 
>>> 
>>> So I did the following coding in principle:
>>> 
>>> 
>>> 
>>> Set matrix A the first time with preallocation
>>> 
>>> 
>>> 
>>> iteration-loop start
>>> 
>>> MatResetPreallocation(A,...)
>>> 
>>> MatZeroEntries (A)
>>> 
>>> Fill Matrix A
>>> 
>>> MatAssemblyXX(A_WAVE,MAT_FINAL_ASSEMBLY,IERR)
>>> 
>>> 
>>> 
>>> iteration-loop end
>>> 
>>> 
>>> 
>>> With these settings, the code run with 2 CPU.
>>> But with 1 CPU I got an error, which was in MatResetPreallocation.
>>> I could not understand, why the above code works with 2 CPU but not with 1 
>>> CPU.
>>> 
>>> At the moment, I believe the reason for this error seems to be a pre-check, 
>>> that is done in SeqAIJ but not in MPIAIJ fo a valid and present matrix A.
>>> 
>>> (Now, an image is included showing the codings of   :
>>> https://petsc.org/release/src/mat/impls/aij/seq/aij.c.html#MatResetPreallocation_SeqAIJ
>>> https://petsc.org/release/src/mat/impls/aij/mpi/mpiaij.c.html#MatResetPreallocation_MPIAIJ
>>>)
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> So, it seems for me at the moment, that the first MatResetPreallocation 
>>> (when the iteration loop is entered the first time) is done on an 
>>> not-assembled matrix A.
>>> So for one CPU I got an error, while 2 CPUs seem to have been more tolerant.
>>> (I'm not sure, if this interpretation is correct.)
>>> 
>>> 
>>> 
>>> So, I changed the coding in that way, 

Re: [petsc-users] Question about using MatCreateAIJ

[Matthew Knepley]

On Thu, Jun 15, 2023 at 8:32 AM Karsten Lettmann <
karsten.lettm...@uni-oldenburg.de> wrote:

> Dear all,
>
>
> I'm quite new to PETSC. So I hope the following questions are not too
> stupid.
>
>
> 1) We have a (Fortran) code, that we want to update from an older PETSC
> version (petsc.2.3.3-p16) to a newer version.
>
> Inside the old code, for creating matrices A, there are function calls
> of the from:
> MatCreateMPIAIJ
>
> In the reference page for this old version it says:
> When calling this routine with a single process communicator, a matrix
> of type SEQAIJ is returned.
>
> So I assume the following behavior of this old routine:
> - for N_proc == 1:
> a matrix of type SEQAIJ is returned.
>
> - for N_proc > 1:
> a matrix of type MPIAIJ is returned
>
>
>
> 2a) So, in the new code, we want to have a similar behavior.
>
> I found that this function is not present any more in the newer PETSC
> versions.
>
> Instead, one might use: MatCreateAIJ(….)
> ( https://petsc.org/release/manualpages/Mat/MatCreateAIJ/ )
>
> If I understand the reference page of the function correctly, then,
> actually, a similar behavior should be expected:
>
> - for N_proc == 1:
> a matrix of type SEQAIJ is returned.
>
> - for N_proc > 1:
> a matrix of type MPIAIJ is returned
>
>
> 2b) However, on the reference page, there is the note:
>
> It is recommended that one use the MatCreate(), MatSetType() and/or
> MatSetFromOptions(), MatSetPreallocation() paradigm instead of this
> routine directly.
>
> So, if I want the behavior above, it is recommended to code it like
> this, isn't it:
>
> If (N_Proc == 1)
>
>  MatCreate(.. ,A ,...)
>  MatSetType(…,A, MATSEQAIJ,..)
>  MatSetSizes(…,A, ..)
>  MatSeqAIJSetPreallocation(,...A,...)
>
> else
>
>  MatCreate(.. ,A ,...)
>  MatSetType(…,A, MATMPIAIJ,..)
>  MatSetSizes(…,A, ..)
>  MatMPIAIJSetPreallocation(,...A,...)
>

You can use

  MatCreate(comm, );
  MatSetType(A, MATAIJ);
  MatSetSizes(A, ...);
  MatXAIJSetPreallocation(A, ...);

We recommend this because we would like to get rid of the convenience
functions that
wrap up exactly this code.

  Thanks,

 Matt


> end
>
>
>
> 3) So my questions are:
>
> - Is my present understanding correct?
>
> If  yes:
>
> - Why might using MatCreateAIJ(….) for my case not be helpful?
>
> - So, why is it recommended to use the way 2b) instead of this
> MatCreateAIJ(….) ?
>
>
> Best, Karsten
>
>
>
>
> --
> ICBM
> Section: Physical Oceanography
> Universitaet Oldenburg
> Postfach 5634
> D-26046 Oldenburg
> Germany
>
> Tel:+49 (0)441 798 4061
> email:  karsten.lettm...@uni-oldenburg.de
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question on ISLocalToGlobalMappingGetIndices Fortran Interface

[Danyang Su]

Thanks, Mark. Yes, it actually works when I update to 
ISLocalToGlobalMappingGetIndicesF90. I made a mistake reporting this does not 
work.

 

Danyang

 

From: Mark Adams 
Date: Monday, May 8, 2023 at 7:22 PM
To: Danyang Su 
Cc: petsc-users 
Subject: Re: [petsc-users] Question on ISLocalToGlobalMappingGetIndices Fortran 
Interface

 

 

 

On Mon, May 8, 2023 at 6:50 PM Danyang Su  wrote:

Dear PETSc-Users,

 

Is there any changes in ISLocalToGlobalMappingGetIndices function after PETSc 
3.17? 

 

In the previous PETSc version (<= 3.17), the function 
‘ISLocalToGlobalMappingGetIndices(ltogm,ltog,idltog,ierr)’ works fine, even 
though the value of idltog looks out of bound (-11472655627481), 
https://www.mcs.anl.gov/petsc/petsc-3.14/src/ksp/ksp/tutorials/ex14f.F90.html. 
The value of idltog is not clear.

 

In the latest PETSc version,  this function can be called, but due to the 
extreme value of idltog, the code fails. I also tried to use 
‘ISLocalToGlobalMappingGetIndicesF90(ltogm,ltog,ierr)’ () but no success.

 

* You do want the latter:

 

doc/changes/319.rst:- Deprecate ``ISLocalToGlobalMappingGetIndices()`` in favor 
of ``ISLocalToGlobalMappingGetIndicesF90()``

 

* You might look at a test:

 

src/ksp/ksp/tutorials/ex14f.F90:  
PetscCall(ISLocalToGlobalMappingGetIndicesF90(ltogm,ltog,ierr))

 

* If you use 64 bit integers be careful. 

 

* You want to use a memory checker like Valgrind or Sanitize.

 

Mark

 

 

#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR <= 4)

call DMDAGetGlobalIndicesF90(dmda_flow%da,PETSC_NULL_INTEGER,  &

 idx,ierr)

CHKERRQ(ierr)

#else

call DMGetLocalToGlobalMapping(dmda_flow%da,ltogm,ierr)

CHKERRQ(ierr)

call ISLocalToGlobalMappingGetIndices(ltogm,ltog,idltog,ierr)

CHKERRQ(ierr)

#endif



dof = dmda_flow%dof

  

do ivol = 1, nngl

 

#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR <= 4)

node_idx_lg2pg(ivol) = (idx(ivol*dof)+1)/dof

#else

node_idx_lg2pg(ivol) = (ltog(ivol*dof + idltog)+1)/dof

#endif

   end do

 

Any suggestions on that. 

 

Thanks,

 

Danyang

Re: [petsc-users] Question on ISLocalToGlobalMappingGetIndices Fortran Interface

[Mark Adams]

On Mon, May 8, 2023 at 6:50 PM Danyang Su  wrote:

> Dear PETSc-Users,
>
>
>
> Is there any changes in ISLocalToGlobalMappingGetIndices function after
> PETSc 3.17?
>
>
>
> In the previous PETSc version (<= 3.17), the function
> ‘ISLocalToGlobalMappingGetIndices(ltogm,ltog,idltog,ierr)’ works fine, even
> though the value of idltog looks out of bound (-11472655627481),
> https://www.mcs.anl.gov/petsc/petsc-3.14/src/ksp/ksp/tutorials/ex14f.F90.html.
> The value of idltog is not clear.
>
>
>
> In the latest PETSc version,  this function can be called, but due to the
> extreme value of idltog, the code fails. I also tried to use
> ‘ISLocalToGlobalMappingGetIndicesF90(ltogm,ltog,ierr)’ () but no success.
>

* You do want the latter:

doc/changes/319.rst:- Deprecate ``ISLocalToGlobalMappingGetIndices()`` in
favor of ``ISLocalToGlobalMappingGetIndicesF90()``

* You might look at a test:

src/ksp/ksp/tutorials/ex14f.F90:
 PetscCall(ISLocalToGlobalMappingGetIndicesF90(ltogm,ltog,ierr))

* If you use 64 bit integers be careful.

* You want to use a memory checker like Valgrind or Sanitize.

Mark


>
> #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR <= 4)
>
> call DMDAGetGlobalIndicesF90(dmda_flow%da,PETSC_NULL_INTEGER,  &
>
>  idx,ierr)
>
> CHKERRQ(ierr)
>
> #else
>
> call DMGetLocalToGlobalMapping(dmda_flow%da,ltogm,ierr)
>
> CHKERRQ(ierr)
>
> call ISLocalToGlobalMappingGetIndices(ltogm,ltog,idltog,ierr)
>
> CHKERRQ(ierr)
>
> #endif
>
>
>
> dof = dmda_flow%dof
>
>
>
> do ivol = 1, nngl
>
>
>
> #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR <= 4)
>
> node_idx_lg2pg(ivol) = (idx(ivol*dof)+1)/dof
>
> #else
>
> node_idx_lg2pg(ivol) = (ltog(ivol*dof + idltog)+1)/dof
>
> #endif
>
>end do
>
>
>
> Any suggestions on that.
>
>
>
> Thanks,
>
>
>
> Danyang
>

Re: [petsc-users] question about leap-frog for wave equation in petsc

[Matthew Knepley]

On Sat, May 6, 2023 at 11:47 AM Huidong Yang 
wrote:

> Hi Petsc developer.
>
> may I ask if there is any available implementations in petsc
> using leap-frog scheme?
>

I don't think we have leapfrog, but we do have Stormer-Verlet, which is
also a 2nd order symplectic method.

  Thanks,

 Matt


> Thanks.
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about linking LAPACK library

[Barry Smith]

   LAPACK is always linked with PETSc so you can always make direct calls to 
LAPACK routines from PETSc code..

   Barry


> On Apr 25, 2023, at 6:03 AM, ­권승리 / 학생 / 항공우주공학과  wrote:
> 
> Thank you for your reply.
> 
> I think I gave an example of an unrealistic problem.
> 
> I just wanted to know how to compute the inverse matrix, so I was wondering 
> if there is an example of computing the inverse matrix in PETSc.
> 
> Alternatively, I want to know how to link the LAPACK library.
> 
> best,
> 
> Seung Lee Kwon
> 
> 2023년 4월 25일 (화) 오후 6:44, Matthew Knepley  >님이 작성:
>> On Mon, Apr 24, 2023 at 11:47 PM ­권승리 / 학생 / 항공우주공학과 > > wrote:
>>> Dear all
>>> 
>>> It depends on the problem. It can have hundreds of thousands of degrees of 
>>> freedom.
>> 
>> Suppose your matrix was dense and had 1e6 dofs. The work to invert a matrix 
>> is O(N^3) with a small
>> constant, so it would take 1e18 = 1 exaflop to invert this matrix and about 
>> 10 Terabytes of RAM to store
>> it. Is this available to you? PETSc's supports Elemental and SCALAPACK for 
>> this kind of calculation.
>> 
>> If the system is sparse, you could invert it using MUMPS, SuperLU_dist, or 
>> Pardiso. Then the work and
>> storage depend on the density. There are good estimates for connectivity 
>> based on regular grids of given
>> dimension. The limiting resource here is usually memory, which motivates 
>> people to try iterative methods.
>> The convergence of iterative methods depend on detailed properties of your 
>> system, like the operator spectrum.
>> 
>>   Thanks,
>> 
>>  Matt
>>  
>>> best,
>>> 
>>> Seung Lee Kwon
>>> 
>>> 2023년 4월 25일 (화) 오후 12:32, Barry Smith >> >님이 작성:
 
   How large are the dense matrices you would like to invert?
 
> On Apr 24, 2023, at 11:27 PM, ­권승리 / 학생 / 항공우주공학과  > wrote:
> 
> Dear all
> 
> Hello.
> I want to make an inverse matrix like inv(A) in MATLAB.
> 
> Are there some methods to inverse matrix in petsc?
> 
> If not, I want to use the inverse function in the LAPACK library.
> 
> Then, how to use the LAPACK library in petsc? I use the C language.
> 
> Best,
> 
> Seung Lee Kwon
> 
> -- 
> Seung Lee Kwon, Ph.D.Candidate
> Aerospace Structures and Materials Laboratory
> Department of Mechanical and Aerospace Engineering
> Seoul National University
> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
> E-mail : ksl7...@snu.ac.kr 
> Office : +82-2-880-7389
> C. P : +82-10-4695-1062
 
>>> 
>>> 
>>> -- 
>>> Seung Lee Kwon, Ph.D.Candidate
>>> Aerospace Structures and Materials Laboratory
>>> Department of Mechanical and Aerospace Engineering
>>> Seoul National University
>>> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
>>> E-mail : ksl7...@snu.ac.kr 
>>> Office : +82-2-880-7389
>>> C. P : +82-10-4695-1062
>> 
>> 
>> -- 
>> What most experimenters take for granted before they begin their experiments 
>> is infinitely more interesting than any results to which their experiments 
>> lead.
>> -- Norbert Wiener
>> 
>> https://www.cse.buffalo.edu/~knepley/ 
> 
> 
> -- 
> Seung Lee Kwon, Ph.D.Candidate
> Aerospace Structures and Materials Laboratory
> Department of Mechanical and Aerospace Engineering
> Seoul National University
> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
> E-mail : ksl7...@snu.ac.kr 
> Office : +82-2-880-7389
> C. P : +82-10-4695-1062

Re: [petsc-users] Question about linking LAPACK library

[Satish Balay via petsc-users]

On Tue, 25 Apr 2023, Matthew Knepley wrote:

> On Tue, Apr 25, 2023 at 6:03 AM ­권승리 / 학생 / 항공우주공학과 
> wrote:
> 
> > Thank you for your reply.
> >
> > I think I gave an example of an unrealistic problem.
> >
> > I just wanted to know how to compute the inverse matrix, so I was
> > wondering if there is an example of computing the inverse matrix in PETSc.
> >
> 
> You just use
> 
>   KSPCreate(PETSC_COMM_WORLD, );
>   KSPSetOperator(ksp, A, A);
>   KSPSetFromOptions(ksp);
>   KSPSolve(ksp, b, x);
> 
> This solves the system. If you want an exact solution use
> 
>   -pc_type lu
> 
> There is a manual chapter on linear solves such as there.
> 
>   Thanks,
> 
>  Matt
> 
> 
> > Alternatively, I want to know how to link the LAPACK library.

You might want to check:

src/mat/impls/dense/seq/dense.c:  PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(, 
, mat->v, >lda, mat->pivots, ));

Satish


> >
> > best,
> >
> > Seung Lee Kwon
> >
> > 2023년 4월 25일 (화) 오후 6:44, Matthew Knepley 님이 작성:
> >
> >> On Mon, Apr 24, 2023 at 11:47 PM ­권승리 / 학생 / 항공우주공학과 
> >> wrote:
> >>
> >>> Dear all
> >>>
> >>> It depends on the problem. It can have hundreds of thousands of degrees
> >>> of freedom.
> >>>
> >>
> >> Suppose your matrix was dense and had 1e6 dofs. The work to invert a
> >> matrix is O(N^3) with a small
> >> constant, so it would take 1e18 = 1 exaflop to invert this matrix and
> >> about 10 Terabytes of RAM to store
> >> it. Is this available to you? PETSc's supports Elemental and SCALAPACK
> >> for this kind of calculation.
> >>
> >> If the system is sparse, you could invert it using MUMPS, SuperLU_dist,
> >> or Pardiso. Then the work and
> >> storage depend on the density. There are good estimates for connectivity
> >> based on regular grids of given
> >> dimension. The limiting resource here is usually memory, which motivates
> >> people to try iterative methods.
> >> The convergence of iterative methods depend on detailed properties of
> >> your system, like the operator spectrum.
> >>
> >>   Thanks,
> >>
> >>  Matt
> >>
> >>
> >>> best,
> >>>
> >>> Seung Lee Kwon
> >>>
> >>> 2023년 4월 25일 (화) 오후 12:32, Barry Smith 님이 작성:
> >>>
> 
>    How large are the dense matrices you would like to invert?
> 
>  On Apr 24, 2023, at 11:27 PM, ­권승리 / 학생 / 항공우주공학과 
>  wrote:
> 
>  Dear all
> 
>  Hello.
>  I want to make an inverse matrix like inv(A) in MATLAB.
> 
>  Are there some methods to inverse matrix in petsc?
> 
>  If not, I want to use the inverse function in the LAPACK library.
> 
>  Then, how to use the LAPACK library in petsc? I use the C language.
> 
>  Best,
> 
>  Seung Lee Kwon
> 
>  --
>  Seung Lee Kwon, Ph.D.Candidate
>  Aerospace Structures and Materials Laboratory
>  Department of Mechanical and Aerospace Engineering
>  Seoul National University
>  Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
>  E-mail : ksl7...@snu.ac.kr
>  Office : +82-2-880-7389
>  C. P : +82-10-4695-1062
> 
> 
> 
> >>>
> >>> --
> >>> Seung Lee Kwon, Ph.D.Candidate
> >>> Aerospace Structures and Materials Laboratory
> >>> Department of Mechanical and Aerospace Engineering
> >>> Seoul National University
> >>> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
> >>> E-mail : ksl7...@snu.ac.kr
> >>> Office : +82-2-880-7389
> >>> C. P : +82-10-4695-1062
> >>>
> >>
> >>
> >> --
> >> What most experimenters take for granted before they begin their
> >> experiments is infinitely more interesting than any results to which their
> >> experiments lead.
> >> -- Norbert Wiener
> >>
> >> https://www.cse.buffalo.edu/~knepley/
> >> 
> >>
> >
> >
> > --
> > Seung Lee Kwon, Ph.D.Candidate
> > Aerospace Structures and Materials Laboratory
> > Department of Mechanical and Aerospace Engineering
> > Seoul National University
> > Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
> > E-mail : ksl7...@snu.ac.kr
> > Office : +82-2-880-7389
> > C. P : +82-10-4695-1062
> >
> 
> 
> 

Re: [petsc-users] Question about linking LAPACK library

[Matthew Knepley]

On Tue, Apr 25, 2023 at 6:03 AM ­권승리 / 학생 / 항공우주공학과 
wrote:

> Thank you for your reply.
>
> I think I gave an example of an unrealistic problem.
>
> I just wanted to know how to compute the inverse matrix, so I was
> wondering if there is an example of computing the inverse matrix in PETSc.
>

You just use

  KSPCreate(PETSC_COMM_WORLD, );
  KSPSetOperator(ksp, A, A);
  KSPSetFromOptions(ksp);
  KSPSolve(ksp, b, x);

This solves the system. If you want an exact solution use

  -pc_type lu

There is a manual chapter on linear solves such as there.

  Thanks,

 Matt


> Alternatively, I want to know how to link the LAPACK library.
>
> best,
>
> Seung Lee Kwon
>
> 2023년 4월 25일 (화) 오후 6:44, Matthew Knepley 님이 작성:
>
>> On Mon, Apr 24, 2023 at 11:47 PM ­권승리 / 학생 / 항공우주공학과 
>> wrote:
>>
>>> Dear all
>>>
>>> It depends on the problem. It can have hundreds of thousands of degrees
>>> of freedom.
>>>
>>
>> Suppose your matrix was dense and had 1e6 dofs. The work to invert a
>> matrix is O(N^3) with a small
>> constant, so it would take 1e18 = 1 exaflop to invert this matrix and
>> about 10 Terabytes of RAM to store
>> it. Is this available to you? PETSc's supports Elemental and SCALAPACK
>> for this kind of calculation.
>>
>> If the system is sparse, you could invert it using MUMPS, SuperLU_dist,
>> or Pardiso. Then the work and
>> storage depend on the density. There are good estimates for connectivity
>> based on regular grids of given
>> dimension. The limiting resource here is usually memory, which motivates
>> people to try iterative methods.
>> The convergence of iterative methods depend on detailed properties of
>> your system, like the operator spectrum.
>>
>>   Thanks,
>>
>>  Matt
>>
>>
>>> best,
>>>
>>> Seung Lee Kwon
>>>
>>> 2023년 4월 25일 (화) 오후 12:32, Barry Smith 님이 작성:
>>>

   How large are the dense matrices you would like to invert?

 On Apr 24, 2023, at 11:27 PM, ­권승리 / 학생 / 항공우주공학과 
 wrote:

 Dear all

 Hello.
 I want to make an inverse matrix like inv(A) in MATLAB.

 Are there some methods to inverse matrix in petsc?

 If not, I want to use the inverse function in the LAPACK library.

 Then, how to use the LAPACK library in petsc? I use the C language.

 Best,

 Seung Lee Kwon

 --
 Seung Lee Kwon, Ph.D.Candidate
 Aerospace Structures and Materials Laboratory
 Department of Mechanical and Aerospace Engineering
 Seoul National University
 Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
 E-mail : ksl7...@snu.ac.kr
 Office : +82-2-880-7389
 C. P : +82-10-4695-1062



>>>
>>> --
>>> Seung Lee Kwon, Ph.D.Candidate
>>> Aerospace Structures and Materials Laboratory
>>> Department of Mechanical and Aerospace Engineering
>>> Seoul National University
>>> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
>>> E-mail : ksl7...@snu.ac.kr
>>> Office : +82-2-880-7389
>>> C. P : +82-10-4695-1062
>>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> 
>>
>
>
> --
> Seung Lee Kwon, Ph.D.Candidate
> Aerospace Structures and Materials Laboratory
> Department of Mechanical and Aerospace Engineering
> Seoul National University
> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
> E-mail : ksl7...@snu.ac.kr
> Office : +82-2-880-7389
> C. P : +82-10-4695-1062
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about linking LAPACK library

[­권승리 / 학생 / 항공우주공학과]

Thank you for your reply.

I think I gave an example of an unrealistic problem.

I just wanted to know how to compute the inverse matrix, so I was wondering
if there is an example of computing the inverse matrix in PETSc.

Alternatively, I want to know how to link the LAPACK library.

best,

Seung Lee Kwon

2023년 4월 25일 (화) 오후 6:44, Matthew Knepley 님이 작성:

> On Mon, Apr 24, 2023 at 11:47 PM ­권승리 / 학생 / 항공우주공학과 
> wrote:
>
>> Dear all
>>
>> It depends on the problem. It can have hundreds of thousands of degrees
>> of freedom.
>>
>
> Suppose your matrix was dense and had 1e6 dofs. The work to invert a
> matrix is O(N^3) with a small
> constant, so it would take 1e18 = 1 exaflop to invert this matrix and
> about 10 Terabytes of RAM to store
> it. Is this available to you? PETSc's supports Elemental and SCALAPACK for
> this kind of calculation.
>
> If the system is sparse, you could invert it using MUMPS, SuperLU_dist, or
> Pardiso. Then the work and
> storage depend on the density. There are good estimates for connectivity
> based on regular grids of given
> dimension. The limiting resource here is usually memory, which motivates
> people to try iterative methods.
> The convergence of iterative methods depend on detailed properties of your
> system, like the operator spectrum.
>
>   Thanks,
>
>  Matt
>
>
>> best,
>>
>> Seung Lee Kwon
>>
>> 2023년 4월 25일 (화) 오후 12:32, Barry Smith 님이 작성:
>>
>>>
>>>   How large are the dense matrices you would like to invert?
>>>
>>> On Apr 24, 2023, at 11:27 PM, ­권승리 / 학생 / 항공우주공학과 
>>> wrote:
>>>
>>> Dear all
>>>
>>> Hello.
>>> I want to make an inverse matrix like inv(A) in MATLAB.
>>>
>>> Are there some methods to inverse matrix in petsc?
>>>
>>> If not, I want to use the inverse function in the LAPACK library.
>>>
>>> Then, how to use the LAPACK library in petsc? I use the C language.
>>>
>>> Best,
>>>
>>> Seung Lee Kwon
>>>
>>> --
>>> Seung Lee Kwon, Ph.D.Candidate
>>> Aerospace Structures and Materials Laboratory
>>> Department of Mechanical and Aerospace Engineering
>>> Seoul National University
>>> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
>>> E-mail : ksl7...@snu.ac.kr
>>> Office : +82-2-880-7389
>>> C. P : +82-10-4695-1062
>>>
>>>
>>>
>>
>> --
>> Seung Lee Kwon, Ph.D.Candidate
>> Aerospace Structures and Materials Laboratory
>> Department of Mechanical and Aerospace Engineering
>> Seoul National University
>> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
>> E-mail : ksl7...@snu.ac.kr
>> Office : +82-2-880-7389
>> C. P : +82-10-4695-1062
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>


-- 
Seung Lee Kwon, Ph.D.Candidate
Aerospace Structures and Materials Laboratory
Department of Mechanical and Aerospace Engineering
Seoul National University
Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
E-mail : ksl7...@snu.ac.kr
Office : +82-2-880-7389
C. P : +82-10-4695-1062

Re: [petsc-users] Question about linking LAPACK library

[Pierre Jolivet]

> On 25 Apr 2023, at 11:43 AM, Matthew Knepley  wrote:
> 
> On Mon, Apr 24, 2023 at 11:47 PM ­권승리 / 학생 / 항공우주공학과  > wrote:
>> Dear all
>> 
>> It depends on the problem. It can have hundreds of thousands of degrees of 
>> freedom.
> 
> Suppose your matrix was dense and had 1e6 dofs. The work to invert a matrix 
> is O(N^3) with a small
> constant, so it would take 1e18 = 1 exaflop to invert this matrix and about 
> 10 Terabytes of RAM to store
> it. Is this available to you? PETSc's supports Elemental and SCALAPACK for 
> this kind of calculation.
> 
> If the system is sparse, you could invert it using MUMPS, SuperLU_dist, or 
> Pardiso. Then the work and
> storage depend on the density. There are good estimates for connectivity 
> based on regular grids of given
> dimension. The limiting resource here is usually memory, which motivates 
> people to try iterative methods.
> The convergence of iterative methods depend on detailed properties of your 
> system, like the operator spectrum.

And to wrap this up, if your operator is truly dense, e.g., BEM or non-local 
discretizations, their are available hierarchical formats such as MatH2Opus and 
MatHtool.
They have efficient matrix-vector product implementations such that you can 
solve linear systems without having to invert (or even store) the coefficient 
matrix explicitly.

Thanks,
Pierre

>   Thanks,
> 
>  Matt
>  
>> best,
>> 
>> Seung Lee Kwon
>> 
>> 2023년 4월 25일 (화) 오후 12:32, Barry Smith > >님이 작성:
>>> 
>>>   How large are the dense matrices you would like to invert?
>>> 
 On Apr 24, 2023, at 11:27 PM, ­권승리 / 학생 / 항공우주공학과 >>> > wrote:
 
 Dear all
 
 Hello.
 I want to make an inverse matrix like inv(A) in MATLAB.
 
 Are there some methods to inverse matrix in petsc?
 
 If not, I want to use the inverse function in the LAPACK library.
 
 Then, how to use the LAPACK library in petsc? I use the C language.
 
 Best,
 
 Seung Lee Kwon
 
 -- 
 Seung Lee Kwon, Ph.D.Candidate
 Aerospace Structures and Materials Laboratory
 Department of Mechanical and Aerospace Engineering
 Seoul National University
 Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
 E-mail : ksl7...@snu.ac.kr 
 Office : +82-2-880-7389
 C. P : +82-10-4695-1062
>>> 
>> 
>> 
>> -- 
>> Seung Lee Kwon, Ph.D.Candidate
>> Aerospace Structures and Materials Laboratory
>> Department of Mechanical and Aerospace Engineering
>> Seoul National University
>> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
>> E-mail : ksl7...@snu.ac.kr 
>> Office : +82-2-880-7389
>> C. P : +82-10-4695-1062
> 
> 
> -- 
> What most experimenters take for granted before they begin their experiments 
> is infinitely more interesting than any results to which their experiments 
> lead.
> -- Norbert Wiener
> 
> https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about linking LAPACK library

[Matthew Knepley]

On Mon, Apr 24, 2023 at 11:47 PM ­권승리 / 학생 / 항공우주공학과 
wrote:

> Dear all
>
> It depends on the problem. It can have hundreds of thousands of degrees of
> freedom.
>

Suppose your matrix was dense and had 1e6 dofs. The work to invert a matrix
is O(N^3) with a small
constant, so it would take 1e18 = 1 exaflop to invert this matrix and about
10 Terabytes of RAM to store
it. Is this available to you? PETSc's supports Elemental and SCALAPACK for
this kind of calculation.

If the system is sparse, you could invert it using MUMPS, SuperLU_dist, or
Pardiso. Then the work and
storage depend on the density. There are good estimates for connectivity
based on regular grids of given
dimension. The limiting resource here is usually memory, which motivates
people to try iterative methods.
The convergence of iterative methods depend on detailed properties of your
system, like the operator spectrum.

  Thanks,

 Matt


> best,
>
> Seung Lee Kwon
>
> 2023년 4월 25일 (화) 오후 12:32, Barry Smith 님이 작성:
>
>>
>>   How large are the dense matrices you would like to invert?
>>
>> On Apr 24, 2023, at 11:27 PM, ­권승리 / 학생 / 항공우주공학과 
>> wrote:
>>
>> Dear all
>>
>> Hello.
>> I want to make an inverse matrix like inv(A) in MATLAB.
>>
>> Are there some methods to inverse matrix in petsc?
>>
>> If not, I want to use the inverse function in the LAPACK library.
>>
>> Then, how to use the LAPACK library in petsc? I use the C language.
>>
>> Best,
>>
>> Seung Lee Kwon
>>
>> --
>> Seung Lee Kwon, Ph.D.Candidate
>> Aerospace Structures and Materials Laboratory
>> Department of Mechanical and Aerospace Engineering
>> Seoul National University
>> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
>> E-mail : ksl7...@snu.ac.kr
>> Office : +82-2-880-7389
>> C. P : +82-10-4695-1062
>>
>>
>>
>
> --
> Seung Lee Kwon, Ph.D.Candidate
> Aerospace Structures and Materials Laboratory
> Department of Mechanical and Aerospace Engineering
> Seoul National University
> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
> E-mail : ksl7...@snu.ac.kr
> Office : +82-2-880-7389
> C. P : +82-10-4695-1062
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about linking LAPACK library

[­권승리 / 학생 / 항공우주공학과]

Dear all

It depends on the problem. It can have hundreds of thousands of degrees of
freedom.

best,

Seung Lee Kwon

2023년 4월 25일 (화) 오후 12:32, Barry Smith 님이 작성:

>
>   How large are the dense matrices you would like to invert?
>
> On Apr 24, 2023, at 11:27 PM, ­권승리 / 학생 / 항공우주공학과 
> wrote:
>
> Dear all
>
> Hello.
> I want to make an inverse matrix like inv(A) in MATLAB.
>
> Are there some methods to inverse matrix in petsc?
>
> If not, I want to use the inverse function in the LAPACK library.
>
> Then, how to use the LAPACK library in petsc? I use the C language.
>
> Best,
>
> Seung Lee Kwon
>
> --
> Seung Lee Kwon, Ph.D.Candidate
> Aerospace Structures and Materials Laboratory
> Department of Mechanical and Aerospace Engineering
> Seoul National University
> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
> E-mail : ksl7...@snu.ac.kr
> Office : +82-2-880-7389
> C. P : +82-10-4695-1062
>
>
>

-- 
Seung Lee Kwon, Ph.D.Candidate
Aerospace Structures and Materials Laboratory
Department of Mechanical and Aerospace Engineering
Seoul National University
Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
E-mail : ksl7...@snu.ac.kr
Office : +82-2-880-7389
C. P : +82-10-4695-1062

Re: [petsc-users] Question about linking LAPACK library

[Barry Smith]

  How large are the dense matrices you would like to invert?

> On Apr 24, 2023, at 11:27 PM, ­권승리 / 학생 / 항공우주공학과  wrote:
> 
> Dear all
> 
> Hello.
> I want to make an inverse matrix like inv(A) in MATLAB.
> 
> Are there some methods to inverse matrix in petsc?
> 
> If not, I want to use the inverse function in the LAPACK library.
> 
> Then, how to use the LAPACK library in petsc? I use the C language.
> 
> Best,
> 
> Seung Lee Kwon
> 
> -- 
> Seung Lee Kwon, Ph.D.Candidate
> Aerospace Structures and Materials Laboratory
> Department of Mechanical and Aerospace Engineering
> Seoul National University
> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
> E-mail : ksl7...@snu.ac.kr 
> Office : +82-2-880-7389
> C. P : +82-10-4695-1062

Re: [petsc-users] question about MatSetLocalToGlobalMapping

[Zhang, Hong via petsc-users]

Pierre,


1) Is there any hope to get PDIPDM to use a MatNest?

KKT matrix is indefinite and ill-conditioned, which must be solved using a 
direct matrix factorization method.
But you are using PCBJACOBI in the paper you attached?
In any case, there are many such systems, e.g., a discretization of Stokes 
equations, that can be solved with something else than a direct factorization.

The KKT in this paper is very special, each block represents a time step with 
loose and weak couplings. We tested larger and slightly varying physical 
parameters, PCZBJACOBI fails convergence while CHOLESKY encounters out of 
memory.
For the current implementation, we use MUMPS Cholesky as default. To use 
MatNest, what direct solver to use, SCHUR_FACTOR? I do not know how to get it 
work.
On the one hand, MatNest can efficiently convert to AIJ or SBAIJ if you want to 
stick to PCLU or PCCHOLESKY.
On the other hand, it allows to easily switch to PCFIELDSPLIT which can be used 
to solve saddle-point problems.

MatNest would be a natural way for these type of applications. PETSc has 
MatConvert_Nest_AIJ(), which could enable users to assemble their matrices in 
MatNest, and apply mumps/superlu direct solvers. I'm not concerned about the 
overhead of matrix conversion, because matrix factorization dominates 
computation in general.
2) Is this fixed 
https://lists.mcs.anl.gov/pipermail/petsc-dev/2020-September/026398.html ?
I cannot get users to transition away from Ipopt because of these two missing 
features.

The existing pdipm is the result of a MS student intern project. None of us 
involved are experts on the optimization solvers. We made a straightforward 
parallelization of Ipopt. It indeed needs further work, e.g., more features, 
better matrix storage, convergence criteria... To our knowledge, parallel pdipm 
is not available other than our pdipm.
Ipopt can use MUMPS and PARDISO internally, so it’s in some sense parallel 
(using shared memory).
Also, this is not a very potent selling point.
My users that are satisfied with Ipopt as a "non-parallel" black box don’t want 
to have to touch part of their code just to stick it in a parallel black box 
which is limited to the same kind of linear solver and which has severe 
limitations with respect to Hessian/Jacobian/constraint distributions.

We saw exiting 'parallel' pdipm papers, which conduct sequential computation 
and send KKT matrices to multiple processors, then they show scaling results on 
matrix computation only. Again, I wish someone can help to improve tao/pdipm. 
This is a useful solver.
We should improve our pdipm.
Hong

On 20 Apr 2023, at 5:47 PM, Zhang, Hong via petsc-users 
mailto:petsc-users@mcs.anl.gov>> wrote:

Karthik,
We built a KKT matrix in  TaoSNESJacobian_PDIPM() (see 
petsc/src/tao/constrained/impls/ipm/pdipm.c) which assembles several small 
matrices into a large KKT matrix in mpiaij format. You could take the same 
approach to insert P and P^T into your K.
FYI, I attached our paper.
Hong


From: petsc-users 
mailto:petsc-users-boun...@mcs.anl.gov>> on 
behalf of Matthew Knepley mailto:knep...@gmail.com>>
Sent: Thursday, April 20, 2023 5:37 AM
To: Karthikeyan Chockalingam - STFC UKRI 
mailto:karthikeyan.chockalin...@stfc.ac.uk>>
Cc: petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov> 
mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] question about MatSetLocalToGlobalMapping

On Thu, Apr 20, 2023 at 6:13 AM Karthikeyan Chockalingam - STFC UKRI via 
petsc-users mailto:petsc-users@mcs.anl.gov>> wrote:
Hello,

I created a new thread, thought would it be more appropriate (and is a 
continuation of my previous post). I want to construct the below K matrix 
(which is composed of submatrices)

K = [A P^T
   P   0]

Where K is of type MatMPIAIJ. I first constructed the top left [A] using 
MatSetValues().

Now, I would like to construct the bottom left [p] and top right [p^T] using 
MatSetValuesLocal().

To use  MatSetValuesLocal(),  I first have to create a local-to-global mapping 
using ISLocalToGlobalMappingCreate. I have created two mapping row_mapping and 
column_mapping.

I do not understand why they are not the same map. Maybe I was unclear before. 
It looks like you have two fields, say phi and lambda, where lambda is a 
Lagrange multiplier imposing some constraint. Then you get a saddle point like 
this. You can imagine matrices

  (phi, phi)--> A
  (phi, lambda) --> P^T
  (lambda, phi) --> P

So you make a L2G map for the phi field and the lambda field. Oh, you are 
calling them row and col map, but they are my phi and lambda
maps. I do not like the row and col names since in P they reverse.

Q1) At what point should I declare MatSetLocalToGlobalMapping – is it just 
before I use MatSetValuesLocal()?

Okay, it is good you are asking this because my thinking was somewhat confused. 
I think the pr

Re: [petsc-users] question about MatSetLocalToGlobalMapping

[Karthikeyan Chockalingam - STFC UKRI via petsc-users]

Thank you Hong. I will look into it.

Best,
Karthik.

From: Zhang, Hong 
Date: Thursday, 20 April 2023 at 16:47
To: Matthew Knepley , Chockalingam, Karthikeyan (STFC,DL,HC) 

Cc: petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] question about MatSetLocalToGlobalMapping
Karthik,
We built a KKT matrix in  TaoSNESJacobian_PDIPM() (see 
petsc/src/tao/constrained/impls/ipm/pdipm.c) which assembles several small 
matrices into a large KKT matrix in mpiaij format. You could take the same 
approach to insert P and P^T into your K.
FYI, I attached our paper.
Hong

From: petsc-users  on behalf of Matthew 
Knepley 
Sent: Thursday, April 20, 2023 5:37 AM
To: Karthikeyan Chockalingam - STFC UKRI 
Cc: petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] question about MatSetLocalToGlobalMapping

On Thu, Apr 20, 2023 at 6:13 AM Karthikeyan Chockalingam - STFC UKRI via 
petsc-users mailto:petsc-users@mcs.anl.gov>> wrote:

Hello,



I created a new thread, thought would it be more appropriate (and is a 
continuation of my previous post). I want to construct the below K matrix 
(which is composed of submatrices)



K = [A P^T

   P   0]



Where K is of type MatMPIAIJ. I first constructed the top left [A] using 
MatSetValues().



Now, I would like to construct the bottom left [p] and top right [p^T] using 
MatSetValuesLocal().



To use  MatSetValuesLocal(),  I first have to create a local-to-global mapping 
using ISLocalToGlobalMappingCreate. I have created two mapping row_mapping and 
column_mapping.

I do not understand why they are not the same map. Maybe I was unclear before. 
It looks like you have two fields, say phi and lambda, where lambda is a 
Lagrange multiplier imposing some constraint. Then you get a saddle point like 
this. You can imagine matrices

  (phi, phi)--> A
  (phi, lambda) --> P^T
  (lambda, phi) --> P

So you make a L2G map for the phi field and the lambda field. Oh, you are 
calling them row and col map, but they are my phi and lambda
maps. I do not like the row and col names since in P they reverse.


Q1) At what point should I declare MatSetLocalToGlobalMapping – is it just 
before I use MatSetValuesLocal()?

Okay, it is good you are asking this because my thinking was somewhat confused. 
I think the precise steps are:

  1) Create the large saddle point matrix K

1a) We must call 
https://petsc.org/main/manualpages/Mat/MatSetLocalToGlobalMapping/ on it. In 
the simplest case, this just maps
   the local rows numbers [0, Nrows) to the global rows numbers 
[rowStart, rowStart + Nrows).

  2) To form each piece:

2a) Extract that block using 
https://petsc.org/main/manualpages/Mat/MatGetLocalSubMatrix/

   This gives back a Mat object that you subsequently restore using 
https://petsc.org/main/manualpages/Mat/MatRestoreLocalSubMatrix/

 2b) Insert values using 
https://petsc.org/main/manualpages/Mat/MatSetValuesLocal/

The local indices used for insertion here are indices relative to 
the block itself, and the L2G map for this matrix
has been rewritten to insert into that block in the larger matrix. 
Thus this looks like just calling MatSetValuesLocal()
on the smaller matrix block, but inserts correctly into the larger 
matrix.

Therefore, the code you write code in 2) could work equally well making the 
large matrix from 1), or independent smaller matrix blocks.

Does this make sense?

  Thanks,

 Matt


I will use MatSetLocalToGlobalMapping(K, row_mapping, column_mapping) to build 
the bottom left [P].





Q2) Can now I reset the mapping as MatSetLocalToGlobalMapping(K, 
column_mapping, row_mapping) to build the top right [P^T]?





Many thanks!



Kind regards,

Karthik.




















--
What most experimenters take for granted before they begin their experiments is 
infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>

Re: [petsc-users] question about MatSetLocalToGlobalMapping

[Karthikeyan Chockalingam - STFC UKRI via petsc-users]

Thank you, Matt. It took me a while but I understand how to work with 
submatrices.

(Below is simple stand-alone code which assembles two submatrices into a matrix 
of 5 x 5, which might be of use to someone else)

Q1) Looks like I don’t need to call MatrixAssembleBegin/End of on submatrices?
Q2) Though I have preallocated A for 5 x 5 – why do I still need to call

  MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
  else I get an error while assembling using MatSetValuesLocal()?
Q3) Am I calling MatRestoreLocalSubmatrix() at the right place and what does it 
do?
Q4) When assembling a large submatrix, can I call MatGetOwershipRange on the 
submatrix?

Best,
Karthik.





  #include 

  MatA;/* linear system matrix */

  PetscInt   i,j,m = 3,n = 3, P = 5, Q = 5;

  PetscErrorCode ierr;

  PetscScalarv;



  // create matrix A

  ierr = MatCreate(PETSC_COMM_WORLD,);CHKERRQ(ierr);

  ierr = MatSetType(A, MATMPIAIJ); CHKERRQ(ierr);

  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,P,Q);CHKERRQ(ierr);

  ierr = MatSetFromOptions(A);CHKERRQ(ierr);

  ierr = MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);CHKERRQ(ierr);

  //ierr = MatSeqAIJSetPreallocation(A,5,NULL);CHKERRQ(ierr);



  // local indices is always 0  1  2

 PetscInt  c = 3, global_col_indices[] = {0, 1, 2};

 PetscInt  r = 2, global_row_indices[] = {3, 4};



  ISLocalToGlobalMapping col_mapping, row_mapping;



  ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, c, global_col_indices, 
PETSC_COPY_VALUES, _mapping);

  ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, r, global_row_indices, 
PETSC_COPY_VALUES, _mapping);

  MatSetLocalToGlobalMapping(A,row_mapping, col_mapping);



  for (i=0; i
Date: Thursday, 20 April 2023 at 11:37
To: Chockalingam, Karthikeyan (STFC,DL,HC) 
Cc: petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] question about MatSetLocalToGlobalMapping
On Thu, Apr 20, 2023 at 6:13 AM Karthikeyan Chockalingam - STFC UKRI via 
petsc-users mailto:petsc-users@mcs.anl.gov>> wrote:
Hello,

I created a new thread, thought would it be more appropriate (and is a 
continuation of my previous post). I want to construct the below K matrix 
(which is composed of submatrices)

K = [A P^T
   P   0]

Where K is of type MatMPIAIJ. I first constructed the top left [A] using 
MatSetValues().

Now, I would like to construct the bottom left [p] and top right [p^T] using 
MatSetValuesLocal().

To use  MatSetValuesLocal(),  I first have to create a local-to-global mapping 
using ISLocalToGlobalMappingCreate. I have created two mapping row_mapping and 
column_mapping.

I do not understand why they are not the same map. Maybe I was unclear before. 
It looks like you have two fields, say phi and lambda, where lambda is a 
Lagrange multiplier imposing some constraint. Then you get a saddle point like 
this. You can imagine matrices

  (phi, phi)--> A
  (phi, lambda) --> P^T
  (lambda, phi) --> P

So you make a L2G map for the phi field and the lambda field. Oh, you are 
calling them row and col map, but they are my phi and lambda
maps. I do not like the row and col names since in P they reverse.

Q1) At what point should I declare MatSetLocalToGlobalMapping – is it just 
before I use MatSetValuesLocal()?

Okay, it is good you are asking this because my thinking was somewhat confused. 
I think the precise steps are:

  1) Create the large saddle point matrix K

1a) We must call 
https://petsc.org/main/manualpages/Mat/MatSetLocalToGlobalMapping/ on it. In 
the simplest case, this just maps
   the local rows numbers [0, Nrows) to the global rows numbers 
[rowStart, rowStart + Nrows).

  2) To form each piece:

2a) Extract that block using 
https://petsc.org/main/manualpages/Mat/MatGetLocalSubMatrix/

   This gives back a Mat object that you subsequently restore using 
https://petsc.org/main/manualpages/Mat/MatRestoreLocalSubMatrix/

 2b) Insert values using 
https://petsc.org/main/manualpages/Mat/MatSetValuesLocal/

The local indices used for insertion here are indices relative to 
the block itself, and the L2G map for this matrix
has been rewritten to insert into that block in the larger matrix. 
Thus this looks like just calling MatSetValuesLocal()
on the smaller matrix block, but inserts correctly into the larger 
matrix.

Therefore, the code you write code in 2) could work equally well making the 
large matrix from 1), or independent smaller matrix blocks.

Does this make sense?

  Thanks,

 Matt

I will use MatSetLocalToGlobalMapping(K, row_mapping, column_mapping) to build 
the bottom left [P].


Q2) Can now I reset the mapping as MatSetLocalToGlobalMapping(K, 
column_mapping, row_mapping) to build the top right [P^T]?


Many thanks!

Kind regards,
Karthik.











--
What most experimenters take for granted before they begin th

Re: [petsc-users] question about MatSetLocalToGlobalMapping

[Pierre Jolivet]

> On 20 Apr 2023, at 10:28 PM, Zhang, Hong  wrote:
> 
> Pierre,
> 1) Is there any hope to get PDIPDM to use a MatNest?
> 
> KKT matrix is indefinite and ill-conditioned, which must be solved using a 
> direct matrix factorization method.

But you are using PCBJACOBI in the paper you attached?
In any case, there are many such systems, e.g., a discretization of Stokes 
equations, that can be solved with something else than a direct factorization.

> For the current implementation, we use MUMPS Cholesky as default. To use 
> MatNest, what direct solver to use, SCHUR_FACTOR? I do not know how to get it 
> work. 

On the one hand, MatNest can efficiently convert to AIJ or SBAIJ if you want to 
stick to PCLU or PCCHOLESKY.
On the other hand, it allows to easily switch to PCFIELDSPLIT which can be used 
to solve saddle-point problems.

> 2) Is this fixed 
> https://lists.mcs.anl.gov/pipermail/petsc-dev/2020-September/026398.html ?
> I cannot get users to transition away from Ipopt because of these two missing 
> features.
> 
> The existing pdipm is the result of a MS student intern project. None of us 
> involved are experts on the optimization solvers. We made a straightforward 
> parallelization of Ipopt. It indeed needs further work, e.g., more features, 
> better matrix storage, convergence criteria... To our knowledge, parallel 
> pdipm is not available other than our pdipm.

Ipopt can use MUMPS and PARDISO internally, so it’s in some sense parallel 
(using shared memory).
Also, this is not a very potent selling point.
My users that are satisfied with Ipopt as a "non-parallel" black box don’t want 
to have to touch part of their code just to stick it in a parallel black box 
which is limited to the same kind of linear solver and which has severe 
limitations with respect to Hessian/Jacobian/constraint distributions.

Thanks,
Pierre

> We should improve our pdipm. 
> Hong
> 
>> On 20 Apr 2023, at 5:47 PM, Zhang, Hong via petsc-users 
>> mailto:petsc-users@mcs.anl.gov>> wrote:
>> 
>> Karthik,
>> We built a KKT matrix in  TaoSNESJacobian_PDIPM() (see 
>> petsc/src/tao/constrained/impls/ipm/pdipm.c) which assembles several small 
>> matrices into a large KKT matrix in mpiaij format. You could take the same 
>> approach to insert P and P^T into your K.
>> FYI, I attached our paper.
>> Hong
>>  
>> From: petsc-users > <mailto:petsc-users-boun...@mcs.anl.gov>> on behalf of Matthew Knepley 
>> mailto:knep...@gmail.com>>
>> Sent: Thursday, April 20, 2023 5:37 AM
>> To: Karthikeyan Chockalingam - STFC UKRI 
>> > <mailto:karthikeyan.chockalin...@stfc.ac.uk>>
>> Cc: petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov> 
>> mailto:petsc-users@mcs.anl.gov>>
>> Subject: Re: [petsc-users] question about MatSetLocalToGlobalMapping
>>  
>> On Thu, Apr 20, 2023 at 6:13 AM Karthikeyan Chockalingam - STFC UKRI via 
>> petsc-users mailto:petsc-users@mcs.anl.gov>> wrote:
>> Hello,
>>  
>> I created a new thread, thought would it be more appropriate (and is a 
>> continuation of my previous post). I want to construct the below K matrix 
>> (which is composed of submatrices)
>>  
>> K = [A P^T
>>P   0]
>>  
>> Where K is of type MatMPIAIJ. I first constructed the top left [A] using 
>> MatSetValues().
>>  
>> Now, I would like to construct the bottom left [p] and top right [p^T] using 
>> MatSetValuesLocal().
>>  
>> To use  MatSetValuesLocal(),  I first have to create a local-to-global 
>> mapping using ISLocalToGlobalMappingCreate. I have created two mapping 
>> row_mapping and column_mapping.
>> 
>> I do not understand why they are not the same map. Maybe I was unclear 
>> before. It looks like you have two fields, say phi and lambda, where lambda 
>> is a Lagrange multiplier imposing some constraint. Then you get a saddle 
>> point like this. You can imagine matrices
>> 
>>   (phi, phi)--> A
>>   (phi, lambda) --> P^T
>>   (lambda, phi) --> P
>> 
>> So you make a L2G map for the phi field and the lambda field. Oh, you are 
>> calling them row and col map, but they are my phi and lambda
>> maps. I do not like the row and col names since in P they reverse.
>>  
>> Q1) At what point should I declare MatSetLocalToGlobalMapping – is it just 
>> before I use MatSetValuesLocal()?
>> 
>> Okay, it is good you are asking this because my thinking was somewhat 
>> confused. I think the precise steps are:
>> 
>>   1) Create the large saddle point matrix K
>> 
>> 1a) We must call 
>> https://petsc.or

Re: [petsc-users] question about MatSetLocalToGlobalMapping

[Zhang, Hong via petsc-users]

Pierre,
1) Is there any hope to get PDIPDM to use a MatNest?

KKT matrix is indefinite and ill-conditioned, which must be solved using a 
direct matrix factorization method. For the current implementation, we use 
MUMPS Cholesky as default. To use MatNest, what direct solver to use, 
SCHUR_FACTOR? I do not know how to get it work.

2) Is this fixed 
https://lists.mcs.anl.gov/pipermail/petsc-dev/2020-September/026398.html ?
I cannot get users to transition away from Ipopt because of these two missing 
features.

The existing pdipm is the result of a MS student intern project. None of us 
involved are experts on the optimization solvers. We made a straightforward 
parallelization of Ipopt. It indeed needs further work, e.g., more features, 
better matrix storage, convergence criteria... To our knowledge, parallel pdipm 
is not available other than our pdipm. We should improve our pdipm.
Hong

On 20 Apr 2023, at 5:47 PM, Zhang, Hong via petsc-users 
 wrote:

Karthik,
We built a KKT matrix in  TaoSNESJacobian_PDIPM() (see 
petsc/src/tao/constrained/impls/ipm/pdipm.c) which assembles several small 
matrices into a large KKT matrix in mpiaij format. You could take the same 
approach to insert P and P^T into your K.
FYI, I attached our paper.
Hong

From: petsc-users  on behalf of Matthew 
Knepley 
Sent: Thursday, April 20, 2023 5:37 AM
To: Karthikeyan Chockalingam - STFC UKRI 
Cc: petsc-users@mcs.anl.gov 
Subject: Re: [petsc-users] question about MatSetLocalToGlobalMapping

On Thu, Apr 20, 2023 at 6:13 AM Karthikeyan Chockalingam - STFC UKRI via 
petsc-users mailto:petsc-users@mcs.anl.gov>> wrote:
Hello,

I created a new thread, thought would it be more appropriate (and is a 
continuation of my previous post). I want to construct the below K matrix 
(which is composed of submatrices)

K = [A P^T
   P   0]

Where K is of type MatMPIAIJ. I first constructed the top left [A] using 
MatSetValues().

Now, I would like to construct the bottom left [p] and top right [p^T] using 
MatSetValuesLocal().

To use  MatSetValuesLocal(),  I first have to create a local-to-global mapping 
using ISLocalToGlobalMappingCreate. I have created two mapping row_mapping and 
column_mapping.

I do not understand why they are not the same map. Maybe I was unclear before. 
It looks like you have two fields, say phi and lambda, where lambda is a 
Lagrange multiplier imposing some constraint. Then you get a saddle point like 
this. You can imagine matrices

  (phi, phi)--> A
  (phi, lambda) --> P^T
  (lambda, phi) --> P

So you make a L2G map for the phi field and the lambda field. Oh, you are 
calling them row and col map, but they are my phi and lambda
maps. I do not like the row and col names since in P they reverse.

Q1) At what point should I declare MatSetLocalToGlobalMapping – is it just 
before I use MatSetValuesLocal()?

Okay, it is good you are asking this because my thinking was somewhat confused. 
I think the precise steps are:

  1) Create the large saddle point matrix K

1a) We must call 
https://petsc.org/main/manualpages/Mat/MatSetLocalToGlobalMapping/ on it. In 
the simplest case, this just maps
   the local rows numbers [0, Nrows) to the global rows numbers 
[rowStart, rowStart + Nrows).

  2) To form each piece:

2a) Extract that block using 
https://petsc.org/main/manualpages/Mat/MatGetLocalSubMatrix/

   This gives back a Mat object that you subsequently restore using 
https://petsc.org/main/manualpages/Mat/MatRestoreLocalSubMatrix/

 2b) Insert values using 
https://petsc.org/main/manualpages/Mat/MatSetValuesLocal/

The local indices used for insertion here are indices relative to 
the block itself, and the L2G map for this matrix
has been rewritten to insert into that block in the larger matrix. 
Thus this looks like just calling MatSetValuesLocal()
on the smaller matrix block, but inserts correctly into the larger 
matrix.

Therefore, the code you write code in 2) could work equally well making the 
large matrix from 1), or independent smaller matrix blocks.

Does this make sense?

  Thanks,

 Matt

I will use MatSetLocalToGlobalMapping(K, row_mapping, column_mapping) to build 
the bottom left [P].


Q2) Can now I reset the mapping as MatSetLocalToGlobalMapping(K, 
column_mapping, row_mapping) to build the top right [P^T]?


Many thanks!

Kind regards,
Karthik.











--
What most experimenters take for granted before they begin their experiments is 
infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>

Re: [petsc-users] question about MatSetLocalToGlobalMapping

[Pierre Jolivet]

Hong,
1) Is there any hope to get PDIPDM to use a MatNest?
2) Is this fixed 
https://lists.mcs.anl.gov/pipermail/petsc-dev/2020-September/026398.html ?
I cannot get users to transition away from Ipopt because of these two missing 
features.

Thanks,
Pierre

> On 20 Apr 2023, at 5:47 PM, Zhang, Hong via petsc-users 
>  wrote:
> 
> Karthik,
> We built a KKT matrix in  TaoSNESJacobian_PDIPM() (see 
> petsc/src/tao/constrained/impls/ipm/pdipm.c) which assembles several small 
> matrices into a large KKT matrix in mpiaij format. You could take the same 
> approach to insert P and P^T into your K.
> FYI, I attached our paper.
> Hong
> From: petsc-users  on behalf of Matthew 
> Knepley 
> Sent: Thursday, April 20, 2023 5:37 AM
> To: Karthikeyan Chockalingam - STFC UKRI 
> Cc: petsc-users@mcs.anl.gov 
> Subject: Re: [petsc-users] question about MatSetLocalToGlobalMapping
>  
> On Thu, Apr 20, 2023 at 6:13 AM Karthikeyan Chockalingam - STFC UKRI via 
> petsc-users mailto:petsc-users@mcs.anl.gov>> wrote:
> Hello,
>  
> I created a new thread, thought would it be more appropriate (and is a 
> continuation of my previous post). I want to construct the below K matrix 
> (which is composed of submatrices)
>  
> K = [A P^T
>P   0]
>  
> Where K is of type MatMPIAIJ. I first constructed the top left [A] using 
> MatSetValues().
>  
> Now, I would like to construct the bottom left [p] and top right [p^T] using 
> MatSetValuesLocal().
>  
> To use  MatSetValuesLocal(),  I first have to create a local-to-global 
> mapping using ISLocalToGlobalMappingCreate. I have created two mapping 
> row_mapping and column_mapping.
> 
> I do not understand why they are not the same map. Maybe I was unclear 
> before. It looks like you have two fields, say phi and lambda, where lambda 
> is a Lagrange multiplier imposing some constraint. Then you get a saddle 
> point like this. You can imagine matrices
> 
>   (phi, phi)--> A
>   (phi, lambda) --> P^T
>   (lambda, phi) --> P
> 
> So you make a L2G map for the phi field and the lambda field. Oh, you are 
> calling them row and col map, but they are my phi and lambda
> maps. I do not like the row and col names since in P they reverse.
>  
> Q1) At what point should I declare MatSetLocalToGlobalMapping – is it just 
> before I use MatSetValuesLocal()?
> 
> Okay, it is good you are asking this because my thinking was somewhat 
> confused. I think the precise steps are:
> 
>   1) Create the large saddle point matrix K
> 
> 1a) We must call 
> https://petsc.org/main/manualpages/Mat/MatSetLocalToGlobalMapping/ on it. In 
> the simplest case, this just maps
>the local rows numbers [0, Nrows) to the global rows numbers 
> [rowStart, rowStart + Nrows).
> 
>   2) To form each piece:
> 
> 2a) Extract that block using 
> https://petsc.org/main/manualpages/Mat/MatGetLocalSubMatrix/
> 
>This gives back a Mat object that you subsequently restore using 
> https://petsc.org/main/manualpages/Mat/MatRestoreLocalSubMatrix/
> 
>  2b) Insert values using 
> https://petsc.org/main/manualpages/Mat/MatSetValuesLocal/
> 
> The local indices used for insertion here are indices relative to 
> the block itself, and the L2G map for this matrix
> has been rewritten to insert into that block in the larger 
> matrix. Thus this looks like just calling MatSetValuesLocal()
> on the smaller matrix block, but inserts correctly into the 
> larger matrix.
> 
> Therefore, the code you write code in 2) could work equally well making the 
> large matrix from 1), or independent smaller matrix blocks.
> 
> Does this make sense?
> 
>   Thanks,
> 
>  Matt
>  
> I will use MatSetLocalToGlobalMapping(K, row_mapping, column_mapping) to 
> build the bottom left [P].
>  
>  
> Q2) Can now I reset the mapping as MatSetLocalToGlobalMapping(K, 
> column_mapping, row_mapping) to build the top right [P^T]? 
>  
>  
> Many thanks!
>  
> Kind regards,
> Karthik.
>  
>  
>  
>  
>  
>  
>  
>  
>  
> 
> 
> -- 
> What most experimenters take for granted before they begin their experiments 
> is infinitely more interesting than any results to which their experiments 
> lead.
> -- Norbert Wiener
> 
> https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
>  interior point method for the solution of dynamic.pdf>

Re: [petsc-users] question about MatSetLocalToGlobalMapping

[Matthew Knepley]

On Thu, Apr 20, 2023 at 6:13 AM Karthikeyan Chockalingam - STFC UKRI via
petsc-users  wrote:

> Hello,
>
>
>
> I created a new thread, thought would it be more appropriate (and is a
> continuation of my previous post). I want to construct the below K matrix
> (which is composed of submatrices)
>
>
>
> K = [A P^T
>
>P   0]
>
>
>
> Where K is of type MatMPIAIJ. I first constructed the top left [A] using
> MatSetValues().
>
>
>
> Now, I would like to construct the bottom left [p] and top right [p^T] using
> MatSetValuesLocal().
>
>
>
> To use  MatSetValuesLocal(),  I first have to create a local-to-global
> mapping using ISLocalToGlobalMappingCreate. I have created two mapping
> row_mapping and column_mapping.
>

I do not understand why they are not the same map. Maybe I was unclear
before. It looks like you have two fields, say phi and lambda, where lambda
is a Lagrange multiplier imposing some constraint. Then you get a saddle
point like this. You can imagine matrices

  (phi, phi)--> A
  (phi, lambda) --> P^T
  (lambda, phi) --> P

So you make a L2G map for the phi field and the lambda field. Oh, you are
calling them row and col map, but they are my phi and lambda
maps. I do not like the row and col names since in P they reverse.


> Q1) At what point should I declare MatSetLocalToGlobalMapping – is it
> just before I use MatSetValuesLocal()?
>

Okay, it is good you are asking this because my thinking was somewhat
confused. I think the precise steps are:

  1) Create the large saddle point matrix K

1a) We must call
https://petsc.org/main/manualpages/Mat/MatSetLocalToGlobalMapping/ on it.
In the simplest case, this just maps
   the local rows numbers [0, Nrows) to the global rows numbers
[rowStart, rowStart + Nrows).

  2) To form each piece:

2a) Extract that block using
https://petsc.org/main/manualpages/Mat/MatGetLocalSubMatrix/

   This gives back a Mat object that you subsequently restore using
https://petsc.org/main/manualpages/Mat/MatRestoreLocalSubMatrix/

 2b) Insert values using
https://petsc.org/main/manualpages/Mat/MatSetValuesLocal/

The local indices used for insertion here are indices relative
to the block itself, and the L2G map for this matrix
has been rewritten to insert into that block in the larger
matrix. Thus this looks like just calling MatSetValuesLocal()
on the smaller matrix block, but inserts correctly into the
larger matrix.

Therefore, the code you write code in 2) could work equally well making the
large matrix from 1), or independent smaller matrix blocks.

Does this make sense?

  Thanks,

 Matt


> I will use MatSetLocalToGlobalMapping(K, row_mapping, column_mapping) to
> build the bottom left [P].
>
>
>
>
>
> Q2) Can now I reset the mapping as MatSetLocalToGlobalMapping(K,
> column_mapping, row_mapping) to build the top right [P^T]?
>
>
>
>
>
> Many thanks!
>
>
>
> Kind regards,
>
> Karthik.
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about -memory_view

[Barry Smith]

> On Apr 12, 2023, at 4:08 PM, Jorti, Zakariae via petsc-users 
>  wrote:
> 
> Hello,
> 
> I am running some matrix computations on 64 nodes, using 640 MPI tasks. 
> And I wanted to check the memory usage with the -memory_view flag.
> I get the following output: 
> 
> Summary of Memory Usage in PETSc
> Maximum (over computational time) process memory:total 3.1056e+11 max 
> 5.9918e+08 min 4.2213e+08
> Current process memory:  total 1.9194e+11 max 
> 3.9960e+08 min 2.2761e+08
> 
> 
> What is the difference between maximum process memory and current process 
> memory? 

   Over computational time means the maximum it ever was (the high water mark) 
and current means what it is right now. For memory usage obtained from the OS 
(what we call  "process memory" in the output as opposed to PetscMalloc()ed 
memory)  often the current does not ever go below the maximum it ever was 
because the "extra now unneeded memory" is not returned to the OS. 


> What does total mean here? (in each node or total of all the nodes)

 Total is sum over all MPI ranks, max is maximum over all ranks, min is minimum 
over all ranks

> Also, if the job fails because of memory shortage, will this -memory_view 
> still output some information?

   Generally not if the job fails before the memory information is printed. 
Usually one runs with smaller memory usage increasing the problem size several 
times to see how the memory usage scales with the problem size (linearly, 
quadratically, etc) and this guides understanding the memory usage and if it 
can be improved. Just running with a large memory usage alone is not that 
useful in providing information.


> 
> Thank you.
> 
> Zakariae

Re: [petsc-users] Question about NASM initialization

[Takahashi, Tadanaga]

Ok, thanks for the clarification.

On Thu, Apr 6, 2023 at 10:25 AM Matthew Knepley  wrote:

> On Thu, Apr 6, 2023 at 10:21 AM Takahashi, Tadanaga  wrote:
>
>> I am following up from the last inquiry. I read the source code nasm.c
>>  and it looks
>> like sub-snes iteration is being initialized with a scatter call from the
>> previous solution. In other words, if I use Newton's method for the local
>> solver, then in each NASM iteration the Newton's method uses the previous
>> local solution as the initial guess. Can anyone confirm this?
>>
>
> This is the intention. There are not many tests, so it is possible there
> is a bug, but it is supposed to use the existing solution.
>
>   Thanks,
>
> Matt
>
>
>> On Sun, Apr 2, 2023 at 6:14 PM Takahashi, Tadanaga  wrote:
>>
>>> Hello PETSc devs,
>>>
>>> I am using SNES NASM with Newton LS on the sub-SNES. I was wondering how
>>> the sub-SNES chooses the initial guess during each NASM iteration. Is it
>>> using the previously computed solution or is it restarting from zero?
>>>
>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>

Re: [petsc-users] Question about NASM initialization

[Matthew Knepley]

On Thu, Apr 6, 2023 at 10:21 AM Takahashi, Tadanaga  wrote:

> I am following up from the last inquiry. I read the source code nasm.c
>  and it looks
> like sub-snes iteration is being initialized with a scatter call from the
> previous solution. In other words, if I use Newton's method for the local
> solver, then in each NASM iteration the Newton's method uses the previous
> local solution as the initial guess. Can anyone confirm this?
>

This is the intention. There are not many tests, so it is possible there is
a bug, but it is supposed to use the existing solution.

  Thanks,

Matt


> On Sun, Apr 2, 2023 at 6:14 PM Takahashi, Tadanaga  wrote:
>
>> Hello PETSc devs,
>>
>> I am using SNES NASM with Newton LS on the sub-SNES. I was wondering how
>> the sub-SNES chooses the initial guess during each NASM iteration. Is it
>> using the previously computed solution or is it restarting from zero?
>>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about NASM initialization

[Takahashi, Tadanaga]

I am following up from the last inquiry. I read the source code nasm.c
 and it looks
like sub-snes iteration is being initialized with a scatter call from the
previous solution. In other words, if I use Newton's method for the local
solver, then in each NASM iteration the Newton's method uses the previous
local solution as the initial guess. Can anyone confirm this?

On Sun, Apr 2, 2023 at 6:14 PM Takahashi, Tadanaga  wrote:

> Hello PETSc devs,
>
> I am using SNES NASM with Newton LS on the sub-SNES. I was wondering how
> the sub-SNES chooses the initial guess during each NASM iteration. Is it
> using the previously computed solution or is it restarting from zero?
>

Re: [petsc-users] Question about MatView

[Barry Smith]

src/mat/tests/ex9.c and other examples in that directory use it.



> On Mar 17, 2023, at 11:05 AM, user_gong Kim  wrote:
> 
> PETSC_COMM_SELF generates an error in more than 2 processes.

It should not


> 
> I would like to use the other method you said, matcreateredundantmatrix. 
> However, there is no example in the manual. Can you give an example using 
> this function?
> 
> 
> 
> 2023년 3월 17일 (금) 오후 11:53, Barry Smith  >님이 작성:
>> 
>>Use 
>> 
 MatCreate(PETSC_COMM_SELF, );
  PetscViewerBinaryOpen(PETSC_COMM_SELF, "mat.bin", FILE_MODE_READ, 
 );
>> 
>> 
>> If it one program running that both views and loads the matrix you can 
>> use MatCreateRedundantMatrix() to reproduce the entire matrix on each MPI 
>> rank. It is better than using the filesystem to do it.
>> 
>> 
>>> On Mar 17, 2023, at 9:45 AM, user_gong Kim >> > wrote:
>>> 
>>> Following your comments,  I did an test.
>>> However, if I run the application in parallel. 
>>> In all processes, it is not possible to obtain values at all positions in 
>>> the matrix through MatGetValue.
>>> As in the previous case of saving in binary, it is read in parallel divided 
>>> form.
>>> Is it impossible to want to get the all value in the whole process?
>>> 
>>> 
>>> Thanks,
>>> Hyung Kim
>>> 
>>> 2023년 3월 17일 (금) 오후 7:35, Matthew Knepley >> >님이 작성:
 On Fri, Mar 17, 2023 at 5:51 AM user_gong Kim >>> > wrote:
> Hello,
> 
>  
> I have 2 questions about MatView.
> 
>  
> 1.I would like to ask if the process below is possible.
> When running in parallel, is it possible to make the matrix of the mpiaij 
> format into a txt file, output it, and read it again so that the entire 
> process has the same matrix?
> 
 No. However, you can do this with a binary viewer. I suggest using
 
   MatViewFromOptions(mat, NULL, "-my_view");
 
 and then the command line argument
 
   -my_view binary:mat.bin
 
 and then you can read this in using
 
   MatCreate(PETSC_COMM_WORLD, );
   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mat.bin", FILE_MODE_READ, 
 );
   MatLoad(mat, viewer);
   ViewerDestroy();
 
   THanks,
 
  Matt
 
  
> 2.If possible, please let me know which function can be used to 
> create a txt file and how to read the txt file.
> 
>  
> Thanks,
> 
> Hyung Kim
> 
 
 
 -- 
 What most experimenters take for granted before they begin their 
 experiments is infinitely more interesting than any results to which their 
 experiments lead.
 -- Norbert Wiener
 
 https://www.cse.buffalo.edu/~knepley/ 
 
>> 

Re: [petsc-users] Question about MatView

[user_gong Kim]

PETSC_COMM_SELF generates an error in more than 2 processes.

I would like to use the other method you said, matcreateredundantmatrix.
However, there is no example in the manual. Can you give an example using
this function?



2023년 3월 17일 (금) 오후 11:53, Barry Smith 님이 작성:

>
>Use
>
> MatCreate(PETSC_COMM_SELF, );
>>  PetscViewerBinaryOpen(PETSC_COMM_SELF, "mat.bin", FILE_MODE_READ,
>> );
>>
>
>
> If it one program running that both views and loads the matrix you can
> use MatCreateRedundantMatrix() to reproduce the entire matrix on each MPI
> rank. It is better than using the filesystem to do it.
>
>
> On Mar 17, 2023, at 9:45 AM, user_gong Kim  wrote:
>
> Following your comments,  I did an test.
> However, if I run the application in parallel.
> In all processes, it is not possible to obtain values at all positions in
> the matrix through MatGetValue.
> As in the previous case of saving in binary, it is read in parallel
> divided form.
> Is it impossible to want to get the all value in the whole process?
>
>
> Thanks,
> Hyung Kim
>
> 2023년 3월 17일 (금) 오후 7:35, Matthew Knepley 님이 작성:
>
>> On Fri, Mar 17, 2023 at 5:51 AM user_gong Kim  wrote:
>>
>>> Hello,
>>>
>>>
>>> I have 2 questions about MatView.
>>>
>>>
>>> 1.I would like to ask if the process below is possible.
>>> When running in parallel, is it possible to make the matrix of the
>>> mpiaij format into a txt file, output it, and read it again so that the
>>> entire process has the same matrix?
>>>
>> No. However, you can do this with a binary viewer. I suggest using
>>
>>   MatViewFromOptions(mat, NULL, "-my_view");
>>
>> and then the command line argument
>>
>>   -my_view binary:mat.bin
>>
>> and then you can read this in using
>>
>>   MatCreate(PETSC_COMM_WORLD, );
>>   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mat.bin", FILE_MODE_READ,
>> );
>>   MatLoad(mat, viewer);
>>   ViewerDestroy();
>>
>>   THanks,
>>
>>  Matt
>>
>>
>>
>>> 2.If possible, please let me know which function can be used to
>>> create a txt file and how to read the txt file.
>>>
>>>
>>> Thanks,
>>>
>>> Hyung Kim
>>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> 
>>
>
>

Re: [petsc-users] Question about MatView

[Barry Smith]

   Use 

>> MatCreate(PETSC_COMM_SELF, );
>>  PetscViewerBinaryOpen(PETSC_COMM_SELF, "mat.bin", FILE_MODE_READ, );


If it one program running that both views and loads the matrix you can use 
MatCreateRedundantMatrix() to reproduce the entire matrix on each MPI rank. It 
is better than using the filesystem to do it.


> On Mar 17, 2023, at 9:45 AM, user_gong Kim  wrote:
> 
> Following your comments,  I did an test.
> However, if I run the application in parallel. 
> In all processes, it is not possible to obtain values at all positions in the 
> matrix through MatGetValue.
> As in the previous case of saving in binary, it is read in parallel divided 
> form.
> Is it impossible to want to get the all value in the whole process?
> 
> 
> Thanks,
> Hyung Kim
> 
> 2023년 3월 17일 (금) 오후 7:35, Matthew Knepley  >님이 작성:
>> On Fri, Mar 17, 2023 at 5:51 AM user_gong Kim > > wrote:
>>> Hello,
>>> 
>>>  
>>> I have 2 questions about MatView.
>>> 
>>>  
>>> 1.I would like to ask if the process below is possible.
>>> When running in parallel, is it possible to make the matrix of the mpiaij 
>>> format into a txt file, output it, and read it again so that the entire 
>>> process has the same matrix?
>>> 
>> No. However, you can do this with a binary viewer. I suggest using
>> 
>>   MatViewFromOptions(mat, NULL, "-my_view");
>> 
>> and then the command line argument
>> 
>>   -my_view binary:mat.bin
>> 
>> and then you can read this in using
>> 
>>   MatCreate(PETSC_COMM_WORLD, );
>>   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mat.bin", FILE_MODE_READ, 
>> );
>>   MatLoad(mat, viewer);
>>   ViewerDestroy();
>> 
>>   THanks,
>> 
>>  Matt
>> 
>>  
>>> 2.If possible, please let me know which function can be used to create 
>>> a txt file and how to read the txt file.
>>> 
>>>  
>>> Thanks,
>>> 
>>> Hyung Kim
>>> 
>> 
>> 
>> -- 
>> What most experimenters take for granted before they begin their experiments 
>> is infinitely more interesting than any results to which their experiments 
>> lead.
>> -- Norbert Wiener
>> 
>> https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about MatView

[user_gong Kim]

Following your comments,  I did an test.
However, if I run the application in parallel.
In all processes, it is not possible to obtain values at all positions in
the matrix through MatGetValue.
As in the previous case of saving in binary, it is read in parallel divided
form.
Is it impossible to want to get the all value in the whole process?


Thanks,
Hyung Kim

2023년 3월 17일 (금) 오후 7:35, Matthew Knepley 님이 작성:

> On Fri, Mar 17, 2023 at 5:51 AM user_gong Kim  wrote:
>
>> Hello,
>>
>>
>>
>> I have 2 questions about MatView.
>>
>>
>>
>> 1.I would like to ask if the process below is possible.
>> When running in parallel, is it possible to make the matrix of the mpiaij
>> format into a txt file, output it, and read it again so that the entire
>> process has the same matrix?
>>
> No. However, you can do this with a binary viewer. I suggest using
>
>   MatViewFromOptions(mat, NULL, "-my_view");
>
> and then the command line argument
>
>   -my_view binary:mat.bin
>
> and then you can read this in using
>
>   MatCreate(PETSC_COMM_WORLD, );
>   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mat.bin", FILE_MODE_READ,
> );
>   MatLoad(mat, viewer);
>   ViewerDestroy();
>
>   THanks,
>
>  Matt
>
>
>
>> 2.If possible, please let me know which function can be used to
>> create a txt file and how to read the txt file.
>>
>>
>>
>> Thanks,
>>
>> Hyung Kim
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>

Re: [petsc-users] Question about MatView

[Matthew Knepley]

On Fri, Mar 17, 2023 at 5:51 AM user_gong Kim  wrote:

> Hello,
>
>
>
> I have 2 questions about MatView.
>
>
>
> 1.I would like to ask if the process below is possible.
> When running in parallel, is it possible to make the matrix of the mpiaij
> format into a txt file, output it, and read it again so that the entire
> process has the same matrix?
>
No. However, you can do this with a binary viewer. I suggest using

  MatViewFromOptions(mat, NULL, "-my_view");

and then the command line argument

  -my_view binary:mat.bin

and then you can read this in using

  MatCreate(PETSC_COMM_WORLD, );
  PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mat.bin", FILE_MODE_READ,
);
  MatLoad(mat, viewer);
  ViewerDestroy();

  THanks,

 Matt



> 2.If possible, please let me know which function can be used to
> create a txt file and how to read the txt file.
>
>
>
> Thanks,
>
> Hyung Kim
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about time issues in parallel computing

[Mark Adams]

Speed up to 4 processors should have at least 40,000 equations for 3D
problems and more for 2D.
At least for iterative solvers. This is probably a good place to start with
direct solvers but you might see benefit with a little less.

Mark

On Wed, Mar 15, 2023 at 9:09 PM ­권승리 / 학생 / 항공우주공학과 
wrote:

> Thank you for your reply.
>
> It was a simple problem, but it has more than 1000 degrees of freedom.
>
> Is this not enough to check speedup?
>
> Best regards
> Seung Lee Kwon
>
> 2023년 3월 15일 (수) 오후 8:50, Matthew Knepley 님이 작성:
>
>> On Wed, Mar 15, 2023 at 3:38 AM ­권승리 / 학생 / 항공우주공학과 
>> wrote:
>>
>>> Dear petsc developers.
>>>
>>> Hello.
>>> I am trying to solve the structural problem with FEM and test parallel
>>> computing works well.
>>>
>>> However, even if I change the number of cores, the total time is
>>> calculated the same.
>>>
>>> I have tested on a simple problem using a MUMPS solver using:
>>> mpiexec -n 1
>>> mpiexec -n 2
>>> mpiexec -n 4
>>> ...
>>>
>>> Could you give me some advice if you have experienced this problem?
>>>
>>
>> If your problem is small, you could very well see no speedup:
>>
>>
>> https://petsc.org/main/faq/#what-kind-of-parallel-computers-or-clusters-are-needed-to-use-petsc-or-why-do-i-get-little-speedup
>>
>>   Thanks,
>>
>>  Matt
>>
>>
>>> Best regards
>>> Seung Lee Kwon
>>> --
>>> Seung Lee Kwon, Ph.D.Candidate
>>> Aerospace Structures and Materials Laboratory
>>> Department of Mechanical and Aerospace Engineering
>>> Seoul National University
>>> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
>>> E-mail : ksl7...@snu.ac.kr
>>> Office : +82-2-880-7389
>>> C. P : +82-10-4695-1062
>>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> 
>>
>
>
> --
> Seung Lee Kwon, Ph.D.Candidate
> Aerospace Structures and Materials Laboratory
> Department of Mechanical and Aerospace Engineering
> Seoul National University
> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
> E-mail : ksl7...@snu.ac.kr
> Office : +82-2-880-7389
> C. P : +82-10-4695-1062
>

Re: [petsc-users] Question about time issues in parallel computing

[­권승리 / 학생 / 항공우주공학과]

Thank you for your reply.

It was a simple problem, but it has more than 1000 degrees of freedom.

Is this not enough to check speedup?

Best regards
Seung Lee Kwon

2023년 3월 15일 (수) 오후 8:50, Matthew Knepley 님이 작성:

> On Wed, Mar 15, 2023 at 3:38 AM ­권승리 / 학생 / 항공우주공학과 
> wrote:
>
>> Dear petsc developers.
>>
>> Hello.
>> I am trying to solve the structural problem with FEM and test parallel
>> computing works well.
>>
>> However, even if I change the number of cores, the total time is
>> calculated the same.
>>
>> I have tested on a simple problem using a MUMPS solver using:
>> mpiexec -n 1
>> mpiexec -n 2
>> mpiexec -n 4
>> ...
>>
>> Could you give me some advice if you have experienced this problem?
>>
>
> If your problem is small, you could very well see no speedup:
>
>
> https://petsc.org/main/faq/#what-kind-of-parallel-computers-or-clusters-are-needed-to-use-petsc-or-why-do-i-get-little-speedup
>
>   Thanks,
>
>  Matt
>
>
>> Best regards
>> Seung Lee Kwon
>> --
>> Seung Lee Kwon, Ph.D.Candidate
>> Aerospace Structures and Materials Laboratory
>> Department of Mechanical and Aerospace Engineering
>> Seoul National University
>> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
>> E-mail : ksl7...@snu.ac.kr
>> Office : +82-2-880-7389
>> C. P : +82-10-4695-1062
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>


-- 
Seung Lee Kwon, Ph.D.Candidate
Aerospace Structures and Materials Laboratory
Department of Mechanical and Aerospace Engineering
Seoul National University
Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
E-mail : ksl7...@snu.ac.kr
Office : +82-2-880-7389
C. P : +82-10-4695-1062

Re: [petsc-users] Question about time issues in parallel computing

[Matthew Knepley]

On Wed, Mar 15, 2023 at 3:38 AM ­권승리 / 학생 / 항공우주공학과 
wrote:

> Dear petsc developers.
>
> Hello.
> I am trying to solve the structural problem with FEM and test parallel
> computing works well.
>
> However, even if I change the number of cores, the total time is
> calculated the same.
>
> I have tested on a simple problem using a MUMPS solver using:
> mpiexec -n 1
> mpiexec -n 2
> mpiexec -n 4
> ...
>
> Could you give me some advice if you have experienced this problem?
>

If your problem is small, you could very well see no speedup:


https://petsc.org/main/faq/#what-kind-of-parallel-computers-or-clusters-are-needed-to-use-petsc-or-why-do-i-get-little-speedup

  Thanks,

 Matt


> Best regards
> Seung Lee Kwon
> --
> Seung Lee Kwon, Ph.D.Candidate
> Aerospace Structures and Materials Laboratory
> Department of Mechanical and Aerospace Engineering
> Seoul National University
> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
> E-mail : ksl7...@snu.ac.kr
> Office : +82-2-880-7389
> C. P : +82-10-4695-1062
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about rank of matrix

[Barry Smith]

> On Feb 17, 2023, at 2:43 AM, user_gong Kim  wrote:
> 
> Hello,
> 
> I have a question about rank of matrix.
> At the problem 
> Au = b, 
> 
> In my case, sometimes global matrix A is not full rank.
> In this case, the global matrix A is more likely to be singular, and if it 
> becomes singular, the problem cannot be solved even in the case of the direct 
> solver.
> I haven't solved the problem with an iterative solver yet, but I would like 
> to ask someone who has experienced this kind of problem.
> 
> 1. If it is not full rank, is there a numerical technique to solve it by 
> catching rows and columns with empty ranks in advance?

   Matrices with completely empty (zero) rows and their corresponding columns 
are a simple type of singularity (so long as the corresponding right hand side 
values are zero) that are not difficult to solve.

   Essentially iterative solvers never "see" the empty rows and columns and 
successfully solve the linear systems. 

   You can also use PCREDISTRIBUTE (even on one MPI rank) in the main git 
branch of PETSc to automatically extract the nontrivial rows and columns from 
the matrix and solve that smaller system with any solver including direct 
solvers.



> 
> 2.If anyone has solved it in a different way than the above numerical 
> analysis method, please tell me your experience.
> 
> Thanks,
> Hyung Kim
> 
> 

Re: [petsc-users] Question about rank of matrix

[Matthew Knepley]

On Fri, Feb 17, 2023 at 2:43 AM user_gong Kim  wrote:

> Hello,
>
> I have a question about rank of matrix.
> At the problem
> Au = b,
>
> In my case, sometimes global matrix A is not full rank.
> In this case, the global matrix A is more likely to be singular, and if it
> becomes singular, the problem cannot be solved even in the case of the
> direct solver.
> I haven't solved the problem with an iterative solver yet, but I would
> like to ask someone who has experienced this kind of problem.
>
> 1. If it is not full rank, is there a numerical technique to solve it by
> catching rows and columns with empty ranks in advance?
>
> 2.If anyone has solved it in a different way than the above numerical
> analysis method, please tell me your experience.
>

As Pierre points out,  MUMPS can solve singular systems.

If you have an explicit characterization of the null space, then many
iterative methods can also solve it by projecting out
the null space. You call MatSetNullSpace() on the system matrix.

  Thanks,

 Matt


> Thanks,
> Hyung Kim
>
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about rank of matrix

[Pierre Jolivet]

> On 17 Feb 2023, at 8:56 AM, Stefano Zampini  wrote:
> 
> 
> On Fri, Feb 17, 2023, 10:43 user_gong Kim  > wrote:
>> Hello,
>> 
>> I have a question about rank of matrix.
>> At the problem 
>> Au = b, 
>> 
>> In my case, sometimes global matrix A is not full rank.
>> In this case, the global matrix A is more likely to be singular, and if it 
>> becomes singular, the problem cannot be solved even in the case of the 
>> direct solver.
>> I haven't solved the problem with an iterative solver yet, but I would like 
>> to ask someone who has experienced this kind of problem.
>> 
>> 1. If it is not full rank, is there a numerical technique to solve it by 
>> catching rows and columns with empty ranks in advance?
>> 
>> 2.If anyone has solved it in a different way than the above numerical 
>> analysis method, please tell me your experience.
>> 
>> Thanks,
>> Hyung Kim
> 
> 
> My experience with this is usually associated to reading a book and find the 
> solution I'm looking for. 

On top of that, some exact factorization packages can solve singular systems, 
unlike what you are stating.
E.g., MUMPS, together with the option -mat_mumps_icntl_24, see 
https://mumps-solver.org/doc/userguide_5.5.1.pdf

Thanks,
Pierre

Re: [petsc-users] Question about rank of matrix

[Stefano Zampini]

On Fri, Feb 17, 2023, 10:43 user_gong Kim  wrote:

> Hello,
>
> I have a question about rank of matrix.
> At the problem
> Au = b,
>
> In my case, sometimes global matrix A is not full rank.
> In this case, the global matrix A is more likely to be singular, and if it
> becomes singular, the problem cannot be solved even in the case of the
> direct solver.
> I haven't solved the problem with an iterative solver yet, but I would
> like to ask someone who has experienced this kind of problem.
>
> 1. If it is not full rank, is there a numerical technique to solve it by
> catching rows and columns with empty ranks in advance?
>
> 2.If anyone has solved it in a different way than the above numerical
> analysis method, please tell me your experience.
>
> Thanks,
> Hyung Kim
>

My experience with this is usually associated to reading a book and find
the solution I'm looking for.

>
>
>

Re: [petsc-users] Question for Petsc

[Stefano Zampini]

For bddc, you can also take a look at
https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tutorials/ex71.c

On Thu, Feb 16, 2023, 19:41 Matthew Knepley  wrote:

> On Thu, Feb 16, 2023 at 9:14 AM ziming xiong 
> wrote:
>
>> Hello，
>> I want to use Petsc to implement high performance computing, and I mainly
>> want to apply DDM methods to parallel computing. I have implemented some of
>> the DDM methods (such as ASM, Bjacobi, etc.), but I don't understand the
>> PCBDDC method. The official example (src/ksp/ksp/tutorials/ex59.c.html) is
>> too complicated, so I have not been able to figure out the setting process.
>> I would like to ask you if you have other simple and clearer examples for
>> reference.
>>
>
> You could look at the paper:
> https://epubs.siam.org/doi/abs/10.1137/15M1025785
>
>
>> Secondly, I would like to apply mklPardiso to Petsc. But not work, can u
>> help me figure out the problem? i use oneAPI for the mklpardiso, and when i
>> configure, i give the blaslapack lib.
>>
>
> You should reconfigure with --download-mkl_pardiso
>
>   Thanks,
>
>  Matt
>
>
>> there are the errors:
>>
>> [0]PETSC ERROR: See
>> https://petsc.org/release/overview/linear_solve_table/ for possible LU
>> and Cholesky solvers
>> [0]PETSC ERROR: Could not locate solver type mkl_pardiso for
>> factorization type LU and matrix type seqaij. Perhaps you must ./configure
>> with --download-mkl_pardiso
>> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
>> [0]PETSC ERROR: Petsc Release Version 3.18.3, Dec 28, 2022
>> [0]PETSC ERROR:
>> C:\Users\XiongZiming\Desktop\test_petsc_FEM\test_petsc_fem\x64\Release\test_petsc_fem.exe
>> on a arch-mswin-c-debug named lmeep-329 by XiongZiming Thu Feb 16 15:05:14
>> 2023
>> [0]PETSC ERROR: Configure options --with-cc="win32fe cl" --with-fc=0
>> --with-cxx="win32fe cl" --with-shared-libraries=0
>> --with-mpi-include="[/cygdrive/c/PROGRA~2/Intel/MPI/Include,/cygdrive/c/PROGRA~2/Intel/MPI/Include/x64]"
>> --with-mpi-lib="-L/cygdrive/c/PROGRA~2/Intel/MPI/lib/x64 msmpifec.lib
>> msmpi.lib" --with-mpiexec=/cygdrive/c/PROGRA~1/Microsoft_MPI/Bin/mpiexec
>> --with-blaslapack-lib="-L/cygdrive/c/PROGRA~2/Intel/oneAPI/mkl/2023.0.0/lib/intel64
>> mkl_intel_lp64_dll.lib mkl_sequential_dll.lib mkl_core_dll.lib"
>>
>> Thanks,
>> Ziming XIONG
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>

Re: [petsc-users] Question for Petsc

[Matthew Knepley]

On Thu, Feb 16, 2023 at 9:14 AM ziming xiong 
wrote:

> Hello，
> I want to use Petsc to implement high performance computing, and I mainly
> want to apply DDM methods to parallel computing. I have implemented some of
> the DDM methods (such as ASM, Bjacobi, etc.), but I don't understand the
> PCBDDC method. The official example (src/ksp/ksp/tutorials/ex59.c.html) is
> too complicated, so I have not been able to figure out the setting process.
> I would like to ask you if you have other simple and clearer examples for
> reference.
>

You could look at the paper:
https://epubs.siam.org/doi/abs/10.1137/15M1025785


> Secondly, I would like to apply mklPardiso to Petsc. But not work, can u
> help me figure out the problem? i use oneAPI for the mklpardiso, and when i
> configure, i give the blaslapack lib.
>

You should reconfigure with --download-mkl_pardiso

  Thanks,

 Matt


> there are the errors:
>
> [0]PETSC ERROR: See https://petsc.org/release/overview/linear_solve_table/
> for possible LU and Cholesky solvers
> [0]PETSC ERROR: Could not locate solver type mkl_pardiso for factorization
> type LU and matrix type seqaij. Perhaps you must ./configure with
> --download-mkl_pardiso
> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.18.3, Dec 28, 2022
> [0]PETSC ERROR:
> C:\Users\XiongZiming\Desktop\test_petsc_FEM\test_petsc_fem\x64\Release\test_petsc_fem.exe
> on a arch-mswin-c-debug named lmeep-329 by XiongZiming Thu Feb 16 15:05:14
> 2023
> [0]PETSC ERROR: Configure options --with-cc="win32fe cl" --with-fc=0
> --with-cxx="win32fe cl" --with-shared-libraries=0
> --with-mpi-include="[/cygdrive/c/PROGRA~2/Intel/MPI/Include,/cygdrive/c/PROGRA~2/Intel/MPI/Include/x64]"
> --with-mpi-lib="-L/cygdrive/c/PROGRA~2/Intel/MPI/lib/x64 msmpifec.lib
> msmpi.lib" --with-mpiexec=/cygdrive/c/PROGRA~1/Microsoft_MPI/Bin/mpiexec
> --with-blaslapack-lib="-L/cygdrive/c/PROGRA~2/Intel/oneAPI/mkl/2023.0.0/lib/intel64
> mkl_intel_lp64_dll.lib mkl_sequential_dll.lib mkl_core_dll.lib"
>
> Thanks,
> Ziming XIONG
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about preconditioner

[Barry Smith]

   If your  matrix has the form 

(   A   B )
(   C   0 )

then often PCFIELDSPLIT can be a useful preconditioner with the option 
-pc_fieldsplit_detect_saddle_point



> On Feb 16, 2023, at 2:42 AM, user_gong Kim  wrote:
> 
>  
> Hello,
> 
>  
> There are some questions about some preconditioners.
> 
> The questions are from problem Au=b. The global matrix A has zero value 
> diagonal terms.
> 
> 1. Which preconditioner is preferred for matrix A which has zero value in 
> diagonal terms?
> The most frequently used basic 2 preconditioners are jacobi and SOR (gauss 
> seidel). As people knows both methods should have non zero diagonal terms. 
> Although the improved method is applied in PETSc, jacobi can also solve the 
> case with zero diagonal term, but I ask because I know that it is not 
> recommended.
> 
> 2. Second question is about running code with the two command options 
> below in a single process.
> 1st command : -ksp_type gmres -pc_type bjacobi -sub_pc_type jacobi
> 2nd command : -ksp_type gmres -pc_type hpddm -sub_pc_type jacobi
> When domain decomposition methods such as bjacobi or hpddm are parallel, the 
> global matrix is divided for each process. As far as I know, running it in a 
> single process should eventually produce the same result if the sub pc type 
> is the same. However, in the second option, ksp did not converge.
> In this case, I wonder how to analyze the situation.
> How can I monitor and see the difference between the two?
> 
>  
>  
> Thanks,
> 
> Hyung Kim
> 

Re: [petsc-users] Question about preconditioner

[Matthew Knepley]

On Thu, Feb 16, 2023 at 2:43 AM user_gong Kim  wrote:

>
>
> Hello,
>
>
>
> There are some questions about some preconditioners.
>
> The questions are from problem Au=b. The global matrix A has zero value
> diagonal terms.
>
> 1. Which preconditioner is preferred for matrix A which has zero
> value in diagonal terms?
> The most frequently used basic 2 preconditioners are jacobi and SOR (gauss
> seidel). As people knows both methods should have non zero diagonal terms.
> Although the improved method is applied in PETSc, jacobi can also solve the
> case with zero diagonal term, but I ask because I know that it is not
> recommended.
>
> 2. Second question is about running code with the two command options
> below in a single process.
> 1st command : -ksp_type gmres -pc_type bjacobi -sub_pc_type jacobi
> 2nd command : -ksp_type gmres -pc_type hpddm -sub_pc_type jacobi
> When domain decomposition methods such as bjacobi or hpddm are parallel,
> the global matrix is divided for each process. As far as I know, running it
> in a single process should eventually produce the same result if the sub pc
> type is the same. However, in the second option, ksp did not converge.
> In this case, I wonder how to analyze the situation.
> How can I monitor and see the difference between the two?
>
> Pierre is correct. I will just note that the best way to use PETSc is
generally to find the preconditioner you need
by looking in literature for what has been successful for other people, and
then reproducing it in PETSc, which
should be easy.

  Thanks,

 Matt


>
>
> Thanks,
>
> Hyung Kim
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about preconditioner

[Pierre Jolivet]

> On 16 Feb 2023, at 8:43 AM, user_gong Kim  wrote:
> 
> 
>  
> 
> Hello,
> 
>  
> 
> There are some questions about some preconditioners.
> 
> The questions are from problem Au=b. The global matrix A has zero value 
> diagonal terms.
> 
> 1. Which preconditioner is preferred for matrix A which has zero value in 
> diagonal terms?
> 

This question has not a single answer. It all depends on where your A and b are 
coming from.

> The most frequently used basic 2 preconditioners are jacobi and SOR (gauss 
> seidel).
> 
They are not the most frequently used. And rightfully so, as they very often 
can’t handle non-trivial systems.

> As people knows both methods should have non zero diagonal terms. Although 
> the improved method is applied in PETSc, jacobi can also solve the case with 
> zero diagonal term, but I ask because I know that it is not recommended.
> 
> 2. Second question is about running code with the two command options 
> below in a single process.
> 1st command : -ksp_type gmres -pc_type bjacobi -sub_pc_type jacobi
> 2nd command : -ksp_type gmres -pc_type hpddm -sub_pc_type jacobi
> When domain decomposition methods such as bjacobi or hpddm are parallel, the 
> global matrix is divided for each process. As far as I know, running it in a 
> single process should eventually produce the same result if the sub pc type 
> is the same. However, in the second option, ksp did not converge.
> 
1st command: it’s pointless to couple PCBJACOBI with PCJABOCI, it’s equivalent 
to only using PCJABOBI.
2nd command: it’s pointless to use PCHPDDM if you don’t specify in some way how 
to coarsen your problem (either algebraically or via an auxiliary operator). 
You just have a single level (equivalent to PCBJACOBI), but its options are 
prefixed by -pc_hpddm_coarse_ instead of -sub_
Again, both sets of options do not make sense.
If you want, you could share your A and b (or tell us what you are 
discretizing) and we will be able to provide a better feedback.

Thanks,
Pierre

> In this case, I wonder how to analyze the situation.
> How can I monitor and see the difference between the two?
> 
>  
> 
>  
> 
> Thanks,
> 
> Hyung Kim

Re: [petsc-users] Question about applying algebraic multigrid (AMG) PC and multilevel Schwarz PC

[Mark Adams]

On Fri, Feb 10, 2023 at 9:17 AM ­권승리  wrote:

> Dear petsc developers.
>
> Hello.
> While I was applying the preconditioner in contact with the FEM code, I
> had some questions.
>
> How can I apply the AMG PC and the multilevel Schwarz PC in the FEM code?
> Purpose is to compare which PC is better in the FEM code with contact
> situation.
>
> In my case, the Jacobi PC could be easily applied as shown in the example
> below. However, AMG and multilevel have too many options.
>
> ex) ./app -ksp_type gmres -pc_type jacobi -ksp_view
>

You can use the built-in AMG PC with: -pc_type gamg

Technically AMG is a  multilevel Schwarz PC, but I think you may be
referring to: -pc_type bddc

All solvers are equipped with default parameters that tend to be
optimized for robustness rather than peak performance and they are a good
place to start.

Thanks,
Mark


> Best regards
>
> Seung Lee Kwon
> --
> Seung Lee Kwon, Ph.D.Candidate
> Aerospace Structures and Materials Laboratory
> Department of Mechanical and Aerospace Engineering
> Seoul National University
> Building 300 Rm 503, Gwanak-ro 1, Gwanak-gu, Seoul, South Korea, 08826
> E-mail : ksl7...@snu.ac.kr
> Office : +82-2-880-7389
> C. P : +82-10-4695-1062
>

Re: [petsc-users] Question about MatGetRow

[Barry Smith]

> PETSC_OWN_POINTER,

indicates that the IS is to take ownership of the memory in indices, hence you 
are no longer reasonable for that memory and should not call PetscFree() on it. 
So the code is running correctly.


> On Feb 4, 2023, at 3:24 AM, 김성익  wrote:
> 
> Following your comments, that works in multi processes.
> After extracting procedure from 'newmat' which is submatrix.
> 
> PetscInt *indices;
> PetscMalloc1(1, );
> Indices[0] = 0;
> ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_OWN_POINTER, );
> MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);
> (extract from newmat)
> 
> I did 'ISDestroy(); PetscFree(indices);'
> However I got an error 'double free'.
> So I deleted PetscFree. 
> Is this correct response of that error?
> If not, how should I deal with that error??
> 
> Thanks,
> Hyung Kim 
> 
> 
> 
> 
> 
> 
> 
> 2023년 2월 3일 (금) 오후 11:33, Matthew Knepley  >님이 작성:
>> On Fri, Feb 3, 2023 at 9:04 AM 김성익 > > wrote:
>>> Actually in the last mail, below scripts are running in all processes.
>>> 
>>> IS isrow;
>>> PetscInt *indices;
>>> PetscMalloc1(1, );
>>> Indices[0] = 0;
>>> ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_OWN_POINTER, );
>>> MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);
>>> (extract from newmat)
>>> 
>>> However, you said it cannot get the values of first row of global matrix.
>>> Please let me know how can I fix this scripts for getting the 1st row of 
>>> global matrix in all processes.
>> 
>> Did you run and see what you get? If it is on all processes, it should work.
>> 
>>   Thanks,
>> 
>>  Matt
>>  
>>> Hyung Kim
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 2023년 2월 3일 (금) 오후 10:54, Matthew Knepley >> >님이 작성:
 On Fri, Feb 3, 2023 at 8:52 AM 김성익 >>> > wrote:
> I want to extract same row values of global matrix in all processes.
> Then how can I do this??
 
 Create the same IS on each process.
 
   THanks,
 
 Matt
  
> The case of same problem of vector, I just use vecscattertoall.
> However, I can't find same function for matrix.
> 
> Hyung Kim
> 
> 2023년 2월 3일 (금) 오후 10:47, Matthew Knepley  >님이 작성:
>> On Fri, Feb 3, 2023 at 8:45 AM 김성익 > > wrote:
>>> Following your comments,
>>> If I extract first row of below matrix.
>>> 
>>> IS isrow;
>>> PetscInt *indices;
>>> PetscMalloc1(1, *indices);
>> 
>> That should be 
>>  
>>> Indices[0] = 0;
>>> ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_COPY_VALUES, 
>>> );
>> 
>> You should use PETSC_OWN_POINTER.
>>  
>>> MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);
>>> 
>>> Then can I get the array about first row of global matrix in all 
>>> processes?
>> 
>> No, just on the process which gives 0. If you do that on every process, 
>> every rank with get row 0.
>> 
>>   Thanks,
>> 
>>  Matt
>>  
>>> Hyung Kim
>>> 
>>> 2023년 2월 3일 (금) 오후 10:26, Matthew Knepley >> >님이 작성:
 On Fri, Feb 3, 2023 at 8:06 AM 김성익 >>> > wrote:
> Following your comments,
> I want to check below things.
> For example, the global dense matrix are as below.
> 
> If I want to get first row ('1 2 0 0 3 0 0 4') in Proc 1.
> Then I should put 'MatCreateSubMatrix 
> (mat,
>  isrow, NULL, MAT_INITIAL_MATRIX, )'
> and isrow will be [0 1 2 3 4 5 6 7].
> 
> In this case, How can I make isrow?
> Actually I can't understand the procedure of handling isrow.
 
 You create an IS object of type ISGENERAL and give it the array of 
 global indices that you want to extract.
 
   Thanks,
 
  Matt
  
> Hyung Kim
> 
> 2023년 2월 3일 (금) 오후 9:03, Mark Adams  >님이 작성:
>> https://petsc.org/main/docs/manualpages/Mat/MatCreateSubMatrix/
>> 
>> Note, PETSc lets you give NULL arguments if there is a reasonable 
>> default.
>> In this case give NULL for the column IS and you will get the whole 
>> columns.
>> 
>> Mark
>> 
>> On Fri, Feb 3, 2023 at 4:05 AM 김성익 > > wrote:
>>> Hello,
>>> 
>>> 
>>> By using MatGetRow, user can get vectors from local matrix (at each 
>>> process).
>>> 
>>> However, I need other process's row values.
>>> So I have 2 questions.
>>> 
>>> 1. Is there any function for getting arrays from other process's??
>>> 
>>> 2. Or is 

Re: [petsc-users] Question about MatGetRow

[김성익]

Following your comments, that works in multi processes.
After extracting procedure from 'newmat' which is submatrix.

PetscInt *indices;
PetscMalloc1(1, );
Indices[0] = 0;
ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_OWN_POINTER, );
MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);
(extract from newmat)

I did 'ISDestroy(); PetscFree(indices);'
However I got an error 'double free'.
So I deleted PetscFree.
Is this correct response of that error?
If not, how should I deal with that error??

Thanks,
Hyung Kim







2023년 2월 3일 (금) 오후 11:33, Matthew Knepley 님이 작성:

> On Fri, Feb 3, 2023 at 9:04 AM 김성익  wrote:
>
>> Actually in the last mail, below scripts are running in all processes.
>>
>> IS isrow;
>> PetscInt *indices;
>> PetscMalloc1(1, );
>> Indices[0] = 0;
>> ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_OWN_POINTER, );
>> MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);
>> (extract from newmat)
>>
>> However, you said it cannot get the values of first row of global matrix.
>> Please let me know how can I fix this scripts for getting the 1st row of
>> global matrix in all processes.
>>
>
> Did you run and see what you get? If it is on all processes, it should
> work.
>
>   Thanks,
>
>  Matt
>
>
>> Hyung Kim
>>
>>
>>
>>
>>
>>
>>
>> 2023년 2월 3일 (금) 오후 10:54, Matthew Knepley 님이 작성:
>>
>>> On Fri, Feb 3, 2023 at 8:52 AM 김성익  wrote:
>>>
 I want to extract same row values of global matrix in all processes.
 Then how can I do this??

>>>
>>> Create the same IS on each process.
>>>
>>>   THanks,
>>>
>>> Matt
>>>
>>>
 The case of same problem of vector, I just use vecscattertoall.
 However, I can't find same function for matrix.

 Hyung Kim

 2023년 2월 3일 (금) 오후 10:47, Matthew Knepley 님이 작성:

> On Fri, Feb 3, 2023 at 8:45 AM 김성익  wrote:
>
>> Following your comments,
>> If I extract first row of below matrix.
>> [image: image.png]
>> IS isrow;
>> PetscInt *indices;
>> PetscMalloc1(1, *indices);
>>
>
> That should be 
>
>
>> Indices[0] = 0;
>> ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_COPY_VALUES,
>> );
>>
>
> You should use PETSC_OWN_POINTER.
>
>
>> MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);
>>
>> Then can I get the array about first row of global matrix in all
>> processes?
>>
>
> No, just on the process which gives 0. If you do that on every
> process, every rank with get row 0.
>
>   Thanks,
>
>  Matt
>
>
>> Hyung Kim
>>
>> 2023년 2월 3일 (금) 오후 10:26, Matthew Knepley 님이 작성:
>>
>>> On Fri, Feb 3, 2023 at 8:06 AM 김성익  wrote:
>>>
 Following your comments,
 I want to check below things.
 For example, the global dense matrix are as below.
 [image: image.png]
 If I want to get first row ('1 2 0 0 3 0 0 4') in Proc 1.
 Then I should put 'MatCreateSubMatrix
 (
 mat, isrow, NULL, MAT_INITIAL_MATRIX, *&*newmat)'
 and isrow will be [0 1 2 3 4 5 6 7].

 In this case, How can I make isrow?
 Actually I can't understand the procedure of handling isrow.

>>>
>>> You create an IS object of type ISGENERAL and give it the array of
>>> global indices that you want to extract.
>>>
>>>   Thanks,
>>>
>>>  Matt
>>>
>>>
 Hyung Kim

 2023년 2월 3일 (금) 오후 9:03, Mark Adams 님이 작성:

> https://petsc.org/main/docs/manualpages/Mat/MatCreateSubMatrix/
>
> Note, PETSc lets you give NULL arguments if there is a reasonable
> default.
> In this case give NULL for the column IS and you will get the
> whole columns.
>
> Mark
>
> On Fri, Feb 3, 2023 at 4:05 AM 김성익  wrote:
>
>> Hello,
>>
>>
>> By using MatGetRow, user can get vectors from local matrix (at
>> each process).
>>
>> However, I need other process's row values.
>> So I have 2 questions.
>>
>> 1. Is there any function for getting arrays from other process's??
>>
>> 2. Or is there any function like matrix version of
>> vecscattertoall??
>>
>> Thanks,
>> Hyung Kim
>>
>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which 
>>> their
>>> experiments lead.
>>> -- Norbert Wiener
>>>
>>> https://www.cse.buffalo.edu/~knepley/
>>> 
>>>
>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which 

Re: [petsc-users] Question about MatGetRow

[Matthew Knepley]

On Fri, Feb 3, 2023 at 9:04 AM 김성익  wrote:

> Actually in the last mail, below scripts are running in all processes.
>
> IS isrow;
> PetscInt *indices;
> PetscMalloc1(1, );
> Indices[0] = 0;
> ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_OWN_POINTER, );
> MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);
> (extract from newmat)
>
> However, you said it cannot get the values of first row of global matrix.
> Please let me know how can I fix this scripts for getting the 1st row of
> global matrix in all processes.
>

Did you run and see what you get? If it is on all processes, it should work.

  Thanks,

 Matt


> Hyung Kim
>
>
>
>
>
>
>
> 2023년 2월 3일 (금) 오후 10:54, Matthew Knepley 님이 작성:
>
>> On Fri, Feb 3, 2023 at 8:52 AM 김성익  wrote:
>>
>>> I want to extract same row values of global matrix in all processes.
>>> Then how can I do this??
>>>
>>
>> Create the same IS on each process.
>>
>>   THanks,
>>
>> Matt
>>
>>
>>> The case of same problem of vector, I just use vecscattertoall.
>>> However, I can't find same function for matrix.
>>>
>>> Hyung Kim
>>>
>>> 2023년 2월 3일 (금) 오후 10:47, Matthew Knepley 님이 작성:
>>>
 On Fri, Feb 3, 2023 at 8:45 AM 김성익  wrote:

> Following your comments,
> If I extract first row of below matrix.
> [image: image.png]
> IS isrow;
> PetscInt *indices;
> PetscMalloc1(1, *indices);
>

 That should be 


> Indices[0] = 0;
> ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_COPY_VALUES,
> );
>

 You should use PETSC_OWN_POINTER.


> MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);
>
> Then can I get the array about first row of global matrix in all
> processes?
>

 No, just on the process which gives 0. If you do that on every process,
 every rank with get row 0.

   Thanks,

  Matt


> Hyung Kim
>
> 2023년 2월 3일 (금) 오후 10:26, Matthew Knepley 님이 작성:
>
>> On Fri, Feb 3, 2023 at 8:06 AM 김성익  wrote:
>>
>>> Following your comments,
>>> I want to check below things.
>>> For example, the global dense matrix are as below.
>>> [image: image.png]
>>> If I want to get first row ('1 2 0 0 3 0 0 4') in Proc 1.
>>> Then I should put 'MatCreateSubMatrix
>>> (
>>> mat, isrow, NULL, MAT_INITIAL_MATRIX, *&*newmat)'
>>> and isrow will be [0 1 2 3 4 5 6 7].
>>>
>>> In this case, How can I make isrow?
>>> Actually I can't understand the procedure of handling isrow.
>>>
>>
>> You create an IS object of type ISGENERAL and give it the array of
>> global indices that you want to extract.
>>
>>   Thanks,
>>
>>  Matt
>>
>>
>>> Hyung Kim
>>>
>>> 2023년 2월 3일 (금) 오후 9:03, Mark Adams 님이 작성:
>>>
 https://petsc.org/main/docs/manualpages/Mat/MatCreateSubMatrix/

 Note, PETSc lets you give NULL arguments if there is a reasonable
 default.
 In this case give NULL for the column IS and you will get the whole
 columns.

 Mark

 On Fri, Feb 3, 2023 at 4:05 AM 김성익  wrote:

> Hello,
>
>
> By using MatGetRow, user can get vectors from local matrix (at
> each process).
>
> However, I need other process's row values.
> So I have 2 questions.
>
> 1. Is there any function for getting arrays from other process's??
>
> 2. Or is there any function like matrix version of
> vecscattertoall??
>
> Thanks,
> Hyung Kim
>

>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which 
>> their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> 
>>
>

 --
 What most experimenters take for granted before they begin their
 experiments is infinitely more interesting than any results to which their
 experiments lead.
 -- Norbert Wiener

 https://www.cse.buffalo.edu/~knepley/
 

>>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> 
>>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about MatGetRow

[김성익]

Actually in the last mail, below scripts are running in all processes.

IS isrow;
PetscInt *indices;
PetscMalloc1(1, );
Indices[0] = 0;
ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_OWN_POINTER, );
MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);
(extract from newmat)

However, you said it cannot get the values of first row of global matrix.
Please let me know how can I fix this scripts for getting the 1st row of
global matrix in all processes.

Hyung Kim







2023년 2월 3일 (금) 오후 10:54, Matthew Knepley 님이 작성:

> On Fri, Feb 3, 2023 at 8:52 AM 김성익  wrote:
>
>> I want to extract same row values of global matrix in all processes.
>> Then how can I do this??
>>
>
> Create the same IS on each process.
>
>   THanks,
>
> Matt
>
>
>> The case of same problem of vector, I just use vecscattertoall.
>> However, I can't find same function for matrix.
>>
>> Hyung Kim
>>
>> 2023년 2월 3일 (금) 오후 10:47, Matthew Knepley 님이 작성:
>>
>>> On Fri, Feb 3, 2023 at 8:45 AM 김성익  wrote:
>>>
 Following your comments,
 If I extract first row of below matrix.
 [image: image.png]
 IS isrow;
 PetscInt *indices;
 PetscMalloc1(1, *indices);

>>>
>>> That should be 
>>>
>>>
 Indices[0] = 0;
 ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_COPY_VALUES,
 );

>>>
>>> You should use PETSC_OWN_POINTER.
>>>
>>>
 MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);

 Then can I get the array about first row of global matrix in all
 processes?

>>>
>>> No, just on the process which gives 0. If you do that on every process,
>>> every rank with get row 0.
>>>
>>>   Thanks,
>>>
>>>  Matt
>>>
>>>
 Hyung Kim

 2023년 2월 3일 (금) 오후 10:26, Matthew Knepley 님이 작성:

> On Fri, Feb 3, 2023 at 8:06 AM 김성익  wrote:
>
>> Following your comments,
>> I want to check below things.
>> For example, the global dense matrix are as below.
>> [image: image.png]
>> If I want to get first row ('1 2 0 0 3 0 0 4') in Proc 1.
>> Then I should put 'MatCreateSubMatrix
>> (mat
>> , isrow, NULL, MAT_INITIAL_MATRIX, *&*newmat)'
>> and isrow will be [0 1 2 3 4 5 6 7].
>>
>> In this case, How can I make isrow?
>> Actually I can't understand the procedure of handling isrow.
>>
>
> You create an IS object of type ISGENERAL and give it the array of
> global indices that you want to extract.
>
>   Thanks,
>
>  Matt
>
>
>> Hyung Kim
>>
>> 2023년 2월 3일 (금) 오후 9:03, Mark Adams 님이 작성:
>>
>>> https://petsc.org/main/docs/manualpages/Mat/MatCreateSubMatrix/
>>>
>>> Note, PETSc lets you give NULL arguments if there is a reasonable
>>> default.
>>> In this case give NULL for the column IS and you will get the whole
>>> columns.
>>>
>>> Mark
>>>
>>> On Fri, Feb 3, 2023 at 4:05 AM 김성익  wrote:
>>>
 Hello,


 By using MatGetRow, user can get vectors from local matrix (at each
 process).

 However, I need other process's row values.
 So I have 2 questions.

 1. Is there any function for getting arrays from other process's??

 2. Or is there any function like matrix version of vecscattertoall??

 Thanks,
 Hyung Kim

>>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>

>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> experiments lead.
>>> -- Norbert Wiener
>>>
>>> https://www.cse.buffalo.edu/~knepley/
>>> 
>>>
>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>

Re: [petsc-users] Question about MatGetRow

[Matthew Knepley]

On Fri, Feb 3, 2023 at 8:52 AM 김성익  wrote:

> I want to extract same row values of global matrix in all processes.
> Then how can I do this??
>

Create the same IS on each process.

  THanks,

Matt


> The case of same problem of vector, I just use vecscattertoall.
> However, I can't find same function for matrix.
>
> Hyung Kim
>
> 2023년 2월 3일 (금) 오후 10:47, Matthew Knepley 님이 작성:
>
>> On Fri, Feb 3, 2023 at 8:45 AM 김성익  wrote:
>>
>>> Following your comments,
>>> If I extract first row of below matrix.
>>> [image: image.png]
>>> IS isrow;
>>> PetscInt *indices;
>>> PetscMalloc1(1, *indices);
>>>
>>
>> That should be 
>>
>>
>>> Indices[0] = 0;
>>> ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_COPY_VALUES,
>>> );
>>>
>>
>> You should use PETSC_OWN_POINTER.
>>
>>
>>> MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);
>>>
>>> Then can I get the array about first row of global matrix in all
>>> processes?
>>>
>>
>> No, just on the process which gives 0. If you do that on every process,
>> every rank with get row 0.
>>
>>   Thanks,
>>
>>  Matt
>>
>>
>>> Hyung Kim
>>>
>>> 2023년 2월 3일 (금) 오후 10:26, Matthew Knepley 님이 작성:
>>>
 On Fri, Feb 3, 2023 at 8:06 AM 김성익  wrote:

> Following your comments,
> I want to check below things.
> For example, the global dense matrix are as below.
> [image: image.png]
> If I want to get first row ('1 2 0 0 3 0 0 4') in Proc 1.
> Then I should put 'MatCreateSubMatrix
> (mat,
> isrow, NULL, MAT_INITIAL_MATRIX, *&*newmat)'
> and isrow will be [0 1 2 3 4 5 6 7].
>
> In this case, How can I make isrow?
> Actually I can't understand the procedure of handling isrow.
>

 You create an IS object of type ISGENERAL and give it the array of
 global indices that you want to extract.

   Thanks,

  Matt


> Hyung Kim
>
> 2023년 2월 3일 (금) 오후 9:03, Mark Adams 님이 작성:
>
>> https://petsc.org/main/docs/manualpages/Mat/MatCreateSubMatrix/
>>
>> Note, PETSc lets you give NULL arguments if there is a reasonable
>> default.
>> In this case give NULL for the column IS and you will get the whole
>> columns.
>>
>> Mark
>>
>> On Fri, Feb 3, 2023 at 4:05 AM 김성익  wrote:
>>
>>> Hello,
>>>
>>>
>>> By using MatGetRow, user can get vectors from local matrix (at each
>>> process).
>>>
>>> However, I need other process's row values.
>>> So I have 2 questions.
>>>
>>> 1. Is there any function for getting arrays from other process's??
>>>
>>> 2. Or is there any function like matrix version of vecscattertoall??
>>>
>>> Thanks,
>>> Hyung Kim
>>>
>>

 --
 What most experimenters take for granted before they begin their
 experiments is infinitely more interesting than any results to which their
 experiments lead.
 -- Norbert Wiener

 https://www.cse.buffalo.edu/~knepley/
 

>>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> 
>>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about MatGetRow

[김성익]

I want to extract same row values of global matrix in all processes.
Then how can I do this??

The case of same problem of vector, I just use vecscattertoall.
However, I can't find same function for matrix.

Hyung Kim

2023년 2월 3일 (금) 오후 10:47, Matthew Knepley 님이 작성:

> On Fri, Feb 3, 2023 at 8:45 AM 김성익  wrote:
>
>> Following your comments,
>> If I extract first row of below matrix.
>> [image: image.png]
>> IS isrow;
>> PetscInt *indices;
>> PetscMalloc1(1, *indices);
>>
>
> That should be 
>
>
>> Indices[0] = 0;
>> ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_COPY_VALUES, );
>>
>
> You should use PETSC_OWN_POINTER.
>
>
>> MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);
>>
>> Then can I get the array about first row of global matrix in all
>> processes?
>>
>
> No, just on the process which gives 0. If you do that on every process,
> every rank with get row 0.
>
>   Thanks,
>
>  Matt
>
>
>> Hyung Kim
>>
>> 2023년 2월 3일 (금) 오후 10:26, Matthew Knepley 님이 작성:
>>
>>> On Fri, Feb 3, 2023 at 8:06 AM 김성익  wrote:
>>>
 Following your comments,
 I want to check below things.
 For example, the global dense matrix are as below.
 [image: image.png]
 If I want to get first row ('1 2 0 0 3 0 0 4') in Proc 1.
 Then I should put 'MatCreateSubMatrix
 (mat,
 isrow, NULL, MAT_INITIAL_MATRIX, *&*newmat)'
 and isrow will be [0 1 2 3 4 5 6 7].

 In this case, How can I make isrow?
 Actually I can't understand the procedure of handling isrow.

>>>
>>> You create an IS object of type ISGENERAL and give it the array of
>>> global indices that you want to extract.
>>>
>>>   Thanks,
>>>
>>>  Matt
>>>
>>>
 Hyung Kim

 2023년 2월 3일 (금) 오후 9:03, Mark Adams 님이 작성:

> https://petsc.org/main/docs/manualpages/Mat/MatCreateSubMatrix/
>
> Note, PETSc lets you give NULL arguments if there is a reasonable
> default.
> In this case give NULL for the column IS and you will get the whole
> columns.
>
> Mark
>
> On Fri, Feb 3, 2023 at 4:05 AM 김성익  wrote:
>
>> Hello,
>>
>>
>> By using MatGetRow, user can get vectors from local matrix (at each
>> process).
>>
>> However, I need other process's row values.
>> So I have 2 questions.
>>
>> 1. Is there any function for getting arrays from other process's??
>>
>> 2. Or is there any function like matrix version of vecscattertoall??
>>
>> Thanks,
>> Hyung Kim
>>
>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> experiments lead.
>>> -- Norbert Wiener
>>>
>>> https://www.cse.buffalo.edu/~knepley/
>>> 
>>>
>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>

Re: [petsc-users] Question about MatGetRow

[김성익]

Following your comments,
If I extract first row of below matrix.
[image: image.png]
IS isrow;
PetscInt *indices;
PetscMalloc1(1, *indices);
Indices[0] = 0;
ISCreateGenreral(PETSC_COMM_WORLD, 1, indices, PETSC_COPY_VALUES, );
MatCreateSubMatrix(mat,isrow,NULL,MAT_INITIAL_MATRIX,);

Then can I get the array about first row of global matrix in all processes?

Hyung Kim

2023년 2월 3일 (금) 오후 10:26, Matthew Knepley 님이 작성:

> On Fri, Feb 3, 2023 at 8:06 AM 김성익  wrote:
>
>> Following your comments,
>> I want to check below things.
>> For example, the global dense matrix are as below.
>> [image: image.png]
>> If I want to get first row ('1 2 0 0 3 0 0 4') in Proc 1.
>> Then I should put 'MatCreateSubMatrix
>> (mat,
>> isrow, NULL, MAT_INITIAL_MATRIX, *&*newmat)'
>> and isrow will be [0 1 2 3 4 5 6 7].
>>
>> In this case, How can I make isrow?
>> Actually I can't understand the procedure of handling isrow.
>>
>
> You create an IS object of type ISGENERAL and give it the array of global
> indices that you want to extract.
>
>   Thanks,
>
>  Matt
>
>
>> Hyung Kim
>>
>> 2023년 2월 3일 (금) 오후 9:03, Mark Adams 님이 작성:
>>
>>> https://petsc.org/main/docs/manualpages/Mat/MatCreateSubMatrix/
>>>
>>> Note, PETSc lets you give NULL arguments if there is a reasonable
>>> default.
>>> In this case give NULL for the column IS and you will get the whole
>>> columns.
>>>
>>> Mark
>>>
>>> On Fri, Feb 3, 2023 at 4:05 AM 김성익  wrote:
>>>
 Hello,


 By using MatGetRow, user can get vectors from local matrix (at each
 process).

 However, I need other process's row values.
 So I have 2 questions.

 1. Is there any function for getting arrays from other process's??

 2. Or is there any function like matrix version of vecscattertoall??

 Thanks,
 Hyung Kim

>>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>

Re: [petsc-users] Question about MatGetRow

[Matthew Knepley]

On Fri, Feb 3, 2023 at 8:06 AM 김성익  wrote:

> Following your comments,
> I want to check below things.
> For example, the global dense matrix are as below.
> [image: image.png]
> If I want to get first row ('1 2 0 0 3 0 0 4') in Proc 1.
> Then I should put 'MatCreateSubMatrix
> (mat,
> isrow, NULL, MAT_INITIAL_MATRIX, *&*newmat)'
> and isrow will be [0 1 2 3 4 5 6 7].
>
> In this case, How can I make isrow?
> Actually I can't understand the procedure of handling isrow.
>

You create an IS object of type ISGENERAL and give it the array of global
indices that you want to extract.

  Thanks,

 Matt


> Hyung Kim
>
> 2023년 2월 3일 (금) 오후 9:03, Mark Adams 님이 작성:
>
>> https://petsc.org/main/docs/manualpages/Mat/MatCreateSubMatrix/
>>
>> Note, PETSc lets you give NULL arguments if there is a reasonable default.
>> In this case give NULL for the column IS and you will get the whole
>> columns.
>>
>> Mark
>>
>> On Fri, Feb 3, 2023 at 4:05 AM 김성익  wrote:
>>
>>> Hello,
>>>
>>>
>>> By using MatGetRow, user can get vectors from local matrix (at each
>>> process).
>>>
>>> However, I need other process's row values.
>>> So I have 2 questions.
>>>
>>> 1. Is there any function for getting arrays from other process's??
>>>
>>> 2. Or is there any function like matrix version of vecscattertoall??
>>>
>>> Thanks,
>>> Hyung Kim
>>>
>>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about MatGetRow

[김성익]

Following your comments,
I want to check below things.
For example, the global dense matrix are as below.
[image: image.png]
If I want to get first row ('1 2 0 0 3 0 0 4') in Proc 1.
Then I should put 'MatCreateSubMatrix
(mat, isrow
, NULL, MAT_INITIAL_MATRIX, *&*newmat)'
and isrow will be [0 1 2 3 4 5 6 7].

In this case, How can I make isrow?
Actually I can't understand the procedure of handling isrow.

Hyung Kim

2023년 2월 3일 (금) 오후 9:03, Mark Adams 님이 작성:

> https://petsc.org/main/docs/manualpages/Mat/MatCreateSubMatrix/
>
> Note, PETSc lets you give NULL arguments if there is a reasonable default.
> In this case give NULL for the column IS and you will get the whole
> columns.
>
> Mark
>
> On Fri, Feb 3, 2023 at 4:05 AM 김성익  wrote:
>
>> Hello,
>>
>>
>> By using MatGetRow, user can get vectors from local matrix (at each
>> process).
>>
>> However, I need other process's row values.
>> So I have 2 questions.
>>
>> 1. Is there any function for getting arrays from other process's??
>>
>> 2. Or is there any function like matrix version of vecscattertoall??
>>
>> Thanks,
>> Hyung Kim
>>
>

Re: [petsc-users] Question about MatGetRow

[Mark Adams]

https://petsc.org/main/docs/manualpages/Mat/MatCreateSubMatrix/

Note, PETSc lets you give NULL arguments if there is a reasonable default.
In this case give NULL for the column IS and you will get the whole columns.

Mark

On Fri, Feb 3, 2023 at 4:05 AM 김성익  wrote:

> Hello,
>
>
> By using MatGetRow, user can get vectors from local matrix (at each
> process).
>
> However, I need other process's row values.
> So I have 2 questions.
>
> 1. Is there any function for getting arrays from other process's??
>
> 2. Or is there any function like matrix version of vecscattertoall??
>
> Thanks,
> Hyung Kim
>

Re: [petsc-users] Question about handling matrix

[김성익]

Both small and large matrix are sparse matrix.
I will try following comments.

Thanks,
Hyung Kim

2023년 2월 1일 (수) 오후 11:30, Jed Brown 님이 작성:

> Is the small matrix dense? Then you can use MatSetValues. If the small
> matrix is sparse, you can assemble it with larger dimension (empty rows and
> columns) and use MatAXPY.
>
> 김성익  writes:
>
> > Hello,
> >
> >
> > I want to put small matrix to large matrix.
> > The schematic of operation is as below.
> > [image: image.png]
> > Is there any function for put small matrix to large matrix at once?
> >
> > Thanks,
> > Hyung Kim
>

Re: [petsc-users] Question about handling matrix

[Jed Brown]

Is the small matrix dense? Then you can use MatSetValues. If the small matrix 
is sparse, you can assemble it with larger dimension (empty rows and columns) 
and use MatAXPY.

김성익  writes:

> Hello,
>
>
> I want to put small matrix to large matrix.
> The schematic of operation is as below.
> [image: image.png]
> Is there any function for put small matrix to large matrix at once?
>
> Thanks,
> Hyung Kim

Re: [petsc-users] Question about handling matrix

[Mark Adams]

Maybe create the large matrix and use
https://www.mcs.anl.gov/petsc/petsc-3.7/docs/manualpages/Mat/MatGetSubMatrix.html
with MAT_REUSE_MATRIX.

Or pad the IS arguments to MatGetSubMatrix with -1, so the size is correct
and PETSc should ignore -1, and use MAT_INITIAL_MATRIX.

Others may know what will work.

Good luck,
Mark

On Wed, Feb 1, 2023 at 5:26 AM 김성익  wrote:

> Hello,
>
>
> I want to put small matrix to large matrix.
> The schematic of operation is as below.
> [image: image.png]
> Is there any function for put small matrix to large matrix at once?
>
> Thanks,
> Hyung Kim
>

Re: [petsc-users] Question - about the 'Hint for performance tuning'

[Barry Smith]

> On Jan 2, 2023, at 9:27 AM, 김성익  wrote:
> 
> There are more questions.
> 
> 1. Following your comments, but there is an error as below.
> 
> How can I fix this?

   The previous value of PETSC_ARCH was not arch-main-debug do ls $PETSC_DIR 
and locate the directory whose name begins with arch- that will tell you the 
previous value used for PETSC_ARCH.

> 
> 2. After changing the optimized build, then how can I set the debug mode 
> again?

   export PETSC_ARCH=arch- whatever the previous value was and then recompile 
your code (note you do not need to recompile PETSc just your executable).

> 
> 3.Following your comments, the new makefile is as below. Is it right?
> CFLAGS   = -O3
> FFLAGS   =
> CPPFLAGS =
> FPPFLAGS =
> COPTFLAGS= -march=native
>  
> app : a1.o a2.o a3.o a4.o
>$(LINK.C) -o $@ $^ $(LDLIBS)
>  
> include ${PETSC_DIR}/lib/petsc/conf/rules
> include ${PETSC_DIR}/lib/petsc/conf/test

Best not to set these values in the makefile at all because they will affect 
all compilers. Just set them with ./configure CCOPTFLAGS="-O3 -march=native"

> 
> 
> 4. I have no such processors. Where can I find benchmark information about 
> STREAMS?

  do make mpistreams in PETSC_DIR


> 
> 
> Thanks, 
> Hyung Kim
> 
> 
> 
> 
> 
> 
> 
> 
> 2023년 1월 2일 (월) 오후 11:03, Matthew Knepley  >님이 작성:
>> On Mon, Jan 2, 2023 at 4:16 AM 김성익 > > wrote:
>>> Hello,
>>> 
>>> Happy new year!!
>>> 
>>>  
>>> I have some questions about “Hint for performance tuning” in user guide of 
>>> petsc.
>>> 
>>>  
>>> 1. In the “Performance Pitfalls and Advice” section, there are 2 modes 
>>> “debug” and “optimized builds. My current setup is debug mode. So I want to 
>>> change for test the performance the optimized build mode. However, if I 
>>> configure again, does the existing debug mode disappear? Is there any way 
>>> to coexist the 2 modes and use them in the run the application?
>>> 
>> Suppose your current arch is named "arch-main-debug". Then you can make an 
>> optimized version using
>> 
>>   cd $PETSC_DIR
>>   ./arch-main-debug/lib/petsc/conf/reconfigure-arch-main-debug.py 
>> --with-debugging=0 --PETSC_ARCH=arch-main-opt
>>  
>>> 2. In the guide, there are some paragraphs about optimization level of 
>>> compiler. To control the optimization level of compiler, I put the ‘-O3’ as 
>>> below. Is this right??
>>> 
>>> CFLAGS   = -O3
>>> FFLAGS   =
>>> CPPFLAGS =
>>> FPPFLAGS =
>>>  
>>> app : a1.o a2.o a3.o a4.o
>>>$(LINK.C) -o $@ $^ $(LDLIBS)
>>>  
>>> include ${PETSC_DIR}/lib/petsc/conf/rules
>>> include ${PETSC_DIR}/lib/petsc/conf/test
>> 
>> You could dp this, but that only changes it for that directory. It is best 
>> to do it by reconfiguring.
>>  
>>> 3. In the guide, user should put ‘-march=native’ for using AVX2 or 
>>> AVX-512. Where should I put the ‘-march=native’ for using AVX?
>>> 
>> You can add --COPTFLAGS="" with any flags you want to the 
>> configure.
>> 
>>> 4. After read the “Hint for performance tuning” I understood that for 
>>> good performance and scalability user should use the multiple node and 
>>> multiple socket . However, before composing cluster system, many users just 
>>> can use desktop system. 
>>> In that case, between intel 13th i9 and amd ryzen 7950x, can the 7950x, 
>>> which has an arichitecture similar to the server processor, be more 
>>> suitable for petsc? (Because the architecture of intel desktop cpu is 
>>> big.little.)
>>> 
>> 
>> A good guide is to run the STREAMS benchmark on the processor. PETSc 
>> performance closely tracks that.
>> 
>>   Thanks,
>> 
>>  Matt
>>  
>>>  
>>> Thanks,
>>> 
>>> Hyung Kim
>>> 
>> 
>> 
>> -- 
>> What most experimenters take for granted before they begin their experiments 
>> is infinitely more interesting than any results to which their experiments 
>> lead.
>> -- Norbert Wiener
>> 
>> https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question - about the 'Hint for performance tuning'

[Matthew Knepley]

On Mon, Jan 2, 2023 at 4:16 AM 김성익  wrote:

> Hello,
>
> Happy new year!!
>
>
>
> I have some questions about “Hint for performance tuning” in user guide of
> petsc.
>
>
>
> 1. In the “Performance Pitfalls and Advice” section, there are 2
> modes “debug” and “optimized builds. My current setup is debug mode. So I
> want to change for test the performance the optimized build mode. However,
> if I configure again, does the existing debug mode disappear? Is there any
> way to coexist the 2 modes and use them in the run the application?
>
Suppose your current arch is named "arch-main-debug". Then you can make an
optimized version using

  cd $PETSC_DIR
  ./arch-main-debug/lib/petsc/conf/reconfigure-arch-main-debug.py
--with-debugging=0 --PETSC_ARCH=arch-main-opt


> 2. In the guide, there are some paragraphs about optimization level
> of compiler. To control the optimization level of compiler, I put the ‘-O3’
> as below. Is this right??
>
> CFLAGS   = -O3
>
> FFLAGS   =
>
> CPPFLAGS =
>
> FPPFLAGS =
>
>
>
> app : a1.o a2.o a3.o a4.o
>
>$(LINK.C) -o $@ $^ $(LDLIBS)
>
>
>
> include ${PETSC_DIR}/lib/petsc/conf/rules
>
> include ${PETSC_DIR}/lib/petsc/conf/test
>

You could dp this, but that only changes it for that directory. It is best
to do it by reconfiguring.


> 3. In the guide, user should put ‘-march=native’ for using AVX2 or
> AVX-512. Where should I put the ‘-march=native’ for using AVX?
>
You can add --COPTFLAGS="" with any flags you want to the
configure.

4. After read the “Hint for performance tuning” I understood that for
> good performance and scalability user should use the multiple node and
> multiple socket . However, before composing cluster system, many users just
> can use desktop system.
> In that case, between intel 13th i9 and amd ryzen 7950x, can the 7950x,
> which has an arichitecture similar to the server processor, be more
> suitable for petsc? (Because the architecture of intel desktop cpu is
> big.little.)
>

A good guide is to run the STREAMS benchmark on the processor. PETSc
performance closely tracks that.

  Thanks,

 Matt


>
>
> Thanks,
>
> Hyung Kim
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question-Memory of matsetvalue

[Jed Brown]

This is what I'd expect to observe if you didn't preallocate correctly for the 
second matrix, which has more nonzeros per row.

https://petsc.org/release/docs/manual/mat/#sec-matsparse

김성익  writes:

> Hello,
>
>
>
> I have a question about memory of matsetvalue.
>
> When I assembly the local matrix to global matrix, I’m just using
> matsetvalue.
> Because the connectivity is so complex I can’t use matsetvalues.
>
> I asked this question because I was curious about how ram memory is
> allocated differently for the two simtuations below.
>
>
>
> First situation.
>
> The local matrix size is 24 by 24. And the number of local matrix is
> 125,000.
>
> When assembly procedure by using matsetvalue, memory allocation does not
> increase.
> So I just put Matassemblybegin and matassemblyend after all matsetvalue.
>
>
>
> Second situation.
>
> The local matrix size is 60 by 60. And the number of local matrix is 27,000.
>
> When assembly procedure by using matsetvalue, memory allocation does
> increase.
>
> So I just put Matassemblybegin and matassemblyend after each local matrix
> assembly.
> This did not increase the memory further..
>
>
>
> Why this situation is happen?
>
> And is there any way to prevent the memory allocation from increasing?
>
>
>
>
>
>
>
> Thanks,
>
> Hyung Kim

Re: [petsc-users] Question-Memory of matsetvalue

[Matthew Knepley]

On Fri, Dec 30, 2022 at 4:36 AM 김성익  wrote:

> Hello,
>
>
>
> I have a question about memory of matsetvalue.
>
> When I assembly the local matrix to global matrix, I’m just using
> matsetvalue.
> Because the connectivity is so complex I can’t use matsetvalues.
>
> I asked this question because I was curious about how ram memory is
> allocated differently for the two simtuations below.
>
>
>
> First situation.
>
> The local matrix size is 24 by 24. And the number of local matrix is
> 125,000.
>
> When assembly procedure by using matsetvalue, memory allocation does not
> increase.
> So I just put Matassemblybegin and matassemblyend after all matsetvalue.
>
>
>
> Second situation.
>
> The local matrix size is 60 by 60. And the number of local matrix is
> 27,000.
>
> When assembly procedure by using matsetvalue, memory allocation does
> increase.
>
> So I just put Matassemblybegin and matassemblyend after each local matrix
> assembly.
> This did not increase the memory further..
>
>
>
> Why this situation is happen?
>
> And is there any way to prevent the memory allocation from increasing?
>

Matrix assembly has two stages. First you insert entries into the local
portion of your parallel matrix
using MatSetValue(s). If all values you try to insert are local, this is
the end.

However, if you try to insert values that are local to another process, we
store those values. When you
call MatAssemblyBegin/End(), we communicate them to the correct process and
insert.

For a scalable code, you should insert most values on the correct process.
If not, significant memory
can be consumed storing these values. Anywhere in the assembly process you
can call

  MatAssemblyBegin(A, MAT_ASSEMBLY_FLUSH);
  MatAssemblyEnd(A, MAT_ASSEMBLY_FLUSH);

This will communicate the cache of values, but not end the assembly.

   Thanks,

   Matt



> Thanks,
>
> Hyung Kim
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about eigenvalue and eigenvectors

[Matthew Knepley]

Or just call the LAPACK routine directly.

   Matt

On Sat, Dec 24, 2022 at 7:14 AM Stefano Zampini 
wrote:

> For 3x3 matrices you can use explicit formulas
>
> On Sat, Dec 24, 2022, 11:20 김성익  wrote:
>
>> Hello,
>>
>>
>> I tried to calculate the eigenvalues and eigenvectors in 3 by 3 matrix
>> (real and nonsymmetric).
>> I already checked the kspcomputeeigenvalues and kspcomputeritz.
>>
>> However, the target matrix is just 3 by 3 matrix.
>> So I need another way to calculate the values and vectors.
>> Can anyone recommend other methods that are efficient for such small size
>> problems??
>>
>> Thanks,
>> Hyung Kim
>>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question about eigenvalue and eigenvectors

[Stefano Zampini]

For 3x3 matrices you can use explicit formulas

On Sat, Dec 24, 2022, 11:20 김성익  wrote:

> Hello,
>
>
> I tried to calculate the eigenvalues and eigenvectors in 3 by 3 matrix
> (real and nonsymmetric).
> I already checked the kspcomputeeigenvalues and kspcomputeritz.
>
> However, the target matrix is just 3 by 3 matrix.
> So I need another way to calculate the values and vectors.
> Can anyone recommend other methods that are efficient for such small size
> problems??
>
> Thanks,
> Hyung Kim
>

Re: [petsc-users] Question About Assembly matrix and declaration of KSP & pc

[Mark Adams]

On Wed, Nov 30, 2022 at 8:31 AM Matthew Knepley  wrote:

> On Wed, Nov 30, 2022 at 8:25 AM 김성익  wrote:
>
>>
>> In your comments,
>> KSPSetOperators is called if you want to change the system matrix.
>>
>> "change the system matrix" means the components of matrix are changed?
>> I mean the values of some components of matrix are changed.
>>
>
> If you just change values in the matrix, you do not have to call it again.
>

You do need to call KSPSetOperators before the solve if you want the KSP to
redo the (PC) setup.
In GAMG this is substantial. If the matrix does not change much you can try
not doing this and experiment.


>
>   Thanks,
>
>  Matt
>
>
>> Thanks,
>> Hyung Kim
>>
>>
>> 2022년 11월 30일 (수) 오후 10:00, Matthew Knepley 님이 작성:
>>
>>> On Wed, Nov 30, 2022 at 7:51 AM 김성익  wrote:
>>>
 Thank you for your comments.
 However, I have more questions.

 1. Generally,  (KSPCreate, KSPSetOperators, KSPGetPC, PCSetType,
 PCFactorSetMatSolverType, KSPSetFromOptions )
 above functions are should be called after  each "MatassemblyEnd??"

>>>
>>> KSPCreate is called once.
>>>
>>> You do not need PCSetType, PCFactorSetMatSolverType, KSPSetFromOptions
>>> more than once, unless you want to change the solver type.
>>>
>>> KSPSetOperators is called if you want to change the system matrix.
>>>
>>> KSPSolve is called when you want to change the rhs.
>>>
>>>   Thanks,
>>>
>>>  Matt
>>>
>>>
 2.   Though reading the user guide, I don't fully understand under what
 circumstances the functions mentioned above should be called again. Can you
 explain when each function should be called?

 Thanks,

 Hyung Kim

 2022년 11월 30일 (수) 오후 8:37, Mark Adams 님이 작성:

>
>
> On Wed, Nov 30, 2022 at 5:08 AM 김성익  wrote:
>
>> Hello,
>>
>>
>>
>> I’m working on FEM using PETSc.
>>
>> As everyone knows, it is necessary to repeatedly solve Ax=B.
>>
>> Regarding this, I have 4 questions.
>>
>>
>>
>> 1. There are many steps for preparing KSPSolve. For example
>> KSPcreate, KSPSetOperators, KSPGetPC, PCSetType, 
>> PCFactorSetMatSolverType,
>> KSPSetFromOptions…
>> In Nonlinear FEM, there are repeatedly kspsolve for getting answer
>> vector. Is it correct to do all of the aforementioned processes 
>> (KSPcreate,
>> KSPSetOperators ~~~) for each KSPSolve? Or should I declare it only once 
>> at
>> the beginning and not call it again?
>>
>
> You just do these once at setup but for nonlinear problems
> KSPSetOperators tells the solver that you have a new matrix and so "matrix
> setup" work needs to be done.
>
>
>>
>> 2. If the answer to question 1 is that it must be repeated every
>> time, should this work be done right before kspsolve, that is, when the
>> global matrix assembly is finished, or is it irrelevant to performance at
>> any time?
>>
>
> KSPSetOperators should be set after the new matrix values are set but
> it might work before. It just sets a pointer to the matrix and flags it as
> not setup.
>
>
>>
>>
>> 3. When performing FEM, local matrices are often scattered in
>> global matrices depending on connectivity. In this case, which is better 
>> in
>> terms of performance: adding the values one by one with MatSetValue or
>> adding them all at once with MatSetValues even if they are scattered?
>>
>
> You want to add one element matrix at a time, generally.
>
>
>>
>>
>>
>>
>> 4. I would like to measure the time of each section of the
>> process. Which method is recommended?
>>
>
> PETSc methods are timed separateluy, but setup gets folded into
> KSPSolve unless you call SNESSetUp before the SNES[KSP]Solve.
> You can add you own timers also
> https://petsc.org/release/docs/manualpages/Profiling/PetscLogEventRegister/
>
> Mark
>
>
>
>>
>>
>> Thank you for your help.
>>
>>
>>
>> Hyung Kim
>>
>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> experiments lead.
>>> -- Norbert Wiener
>>>
>>> https://www.cse.buffalo.edu/~knepley/
>>> 
>>>
>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>

Re: [petsc-users] Question About Assembly matrix and declaration of KSP & pc

[Matthew Knepley]

On Wed, Nov 30, 2022 at 8:25 AM 김성익  wrote:

>
> In your comments,
> KSPSetOperators is called if you want to change the system matrix.
>
> "change the system matrix" means the components of matrix are changed?
> I mean the values of some components of matrix are changed.
>

If you just change values in the matrix, you do not have to call it again.

  Thanks,

 Matt


> Thanks,
> Hyung Kim
>
>
> 2022년 11월 30일 (수) 오후 10:00, Matthew Knepley 님이 작성:
>
>> On Wed, Nov 30, 2022 at 7:51 AM 김성익  wrote:
>>
>>> Thank you for your comments.
>>> However, I have more questions.
>>>
>>> 1. Generally,  (KSPCreate, KSPSetOperators, KSPGetPC, PCSetType,
>>> PCFactorSetMatSolverType, KSPSetFromOptions )
>>> above functions are should be called after  each "MatassemblyEnd??"
>>>
>>
>> KSPCreate is called once.
>>
>> You do not need PCSetType, PCFactorSetMatSolverType, KSPSetFromOptions
>> more than once, unless you want to change the solver type.
>>
>> KSPSetOperators is called if you want to change the system matrix.
>>
>> KSPSolve is called when you want to change the rhs.
>>
>>   Thanks,
>>
>>  Matt
>>
>>
>>> 2.   Though reading the user guide, I don't fully understand under what
>>> circumstances the functions mentioned above should be called again. Can you
>>> explain when each function should be called?
>>>
>>> Thanks,
>>>
>>> Hyung Kim
>>>
>>> 2022년 11월 30일 (수) 오후 8:37, Mark Adams 님이 작성:
>>>


 On Wed, Nov 30, 2022 at 5:08 AM 김성익  wrote:

> Hello,
>
>
>
> I’m working on FEM using PETSc.
>
> As everyone knows, it is necessary to repeatedly solve Ax=B.
>
> Regarding this, I have 4 questions.
>
>
>
> 1. There are many steps for preparing KSPSolve. For example
> KSPcreate, KSPSetOperators, KSPGetPC, PCSetType, PCFactorSetMatSolverType,
> KSPSetFromOptions…
> In Nonlinear FEM, there are repeatedly kspsolve for getting answer
> vector. Is it correct to do all of the aforementioned processes 
> (KSPcreate,
> KSPSetOperators ~~~) for each KSPSolve? Or should I declare it only once 
> at
> the beginning and not call it again?
>

 You just do these once at setup but for nonlinear problems
 KSPSetOperators tells the solver that you have a new matrix and so "matrix
 setup" work needs to be done.


>
> 2. If the answer to question 1 is that it must be repeated every
> time, should this work be done right before kspsolve, that is, when the
> global matrix assembly is finished, or is it irrelevant to performance at
> any time?
>

 KSPSetOperators should be set after the new matrix values are set but
 it might work before. It just sets a pointer to the matrix and flags it as
 not setup.


>
>
> 3. When performing FEM, local matrices are often scattered in
> global matrices depending on connectivity. In this case, which is better 
> in
> terms of performance: adding the values one by one with MatSetValue or
> adding them all at once with MatSetValues even if they are scattered?
>

 You want to add one element matrix at a time, generally.


>
>
>
>
> 4. I would like to measure the time of each section of the
> process. Which method is recommended?
>

 PETSc methods are timed separateluy, but setup gets folded into
 KSPSolve unless you call SNESSetUp before the SNES[KSP]Solve.
 You can add you own timers also
 https://petsc.org/release/docs/manualpages/Profiling/PetscLogEventRegister/

 Mark



>
>
> Thank you for your help.
>
>
>
> Hyung Kim
>

>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> 
>>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question About Assembly matrix and declaration of KSP & pc

[김성익]

In your comments,
KSPSetOperators is called if you want to change the system matrix.

"change the system matrix" means the components of matrix are changed?
I mean the values of some components of matrix are changed.


Thanks,
Hyung Kim


2022년 11월 30일 (수) 오후 10:00, Matthew Knepley 님이 작성:

> On Wed, Nov 30, 2022 at 7:51 AM 김성익  wrote:
>
>> Thank you for your comments.
>> However, I have more questions.
>>
>> 1. Generally,  (KSPCreate, KSPSetOperators, KSPGetPC, PCSetType,
>> PCFactorSetMatSolverType, KSPSetFromOptions )
>> above functions are should be called after  each "MatassemblyEnd??"
>>
>
> KSPCreate is called once.
>
> You do not need PCSetType, PCFactorSetMatSolverType, KSPSetFromOptions
> more than once, unless you want to change the solver type.
>
> KSPSetOperators is called if you want to change the system matrix.
>
> KSPSolve is called when you want to change the rhs.
>
>   Thanks,
>
>  Matt
>
>
>> 2.   Though reading the user guide, I don't fully understand under what
>> circumstances the functions mentioned above should be called again. Can you
>> explain when each function should be called?
>>
>> Thanks,
>>
>> Hyung Kim
>>
>> 2022년 11월 30일 (수) 오후 8:37, Mark Adams 님이 작성:
>>
>>>
>>>
>>> On Wed, Nov 30, 2022 at 5:08 AM 김성익  wrote:
>>>
 Hello,



 I’m working on FEM using PETSc.

 As everyone knows, it is necessary to repeatedly solve Ax=B.

 Regarding this, I have 4 questions.



 1. There are many steps for preparing KSPSolve. For example
 KSPcreate, KSPSetOperators, KSPGetPC, PCSetType, PCFactorSetMatSolverType,
 KSPSetFromOptions…
 In Nonlinear FEM, there are repeatedly kspsolve for getting answer
 vector. Is it correct to do all of the aforementioned processes (KSPcreate,
 KSPSetOperators ~~~) for each KSPSolve? Or should I declare it only once at
 the beginning and not call it again?

>>>
>>> You just do these once at setup but for nonlinear problems
>>> KSPSetOperators tells the solver that you have a new matrix and so "matrix
>>> setup" work needs to be done.
>>>
>>>

 2. If the answer to question 1 is that it must be repeated every
 time, should this work be done right before kspsolve, that is, when the
 global matrix assembly is finished, or is it irrelevant to performance at
 any time?

>>>
>>> KSPSetOperators should be set after the new matrix values are set but it
>>> might work before. It just sets a pointer to the matrix and flags it as not
>>> setup.
>>>
>>>


 3. When performing FEM, local matrices are often scattered in
 global matrices depending on connectivity. In this case, which is better in
 terms of performance: adding the values one by one with MatSetValue or
 adding them all at once with MatSetValues even if they are scattered?

>>>
>>> You want to add one element matrix at a time, generally.
>>>
>>>




 4. I would like to measure the time of each section of the
 process. Which method is recommended?

>>>
>>> PETSc methods are timed separateluy, but setup gets folded into KSPSolve
>>> unless you call SNESSetUp before the SNES[KSP]Solve.
>>> You can add you own timers also
>>> https://petsc.org/release/docs/manualpages/Profiling/PetscLogEventRegister/
>>>
>>> Mark
>>>
>>>
>>>


 Thank you for your help.



 Hyung Kim

>>>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> 
>

Re: [petsc-users] Question About Assembly matrix and declaration of KSP & pc

[Matthew Knepley]

On Wed, Nov 30, 2022 at 7:51 AM 김성익  wrote:

> Thank you for your comments.
> However, I have more questions.
>
> 1. Generally,  (KSPCreate, KSPSetOperators, KSPGetPC, PCSetType,
> PCFactorSetMatSolverType, KSPSetFromOptions )
> above functions are should be called after  each "MatassemblyEnd??"
>

KSPCreate is called once.

You do not need PCSetType, PCFactorSetMatSolverType, KSPSetFromOptions more
than once, unless you want to change the solver type.

KSPSetOperators is called if you want to change the system matrix.

KSPSolve is called when you want to change the rhs.

  Thanks,

 Matt


> 2.   Though reading the user guide, I don't fully understand under what
> circumstances the functions mentioned above should be called again. Can you
> explain when each function should be called?
>
> Thanks,
>
> Hyung Kim
>
> 2022년 11월 30일 (수) 오후 8:37, Mark Adams 님이 작성:
>
>>
>>
>> On Wed, Nov 30, 2022 at 5:08 AM 김성익  wrote:
>>
>>> Hello,
>>>
>>>
>>>
>>> I’m working on FEM using PETSc.
>>>
>>> As everyone knows, it is necessary to repeatedly solve Ax=B.
>>>
>>> Regarding this, I have 4 questions.
>>>
>>>
>>>
>>> 1. There are many steps for preparing KSPSolve. For example
>>> KSPcreate, KSPSetOperators, KSPGetPC, PCSetType, PCFactorSetMatSolverType,
>>> KSPSetFromOptions…
>>> In Nonlinear FEM, there are repeatedly kspsolve for getting answer
>>> vector. Is it correct to do all of the aforementioned processes (KSPcreate,
>>> KSPSetOperators ~~~) for each KSPSolve? Or should I declare it only once at
>>> the beginning and not call it again?
>>>
>>
>> You just do these once at setup but for nonlinear problems
>> KSPSetOperators tells the solver that you have a new matrix and so "matrix
>> setup" work needs to be done.
>>
>>
>>>
>>> 2. If the answer to question 1 is that it must be repeated every
>>> time, should this work be done right before kspsolve, that is, when the
>>> global matrix assembly is finished, or is it irrelevant to performance at
>>> any time?
>>>
>>
>> KSPSetOperators should be set after the new matrix values are set but it
>> might work before. It just sets a pointer to the matrix and flags it as not
>> setup.
>>
>>
>>>
>>>
>>> 3. When performing FEM, local matrices are often scattered in
>>> global matrices depending on connectivity. In this case, which is better in
>>> terms of performance: adding the values one by one with MatSetValue or
>>> adding them all at once with MatSetValues even if they are scattered?
>>>
>>
>> You want to add one element matrix at a time, generally.
>>
>>
>>>
>>>
>>>
>>>
>>> 4. I would like to measure the time of each section of the process.
>>> Which method is recommended?
>>>
>>
>> PETSc methods are timed separateluy, but setup gets folded into KSPSolve
>> unless you call SNESSetUp before the SNES[KSP]Solve.
>> You can add you own timers also
>> https://petsc.org/release/docs/manualpages/Profiling/PetscLogEventRegister/
>>
>> Mark
>>
>>
>>
>>>
>>>
>>> Thank you for your help.
>>>
>>>
>>>
>>> Hyung Kim
>>>
>>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Question About Assembly matrix and declaration of KSP & pc

[김성익]

Thank you for your comments.
However, I have more questions.

1. Generally,  (KSPCreate, KSPSetOperators, KSPGetPC, PCSetType,
PCFactorSetMatSolverType, KSPSetFromOptions )
above functions are should be called after  each "MatassemblyEnd??"

2.   Though reading the user guide, I don't fully understand under what
circumstances the functions mentioned above should be called again. Can you
explain when each function should be called?

Thanks,

Hyung Kim

2022년 11월 30일 (수) 오후 8:37, Mark Adams 님이 작성:

>
>
> On Wed, Nov 30, 2022 at 5:08 AM 김성익  wrote:
>
>> Hello,
>>
>>
>>
>> I’m working on FEM using PETSc.
>>
>> As everyone knows, it is necessary to repeatedly solve Ax=B.
>>
>> Regarding this, I have 4 questions.
>>
>>
>>
>> 1. There are many steps for preparing KSPSolve. For example
>> KSPcreate, KSPSetOperators, KSPGetPC, PCSetType, PCFactorSetMatSolverType,
>> KSPSetFromOptions…
>> In Nonlinear FEM, there are repeatedly kspsolve for getting answer
>> vector. Is it correct to do all of the aforementioned processes (KSPcreate,
>> KSPSetOperators ~~~) for each KSPSolve? Or should I declare it only once at
>> the beginning and not call it again?
>>
>
> You just do these once at setup but for nonlinear problems KSPSetOperators
> tells the solver that you have a new matrix and so "matrix setup" work
> needs to be done.
>
>
>>
>> 2. If the answer to question 1 is that it must be repeated every
>> time, should this work be done right before kspsolve, that is, when the
>> global matrix assembly is finished, or is it irrelevant to performance at
>> any time?
>>
>
> KSPSetOperators should be set after the new matrix values are set but it
> might work before. It just sets a pointer to the matrix and flags it as not
> setup.
>
>
>>
>>
>> 3. When performing FEM, local matrices are often scattered in global
>> matrices depending on connectivity. In this case, which is better in terms
>> of performance: adding the values one by one with MatSetValue or adding
>> them all at once with MatSetValues even if they are scattered?
>>
>
> You want to add one element matrix at a time, generally.
>
>
>>
>>
>>
>>
>> 4. I would like to measure the time of each section of the process.
>> Which method is recommended?
>>
>
> PETSc methods are timed separateluy, but setup gets folded into KSPSolve
> unless you call SNESSetUp before the SNES[KSP]Solve.
> You can add you own timers also
> https://petsc.org/release/docs/manualpages/Profiling/PetscLogEventRegister/
>
> Mark
>
>
>
>>
>>
>> Thank you for your help.
>>
>>
>>
>> Hyung Kim
>>
>