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