<html>
<head>
<meta content="text/html; charset=utf-8" http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
Ok, I now realize that I had implemented the boundary conditions in
an unnecessarily complicated way... As you pointed out, I can just
manipulate individual matrix rows to enforce the BC's. In that way,
I never have to call MatPtAP, or do any expensive operations.
Probably that's what is commonly done. I've changed my code, and it
seems to work fine and is much faster.<br>
<br>
Thanks a lot for your help, everyone! This has been very educational
for me.<br>
<br>
Samuel
<p><br>
</p>
<div class="moz-cite-prefix">On 06/01/2018 09:04 PM, Hong wrote:<br>
</div>
<blockquote
cite="mid:CAGCphBudV_MjJqraDAvD8pqz9a9hWamF4DD_h9pQhN=3iaGaBg@mail.gmail.com"
type="cite">
<div dir="ltr"><span
style="color:rgb(80,0,80);font-family:arial,sans-serif;font-size:small;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">Samuel,</span>
<div><font color="#500050">I have following questions:</font></div>
<div>
<div><font color="#500050">1) Why solving <span
style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:12.8px;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">\lambda*Pt*B*P*y=Pt*A*P*y
is better than solving original <span
style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:12.8px;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">\lambda*B*x
= A*x?</span></span></font></div>
<div><font color="#500050"><span
style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:12.8px;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline"><span
style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:12.8px;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline"><br>
</span></span></font></div>
<div><span
style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:12.8px;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline"></span><span
style="font-size:12.8px">2) Does your eigen solver require
matrix factorization? If not, i.e., only uses mat-vec
multiplication, then you may implement</span></div>
<div><span style="font-size:12.8px"><span
style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:12.8px;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">z
= (Pt*(A*(P*y))) in your eigensolver instead of using
mat-mat-mat multi<span
style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:12.8px;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">plication.</span></span></span></div>
<div><span style="font-size:12.8px"><br>
</span></div>
<div><span style="font-size:12.8px"><span
style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:12.8px;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline"></span>3)
petsc PtAP() was implemented for multigrid applications,
in which the product C = <span
style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:12.8px;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">PtAP
is a denser but a much smaller matrix.</span></span></div>
<div><span style="font-size:12.8px">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.</span></div>
<div><span style="font-size:12.8px"><br>
</span></div>
<div><span style="font-size:12.8px">Hong<br>
</span>
<div class="gmail_extra"><br>
<div class="gmail_quote">On Fri, Jun 1, 2018 at 12:35 PM,
Smith, Barry F. <span dir="ltr"><<a
moz-do-not-send="true"
href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span>
wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px
0px 0.8ex;border-left:1px solid
rgb(204,204,204);padding-left:1ex"><br>
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 <a
moz-do-not-send="true"
href="mailto:petsc-maint@mcs.anl.gov">petsc-maint@mcs.anl.gov</a>
? Then we can run the code in the profiler with your
matrices and see if there is any way to speed up the
computation. <br>
<span class="HOEnZb"><font color="#888888"><br>
Barry<br>
</font></span>
<div class="HOEnZb">
<div class="h5"><br>
<br>
> On Jun 1, 2018, at 11:07 AM, Samuel Lanthaler
<<a moz-do-not-send="true"
href="mailto:s.lanthaler@gmail.com">s.lanthaler@gmail.com</a>>
wrote:<br>
> <br>
> On 06/01/2018 03:42 PM, Matthew Knepley
wrote:<br>
>> On Fri, Jun 1, 2018 at 9:21 AM, Samuel
Lanthaler <<a moz-do-not-send="true"
href="mailto:s.lanthaler@gmail.com">s.lanthaler@gmail.com</a>>
wrote:<br>
>> Hi,<br>
>> <br>
>> 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<br>
>> <br>
>> P = [ B, 0, 0, 0, ..., 0, 0 ]<br>
>> [ 0, 1, 0, 0, ..., 0, 0 ]<br>
>> [ ]<br>
>> [ ]<br>
>> [ ..., 1, 0 ]<br>
>> [ 0, 0, 0, 0, ..., 0, C ]<br>
>> <br>
>> with B,C are (k-by-k) matrices.<br>
>> Right now, I'm simply constructing A, P
and calling<br>
>> <br>
>> CALL MatPtAP(petsc_matA,petsc_matP,<wbr>MAT_INITIAL_MATRIX,PETSC_<wbr>DEFAULT_REAL,petsc_matPtAP,<wbr>ierr)<br>
>> <br>
>> where I haven't done anything to
pestc_matPtAP, prior to this call. Is this the way
to do it?<br>
>> <br>
>> 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.<br>
>> <br>
>> Are you sure this is what you want to do?
Usually BC are local, since by definition PDE are
local, and<br>
>> are applied pointwise. What kind of BC do
you have here?<br>
>> <br>
> <br>
> 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. <br>
> <br>
> In reality, I have two matrices A,B, and want
to solve the eigenvalue problem \lambda*B*x = A*x.
<br>
> 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.<br>
> Is this an inherently bad/inefficient way of
enforcing BC's? Thanks.<br>
> <br>
> <br>
> <br>
> <br>
>> Thanks,<br>
>> <br>
>> Matt<br>
>> <br>
>> Thanks in advance for any ideas.<br>
>> <br>
>> Cheers,<br>
>> Samuel<br>
>> <br>
>> <br>
>> <br>
>> -- <br>
>> What most experimenters take for granted
before they begin their experiments is infinitely
more interesting than any results to which their
experiments lead.<br>
>> -- Norbert Wiener<br>
>> <br>
>> <a moz-do-not-send="true"
href="https://www.cse.buffalo.edu/%7Eknepley/"
rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~<wbr>knepley/</a><br>
> <br>
<br>
</div>
</div>
</blockquote>
</div>
<br>
</div>
</div>
</div>
</div>
</blockquote>
<br>
</body>
</html>