I do not want the transpose, but that is what is coming out as the solution. Both my tests show the correct original matrix, and both yield the transpose of the solution verified against Mathematica and the website from which I pulled the 4x4 matrix in test.cpp. Could it be related to the factoring?<br>
<br><div class="gmail_quote">On Tue, Jul 5, 2011 at 4:39 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>
MatSolve() and MatMatSolve() will not return the transpose of the solution, they will return the solution. If you want the transpose solve you can use MatSolveTranspose() for a single Vec or use MatMatSolve() for the matrix case after transposing the original matrix if you want the transpose.<br>
<font color="#888888"><br>
Barry<br>
</font><div><div></div><div class="h5"><br>
<br>
On Jul 5, 2011, at 12:17 PM, Adam Byrd wrote:<br>
<br>
> I have working code that produces the correct answer now, thank you. One (hopefully) final question, though. The solution is actually transposed. What causes this? Presumably I can use MatMatSolveTranspose to get around this, but there's no man page and I want to be sure of what will happen in all cases.<br>
><br>
> Respectfully,<br>
> Adam<br>
><br>
> On Thu, Jun 30, 2011 at 4:18 PM, Hong Zhang <<a href="mailto:hzhang@mcs.anl.gov">hzhang@mcs.anl.gov</a>> wrote:<br>
> Add<br>
> ierr = MatGetOrdering(testMat,MATORDERING_ND,&isrow,&iscol);CHKERRQ(ierr);<br>
> before the line<br>
> ierr = MatFactorInfoInitialize(&luinfo);CHKERRQ(ierr);<br>
><br>
> Somehow, your matrix is numerically singular. With<br>
> MATORDERING_NATURAL<br>
> I get<br>
> [0]PETSC ERROR: Detected zero pivot in LU factorization<br>
><br>
> even using MATDENSE matrix format, which calls lapack.<br>
> With MATORDERING_ND, I get useless inverseMat.<br>
><br>
> The modified code I used is attached.<br>
><br>
> Hong<br>
><br>
> On Thu, Jun 30, 2011 at 2:30 PM, Adam Byrd <<a href="mailto:adam1.byrd@gmail.com">adam1.byrd@gmail.com</a>> wrote:<br>
> > I'm trying to work through what I need to do, again by practicing with a<br>
> > small scale random problem. The general order of events seems to be: create<br>
> > a matrix, fill it, assemble it, factor it, then one can use solvers with it.<br>
> > When I use MatLUFactor on my matrix before using it with a solver I get this<br>
> > error:<br>
> > [0]PETSC ERROR: --------------------- Error Message<br>
> > ------------------------------------<br>
> > [0]PETSC ERROR: Null argument, when expecting valid pointer!<br>
> > [0]PETSC ERROR: Null Object: Parameter # 1!<br>
> > [0]PETSC ERROR:<br>
> > ------------------------------------------------------------------------<br>
> > [0]PETSC ERROR: Petsc Release Version 3.1.0, Patch 8, Thu Mar 17 13:37:48<br>
> > CDT 2011<br>
> > [0]PETSC ERROR: See docs/changes/index.html for recent updates.<br>
> > [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.<br>
> > [0]PETSC ERROR: See docs/index.html for manual pages.<br>
> > [0]PETSC ERROR:<br>
> > ------------------------------------------------------------------------<br>
> > [0]PETSC ERROR: ./test on a osx-gnu-c named Macintosh-3.local by adambyrd<br>
> > Thu Jun 30 15:27:30 2011<br>
> > [0]PETSC ERROR: Libraries linked from<br>
> > /Users/adambyrd/soft/petsc-3.1-p8/osx-gnu-cpp/lib<br>
> > [0]PETSC ERROR: Configure run at Tue Jun 28 12:56:55 2011<br>
> > [0]PETSC ERROR: Configure options PETSC_ARCH=osx-gnu-cpp --with-fc=gfortran<br>
> > -download-f-blas-lapack=1 --download-mpich=1 --with-scalar-type=complex<br>
> > --with-clanguage=c++<br>
> > [0]PETSC ERROR:<br>
> > ------------------------------------------------------------------------<br>
> > [0]PETSC ERROR: ISInvertPermutation() line 209 in<br>
> > src/vec/is/interface/index.c<br>
> > [0]PETSC ERROR: MatLUFactorSymbolic_SeqAIJ() line 306 in<br>
> > src/mat/impls/aij/seq/aijfact.c<br>
> > [0]PETSC ERROR: MatLUFactorSymbolic() line 2534 in<br>
> > src/mat/interface/matrix.c<br>
> > [0]PETSC ERROR: MatLUFactor_SeqAIJ() line 945 in<br>
> > src/mat/impls/aij/seq/aijfact.c<br>
> > [0]PETSC ERROR: MatLUFactor() line 2417 in src/mat/interface/matrix.c<br>
> > [0]PETSC ERROR: main() line 62 in WDtest.cpp<br>
> > application called MPI_Abort(MPI_COMM_WORLD, 85) - process 0[unset]:<br>
> > aborting job:<br>
> > application called MPI_Abort(MPI_COMM_WORLD, 85) - process 0<br>
> > I don't understand what I'm doing wrong.<br>
> > Respectfully,<br>
> > Adam<br>
> > On Wed, Jun 29, 2011 at 1:26 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
> >><br>
> >> On Tue, Jun 28, 2011 at 10:20 PM, Adam Byrd <<a href="mailto:adam1.byrd@gmail.com">adam1.byrd@gmail.com</a>> wrote:<br>
> >>><br>
> >>> Matt,<br>
> >>><br>
> >>> Alright, that means I need to continue learning how to use<br>
> >>> MatSetValues(). With my 6x6 example I tried filling it with four 3x3 sub<br>
> >>> matrices, but when I do that I get the error 'sum of local sizes 12 does not<br>
> >>> equal global size.' I had 4 processors each calling MatSetValues for their<br>
> >>> own 3x3. Graphically, I arranged the nodes 0 1<br>
> >>><br>
> >>> 2 3<br>
> >>> where process 0 had global rows 0-2 and global columns 0-2; process 1 had<br>
> >>> 0-2, 3-5; process 2 had 3-5, 0-2; and process 3 had 3-5, 3-5. >From the<br>
> >>> documentation, I think this should be correct, but I'm not sure. Also, which<br>
> >>> format would you recommend for storing the matrix?<br>
> >><br>
> >> 1) With any error, send the Entire error message.<br>
> >> 2) PETSc matrices are divided by rows, not rows and columns, see the<br>
> >> manual section. Rows & columns only makes sense for dense matrices<br>
> >> 3) You can still set arbitrary blocks no matter how the matrix is divided<br>
> >> 4) The error means you tried to set both local and global dimensions, and<br>
> >> they do not add up correctly. Just set the global dimensions<br>
> >> Matt<br>
> >><br>
> >>><br>
> >>> Jack,<br>
> >>><br>
> >>> I'm a summer intern just getting started with this project, so I don't<br>
> >>> know all the details yet (I can ask though). I know I need to find the<br>
> >>> Green's function which will involve the trace of the inverted Hamiltonian,<br>
> >>> as well as the rest of the matrix. I have inquired about avoiding the<br>
> >>> inversion altogether, but my instructor doesn't believe there is a way<br>
> >>> around it. Once I've worked through the math I want to explore other options<br>
> >>> though.<br>
> >>><br>
> >>> Respectfully,<br>
> >>> Adam<br>
> >>><br>
> >>> On Tue, Jun 28, 2011 at 6:08 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
> >>> wrote:<br>
> >>>><br>
> >>>> On Tue, Jun 28, 2011 at 5:01 PM, Adam Byrd <<a href="mailto:adam1.byrd@gmail.com">adam1.byrd@gmail.com</a>> wrote:<br>
> >>>>><br>
> >>>>> Actually, it's quite sparse. In the 3600x3600 there are only just 4<br>
> >>>>> nonzero entries in each row. This means it's 99.9% empty. My smaller 6x6<br>
> >>>>> example is dense, but it's only practice building and manipulating matrices.<br>
> >>>><br>
> >>>> Ah, then its easy. Just call MatSetValues() with each block. Then use<br>
> >>>> MUMPS to do a sparse direct solve.<br>
> >>>> Matt<br>
> >>>><br>
> >>>>><br>
> >>>>> Respectfully,<br>
> >>>>> Adam<br>
> >>>>><br>
> >>>>> On Tue, Jun 28, 2011 at 5:55 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
> >>>>> wrote:<br>
> >>>>>><br>
> >>>>>> It sounds like you have a dense matrix (from your example). Is this<br>
> >>>>>> true? If so, you should use Elemental (on Google Code).<br>
> >>>>>> Thanks,<br>
> >>>>>> Matt<br>
> >>>>>><br>
> >>>>>> On Tue, Jun 28, 2011 at 8:55 AM, Adam Byrd <<a href="mailto:adam1.byrd@gmail.com">adam1.byrd@gmail.com</a>><br>
> >>>>>> wrote:<br>
> >>>>>>><br>
> >>>>>>> Hi,<br>
> >>>>>>> I'm rather new to PETSc and trying to work out the best way to create<br>
> >>>>>>> and fill a large sparse matrix distributed over many processors. Currently,<br>
> >>>>>>> my goal is to create a 3600x3600 matrix in units of 12x12 blocks with<br>
> >>>>>>> several blocks on any given node. I'd like to create the matrix in such a<br>
> >>>>>>> way that each node only holds the information in it's handful of blocks and<br>
> >>>>>>> not the entire matrix. Eventually, this matrix is to be inverted (I know,<br>
> >>>>>>> inversion should be avoided, but as this is a Hamiltonian matrix from which<br>
> >>>>>>> I need the Green's function, I'm unaware of a way to forgo carrying out the<br>
> >>>>>>> inversion). Additionally, the values will be changed slightly and the matrix<br>
> >>>>>>> will be repeatedly inverted. It's structure will remain the same. In order<br>
> >>>>>>> to learn how to do this is I am starting with a small 6x6 matrix broken into<br>
> >>>>>>> four 3x3 blocks and distributed one block per node. I've been able to create<br>
> >>>>>>> a local 3x3 matrix on each node, with it's own values, and with the global<br>
> >>>>>>> row/column IDs correctly set to [0, 1, 2] or [3, 4, 5] depending on where<br>
> >>>>>>> the block is in the matrix. My problem manifests when I try to create the<br>
> >>>>>>> larger matrix from the individual smaller ones. When the matrix is<br>
> >>>>>>> constructed I'm trying to use MatSetValues and having each node pass in it's<br>
> >>>>>>> 3x3 block. I end up with an error that the sum of local lengths 12x12 does<br>
> >>>>>>> not match the global length 6x6. It appears as though this is from passing<br>
> >>>>>>> in four 3x3s and the program interpreting that as a 12x12 instead of as a<br>
> >>>>>>> 6x6 with the blocks in a grid.<br>
> >>>>>>> My question is then: is it possible to fill a matrix as a grid of<br>
> >>>>>>> blocks, or can I only fill it in groups of rows or columns? Also, am I<br>
> >>>>>>> approaching this problem the correct way, or are there more efficient ways<br>
> >>>>>>> of building this matrix with the ultimate goal of inverting it?<br>
> >>>>>>> I have included my copy of a modified example if it helps. I do<br>
> >>>>>>> apologize if this is answered somewhere in the documentation, I have been<br>
> >>>>>>> unable to find a solution.<br>
> >>>>>>> Respectfully,<br>
> >>>>>>> Adam<br>
> >>>>>><br>
> >>>>>><br>
> >>>>>> --<br>
> >>>>>> What most experimenters take for granted before they begin their<br>
> >>>>>> experiments is infinitely more interesting than any results to which their<br>
> >>>>>> experiments lead.<br>
> >>>>>> -- Norbert Wiener<br>
> >>>>><br>
> >>>><br>
> >>>><br>
> >>>><br>
> >>>> --<br>
> >>>> What most experimenters take for granted before they begin their<br>
> >>>> experiments is infinitely more interesting than any results to which their<br>
> >>>> experiments lead.<br>
> >>>> -- Norbert Wiener<br>
> >>><br>
> >><br>
> >><br>
> >><br>
> >> --<br>
> >> What most experimenters take for granted before they begin their<br>
> >> experiments is infinitely more interesting than any results to which their<br>
> >> experiments lead.<br>
> >> -- Norbert Wiener<br>
> ><br>
> ><br>
><br>
<br>
</div></div></blockquote></div><br>