Barry,<div><br></div><div>Thank you so much for getting back to me so quickly. Just out of curiosity, is there a slightly more in depth version of the documentation that I would be able to look at if I have some other "simple" questions like this in the future? </div>
<div> </div><div><br><div class="gmail_quote">On Wed, May 27, 2009 at 5:59 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<br>
Paul,<br>
<br>
In PETSc for the AIJ sparse matrix storage format we store the inverse of the diagonal entries of U. This is because then when it comes to the triangular solves they only involve floating point multiplies and adds (no time consuming divisions). I think this is pretty common for sparse matrix solvers. You could argue that the MatView() should invert those values so what gets printed out is clearer; I view the MatView() for factored matrices as only a useful tool for debugging for developers so it just prints out the stored values. Yes it is a bit confusing.<br>
<br>
Barry<br>
<br>
For BAIJ matrices it gets even more confusing. With BAIJ the factorizations use the little blocks of the matrix as the basic computational unit of the factorization. Here we actually compute the exact inverse of the diagonal block and store that (a block version of what we do for AIJ). Again this is to make the solves as efficient as possible.<div>
<div></div><div class="h5"><br>
<br>
<br>
<br>
On May 27, 2009, at 6:20 PM, Paul Dostert wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
I'm a beginner with PETSc, so please forgive me if this is obvious, but I<br>
couldn't seem to find any help in the archives.<br>
<br>
I'm trying to just the hang of the software, so I've been messing around<br>
with routines. I'm going to need complex matrices (Maxwell's with PML) so<br>
everything is configured for this. I'm messing around with some very simple<br>
test cases, and have a symmetric (but not Hermitian) complex matrix with<br>
2-2i on the diagonal and -1+i on the upper and lower diagonals. I am reading<br>
this in from a petsc binary file (again, for testing purposes, eventually<br>
I'm going to be just reading in my matrix and RHS). I view the matrix, and<br>
it has been read in correctly. I perform LU factorization by doing the<br>
following (where ISCreateStride(PETSC_COMM_WORLD,m,0,1,&perm); has been<br>
called earlier):<br>
<br>
MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&LU);<br>
MatFactorInfoInitialize(&luinfo);<br>
MatLUFactor(LU,perm,perm,&luinfo);<br>
<br>
I get that LU is:<br>
<br>
(1,1) 0.2500 + 0.2500i<br>
(2,1) -0.5000<br>
(1,2) -1.0000 + 1.0000i<br>
(2,2) 0.3333 + 0.3333i<br>
(3,2) -0.6667<br>
(2,3) -1.0000 + 1.0000i<br>
(3,3) 0.3750 + 0.3750i<br>
(4,3) -0.7500<br>
(3,4) -1.0000 + 1.0000i<br>
(4,4) 0.4000 + 0.4000i<br>
(5,4) -0.8000<br>
(4,5) -1.0000 + 1.0000i<br>
(5,5) 0.4167 + 0.4167i<br>
(6,5) -0.8333<br>
(5,6) -1.0000 + 1.0000i<br>
(6,6) 0.4286 + 0.4286i<br>
<br>
Now, I am interpreting this as L being unit on the diagonal and the<br>
lower diagonal portion of this "LU" matrix, while U being the diagona<br>
+ upper of this "LU" matrix. I can interpret this the other way around<br>
as well, and it doesn't matter.<br>
<br>
However, knowing the LU factorization, it is VERY clear the the proper<br>
LU decomposition would have the inverse of the diagonal elements<br>
presented here. So I believe I should have LU as:<br>
<br>
(1,1) 2.0000 - 2.0000i<br>
(2,1) -0.5000<br>
(1,2) -1.0000 + 1.0000i<br>
(2,2) 1.5000 - 1.5000i<br>
(3,2) -0.6667<br>
(2,3) -1.0000 + 1.0000i<br>
(3,3) 1.3333 - 1.3333i<br>
(4,3) -0.7500<br>
(3,4) -1.0000 + 1.0000i<br>
(4,4) 1.2500 - 1.2500i<br>
(5,4) -0.8000<br>
(4,5) -1.0000 + 1.0000i<br>
(5,5) 1.2000 - 1.2000i<br>
(6,5) -0.8333<br>
(5,6) -1.0000 + 1.0000i<br>
(6,6) 1.1667 - 1.1667i<br>
<br>
Is there some reason this returns the inverse of the diagonal entries,<br>
or am I completely missing something? Is returning the inverse<br>
something standard??<br>
<br>
Since I'm new, I'm not quite sure where to look for actual source<br>
code. Is there a location where the LU factorization code is written<br>
and well commented?<br>
<br>
Thank you very much!<br>
</blockquote>
<br>
</div></div></blockquote></div><br></div>