[petsc-users] using real and complex together
Matthew Knepley
knepley at gmail.com
Mon Sep 28 19:21:42 CDT 2020
On Mon, Sep 28, 2020 at 8:20 PM Sam Guo <sam.guo at cd-adapco.com> wrote:
> Hi Randy and Matt,
> Thanks a lot. I’ll look into it.
>
Another option, with less code for you but some code in PETSc, is to check
for a matrix with 0 imaginary part, and then
copy it to a real matrix and call real MUMPS. However, in this case, we
would have to check for real solves as well.
Thanks,
Matt
> BR,
> Sam
>
> On Monday, September 28, 2020, Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Mon, Sep 28, 2020 at 8:15 PM Randall Mackie <rlmackie862 at gmail.com>
>> wrote:
>>
>>> Sam, you can solve a complex matrix using a real version of PETSc by
>>> doubling the size of your matrix and spitting out real/imaginary parts.
>>>
>>> See this paper:
>>>
>>> https://epubs.siam.org/doi/abs/10.1137/S1064827500372262?mobileUi=0
>>>
>>
>> Thanks Randy. Yes, I meant that one.
>>
>> Matt
>>
>>
>>> Randy M.
>>>
>>>
>>> On Sep 28, 2020, at 5:12 PM, Sam Guo <sam.guo at cd-adapco.com> wrote:
>>>
>>> All I want is to solve both real and complex matrix using either real or
>>> complex version of PETSc. If I load complex version of PETSc, I waste
>>> memory solving real matrix. If I can solve complex matrix using real
>>> version of PETSc, I will accept it.
>>>
>>> On Monday, September 28, 2020, Sam Guo <sam.guo at cd-adapco.com> wrote:
>>>
>>>> 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/>
>>>>>
>>>>
>>>
>>
>> --
>> 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/76706037/attachment-0001.html>
More information about the petsc-users
mailing list