<div dir="ltr">This looks like a false positive or there is some subtle bug here that we are not seeing.<br><div>Could this be the first time parallel PtAP has been used (and reported) in petsc4py?</div><div><br></div><div>Mark</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Oct 10, 2023 at 8:27 PM Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">On Tue, Oct 10, 2023 at 5:34 PM Thanasis Boutsikakis <<a href="mailto:thanasis.boutsikakis@corintis.com" target="_blank">thanasis.boutsikakis@corintis.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div>Hi all,<div><br></div><div>Revisiting my code and the proposed solution from Pierre, I realized this works only in sequential. The reason is that PETSc partitions those matrices only row-wise, which leads to an error due to the mismatch between number of columns of A (non-partitioned) and the number of rows of Phi (partitioned).</div></div></blockquote><div><br></div><div>Are you positive about this? P^T A P is designed to run in this scenario, so either we have a bug or the diagnosis is wrong.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div><div><div style="color:rgb(29,33,38);background-color:rgb(255,255,255);font-family:Menlo,Monaco,"Courier New",monospace;line-height:18px;white-space:pre-wrap"><div><span style="color:rgb(10,66,131)">"""Experimenting with PETSc mat-mat multiplication"""</span></div><br><div><span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">time</span></div><br><div><span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">numpy</span> <span style="color:rgb(168,3,20)">as</span> <span style="color:rgb(72,31,148)">np</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">colorama</span> <span style="color:rgb(168,3,20)">import</span> Fore</div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">firedrake</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(13,92,181)">COMM_SELF</span>, <span style="color:rgb(13,92,181)">COMM_WORLD</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">firedrake</span>.<span style="color:rgb(72,31,148)">petsc</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">PETSc</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">mpi4py</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">MPI</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">numpy</span>.<span style="color:rgb(72,31,148)">testing</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">assert_array_almost_equal</span></div><br><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">utilities</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">Print</span></div><br><div>nproc <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(13,92,181)">COMM_WORLD</span>.size</div><div>rank <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(13,92,181)">COMM_WORLD</span>.rank</div><div><br></div><div><div style="line-height:18px"><div><span style="color:rgb(168,3,20)">def</span> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>(<span style="color:rgb(166,61,0)">input_array,</span> <span style="color:rgb(166,61,0)">sparse</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">True</span>):</div><div> <span style="color:rgb(10,66,131)">"""Create a PETSc matrix from an input_array</span></div><br><div><span style="color:rgb(10,66,131)"> Args:</span></div><div><span style="color:rgb(10,66,131)"> input_array (np array): Input array</span></div><div><span style="color:rgb(10,66,131)"> partition_like (PETSc mat, optional): Petsc matrix. Defaults to None.</span></div><div><span style="color:rgb(10,66,131)"> sparse (bool, optional): Toggle for sparese or dense. Defaults to True.</span></div><br><div><span style="color:rgb(10,66,131)"> Returns:</span></div><div><span style="color:rgb(10,66,131)"> PETSc mat: PETSc mpi matrix</span></div><div><span style="color:rgb(10,66,131)"> """</span></div><div> <span style="color:rgb(164,176,188)"># Check if input_array is 1D and reshape if necessary</span></div><div> <span style="color:rgb(168,3,20)">assert</span> <span style="color:rgb(72,31,148)">len</span>(<span style="color:rgb(166,61,0)">input_array</span>.shape) <span style="color:rgb(168,3,20)">==</span> <span style="color:rgb(13,92,181)">2</span>, <span style="color:rgb(10,66,131)">"Input array should be 2-dimensional"</span></div><div> global_rows, global_cols <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(166,61,0)">input_array</span>.shape</div><div> size <span style="color:rgb(168,3,20)">=</span> ((<span style="color:rgb(13,92,181)">None</span>, global_rows), (global_cols, global_cols))</div><br><div> <span style="color:rgb(164,176,188)"># Create a sparse or dense matrix based on the 'sparse' argument</span></div><div> <span style="color:rgb(168,3,20)">if</span> <span style="color:rgb(166,61,0)">sparse</span>:</div><div> matrix <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">PETSc</span>.Mat().createAIJ(<span style="color:rgb(166,61,0)">size</span><span style="color:rgb(168,3,20)">=</span>size, <span style="color:rgb(166,61,0)">comm</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">COMM_WORLD</span>)</div><div> <span style="color:rgb(168,3,20)">else</span>:</div><div> matrix <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">PETSc</span>.Mat().createDense(<span style="color:rgb(166,61,0)">size</span><span style="color:rgb(168,3,20)">=</span>size, <span style="color:rgb(166,61,0)">comm</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">COMM_WORLD</span>)</div><div> matrix.setUp()</div><br><div> local_rows_start, local_rows_end <span style="color:rgb(168,3,20)">=</span> matrix.getOwnershipRange()</div><br><div> <span style="color:rgb(168,3,20)">for</span> counter, i <span style="color:rgb(168,3,20)">in</span> <span style="color:rgb(72,31,148)">enumerate</span>(<span style="color:rgb(72,31,148)">range</span>(local_rows_start, local_rows_end)):</div><div> <span style="color:rgb(164,176,188)"># Calculate the correct row in the array for the current process</span></div><div> row_in_array <span style="color:rgb(168,3,20)">=</span> counter <span style="color:rgb(168,3,20)">+</span> local_rows_start</div><div> matrix.setValues(</div><div> i, <span style="color:rgb(72,31,148)">range</span>(global_cols), <span style="color:rgb(166,61,0)">input_array</span>[row_in_array, :], <span style="color:rgb(166,61,0)">addv</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">False</span></div><div> )</div><br><div> <span style="color:rgb(164,176,188)"># Assembly the matrix to compute the final structure</span></div><div> matrix.assemblyBegin()</div><div> matrix.assemblyEnd()</div><br><div> <span style="color:rgb(168,3,20)">return</span> matrix</div></div></div><br><div><span style="color:rgb(164,176,188)"># --------------------------------------------</span></div><div><span style="color:rgb(164,176,188)"># EXP: Galerkin projection of an mpi PETSc matrix A with an mpi PETSc matrix Phi</span></div><div><span style="color:rgb(164,176,188)"># A' = Phi.T * A * Phi</span></div><div><span style="color:rgb(164,176,188)"># [k x k] <- [k x m] x [m x m] x [m x k]</span></div><div><span style="color:rgb(164,176,188)"># --------------------------------------------</span></div><br><div>m, k <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(13,92,181)">100</span>, <span style="color:rgb(13,92,181)">7</span></div><div><span style="color:rgb(164,176,188)"># Generate the random numpy matrices</span></div><div><span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">random</span>.seed(<span style="color:rgb(13,92,181)">0</span>) <span style="color:rgb(164,176,188)"># sets the seed to 0</span></div><div>A_np <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">random</span>.randint(<span style="color:rgb(166,61,0)">low</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">0</span>, <span style="color:rgb(166,61,0)">high</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">6</span>, <span style="color:rgb(166,61,0)">size</span><span style="color:rgb(168,3,20)">=</span>(m, m))</div><div>Phi_np <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">random</span>.randint(<span style="color:rgb(166,61,0)">low</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">0</span>, <span style="color:rgb(166,61,0)">high</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">6</span>, <span style="color:rgb(166,61,0)">size</span><span style="color:rgb(168,3,20)">=</span>(m, k))</div><br><div><span style="color:rgb(164,176,188)"># --------------------------------------------</span></div><div><span style="color:rgb(164,176,188)"># TEST: Galerking projection of numpy matrices A_np and Phi_np</span></div><div><span style="color:rgb(164,176,188)"># --------------------------------------------</span></div><div>Aprime_np <span style="color:rgb(168,3,20)">=</span> Phi_np.T <span style="color:rgb(13,92,181)">@</span> A_np <span style="color:rgb(13,92,181)">@</span> Phi_np</div><div><span style="color:rgb(72,31,148)">Print</span>(<span style="color:rgb(168,3,20)">f</span><span style="color:rgb(10,66,131)">"MATRIX Aprime_np [</span><span style="color:rgb(13,92,181)">{</span>Aprime_np.shape[<span style="color:rgb(13,92,181)">0</span>]<span style="color:rgb(13,92,181)">}</span><span style="color:rgb(10,66,131)">x</span><span style="color:rgb(13,92,181)">{</span>Aprime_np.shape[<span style="color:rgb(13,92,181)">1</span>]<span style="color:rgb(13,92,181)">}</span><span style="color:rgb(10,66,131)">]"</span>)</div><div><span style="color:rgb(72,31,148)">Print</span>(<span style="color:rgb(168,3,20)">f</span><span style="color:rgb(10,66,131)">"</span><span style="color:rgb(13,92,181)">{</span>Aprime_np<span style="color:rgb(13,92,181)">}</span><span style="color:rgb(10,66,131)">"</span>)</div><br><div><span style="color:rgb(164,176,188)"># Create A as an mpi matrix distributed on each process</span></div><div><span style="color:rgb(13,92,181)">A</span> <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>(A_np, <span style="color:rgb(166,61,0)">sparse</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">False</span>)</div><br><div><span style="color:rgb(164,176,188)"># Create Phi as an mpi matrix distributed on each process</span></div><div>Phi <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>(Phi_np, <span style="color:rgb(166,61,0)">sparse</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">False</span>)</div><br><div><span style="color:rgb(164,176,188)"># Create an empty PETSc matrix object to store the result of the PtAP operation.</span></div><div><span style="color:rgb(164,176,188)"># This will hold the result A' = Phi.T * A * Phi after the computation.</span></div><div>A_prime <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>(<span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">zeros</span>((k, k)), <span style="color:rgb(166,61,0)">sparse</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">False</span>)</div><br><div><span style="color:rgb(164,176,188)"># Perform the PtAP (Phi Transpose times A times Phi) operation.</span></div><div><span style="color:rgb(164,176,188)"># In mathematical terms, this operation is A' = Phi.T * A * Phi.</span></div><div><span style="color:rgb(164,176,188)"># A_prime will store the result of the operation.</span></div><div>A_prime <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(13,92,181)">A</span>.ptap(Phi)</div></div><div><br></div><div>Here is the error</div><div><br></div></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><div><div><div>MATRIX mpiaij A [100x100]</div></div></div><div><div><div>Assembled</div></div></div><div><div><div><br></div></div></div><div><div><div>Partitioning for A:</div></div></div><div><div><div> Rank 0: Rows [0, 34)</div></div></div><div><div><div> Rank 1: Rows [34, 67)</div></div></div><div><div><div> Rank 2: Rows [67, 100)</div></div></div><div><div><div><br></div></div></div><div><div><div>MATRIX mpiaij Phi [100x7]</div></div></div><div><div><div>Assembled</div></div></div><div><div><div><br></div></div></div><div><div><div>Partitioning for Phi:</div></div></div><div><div><div> Rank 0: Rows [0, 34)</div></div></div><div><div><div> Rank 1: Rows [34, 67)</div></div></div><div><div><div> Rank 2: Rows [67, 100)</div></div></div><div><div><div><br></div></div></div><div><div><div>Traceback (most recent call last):</div></div></div><div><div><div> File "/Users/boutsitron/work/galerkin_projection.py", line 87, in <module></div></div></div><div><div><div> A_prime = A.ptap(Phi)</div></div></div><div><div><div> ^^^^^^^^^^^</div></div></div><div><div><div> File "petsc4py/PETSc/Mat.pyx", line 1525, in petsc4py.PETSc.Mat.ptap</div></div></div><div><div><div>petsc4py.PETSc.Error: error code 60</div></div></div><div><div><div>[0] MatPtAP() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matrix.c:9896</div></div></div><div><div><div>[0] MatProductSetFromOptions() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matproduct.c:541</div></div></div><div><div><div>[0] MatProductSetFromOptions_Private() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matproduct.c:435</div></div></div><div><div><div>[0] MatProductSetFromOptions_MPIAIJ() at /Users/boutsitron/firedrake/src/petsc/src/mat/impls/aij/mpi/mpimatmatmult.c:2372</div></div></div><div><div><div>[0] MatProductSetFromOptions_MPIAIJ_PtAP() at /Users/boutsitron/firedrake/src/petsc/src/mat/impls/aij/mpi/mpimatmatmult.c:2266</div></div></div><div><div><div>[0] Nonconforming object sizes</div></div></div><div><div><div>[0] Matrix local dimensions are incompatible, Acol (0, 100) != Prow (0,34)</div></div></div><div><div><div>Abort(1) on node 0 (rank 0 in comm 496): application called MPI_Abort(PYOP2_COMM_WORLD, 1) - process 0</div></div></div></blockquote><div><div><div><br></div><div>Any thoughts?</div><div><br></div><div>Thanks,</div><div>Thanos</div></div><div><br><blockquote type="cite"><div>On 5 Oct 2023, at 14:23, Thanasis Boutsikakis <<a href="mailto:thanasis.boutsikakis@corintis.com" target="_blank">thanasis.boutsikakis@corintis.com</a>> wrote:</div><br><div><div>This works Pierre. Amazing input, thanks a lot!<br id="m_2797074109726544632m_-3099518781218751369lineBreakAtBeginningOfMessage"><div><br><blockquote type="cite"><div>On 5 Oct 2023, at 14:17, Pierre Jolivet <<a href="mailto:pierre@joliv.et" target="_blank">pierre@joliv.et</a>> wrote:</div><br><div><div><div>Not a petsc4py expert here, but you may to try instead:</div><div>A_prime = A.ptap(Phi)</div><div><br></div><div>Thanks,</div><div>Pierre</div><div><br><blockquote type="cite"><div>On 5 Oct 2023, at 2:02 PM, Thanasis Boutsikakis <<a href="mailto:thanasis.boutsikakis@corintis.com" target="_blank">thanasis.boutsikakis@corintis.com</a>> wrote:</div><br><div><div>Thanks Pierre! So I tried this and got a segmentation fault. Is this supposed to work right off the bat or am I missing sth?<div><br></div><div><div>[0]PETSC ERROR: ------------------------------------------------------------------------</div><div>[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range</div><div>[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger</div><div>[0]PETSC ERROR: or see <a href="https://petsc.org/release/faq/#valgrind" target="_blank">https://petsc.org/release/faq/#valgrind</a> and <a href="https://petsc.org/release/faq/" target="_blank">https://petsc.org/release/faq/</a></div><div>[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run</div><div>[0]PETSC ERROR: to get more information on the crash.</div><div>[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.</div><div>Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0</div><div><div><br></div><div><div style="color:rgb(29,33,38);background-color:rgb(255,255,255);font-family:Menlo,Monaco,"Courier New",monospace;line-height:18px;white-space:pre-wrap"><div><span style="color:rgb(10,66,131)">"""Experimenting with PETSc mat-mat multiplication"""</span></div><br><div><span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">time</span></div><br><div><span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">numpy</span> <span style="color:rgb(168,3,20)">as</span> <span style="color:rgb(72,31,148)">np</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">colorama</span> <span style="color:rgb(168,3,20)">import</span> Fore</div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">firedrake</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(13,92,181)">COMM_SELF</span>, <span style="color:rgb(13,92,181)">COMM_WORLD</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">firedrake</span>.<span style="color:rgb(72,31,148)">petsc</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">PETSc</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">mpi4py</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">MPI</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">numpy</span>.<span style="color:rgb(72,31,148)">testing</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">assert_array_almost_equal</span></div><br><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">utilities</span> <span style="color:rgb(168,3,20)">import</span> (</div><div> <span style="color:rgb(72,31,148)">Print</span>,</div><div> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>,</div><div> <span style="color:rgb(72,31,148)">print_matrix_partitioning</span>,</div><div>)</div><br><div>nproc <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(13,92,181)">COMM_WORLD</span>.size</div><div>rank <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(13,92,181)">COMM_WORLD</span>.rank</div><br><div><span style="color:rgb(164,176,188)"># --------------------------------------------</span></div><div><span style="color:rgb(164,176,188)"># EXP: Galerkin projection of an mpi PETSc matrix A with an mpi PETSc matrix Phi</span></div><div><span style="color:rgb(164,176,188)"># A' = Phi.T * A * Phi</span></div><div><span style="color:rgb(164,176,188)"># [k x k] <- [k x m] x [m x m] x [m x k]</span></div><div><span style="color:rgb(164,176,188)"># --------------------------------------------</span></div><br><div>m, k <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(13,92,181)">11</span>, <span style="color:rgb(13,92,181)">7</span></div><div><span style="color:rgb(164,176,188)"># Generate the random numpy matrices</span></div><div><span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">random</span>.seed(<span style="color:rgb(13,92,181)">0</span>) <span style="color:rgb(164,176,188)"># sets the seed to 0</span></div><div>A_np <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">random</span>.randint(<span style="color:rgb(166,61,0)">low</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">0</span>, <span style="color:rgb(166,61,0)">high</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">6</span>, <span style="color:rgb(166,61,0)">size</span><span style="color:rgb(168,3,20)">=</span>(m, m))</div><div>Phi_np <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">random</span>.randint(<span style="color:rgb(166,61,0)">low</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">0</span>, <span style="color:rgb(166,61,0)">high</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">6</span>, <span style="color:rgb(166,61,0)">size</span><span style="color:rgb(168,3,20)">=</span>(m, k))</div><br><div><span style="color:rgb(164,176,188)"># --------------------------------------------</span></div><div><span style="color:rgb(164,176,188)"># TEST: Galerking projection of numpy matrices A_np and Phi_np</span></div><div><span style="color:rgb(164,176,188)"># --------------------------------------------</span></div><div>Aprime_np <span style="color:rgb(168,3,20)">=</span> Phi_np.T <span style="color:rgb(13,92,181)">@</span> A_np <span style="color:rgb(13,92,181)">@</span> Phi_np</div><div><span style="color:rgb(72,31,148)">Print</span>(<span style="color:rgb(168,3,20)">f</span><span style="color:rgb(10,66,131)">"MATRIX Aprime_np [</span><span style="color:rgb(13,92,181)">{</span>Aprime_np.shape[<span style="color:rgb(13,92,181)">0</span>]<span style="color:rgb(13,92,181)">}</span><span style="color:rgb(10,66,131)">x</span><span style="color:rgb(13,92,181)">{</span>Aprime_np.shape[<span style="color:rgb(13,92,181)">1</span>]<span style="color:rgb(13,92,181)">}</span><span style="color:rgb(10,66,131)">]"</span>)</div><div><span style="color:rgb(72,31,148)">Print</span>(<span style="color:rgb(168,3,20)">f</span><span style="color:rgb(10,66,131)">"</span><span style="color:rgb(13,92,181)">{</span>Aprime_np<span style="color:rgb(13,92,181)">}</span><span style="color:rgb(10,66,131)">"</span>)</div><br><div><span style="color:rgb(164,176,188)"># Create A as an mpi matrix distributed on each process</span></div><div><span style="color:rgb(13,92,181)">A</span> <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>(A_np, <span style="color:rgb(166,61,0)">sparse</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">False</span>)</div><div><br></div><div><span style="color:rgb(164,176,188)"># Create Phi as an mpi matrix distributed on each process</span></div><div>Phi <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>(Phi_np, <span style="color:rgb(166,61,0)">sparse</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">False</span>)</div><div><br></div><div><span style="color:rgb(164,176,188)"># Create an empty PETSc matrix object to store the result of the PtAP operation.</span></div><div><span style="color:rgb(164,176,188)"># This will hold the result A' = Phi.T * A * Phi after the computation.</span></div><div>A_prime <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>(<span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">zeros</span>((k, k)), <span style="color:rgb(166,61,0)">sparse</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">False</span>)</div><div><br></div><div><span style="color:rgb(164,176,188)"># Perform the PtAP (Phi Transpose times A times Phi) operation.</span></div><div><span style="color:rgb(164,176,188)"># In mathematical terms, this operation is A' = Phi.T * A * Phi.</span></div><div><span style="color:rgb(164,176,188)"># A_prime will store the result of the operation.</span></div><div>Phi.PtAP(<span style="color:rgb(13,92,181)">A</span>, A_prime)</div></div><div><br><blockquote type="cite"><div>On 5 Oct 2023, at 13:22, Pierre Jolivet <<a href="mailto:pierre@joliv.et" target="_blank">pierre@joliv.et</a>> wrote:</div><br><div><div>How about using ptap which will use MatPtAP?<div>It will be more efficient (and it will help you bypass the issue).</div><div><br></div><div>Thanks,</div><div>Pierre<br id="m_2797074109726544632m_-3099518781218751369lineBreakAtBeginningOfMessage"><div><br><blockquote type="cite"><div>On 5 Oct 2023, at 1:18 PM, Thanasis Boutsikakis <<a href="mailto:thanasis.boutsikakis@corintis.com" target="_blank">thanasis.boutsikakis@corintis.com</a>> wrote:</div><br><div><div>Sorry, forgot function <span style="color:rgb(72,31,148);font-family:Menlo,Monaco,"Courier New",monospace;white-space:pre-wrap;background-color:rgb(255,255,255)">create_petsc_matrix()</span><div><br></div><div><div style="color:rgb(29,33,38);background-color:rgb(255,255,255);font-family:Menlo,Monaco,"Courier New",monospace;line-height:18px;white-space:pre-wrap"><div><span style="color:rgb(168,3,20)">def</span> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>(<span style="color:rgb(166,61,0)">input_array</span> <span style="color:rgb(166,61,0)">sparse</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">True</span>):</div><div> <span style="color:rgb(10,66,131)">"""Create a PETSc matrix from an input_array</span></div><br><div><span style="color:rgb(10,66,131)"> Args:</span></div><div><span style="color:rgb(10,66,131)"> input_array (np array): Input array</span></div><div><span style="color:rgb(10,66,131)"> partition_like (PETSc mat, optional): Petsc matrix. Defaults to None.</span></div><div><span style="color:rgb(10,66,131)"> sparse (bool, optional): Toggle for sparese or dense. Defaults to True.</span></div><br><div><span style="color:rgb(10,66,131)"> Returns:</span></div><div><span style="color:rgb(10,66,131)"> PETSc mat: PETSc matrix</span></div><div><span style="color:rgb(10,66,131)"> """</span></div><div> <span style="color:rgb(164,176,188)"># Check if input_array is 1D and reshape if necessary</span></div><div> <span style="color:rgb(168,3,20)">assert</span> <span style="color:rgb(72,31,148)">len</span>(<span style="color:rgb(166,61,0)">input_array</span>.shape) <span style="color:rgb(168,3,20)">==</span> <span style="color:rgb(13,92,181)">2</span>, <span style="color:rgb(10,66,131)">"Input array should be 2-dimensional"</span></div><div> global_rows, global_cols <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(166,61,0)">input_array</span>.shape</div><div><br></div><div> size <span style="color:rgb(168,3,20)">=</span> ((<span style="color:rgb(13,92,181)">None</span>, global_rows), (global_cols, global_cols))</div><br><div> <span style="color:rgb(164,176,188)"># Create a sparse or dense matrix based on the 'sparse' argument</span></div><div> <span style="color:rgb(168,3,20)">if</span> <span style="color:rgb(166,61,0)">sparse</span>:</div><div> matrix <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">PETSc</span>.Mat().createAIJ(<span style="color:rgb(166,61,0)">size</span><span style="color:rgb(168,3,20)">=</span>size, <span style="color:rgb(166,61,0)">comm</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">COMM_WORLD</span>)</div><div> <span style="color:rgb(168,3,20)">else</span>:</div><div> matrix <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">PETSc</span>.Mat().createDense(<span style="color:rgb(166,61,0)">size</span><span style="color:rgb(168,3,20)">=</span>size, <span style="color:rgb(166,61,0)">comm</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">COMM_WORLD</span>)</div><div> matrix.setUp()</div><br><div> local_rows_start, local_rows_end <span style="color:rgb(168,3,20)">=</span> matrix.getOwnershipRange()</div><br><div> <span style="color:rgb(168,3,20)">for</span> counter, i <span style="color:rgb(168,3,20)">in</span> <span style="color:rgb(72,31,148)">enumerate</span>(<span style="color:rgb(72,31,148)">range</span>(local_rows_start, local_rows_end)):</div><div> <span style="color:rgb(164,176,188)"># Calculate the correct row in the array for the current process</span></div><div> row_in_array <span style="color:rgb(168,3,20)">=</span> counter <span style="color:rgb(168,3,20)">+</span> local_rows_start</div><div> matrix.setValues(</div><div> i, <span style="color:rgb(72,31,148)">range</span>(global_cols), <span style="color:rgb(166,61,0)">input_array</span>[row_in_array, :], <span style="color:rgb(166,61,0)">addv</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">False</span></div><div> )</div><br><div> <span style="color:rgb(164,176,188)"># Assembly the matrix to compute the final structure</span></div><div> matrix.assemblyBegin()</div><div> matrix.assemblyEnd()</div><br><div> <span style="color:rgb(168,3,20)">return</span> matrix</div></div><div><br><blockquote type="cite"><div>On 5 Oct 2023, at 13:09, Thanasis Boutsikakis <<a href="mailto:thanasis.boutsikakis@corintis.com" target="_blank">thanasis.boutsikakis@corintis.com</a>> wrote:</div><br><div><div>Hi everyone,<div><br></div><div>I am trying a Galerkin projection (see MFE below) and I cannot get the <span style="color:rgb(29,33,38);font-family:Menlo,Monaco,"Courier New",monospace;white-space:pre-wrap;background-color:rgb(255,255,255)">Phi</span><span style="color:rgb(29,33,38);font-family:Menlo,Monaco,"Courier New",monospace;white-space:pre-wrap;background-color:rgb(255,255,255)">.transposeMatMult(A, A1) work. The error is</span></div><div><span style="color:rgb(29,33,38);font-family:Menlo,Monaco,"Courier New",monospace;white-space:pre-wrap;background-color:rgb(255,255,255)"><br></span></div><div><div> Phi.transposeMatMult(A, A1)</div><div> File "petsc4py/PETSc/Mat.pyx", line 1514, in petsc4py.PETSc.Mat.transposeMatMult</div><div>petsc4py.PETSc.Error: error code 56</div><div>[0] MatTransposeMatMult() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matrix.c:10135</div><div>[0] MatProduct_Private() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matrix.c:9989</div><div>[0] No support for this operation for this object type</div><div>[0] Call MatProductCreate() first</div></div><div><br></div><div>Do you know if these exposed to petsc4py or maybe there is another way? I cannot get the MFE to work (neither in sequential nor in parallel)</div><div><br></div><div><div style="color:rgb(29,33,38);background-color:rgb(255,255,255);font-family:Menlo,Monaco,"Courier New",monospace;line-height:18px;white-space:pre-wrap"><div><span style="color:rgb(10,66,131)">"""Experimenting with PETSc mat-mat multiplication"""</span></div><br><div><span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">time</span></div><br><div><span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">numpy</span> <span style="color:rgb(168,3,20)">as</span> <span style="color:rgb(72,31,148)">np</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">colorama</span> <span style="color:rgb(168,3,20)">import</span> Fore</div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">firedrake</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(13,92,181)">COMM_SELF</span>, <span style="color:rgb(13,92,181)">COMM_WORLD</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">firedrake</span>.<span style="color:rgb(72,31,148)">petsc</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">PETSc</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">mpi4py</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">MPI</span></div><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">numpy</span>.<span style="color:rgb(72,31,148)">testing</span> <span style="color:rgb(168,3,20)">import</span> <span style="color:rgb(72,31,148)">assert_array_almost_equal</span></div><br><div><span style="color:rgb(168,3,20)">from</span> <span style="color:rgb(72,31,148)">utilities</span> <span style="color:rgb(168,3,20)">import</span> (</div><div> <span style="color:rgb(72,31,148)">Print</span>,</div><div> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>,</div><div>)</div><br><div>nproc <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(13,92,181)">COMM_WORLD</span>.size</div><div>rank <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(13,92,181)">COMM_WORLD</span>.rank</div><br><div><span style="color:rgb(164,176,188)"># --------------------------------------------</span></div><div><span style="color:rgb(164,176,188)"># EXP: Galerkin projection of an mpi PETSc matrix A with an mpi PETSc matrix Phi</span></div><div><span style="color:rgb(164,176,188)"># A' = Phi.T * A * Phi</span></div><div><span style="color:rgb(164,176,188)"># [k x k] <- [k x m] x [m x m] x [m x k]</span></div><div><span style="color:rgb(164,176,188)"># --------------------------------------------</span></div><br><div>m, k <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(13,92,181)">11</span>, <span style="color:rgb(13,92,181)">7</span></div><div><span style="color:rgb(164,176,188)"># Generate the random numpy matrices</span></div><div><span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">random</span>.seed(<span style="color:rgb(13,92,181)">0</span>) <span style="color:rgb(164,176,188)"># sets the seed to 0</span></div><div>A_np <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">random</span>.randint(<span style="color:rgb(166,61,0)">low</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">0</span>, <span style="color:rgb(166,61,0)">high</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">6</span>, <span style="color:rgb(166,61,0)">size</span><span style="color:rgb(168,3,20)">=</span>(m, m))</div><div>Phi_np <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">random</span>.randint(<span style="color:rgb(166,61,0)">low</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">0</span>, <span style="color:rgb(166,61,0)">high</span><span style="color:rgb(168,3,20)">=</span><span style="color:rgb(13,92,181)">6</span>, <span style="color:rgb(166,61,0)">size</span><span style="color:rgb(168,3,20)">=</span>(m, k))</div><br><div><span style="color:rgb(164,176,188)"># Create A as an mpi matrix distributed on each process</span></div><div><span style="color:rgb(13,92,181)">A</span> <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>(A_np)</div><div><br></div><div><span style="color:rgb(164,176,188)"># Create Phi as an mpi matrix distributed on each process</span></div><div>Phi <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>(Phi_np)</div><div><br></div><div><span style="color:rgb(13,92,181)">A1</span> <span style="color:rgb(168,3,20)">=</span> <span style="color:rgb(72,31,148)">create_petsc_matrix</span>(<span style="color:rgb(72,31,148)">np</span>.<span style="color:rgb(72,31,148)">zeros</span>((k, m)))</div><div><br></div><div><span style="color:rgb(164,176,188)"># Now A1 contains the result of Phi^T * A</span></div><div>Phi.transposeMatMult(<span style="color:rgb(13,92,181)">A</span>, <span style="color:rgb(13,92,181)">A1</span>)</div><br></div></div></div></div></blockquote></div><br></div></div></div></blockquote></div><br></div></div></div></blockquote></div><br></div></div></div></div></div></blockquote></div><br></div></div></blockquote></div><br></div></div></blockquote></div><br></div></div></div></blockquote></div><br clear="all"><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</blockquote></div>