[petsc-users] using real and complex together

Sam Guo sam.guo at cd-adapco.com
Mon Sep 28 19:34:09 CDT 2020 

That sounds like an useful feature for future PETSc.

On Monday, September 28, 2020, Matthew Knepley <knepley at gmail.com> wrote:

> 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/84419227/attachment.html>

More information about the petsc-users mailing list