On Mon, Sep 28, 2020 at 6:29 PM Sam Guo <sam....@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 <knep...@gmail.com> wrote: > >> On Mon, Sep 28, 2020 at 5:01 PM Sam Guo <sam....@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 <knep...@gmail.com> >>> wrote: >>> >>>> On Mon, Sep 28, 2020 at 3:43 PM Sam Guo <sam....@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....@cd-adapco.com> wrote: >>>>> >>>>>> Thanks >>>>>> >>>>>> On Tuesday, May 26, 2020, Stefano Zampini <stefano.zamp...@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....@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 <hzh...@mcs.anl.gov> wrote: >>>>>>> >>>>>>>> You can build PETSc with complex version, and declare some >>>>>>>> variables as 'PETSC_REAL'. >>>>>>>> Hong >>>>>>>> >>>>>>>> ------------------------------ >>>>>>>> *From:* petsc-users <petsc-users-boun...@mcs.anl.gov> on behalf of >>>>>>>> Sam Guo <sam....@cd-adapco.com> >>>>>>>> *Sent:* Tuesday, May 26, 2020 1:00 PM >>>>>>>> *To:* PETSc <petsc-users@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/>