[petsc-users] using real and complex together
Matthew Knepley
knepley at gmail.com
Mon Sep 28 18:07:41 CDT 2020
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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200928/dee2f3af/attachment.html>
More information about the petsc-users
mailing list