Samuel, I have following questions: 1) Why solving \lambda*Pt*B*P*y=Pt*A*P*y is better than solving original \lambda*B*x = A*x?
2) Does your eigen solver require matrix factorization? If not, i.e., only uses mat-vec multiplication, then you may implement z = (Pt*(A*(P*y))) in your eigensolver instead of using mat-mat-mat multi plication. 3) petsc PtAP() was implemented for multigrid applications, in which the product C = PtAP is a denser but a much smaller matrix. I have not seen the use of your case, that P is square with same size as A. If C is much denser than A, then PtAP consumes a large portion of time is anticipated. Hong On Fri, Jun 1, 2018 at 12:35 PM, Smith, Barry F. <[email protected]> wrote: > > Could you save the P matrix with MatView() using a binary viewer and the > A matrix with MatView() and the binary viewer and email them to > [email protected] ? Then we can run the code in the profiler with > your matrices and see if there is any way to speed up the computation. > > Barry > > > > On Jun 1, 2018, at 11:07 AM, Samuel Lanthaler <[email protected]> > wrote: > > > > On 06/01/2018 03:42 PM, Matthew Knepley wrote: > >> On Fri, Jun 1, 2018 at 9:21 AM, Samuel Lanthaler <[email protected]> > wrote: > >> Hi, > >> > >> I was wondering what the most efficient way to use MatPtAP would be in > the following situation: I am discretizing a PDE system. The discretization > yields a matrix A that has a band structure (with k upper and lower bands, > say). In order to implement the boundary conditions, I use a transformation > matrix P which is essentially the unit matrix, except for the entries > P_{ij} where i,j<k and n-i,n-j<k, so > >> > >> P = [ B, 0, 0, 0, ..., 0, 0 ] > >> [ 0, 1, 0, 0, ..., 0, 0 ] > >> [ ] > >> [ ] > >> [ ..., 1, 0 ] > >> [ 0, 0, 0, 0, ..., 0, C ] > >> > >> with B,C are (k-by-k) matrices. > >> Right now, I'm simply constructing A, P and calling > >> > >> CALL MatPtAP(petsc_matA,petsc_matP,MAT_INITIAL_MATRIX,PETSC_ > DEFAULT_REAL,petsc_matPtAP,ierr) > >> > >> where I haven't done anything to pestc_matPtAP, prior to this call. Is > this the way to do it? > >> > >> I'm asking because, currently, setting up the matrices A and P takes > very little time, whereas the operation MatPtAP is taking quite long, which > seems very odd... The matrices are of type MPIAIJ. In my problem, the total > matrix dimension is around 10'000 and the matrix blocks (B,C) are of size > ~100. > >> > >> Are you sure this is what you want to do? Usually BC are local, since > by definition PDE are local, and > >> are applied pointwise. What kind of BC do you have here? > >> > > > > The boundary conditions are a mixture of Dirichlet and Neumann; in my > case, the PDE is a system involving 8 variables on a disk, where the > periodic direction is discretized using a Fourier series expansion, the > radial direction uses B-splines. > > > > In reality, I have two matrices A,B, and want to solve the eigenvalue > problem \lambda*B*x = A*x. > > I found it quite convenient to use a transformation P to a different set > of variables y, such that x=P*y and x satisfies the BC iff certain > components of y are 0. The latter is enforced by inserting spurious > eigenvalues at the relevant components of y in the transformed eigenvalue > problem \lambda*Pt*B*P*y=Pt*A*P*y. After solving the EVP in terms of y, I > get back x=P*y. > > Is this an inherently bad/inefficient way of enforcing BC's? Thanks. > > > > > > > > > >> Thanks, > >> > >> Matt > >> > >> Thanks in advance for any ideas. > >> > >> Cheers, > >> Samuel > >> > >> > >> > >> -- > >> 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/ > > > >
