[petsc-users] using real and complex together
Matthew Knepley
knepley at gmail.com
Mon Sep 28 19:01:17 CDT 2020
On Mon, Sep 28, 2020 at 7:44 PM Sam Guo <sam.guo at cd-adapco.com> wrote:
> I want to make sure I understand you. I load real version of PETSc. If my
> input matrix is complex,
You said that your matrix was real.
> just get real and imagine parts of PETSc Vec
No, the PETSc Vec would be real. You would have two vectors.
Thanks,
Matt
> and do the matrix vector multiplication. Right?
>
> On Monday, September 28, 2020, Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Mon, Sep 28, 2020 at 6:29 PM Sam Guo <sam.guo at cd-adapco.com> wrote:
>>
>>> “ I think it would be much easier to just decompose your complex work
>>> into real and imaginary parts and use PETSc with real scalars to compute
>>> them separately.
>>> Since you know your matrices have 0 imaginary part, this becomes very
>>> straightforward.”
>>>
>>> I think this is what I am trying to do. But since I have to provide
>>> matrix-vector operator(since I am using shell matrix), the Vec I receive is
>>> complex. I need to convert complex vec to real one and then convert it
>>> back(that’s the code I shown before).
>>>
>>
>> No, this is not what I am advocating. Keep your vectors real, just keep
>> one for the real part and one for the imaginary part. Then you can just
>> call MatMult() twice with your matrix.
>>
>> Thanks,
>>
>> Matt
>>
>>
>>> On Monday, September 28, 2020, Matthew Knepley <knepley at gmail.com>
>>> wrote:
>>>
>>>> On Mon, Sep 28, 2020 at 5:01 PM Sam Guo <sam.guo at cd-adapco.com> wrote:
>>>>
>>>>> Hi Matt,
>>>>> Since I use MUMPS as preconditioner, complex uses too much memory
>>>>> if my input matrix is real. Ideally if I can compile real and complex into
>>>>> different symbols (like MUMPS) , I can load both version without conflict.
>>>>>
>>>>
>>>> What I mean to say is that it would be great if it were as simple as
>>>> using two different symbols, but unfortunately the problem is more
>>>> difficult. I was trying to use
>>>> the example of templates. This would be a very intrusive change no
>>>> matter what technology you are using.
>>>>
>>>> So your main memory usage is from the MUMPS factorization, and you
>>>> cannot afford to double that usage?
>>>>
>>>> You could consider writing a version of AIJ that stores real entries,
>>>> but allows complex vector values. It would promote to complex for the row
>>>> dot product.
>>>> However, you would also have to do the same work for all the solves you
>>>> do with MUMPS.
>>>>
>>>> I think it would be much easier to just decompose your complex work
>>>> into real and imaginary parts and use PETSc with real scalars to compute
>>>> them separately.
>>>> Since you know your matrices have 0 imaginary part, this becomes very
>>>> straightforward.
>>>>
>>>> Thanks,
>>>>
>>>> Matt
>>>>
>>>>
>>>>> Thanks,
>>>>> Sam
>>>>>
>>>>> On Mon, Sep 28, 2020 at 12:52 PM Matthew Knepley <knepley at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> On Mon, Sep 28, 2020 at 3:43 PM Sam Guo <sam.guo at cd-adapco.com>
>>>>>> wrote:
>>>>>>
>>>>>>> Hi Stefano and PETSc dev team,
>>>>>>> I want to try your suggestion to always load complex version of
>>>>>>> PETSc but if my input matrix A is real, I want to create shell matrix to
>>>>>>> matrix-vector and factorization using real only.
>>>>>>>
>>>>>>
>>>>>> I do not think that will work as you expect. I will try to explain
>>>>>> below.
>>>>>>
>>>>>>
>>>>>>> I still need to understand how MatRealPart works. Does it just
>>>>>>> zero out the image numerical values or does it delete the image memory?
>>>>>>>
>>>>>>
>>>>>> When we have complex values, we use the "complex" type to allocate
>>>>>> and store them. Thus you cannot talk about just the memory to store
>>>>>> imaginary parts.
>>>>>> MatRealPart sets the imaginary parts of all the matrix elements to
>>>>>> zero.
>>>>>>
>>>>>>
>>>>>>> If my input matrix A is real, how do I create a shell matrix to
>>>>>>> matrix -vector multiplication y=A*x where A is real, PestcScalar = complex,
>>>>>>> x and y are Vec? I notice there is a VecRealPart but it seems it just zeros
>>>>>>> the image numerical values. It seems I still have to create a PetscReal
>>>>>>> pointer to copy the real part of PetacScalar pointers like following. Can
>>>>>>> you comment on it?
>>>>>>>
>>>>>>
>>>>>> What you suggest would mean rewriting the matrix multiplication
>>>>>> algorithm by hand after extracting the values. I am not sure if this
>>>>>> is really what you want to do. Is the matrix memory really your
>>>>>> limiting factor? Even if you tried to do this with templates, the memory
>>>>>> from temporaries would be very hard to control.
>>>>>>
>>>>>> Thanks,
>>>>>>
>>>>>> Matt
>>>>>>
>>>>>>
>>>>>>> Thanks,
>>>>>>> Sam
>>>>>>>
>>>>>>> PetscScalar *px = nullptr;
>>>>>>> VecGetArrayRead(x, &px);
>>>>>>> PetscScalar *py = nullptr;
>>>>>>> VecGetArray(y, &py);
>>>>>>> int localSize = 0;
>>>>>>> VecGetLocalSize(x, &localSize);
>>>>>>> std::vector<PetasReal> realX(localSize); // I am using c++ to call
>>>>>>> PETSc
>>>>>>>
>>>>>>> //retrieve real part
>>>>>>> for(int i = 0; i < localSize; i++) realX[i] = PetscRealPart(px[i]);
>>>>>>>
>>>>>>> // do real matrix-vector multiplication
>>>>>>> // realY=A*realX
>>>>>>> // here where realY is std::vector<PetscReal>
>>>>>>>
>>>>>>> //put real part back to py
>>>>>>> for(int i = 0; i < localSize; i++) pv[i] = realY[i];
>>>>>>> VecRestoreArray(y,&py);
>>>>>>>
>>>>>>> On Tue, May 26, 2020 at 1:49 PM Sam Guo <sam.guo at cd-adapco.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> Thanks
>>>>>>>>
>>>>>>>> On Tuesday, May 26, 2020, Stefano Zampini <
>>>>>>>> stefano.zampini at gmail.com> wrote:
>>>>>>>>
>>>>>>>>> All the solvers/matrices/vectors works for PetscScalar types (i.e.
>>>>>>>>> in your case complex)
>>>>>>>>> If you need to solve for the real part only, you can duplicate the
>>>>>>>>> matrix and call MatRealPart to zero out the imaginary part. But the solve
>>>>>>>>> will always run in the complex space
>>>>>>>>> You should not be worried about doubling the memory for a matrix
>>>>>>>>> (i.e. real and imaginary part)
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> On May 26, 2020, at 11:28 PM, Sam Guo <sam.guo at cd-adapco.com>
>>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>> complex version is needed since matrix sometimes is real and
>>>>>>>>> sometimes is complex. I want to solve real matrix without allocating memory
>>>>>>>>> for imaginary part((except eigen pairs).
>>>>>>>>>
>>>>>>>>> On Tuesday, May 26, 2020, Zhang, Hong <hzhang at mcs.anl.gov> wrote:
>>>>>>>>>
>>>>>>>>>> You can build PETSc with complex version, and declare some
>>>>>>>>>> variables as 'PETSC_REAL'.
>>>>>>>>>> Hong
>>>>>>>>>>
>>>>>>>>>> ------------------------------
>>>>>>>>>> *From:* petsc-users <petsc-users-bounces at mcs.anl.gov> on behalf
>>>>>>>>>> of Sam Guo <sam.guo at cd-adapco.com>
>>>>>>>>>> *Sent:* Tuesday, May 26, 2020 1:00 PM
>>>>>>>>>> *To:* PETSc <petsc-users at mcs.anl.gov>
>>>>>>>>>> *Subject:* [petsc-users] using real and complex together
>>>>>>>>>>
>>>>>>>>>> Dear PETSc dev team,
>>>>>>>>>> Can I use both real and complex versions together?
>>>>>>>>>>
>>>>>>>>>> Thanks,
>>>>>>>>>> Sam
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>
>>>>>> --
>>>>>> 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/>
>>
>
--
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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200928/b34270ec/attachment.html>
More information about the petsc-users
mailing list