[petsc-users] MatPtAP

[Samuel Lanthaler]

On 06/01/2018 03:42 PM, Matthew Knepley wrote:
> On Fri, Jun 1, 2018 at 9:21 AM, Samuel Lanthaler 
> <s.lanthaler at gmail.com <mailto:s.lanthaler at gmail.com>> 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/ <http://www.caam.rice.edu/%7Emk51/>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180601/931de0b8/attachment.html>