<div dir="ltr"><div dir="ltr">On Fri, May 6, 2022 at 9:28 AM Quentin Chevalier <<a href="mailto:quentin.chevalier@polytechnique.edu">quentin.chevalier@polytechnique.edu</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 dir="ltr"><div>Sorry for forgetting the list. Making two matrices was more of a precaution then a carefully thought strategy.</div><div><br></div><div>It would seem the MWE  as I provided it above (with a setDimensions to reduce calculation time) does work in sequential and gives the eigenvalue 118897.88711586884.</div><div><br></div><div>I had a working MUMPS solver config for a similar problem so I simply put it there. This gives the following code.</div><div><br></div><div>However, running it even in sequential gives a "MatSolverType mumps does not support matrix type python" error.<br></div><div>petsc4py.PETSc.Error: error code 92<br>[0] EPSSolve() at /usr/local/slepc/src/eps/interface/epssolve.c:136<br>[0] EPSSetUp() at /usr/local/slepc/src/eps/interface/epssetup.c:350<br>[0] STSetUp() at /usr/local/slepc/src/sys/classes/st/interface/stsolve.c:582<br>[0] STSetUp_Sinvert() at /usr/local/slepc/src/sys/classes/st/impls/sinvert/sinvert.c:123<br>[0] KSPSetUp() at /usr/local/petsc/src/ksp/ksp/interface/itfunc.c:408<br>[0] PCSetUp() at /usr/local/petsc/src/ksp/pc/interface/precon.c:1016<br>[0] PCSetUp_LU() at /usr/local/petsc/src/ksp/pc/impls/factor/lu/lu.c:82<br>[0] MatGetFactor() at /usr/local/petsc/src/mat/interface/matrix.c:4779<br>[0] See <a href="https://petsc.org/release/overview/linear_solve_table/" target="_blank">https://petsc.org/release/overview/linear_solve_table/</a> for possible LU and Cholesky solvers<br>[0] MatSolverType mumps does not support matrix type python</div></div></blockquote><div><br></div><div>You would need type AIJ for MUMPS.</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 dir="ltr"><div>Did I do something wrong ? I'm a bit confused by your FAQ entry, I never had a problem with this configuration in parallel.<br></div><div> <br><div style="color:rgb(212,212,212);background-color:rgb(30,30,30);font-family:"Droid Sans Mono","monospace",monospace,"Droid Sans Fallback";font-weight:normal;font-size:14px;line-height:19px;white-space:pre-wrap"><div><span style="color:rgb(197,134,192)">from</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(78,201,176)">petsc4py</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(197,134,192)">import</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(78,201,176)">PETSc</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(197,134,192)">as</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(78,201,176)">pet</span></div><div><span style="color:rgb(197,134,192)">from</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(78,201,176)">slepc4py</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(197,134,192)">import</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(78,201,176)">SLEPc</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(197,134,192)">as</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(78,201,176)">slp</span></div><div><span style="color:rgb(197,134,192)">from</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(78,201,176)">mpi4py</span><span style="color:rgb(212,212,212)">.</span><span style="color:rgb(78,201,176)">MPI</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(197,134,192)">import</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(79,193,255)">COMM_WORLD</span></div><br><div><span style="color:rgb(156,220,254)">dir</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(206,145,120)">"./sanity_check_mats/"</span></div><br><div><span style="color:rgb(106,153,85)"># Global sizes</span></div><div><span style="color:rgb(156,220,254)">m</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(181,206,168)">980143</span><span style="color:rgb(212,212,212)">; </span><span style="color:rgb(156,220,254)">m_local</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(156,220,254)">m</span></div><div><span style="color:rgb(156,220,254)">n</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(181,206,168)">904002</span><span style="color:rgb(212,212,212)">; </span><span style="color:rgb(156,220,254)">n_local</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(156,220,254)">n</span></div><div><span style="color:rgb(197,134,192)">if</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(79,193,255)">COMM_WORLD</span><span style="color:rgb(212,212,212)">.</span><span style="color:rgb(156,220,254)">size</span><span style="color:rgb(212,212,212)">></span><span style="color:rgb(181,206,168)">1</span><span style="color:rgb(212,212,212)">:</span></div><div><span style="color:rgb(212,212,212)">    </span><span style="color:rgb(156,220,254)">m_local</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(79,193,255)">COMM_WORLD</span><span style="color:rgb(212,212,212)">.</span><span style="color:rgb(156,220,254)">rank</span><span style="color:rgb(212,212,212)">*(</span><span style="color:rgb(181,206,168)">490363</span><span style="color:rgb(212,212,212)">-</span><span style="color:rgb(181,206,168)">489780</span><span style="color:rgb(212,212,212)">)+</span><span style="color:rgb(181,206,168)">489780</span></div><div><span style="color:rgb(212,212,212)">    </span><span style="color:rgb(156,220,254)">n_local</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(79,193,255)">COMM_WORLD</span><span style="color:rgb(212,212,212)">.</span><span style="color:rgb(156,220,254)">rank</span><span style="color:rgb(212,212,212)">*(</span><span style="color:rgb(181,206,168)">452259</span><span style="color:rgb(212,212,212)">-</span><span style="color:rgb(181,206,168)">451743</span><span style="color:rgb(212,212,212)">)+</span><span style="color:rgb(181,206,168)">451743</span></div><div><span style="color:rgb(212,212,212)">    </span><span style="color:rgb(156,220,254)">dir</span><span style="color:rgb(212,212,212)">+=</span><span style="color:rgb(206,145,120)">"par/"</span></div><div><span style="color:rgb(197,134,192)">else</span><span style="color:rgb(212,212,212)">:</span></div><div><span style="color:rgb(212,212,212)">    </span><span style="color:rgb(156,220,254)">m_local</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(156,220,254)">m</span></div><div><span style="color:rgb(212,212,212)">    </span><span style="color:rgb(156,220,254)">n_local</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(156,220,254)">n</span></div><div><span style="color:rgb(212,212,212)">    </span><span style="color:rgb(156,220,254)">dir</span><span style="color:rgb(212,212,212)">+=</span><span style="color:rgb(206,145,120)">"seq/"</span></div><br><div><span style="color:rgb(79,193,255)">QE</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(78,201,176)">pet</span><span style="color:rgb(212,212,212)">.Mat().createAIJ([[</span><span style="color:rgb(156,220,254)">m_local</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">m</span><span style="color:rgb(212,212,212)">],[</span><span style="color:rgb(156,220,254)">n_local</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">n</span><span style="color:rgb(212,212,212)">]])</span></div><div><span style="color:rgb(79,193,255)">L</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(78,201,176)">pet</span><span style="color:rgb(212,212,212)">.Mat().createAIJ([[</span><span style="color:rgb(156,220,254)">m_local</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">m</span><span style="color:rgb(212,212,212)">],[</span><span style="color:rgb(156,220,254)">m_local</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">m</span><span style="color:rgb(212,212,212)">]])</span></div><div><span style="color:rgb(156,220,254)">Mf</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(78,201,176)">pet</span><span style="color:rgb(212,212,212)">.Mat().createAIJ([[</span><span style="color:rgb(156,220,254)">n_local</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">n</span><span style="color:rgb(212,212,212)">],[</span><span style="color:rgb(156,220,254)">n_local</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">n</span><span style="color:rgb(212,212,212)">]])</span></div><br><div><span style="color:rgb(156,220,254)">viewerQE</span><span style="color:rgb(212,212,212)"> = </span><span style="color:rgb(78,201,176)">pet</span><span style="color:rgb(212,212,212)">.Viewer().createMPIIO(</span><span style="color:rgb(156,220,254)">dir</span><span style="color:rgb(212,212,212)">+</span><span style="color:rgb(206,145,120)">"QE.dat"</span><span style="color:rgb(212,212,212)">, </span><span style="color:rgb(206,145,120)">'r'</span><span style="color:rgb(212,212,212)">, </span><span style="color:rgb(79,193,255)">COMM_WORLD</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(79,193,255)">QE</span><span style="color:rgb(212,212,212)">.load(</span><span style="color:rgb(156,220,254)">viewerQE</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(156,220,254)">viewerL</span><span style="color:rgb(212,212,212)"> = </span><span style="color:rgb(78,201,176)">pet</span><span style="color:rgb(212,212,212)">.Viewer().createMPIIO(</span><span style="color:rgb(156,220,254)">dir</span><span style="color:rgb(212,212,212)">+</span><span style="color:rgb(206,145,120)">"L.dat"</span><span style="color:rgb(212,212,212)">, </span><span style="color:rgb(206,145,120)">'r'</span><span style="color:rgb(212,212,212)">, </span><span style="color:rgb(79,193,255)">COMM_WORLD</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(79,193,255)">L</span><span style="color:rgb(212,212,212)">.load(</span><span style="color:rgb(156,220,254)">viewerL</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(156,220,254)">viewerMf</span><span style="color:rgb(212,212,212)"> = </span><span style="color:rgb(78,201,176)">pet</span><span style="color:rgb(212,212,212)">.Viewer().createMPIIO(</span><span style="color:rgb(156,220,254)">dir</span><span style="color:rgb(212,212,212)">+</span><span style="color:rgb(206,145,120)">"Mf.dat"</span><span style="color:rgb(212,212,212)">, </span><span style="color:rgb(206,145,120)">'r'</span><span style="color:rgb(212,212,212)">, </span><span style="color:rgb(79,193,255)">COMM_WORLD</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(156,220,254)">Mf</span><span style="color:rgb(212,212,212)">.load(</span><span style="color:rgb(156,220,254)">viewerMf</span><span style="color:rgb(212,212,212)">)</span></div><br><div><span style="color:rgb(79,193,255)">QE</span><span style="color:rgb(212,212,212)">.assemble()</span></div><div><span style="color:rgb(79,193,255)">L</span><span style="color:rgb(212,212,212)">.assemble()</span></div><div><span style="color:rgb(156,220,254)">Mf</span><span style="color:rgb(212,212,212)">.assemble()</span></div><br><div><span style="color:rgb(156,220,254)">KSPs</span><span style="color:rgb(212,212,212)"> = []</span></div><div><span style="color:rgb(106,153,85)"># Useful solvers (here to put options for computing a smart R)</span></div><div><span style="color:rgb(197,134,192)">for</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(156,220,254)">Mat</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(197,134,192)">in</span><span style="color:rgb(212,212,212)"> [</span><span style="color:rgb(79,193,255)">L</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(79,193,255)">L</span><span style="color:rgb(212,212,212)">.hermitianTranspose()]:</span></div><div><span style="color:rgb(212,212,212)">    </span><span style="color:rgb(79,193,255)">KSP</span><span style="color:rgb(212,212,212)"> = </span><span style="color:rgb(78,201,176)">pet</span><span style="color:rgb(212,212,212)">.KSP().create()</span></div><div><span style="color:rgb(212,212,212)">    </span><span style="color:rgb(79,193,255)">KSP</span><span style="color:rgb(212,212,212)">.setOperators(</span><span style="color:rgb(156,220,254)">Mat</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(212,212,212)">    </span><span style="color:rgb(79,193,255)">KSP</span><span style="color:rgb(212,212,212)">.setFromOptions()</span></div><div><span style="color:rgb(212,212,212)">    </span><span style="color:rgb(156,220,254)">KSPs</span><span style="color:rgb(212,212,212)">.</span><span style="color:rgb(220,220,170)">append</span><span style="color:rgb(212,212,212)">(</span><span style="color:rgb(79,193,255)">KSP</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(212,212,212)">    </span></div><div><span style="color:rgb(86,156,214)">class</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(78,201,176)">LHS_class</span><span style="color:rgb(212,212,212)">:            </span></div><div><span style="color:rgb(212,212,212)">    </span><span style="color:rgb(86,156,214)">def</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(220,220,170)">mult</span><span style="color:rgb(212,212,212)">(</span><span style="color:rgb(156,220,254)">self</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">A</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">x</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">y</span><span style="color:rgb(212,212,212)">):</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(156,220,254)">w</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(78,201,176)">pet</span><span style="color:rgb(212,212,212)">.Vec().createMPI([</span><span style="color:rgb(156,220,254)">m_local</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">m</span><span style="color:rgb(212,212,212)">],</span><span style="color:rgb(156,220,254)">comm</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(79,193,255)">COMM_WORLD</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(156,220,254)">z</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(78,201,176)">pet</span><span style="color:rgb(212,212,212)">.Vec().createMPI([</span><span style="color:rgb(156,220,254)">m_local</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">m</span><span style="color:rgb(212,212,212)">],</span><span style="color:rgb(156,220,254)">comm</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(79,193,255)">COMM_WORLD</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(79,193,255)">QE</span><span style="color:rgb(212,212,212)">.mult(</span><span style="color:rgb(156,220,254)">x</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">w</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(156,220,254)">KSPs</span><span style="color:rgb(212,212,212)">[</span><span style="color:rgb(181,206,168)">0</span><span style="color:rgb(212,212,212)">].solve(</span><span style="color:rgb(156,220,254)">w</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">z</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(156,220,254)">KSPs</span><span style="color:rgb(212,212,212)">[</span><span style="color:rgb(181,206,168)">1</span><span style="color:rgb(212,212,212)">].solve(</span><span style="color:rgb(156,220,254)">z</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">w</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(79,193,255)">QE</span><span style="color:rgb(212,212,212)">.multTranspose(</span><span style="color:rgb(156,220,254)">w</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">y</span><span style="color:rgb(212,212,212)">)</span></div><br><div><span style="color:rgb(106,153,85)"># Matrix free operator</span></div><div><span style="color:rgb(79,193,255)">LHS</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(78,201,176)">pet</span><span style="color:rgb(212,212,212)">.Mat()</span></div><div><span style="color:rgb(79,193,255)">LHS</span><span style="color:rgb(212,212,212)">.create(</span><span style="color:rgb(156,220,254)">comm</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(79,193,255)">COMM_WORLD</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(79,193,255)">LHS</span><span style="color:rgb(212,212,212)">.setSizes([[</span><span style="color:rgb(156,220,254)">n_local</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">n</span><span style="color:rgb(212,212,212)">],[</span><span style="color:rgb(156,220,254)">n_local</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">n</span><span style="color:rgb(212,212,212)">]])</span></div><div><span style="color:rgb(79,193,255)">LHS</span><span style="color:rgb(212,212,212)">.setType(</span><span style="color:rgb(78,201,176)">pet</span><span style="color:rgb(212,212,212)">.Mat.Type.PYTHON)</span></div><div><span style="color:rgb(79,193,255)">LHS</span><span style="color:rgb(212,212,212)">.setPythonContext(</span><span style="color:rgb(78,201,176)">LHS_class</span><span style="color:rgb(212,212,212)">())</span></div><div><span style="color:rgb(79,193,255)">LHS</span><span style="color:rgb(212,212,212)">.setUp()</span></div><br><div><span style="color:rgb(156,220,254)">x</span><span style="color:rgb(212,212,212)">, </span><span style="color:rgb(156,220,254)">y</span><span style="color:rgb(212,212,212)"> = </span><span style="color:rgb(79,193,255)">LHS</span><span style="color:rgb(212,212,212)">.createVecs()</span></div><div><span style="color:rgb(156,220,254)">x</span><span style="color:rgb(212,212,212)">.set(</span><span style="color:rgb(181,206,168)">1</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(79,193,255)">LHS</span><span style="color:rgb(212,212,212)">.mult(</span><span style="color:rgb(156,220,254)">x</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">y</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(220,220,170)">print</span><span style="color:rgb(212,212,212)">(</span><span style="color:rgb(156,220,254)">y</span><span style="color:rgb(212,212,212)">.norm())</span></div><br><div><span style="color:rgb(106,153,85)"># Eigensolver</span></div><div><span style="color:rgb(79,193,255)">EPS</span><span style="color:rgb(212,212,212)"> = </span><span style="color:rgb(78,201,176)">slp</span><span style="color:rgb(212,212,212)">.EPS(); </span><span style="color:rgb(79,193,255)">EPS</span><span style="color:rgb(212,212,212)">.create()</span></div><div><span style="color:rgb(79,193,255)">EPS</span><span style="color:rgb(212,212,212)">.setOperators(</span><span style="color:rgb(79,193,255)">LHS</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">Mf</span><span style="color:rgb(212,212,212)">) </span><span style="color:rgb(106,153,85)"># Solve QE^T*L^-1H*L^-1*QE*f=sigma^2*Mf*f (cheaper than a proper SVD)</span></div><div><span style="color:rgb(79,193,255)">EPS</span><span style="color:rgb(212,212,212)">.setProblemType(</span><span style="color:rgb(78,201,176)">slp</span><span style="color:rgb(212,212,212)">.EPS.ProblemType.GHEP) </span><span style="color:rgb(106,153,85)"># Specify that A is hermitian (by construction), but M is semi-positive</span></div><div><span style="color:rgb(79,193,255)">EPS</span><span style="color:rgb(212,212,212)">.setWhichEigenpairs(</span><span style="color:rgb(79,193,255)">EPS</span><span style="color:rgb(212,212,212)">.Which.LARGEST_MAGNITUDE) </span><span style="color:rgb(106,153,85)"># Find largest eigenvalues</span></div><div><span style="color:rgb(206,145,120)"># Spectral transform</span></div><div><span style="color:rgb(206,145,120)">ST  = EPS.getST();  ST.setType('sinvert')</span></div><div><span style="color:rgb(206,145,120)"># Krylov subspace</span></div><div><span style="color:rgb(206,145,120)">KSP = ST.getKSP(); KSP.setType('preonly')</span></div><div><span style="color:rgb(206,145,120)"># Preconditioner</span></div><div><span style="color:rgb(206,145,120)">PC  = KSP.getPC();  PC.setType('lu')</span></div><div><span style="color:rgb(206,145,120)">PC.setFactorSolverType('mumps')</span></div><div><span style="color:rgb(206,145,120)"></span></div><div><span style="color:rgb(79,193,255)">EPS</span><span style="color:rgb(212,212,212)">.setDimensions(</span><span style="color:rgb(181,206,168)">1</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(181,206,168)">5</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(79,193,255)">EPS</span><span style="color:rgb(212,212,212)">.setFromOptions()</span></div><div><span style="color:rgb(79,193,255)">EPS</span><span style="color:rgb(212,212,212)">.solve()</span></div><div><span style="color:rgb(220,220,170)">print</span><span style="color:rgb(212,212,212)">(</span><span style="color:rgb(79,193,255)">EPS</span><span style="color:rgb(212,212,212)">.getEigenvalue(</span><span style="color:rgb(181,206,168)">0</span><span style="color:rgb(212,212,212)">))</span></div></div><br>Quentin CHEVALIER – IA parcours recherche<br>LadHyX - Ecole polytechnique<br>__________<br><br><br><br>On Fri, 6 May 2022 at 14:47, Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>> wrote:<br>><br>> Please respond to the list.<br>><br>> Why would you need different matrices for sequential and parallel? PETSc takes care of loading matrices from binary files in parallel codes. If I use the 'seq' matrices for both sequential and parallel runs it works. I get the eigenvalue 128662.745858<br>><br>> Take into account that in parallel you need to use a parallel linear solver, e.g., configure PETSc with MUMPS. See the FAQ #10 <a href="https://slepc.upv.es/documentation/faq.htm#faq10" target="_blank">https://slepc.upv.es/documentation/faq.htm#faq10</a><br>><br>> Jose<br>><br>> > El 5 may 2022, a las 14:57, Quentin Chevalier <<a href="mailto:quentin.chevalier@polytechnique.edu" target="_blank">quentin.chevalier@polytechnique.edu</a>> escribió:<br>> ><br>> > Thank you for your answer, Jose.<br>> ><br>> > It would appear my problem is slightly more complicated than it appears. With slight modification of the original MWE to account for functioning in serial or parallel, and computing the matrices in either sequential or parallel (apologies, it's a large file), and including a slightly modified version of your test case, I obtain two different results in serial and parallel (commands python3 MWE.py and mpirun -n 2 python3 MWE.py).<br>> ><br>> > Serial gives my a finite norm (around 28) and parallel gives me two norms of 0 and a code 77.<br>> ><br>> > This is still a problem as I would really like my code to be parallel. I saved my matrices using :<br>> > viewerQE = pet.Viewer().createMPIIO("QE.dat", 'w', COMM_WORLD)<br>> > QE.view(viewerQE)<br>> > viewerL = pet.Viewer().createMPIIO("L.dat", 'w', COMM_WORLD)<br>> > L.view(viewerL)<br>> > viewerMq = pet.Viewer().createMPIIO("Mq.dat", 'w', COMM_WORLD)<br>> > Mq.view(viewerMq)<br>> > viewerMf = pet.Viewer().createMPIIO("Mf.dat", 'w', COMM_WORLD)<br>> > Mf.view(viewerMf)<br>> ><br>> > New MWE (also attached for convenience) :<br>> > from petsc4py import PETSc as pet<br>> > from slepc4py import SLEPc as slp<br>> > from mpi4py.MPI import COMM_WORLD<br>> ><br>> > dir="./mats/"<br>> ><br>> > # Global sizes<br>> > m=980143; m_local=m<br>> > n=904002; n_local=n<br>> > if COMM_WORLD.size>1:<br>> >     m_local=COMM_WORLD.rank*(490363-489780)+489780<br>> >     n_local=COMM_WORLD.rank*(452259-451743)+451743<br>> >     dir+="par/"<br>> > else:<br>> >     m_local=m<br>> >     n_local=n<br>> >     dir+="seq/"<br>> ><br>> > QE=pet.Mat().createAIJ([[m_local,m],[n_local,n]])<br>> > L=pet.Mat().createAIJ([[m_local,m],[m_local,m]])<br>> > Mf=pet.Mat().createAIJ([[n_local,n],[n_local,n]])<br>> ><br>> > viewerQE = pet.Viewer().createMPIIO(dir+"QE.dat", 'r', COMM_WORLD)<br>> > QE.load(viewerQE)<br>> > viewerL = pet.Viewer().createMPIIO(dir+"L.dat", 'r', COMM_WORLD)<br>> > L.load(viewerL)<br>> > viewerMf = pet.Viewer().createMPIIO(dir+"Mf.dat", 'r', COMM_WORLD)<br>> > Mf.load(viewerMf)<br>> ><br>> > QE.assemble()<br>> > L.assemble()<br>> > Mf.assemble()<br>> ><br>> > KSPs = []<br>> > # Useful solvers (here to put options for computing a smart R)<br>> > for Mat in [L,L.hermitianTranspose()]:<br>> >     KSP = pet.KSP().create()<br>> >     KSP.setOperators(Mat)<br>> >     KSP.setFromOptions()<br>> >     KSPs.append(KSP)<br>> >     <br>> > class LHS_class:           <br>> >     def mult(self,A,x,y):<br>> >         w=pet.Vec().createMPI([m_local,m],comm=COMM_WORLD)<br>> >         z=pet.Vec().createMPI([m_local,m],comm=COMM_WORLD)<br>> >         QE.mult(x,w)<br>> >         KSPs[0].solve(w,z)<br>> >         KSPs[1].solve(z,w)<br>> >         QE.multTranspose(w,y)<br>> ><br>> > # Matrix free operator<br>> > LHS=pet.Mat()<br>> > LHS.create(comm=COMM_WORLD)<br>> > LHS.setSizes([[n_local,n],[n_local,n]])<br>> > LHS.setType(pet.Mat.Type.PYTHON)<br>> > LHS.setPythonContext(LHS_class())<br>> > LHS.setUp()<br>> ><br>> > x, y = LHS.createVecs()<br>> > x.set(1)<br>> > LHS.mult(x,y)<br>> > print(y.norm())<br>> ><br>> > # Eigensolver<br>> > EPS = slp.EPS(); EPS.create()<br>> > EPS.setOperators(LHS,Mf) # Solve QE^T*L^-1H*L^-1*QE*f=sigma^2*Mf*f (cheaper than a proper SVD)<br>> > EPS.setProblemType(slp.EPS.ProblemType.GHEP) # Specify that A is hermitian (by construction), but M is semi-positive<br>> > EPS.setWhichEigenpairs(EPS.Which.LARGEST_MAGNITUDE) # Find largest eigenvalues<br>> > EPS.setFromOptions()<br>> > EPS.solve()<br>> ><br>> ><br>> > Quentin CHEVALIER – IA parcours recherche<br>> > LadHyX - Ecole polytechnique<br>> > __________<br>> ><br>> ><br>> > <br>> ><br>> ><br>> > Quentin CHEVALIER – IA parcours recherche<br>> > LadHyX - Ecole polytechnique<br>> ><br>> > __________<br>> ><br>> ><br>> ><br>> > On Thu, 5 May 2022 at 12:05, Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>> wrote:<br>> > Your operator is not well formed. If you do this:<br>> ><br>> > x, y = LHS.createVecs()<br>> > x.set(1)<br>> > LHS.mult(x,y)<br>> > y.view()<br>> ><br>> > you will see that the output is all zeros. That is why SLEPc complains that "Initial vector is zero or belongs to the deflation space".<br>> ><br>> > Jose<br>> ><br>> ><br>> > > El 5 may 2022, a las 10:46, Quentin Chevalier <<a href="mailto:quentin.chevalier@polytechnique.edu" target="_blank">quentin.chevalier@polytechnique.edu</a>> escribió:<br>> > ><br>> > > Just a quick amend on the previous statement ; the problem arises in<br>> > > sequential and parallel. The MWE as is is provided for the parallel<br>> > > case, but imposing m_local=m makes it go sequential.<br>> > ><br>> > > Cheers,<br>> > ><br>> > > Quentin CHEVALIER – IA parcours recherche<br>> > > LadHyX - Ecole polytechnique<br>> > > __________<br>> > ><br>> > ><br>> > ><br>> > > On Thu, 5 May 2022 at 10:34, Quentin Chevalier<br>> > > <<a href="mailto:quentin.chevalier@polytechnique.edu" target="_blank">quentin.chevalier@polytechnique.edu</a>> wrote:<br>> > >><br>> > >> Hello all and thanks for your great work in bringing this very helpful package to the community !<br>> > >><br>> > >> That said, I wouldn't need this mailing list if everything was running smoothly. I have a rather involved eigenvalue problem that I've been working on that's been throwing a mysterious error :<br>> > >>><br>> > >>> petsc4py.PETSc.Error: error code 77<br>> > >>><br>> > >>> [1] EPSSolve() at /usr/local/slepc/src/eps/interface/epssolve.c:149<br>> > >>> [1] EPSSolve_KrylovSchur_Default() at /usr/local/slepc/src/eps/impls/krylov/krylovschur/krylovschur.c:289<br>> > >>> [1] EPSGetStartVector() at /usr/local/slepc/src/eps/interface/epssolve.c:824<br>> > >>> [1] Petsc has generated inconsistent data<br>> > >>> [1] Initial vector is zero or belongs to the deflation space<br>> > >><br>> > >><br>> > >> This problem occurs in parallel with two processors, using the petsc4py library using the dolfinx/dolfinx docker container. I have PETSc version 3.16.0, in complex mode, python 3, and I'm running all of that on a OpenSUSE Leap 15.2 machine (but I think the docker container has a Ubuntu OS).<br>> > >><br>> > >> I wrote a minimal working example below, but I'm afraid the process for building the matrices is involved, so I decided to directly share the matrices instead : <a href="https://seminaris.polytechnique.fr/share/s/ryJ6L2nR4ketDwP" target="_blank">https://seminaris.polytechnique.fr/share/s/ryJ6L2nR4ketDwP</a><br>> > >><br>> > >> They are in binary format, but inside the container I hope someone could reproduce my issue. A word on the structure and intent behind these matrices :<br>> > >><br>> > >> QE is a diagonal rectangular real matrix. Think of it as some sort of preconditioner<br>> > >> L is the least dense of them all, the only one that is complex, and in order to avoid inverting it I'm using two KSPs to compute solve problems on the fly<br>> > >> Mf is a diagonal square real matrix, its on the right-hand side of the Generalised Hermitian Eigenvalue problem (I'm solving QE^H*L^-1H*L^-1*QE*x=lambda*Mf*x<br>> > >><br>> > >> Full MWE is below :<br>> > >><br>> > >> from petsc4py import PETSc as pet<br>> > >> from slepc4py import SLEPc as slp<br>> > >> from mpi4py.MPI import COMM_WORLD<br>> > >><br>> > >> # Global sizes<br>> > >> m_local=COMM_WORLD.rank*(490363-489780)+489780<br>> > >> n_local=COMM_WORLD.rank*(452259-451743)+451743<br>> > >> m=980143<br>> > >> n=904002<br>> > >><br>> > >> QE=pet.Mat().createAIJ([[m_local,m],[n_local,n]])<br>> > >> L=pet.Mat().createAIJ([[m_local,m],[m_local,m]])<br>> > >> Mf=pet.Mat().createAIJ([[n_local,n],[n_local,n]])<br>> > >><br>> > >> viewerQE = pet.Viewer().createMPIIO("QE.dat", 'r', COMM_WORLD)<br>> > >> QE.load(viewerQE)<br>> > >> viewerL = pet.Viewer().createMPIIO("L.dat", 'r', COMM_WORLD)<br>> > >> L.load(viewerL)<br>> > >> viewerMf = pet.Viewer().createMPIIO("Mf.dat", 'r', COMM_WORLD)<br>> > >> Mf.load(viewerMf)<br>> > >><br>> > >> QE.assemble()<br>> > >> L.assemble()<br>> > >><br>> > >> KSPs = []<br>> > >> # Useful solvers (here to put options for computing a smart R)<br>> > >> for Mat in [L,L.hermitianTranspose()]:<br>> > >> KSP = pet.KSP().create()<br>> > >> KSP.setOperators(Mat)<br>> > >> KSP.setFromOptions()<br>> > >> KSPs.append(KSP)<br>> > >> class LHS_class:<br>> > >> def mult(self,A,x,y):<br>> > >> w=pet.Vec().createMPI([m_local,m],comm=COMM_WORLD)<br>> > >> z=pet.Vec().createMPI([m_local,m],comm=COMM_WORLD)<br>> > >> QE.mult(x,w)<br>> > >> KSPs[0].solve(w,z)<br>> > >> KSPs[1].solve(z,w)<br>> > >> QE.multTranspose(w,y)<br>> > >><br>> > >> # Matrix free operator<br>> > >> LHS=pet.Mat()<br>> > >> LHS.create(comm=COMM_WORLD)<br>> > >> LHS.setSizes([[n_local,n],[n_local,n]])<br>> > >> LHS.setType(pet.Mat.Type.PYTHON)<br>> > >> LHS.setPythonContext(LHS_class())<br>> > >> LHS.setUp()<br>> > >><br>> > >> # Eigensolver<br>> > >> EPS = slp.EPS(); EPS.create()<br>> > >> EPS.setOperators(LHS,Mf) # Solve QE^T*L^-1H*L^-1*QE*x=lambda*Mf*x<br>> > >> EPS.setProblemType(slp.EPS.ProblemType.GHEP) # Specify that A is hermitian (by construction), and B is semi-positive<br>> > >> EPS.setWhichEigenpairs(EPS.Which.LARGEST_MAGNITUDE) # Find largest eigenvalues<br>> > >> EPS.setFromOptions()<br>> > >> EPS.solve()<br>> > >><br>> > >> Quentin CHEVALIER – IA parcours recherche<br>> > >><br>> > >> LadHyX - Ecole polytechnique<br>> > >><br>> > >> __________<br>> ><br>> > <image003.jpg><MWE.py><br>><br></div></div>
</blockquote></div><br clear="all"><div><br></div>-- <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>