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/ >
