<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 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 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 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 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 href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~<wbr>knepley/</a><br>
> <br>
<br>
</div></div></blockquote></div><br></div></div></div></div>