Barry, 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?
On Wed, May 27, 2009 at 5:59 PM, Barry Smith <bsmith at mcs.anl.gov> wrote: > > Paul, > > 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. > > Barry > > 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. > > > > > On May 27, 2009, at 6:20 PM, Paul Dostert wrote: > > I'm a beginner with PETSc, so please forgive me if this is obvious, but I >> couldn't seem to find any help in the archives. >> >> I'm trying to just the hang of the software, so I've been messing around >> with routines. I'm going to need complex matrices (Maxwell's with PML) so >> everything is configured for this. I'm messing around with some very >> simple >> test cases, and have a symmetric (but not Hermitian) complex matrix with >> 2-2i on the diagonal and -1+i on the upper and lower diagonals. I am >> reading >> this in from a petsc binary file (again, for testing purposes, eventually >> I'm going to be just reading in my matrix and RHS). I view the matrix, and >> it has been read in correctly. I perform LU factorization by doing the >> following (where ISCreateStride(PETSC_COMM_WORLD,m,0,1,&perm); has been >> called earlier): >> >> MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&LU); >> MatFactorInfoInitialize(&luinfo); >> MatLUFactor(LU,perm,perm,&luinfo); >> >> I get that LU is: >> >> (1,1) 0.2500 + 0.2500i >> (2,1) -0.5000 >> (1,2) -1.0000 + 1.0000i >> (2,2) 0.3333 + 0.3333i >> (3,2) -0.6667 >> (2,3) -1.0000 + 1.0000i >> (3,3) 0.3750 + 0.3750i >> (4,3) -0.7500 >> (3,4) -1.0000 + 1.0000i >> (4,4) 0.4000 + 0.4000i >> (5,4) -0.8000 >> (4,5) -1.0000 + 1.0000i >> (5,5) 0.4167 + 0.4167i >> (6,5) -0.8333 >> (5,6) -1.0000 + 1.0000i >> (6,6) 0.4286 + 0.4286i >> >> Now, I am interpreting this as L being unit on the diagonal and the >> lower diagonal portion of this "LU" matrix, while U being the diagona >> + upper of this "LU" matrix. I can interpret this the other way around >> as well, and it doesn't matter. >> >> However, knowing the LU factorization, it is VERY clear the the proper >> LU decomposition would have the inverse of the diagonal elements >> presented here. So I believe I should have LU as: >> >> (1,1) 2.0000 - 2.0000i >> (2,1) -0.5000 >> (1,2) -1.0000 + 1.0000i >> (2,2) 1.5000 - 1.5000i >> (3,2) -0.6667 >> (2,3) -1.0000 + 1.0000i >> (3,3) 1.3333 - 1.3333i >> (4,3) -0.7500 >> (3,4) -1.0000 + 1.0000i >> (4,4) 1.2500 - 1.2500i >> (5,4) -0.8000 >> (4,5) -1.0000 + 1.0000i >> (5,5) 1.2000 - 1.2000i >> (6,5) -0.8333 >> (5,6) -1.0000 + 1.0000i >> (6,6) 1.1667 - 1.1667i >> >> Is there some reason this returns the inverse of the diagonal entries, >> or am I completely missing something? Is returning the inverse >> something standard?? >> >> Since I'm new, I'm not quite sure where to look for actual source >> code. Is there a location where the LU factorization code is written >> and well commented? >> >> Thank you very much! >> > > -------------- next part -------------- An HTML attachment was scrubbed... URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20090527/ae96e33e/attachment.htm>
