That sounds like an useful feature for future PETSc. On Monday, September 28, 2020, Matthew Knepley <[email protected]> wrote:
> On Mon, Sep 28, 2020 at 8:20 PM Sam Guo <[email protected]> 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 <[email protected]> wrote: >> >>> On Mon, Sep 28, 2020 at 8:15 PM Randall Mackie <[email protected]> >>> 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 <[email protected]> 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 <[email protected]> wrote: >>>> >>>>> If I load complex version of PETSc, how can Vec be real? >>>>> >>>>> On Monday, September 28, 2020, Matthew Knepley <[email protected]> >>>>> wrote: >>>>> >>>>>> On Mon, Sep 28, 2020 at 7:44 PM Sam Guo <[email protected]> >>>>>> 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 <[email protected]> >>>>>>> wrote: >>>>>>> >>>>>>>> On Mon, Sep 28, 2020 at 6:29 PM Sam Guo <[email protected]> >>>>>>>> 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 <[email protected]> >>>>>>>>> wrote: >>>>>>>>> >>>>>>>>>> On Mon, Sep 28, 2020 at 5:01 PM Sam Guo <[email protected]> >>>>>>>>>> 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 < >>>>>>>>>>> [email protected]> wrote: >>>>>>>>>>> >>>>>>>>>>>> On Mon, Sep 28, 2020 at 3:43 PM Sam Guo <[email protected]> >>>>>>>>>>>> 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 <[email protected]> >>>>>>>>>>>>> wrote: >>>>>>>>>>>>> >>>>>>>>>>>>>> Thanks >>>>>>>>>>>>>> >>>>>>>>>>>>>> On Tuesday, May 26, 2020, Stefano Zampini < >>>>>>>>>>>>>> [email protected]> 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 <[email protected]> >>>>>>>>>>>>>>> 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 <[email protected]> >>>>>>>>>>>>>>> wrote: >>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> You can build PETSc with complex version, and declare some >>>>>>>>>>>>>>>> variables as 'PETSC_REAL'. >>>>>>>>>>>>>>>> Hong >>>>>>>>>>>>>>>> >>>>>>>>>>>>>>>> ------------------------------ >>>>>>>>>>>>>>>> *From:* petsc-users <[email protected]> on >>>>>>>>>>>>>>>> behalf of Sam Guo <[email protected]> >>>>>>>>>>>>>>>> *Sent:* Tuesday, May 26, 2020 1:00 PM >>>>>>>>>>>>>>>> *To:* PETSc <[email protected]> >>>>>>>>>>>>>>>> *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/> >
