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 <https://epubs.siam.org/doi/abs/10.1137/S1064827500372262?mobileUi=0> Randy M. > On Sep 28, 2020, at 5:12 PM, Sam Guo <sam....@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....@cd-adapco.com > <mailto:sam....@cd-adapco.com>> wrote: > If I load complex version of PETSc, how can Vec be real? > > On Monday, September 28, 2020, Matthew Knepley <knep...@gmail.com > <mailto:knep...@gmail.com>> wrote: > On Mon, Sep 28, 2020 at 7:44 PM Sam Guo <sam....@cd-adapco.com > <mailto:sam....@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 <knep...@gmail.com > <mailto:knep...@gmail.com>> wrote: > On Mon, Sep 28, 2020 at 6:29 PM Sam Guo <sam....@cd-adapco.com > <mailto: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 > <mailto:knep...@gmail.com>> wrote: > On Mon, Sep 28, 2020 at 5:01 PM Sam Guo <sam....@cd-adapco.com > <mailto: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 > <mailto:knep...@gmail.com>> wrote: > On Mon, Sep 28, 2020 at 3:43 PM Sam Guo <sam....@cd-adapco.com > <mailto: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 > <mailto:sam....@cd-adapco.com>> wrote: > Thanks > > On Tuesday, May 26, 2020, Stefano Zampini <stefano.zamp...@gmail.com > <mailto: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 >> <mailto: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 >> <mailto: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 >> <mailto:petsc-users-boun...@mcs.anl.gov>> on behalf of Sam Guo >> <sam....@cd-adapco.com <mailto:sam....@cd-adapco.com>> >> Sent: Tuesday, May 26, 2020 1:00 PM >> To: PETSc <petsc-users@mcs.anl.gov <mailto: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/> > > > -- > 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/>
