[petsc-users] using real and complex together
Sam Guo
sam.guo at cd-adapco.com
Mon Sep 28 19:06:22 CDT 2020
If I load complex version of PETSc, how can Vec be real?
On Monday, September 28, 2020, Matthew Knepley <knepley at gmail.com> wrote:
> 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/815fa4d4/attachment-0001.html>
More information about the petsc-users
mailing list