<div dir="ltr">Matrix to matrix products are taking much longer than expected... My snippet is below. m and n are quite large, >1M each. I'm running this on 35 procs. As you can see U, S and V are quite sparse SVD matrices (only their first 4 columns are dense, plus a chop). I expected therefore approximate R to have only rank 4 and computations to run smoothly once the cross products between U and V are computed... Right now my bottle neck is not preconditioner but getting that approximation of M2. What do you think ?<br><div><br></div><div><div style="color:rgb(212,212,212);background-color:rgb(30,30,30);font-family:"Droid Sans Mono","monospace",monospace;font-weight:normal;font-size:14px;line-height:19px;white-space:pre"><div></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)">approximate</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)">k</span><span style="color:rgb(212,212,212)">:</span><span style="color:rgb(78,201,176)">int</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</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)">self</span><span style="color:rgb(212,212,212)">.</span><span style="color:rgb(156,220,254)">R</span><span style="color:rgb(212,212,212)">.getSize()</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)">n_local</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)">R</span><span style="color:rgb(212,212,212)">.getLocalSize()</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(79,193,255)">I</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(79,193,255)">J</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)">tmp3</span><span style="color:rgb(212,212,212)">.</span><span style="color:rgb(156,220,254)">vector</span><span style="color:rgb(212,212,212)">.getOwnershipRange()</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(79,193,255)">S</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().create(</span><span style="color:rgb(156,220,254)">comm</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)">S</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(212,212,212)">        </span><span style="color:rgb(79,193,255)">U</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().create(</span><span style="color:rgb(156,220,254)">comm</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)">U</span><span style="color:rgb(212,212,212)">.setSizes([[</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(212,212,212)">        </span><span style="color:rgb(79,193,255)">V</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().create(</span><span style="color:rgb(156,220,254)">comm</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)">V</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(212,212,212)">        </span><span style="color:rgb(197,134,192)">for</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(79,193,255)">A</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)">U</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(79,193,255)">S</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(79,193,255)">V</span><span style="color:rgb(212,212,212)">): </span><span style="color:rgb(79,193,255)">A</span><span style="color:rgb(212,212,212)">.setUp()</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(156,220,254)">g</span><span style="color:rgb(212,212,212)">=</span><span style="color:rgb(78,201,176)">np</span><span style="color:rgb(212,212,212)">.</span><span style="color:rgb(220,220,170)">loadtxt</span><span style="color:rgb(212,212,212)">(</span><span style="color:rgb(206,145,120)">"gains/txt/dummy.txt"</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(197,134,192)">for</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(156,220,254)">i</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(78,201,176)">range</span><span style="color:rgb(212,212,212)">(</span><span style="color:rgb(156,220,254)">k</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)">S</span><span style="color:rgb(212,212,212)">.setValue(</span><span style="color:rgb(156,220,254)">i</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">i</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(156,220,254)">g</span><span style="color:rgb(212,212,212)">[</span><span style="color:rgb(156,220,254)">i</span><span style="color:rgb(212,212,212)">])</span></div><div><span style="color:rgb(212,212,212)">            </span><span style="color:rgb(220,220,170)">loadStuff</span><span style="color:rgb(212,212,212)">(</span><span style="color:rgb(206,145,120)">"left-vecs/"</span><span style="color:rgb(212,212,212)">, params</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)">tmp3</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)">U</span><span style="color:rgb(212,212,212)">.setValues(</span><span style="color:rgb(156,220,254)">self</span><span style="color:rgb(212,212,212)">.reindex</span><span style="color:rgb(212,212,212)">,[</span><span style="color:rgb(156,220,254)">i</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)">tmp3</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)">.array)</span></div><div><span style="color:rgb(212,212,212)">            </span><span style="color:rgb(220,220,170)">loadStuff</span><span style="color:rgb(212,212,212)">(</span><span style="color:rgb(206,145,120)">"right-vecs/"</span><span style="color:rgb(212,212,212)">,params</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)">tmp3</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)">V</span><span style="color:rgb(212,212,212)">.setValues(</span><span style="color:rgb(78,201,176)">range</span><span style="color:rgb(212,212,212)">(</span><span style="color:rgb(79,193,255)">I</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(79,193,255)">J</span><span style="color:rgb(212,212,212)">),  [</span><span style="color:rgb(156,220,254)">i</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)">tmp3</span><span style="color:rgb(212,212,212)">.</span><span style="color:rgb(156,220,254)">vector</span><span style="color:rgb(212,212,212)">)</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(197,134,192)">for</span><span style="color:rgb(212,212,212)"> </span><span style="color:rgb(79,193,255)">A</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)">U</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(79,193,255)">S</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(79,193,255)">V</span><span style="color:rgb(212,212,212)">): </span><span style="color:rgb(79,193,255)">A</span><span style="color:rgb(212,212,212)">.assemble()</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(79,193,255)">U</span><span style="color:rgb(212,212,212)">.chop(</span><span style="color:rgb(181,206,168)">1e-6</span><span style="color:rgb(212,212,212)">); </span><span style="color:rgb(79,193,255)">V</span><span style="color:rgb(212,212,212)">.chop(</span><span style="color:rgb(181,206,168)">1e-6</span><span style="color:rgb(212,212,212)">)</span></div><br><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(79,193,255)">V</span><span style="color:rgb(212,212,212)">.hermitianTranspose(</span><span style="color:rgb(79,193,255)">V</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)">self</span><span style="color:rgb(212,212,212)">.</span><span style="color:rgb(156,220,254)">P</span><span style="color:rgb(212,212,212)">.matMult(</span><span style="color:rgb(79,193,255)">U</span><span style="color:rgb(212,212,212)">.matMult(</span><span style="color:rgb(79,193,255)">S</span><span style="color:rgb(212,212,212)">.matMult(</span><span style="color:rgb(79,193,255)">V</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(79,193,255)">S</span><span style="color:rgb(212,212,212)">),</span><span style="color:rgb(79,193,255)">U</span><span style="color:rgb(212,212,212)">),</span><span style="color:rgb(79,193,255)">U</span><span style="color:rgb(212,212,212)">) </span><span style="color:rgb(106,153,85)"># Multiplications in place (everyone is square except U)</span></div><div><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)">M2</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)">N</span><span style="color:rgb(212,212,212)">.duplicate()</span></div><div><span style="color:rgb(212,212,212)">        </span><span style="color:rgb(79,193,255)">M2</span><span style="color:rgb(212,212,212)">.matMult(</span><span style="color:rgb(79,193,255)">U</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(79,193,255)">U</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)">U</span><span style="color:rgb(212,212,212)">.hermitianTranspose().matMult(</span><span style="color:rgb(79,193,255)">U</span><span style="color:rgb(212,212,212)">,</span><span style="color:rgb(79,193,255)">M2</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(212,212,212)">        </span><span style="color:rgb(197,134,192)">return</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)">M</span><span style="color:rgb(212,212,212)">+</span><span style="color:rgb(79,193,255)">M2</span></div></div></div><br>Quentin<br><br><br>On Tue, 18 Jul 2023 at 07:48, Quentin Chevalier <<a href="mailto:quentin.chevalier@polytechnique.edu">quentin.chevalier@polytechnique.edu</a>> wrote:<br>><br>> Thank you for that pointer ! I have on hand a partial SVD of R, so I used that to build the approximate matrix instead.<br>><br>> It's great that so many nice features of PETSc like STSetPreconditionerMat are accessible through petsc4py !<br>><br>> Good day,<br>><br>> Quentin<br>><br>><br>><br>> On Mon, 17 Jul 2023 at 17:29, Jose E. Roman <<a href="mailto:jroman@dsic.upv.es">jroman@dsic.upv.es</a>> wrote:<br>> ><br>> > It is possible to pass a different matrix to build the preconditioner. That is, the shell matrix for B (EPSSetOperators) and an explicit matrix (that approximates B) for the preconditioner. For instance, you can try passing M for building the preconditioner. Since M is an explicit matrix, you can try the default preconditioner (block Jacobi with ILU as local solver) or even a full LU decomposition. The effectiveness of the preconditioner will depend on how the update M+R^H P M P R moves the eigenvalues around.<br>> ><br>> > You can do this with STSetSplitPreconditioner() or STSetPreconditionerMat(). In your case any of them will do.<br>> ><br>> > Jose<br>> ><br>> ><br>> > > El 17 jul 2023, a las 15:50, Quentin Chevalier <<a href="mailto:quentin.chevalier@polytechnique.edu">quentin.chevalier@polytechnique.edu</a>> escribió:<br>> > ><br>> > > Thank you for this suggestion, I tried to implement that but it's<br>> > > proven pretty hard to implement MATOP_GET_DIAGONAL without completely<br>> > > tanking performance. After all, B is a shell matrix for a reason : it<br>> > > looks like M+R^H P M P R with R itself a shell matrix.<br>> > ><br>> > > Allow me to point out that I have no shift. My eigenvalue problem is<br>> > > purely about the largest ones out there. Section 8.2 and 3.4.3 led me<br>> > > to think that there was a way to avoid computing (or writing a shell<br>> > > matrix about it) B^-1... But you seem to stress that there's no way<br>> > > around it.<br>> > ><br>> > > Quentin<br>> > ><br>> > ><br>> > ><br>> > > On Mon, 17 Jul 2023 at 11:56, Jose E. Roman <<a href="mailto:jroman@dsic.upv.es">jroman@dsic.upv.es</a>> wrote:<br>> > >><br>> > >> The B-inner product is independent of the ST operator. See Table 3.2. In generalized eigenproblems you always have an inverse.<br>> > >><br>> > >> If your matrix is diagonally dominant, try implementing the MATOP_GET_DIAGONAL operation and using PCJACOBI. Apart from this, you have to build your own preconditioner.<br>> > >><br>> > >> Jose<br>> > >><br>> > >><br>> > >>> El 17 jul 2023, a las 11:48, Quentin Chevalier <<a href="mailto:quentin.chevalier@polytechnique.edu">quentin.chevalier@polytechnique.edu</a>> escribió:<br>> > >>><br>> > >>> Hello Jose,<br>> > >>><br>> > >>> I guess I expected B to not be inverted but instead used as a mass for a problem-specific inner product since I specified GHEP as a problem type. p50 of the same user manual seems to imply that that would indeed be the case. I don't see what problem there would be with using a shell B matrix as a weighting matrix, as long as a mat utility is provided of course.<br>> > >>><br>> > >>> I tried the first approach - I set up my KSP as CG since B is hermitian positive-definite (I made a mistake in my first email), but I'm getting a KSPSolve has not converged, reason DIVERGED_ITS error. I'm letting it run for 1000 iterations already so it seems suspiciously slow for a CG solver.<br>> > >>><br>> > >>> I'm grappling with a shell preconditioner now to try and speed it up, but I'm unsure which one allows for shell matrices.<br>> > >>><br>> > >>> Thank you for your time,<br>> > >>><br>> > >>> Quentin<br>> > >>><br>> > >>><br>> > >>> On Wed, 12 Jul 2023 at 19:24, Jose E. Roman <<a href="mailto:jroman@dsic.upv.es">jroman@dsic.upv.es</a>> wrote:<br>> > >>>><br>> > >>>> By default, it is solving the problem as B^{-1}*A*x=lambda*x (see chapter on Spectral Transformation). That is why A can be a shell matrix without problem. But B needs to be an explicit matrix in order to compute an LU factorization. If B is also a shell matrix then you should set an iterative solver for the associated KSP (see examples in the chapter).<br>> > >>>><br>> > >>>> An alternative is to create a shell matrix M that computes the action of B^{-1}*A, then pass M to the EPS solver as a standard eigenproblem.<br>> > >>>><br>> > >>>> Jose<br>> > >>>><br>> > >>>><br>> > >>>>> El 12 jul 2023, a las 19:04, Quentin Chevalier <<a href="mailto:quentin.chevalier@polytechnique.edu">quentin.chevalier@polytechnique.edu</a>> escribió:<br>> > >>>>><br>> > >>>>> Hello PETSc Users,<br>> > >>>>><br>> > >>>>> I have a generalised eigenvalue problem : Ax= lambda Bx<br>> > >>>>> I used to have only A as a matrix-free method, I used mumps and an LU preconditioner, everything worked fine.<br>> > >>>>><br>> > >>>>> Now B is matrix-free as well, and my solver is returning an error : "MatSolverType mumps does not support matrix type python", which is ironic given it seem to handle A quite fine.<br>> > >>>>><br>> > >>>>> I have read in the user manual here that there some methods may require additional methods to be supplied for B like MATOP_GET_DIAGONAL but it's unclear to me exactly what I should be implementing and what is the best solver for my case.<br>> > >>>>><br>> > >>>>> A is hermitian, B is hermitian positive but not positive-definite or real. Therefore I have specified a GHEP problem type to the EPS object.<br>> > >>>>><br>> > >>>>> I use PETSc in complex mode through the petsc4py bridge.<br>> > >>>>><br>> > >>>>> Any help on how to get EPS to work for a generalised matrix-free case would be welcome. Performance is not a key issue here - I have a tractable high value case on hand.<br>> > >>>>><br>> > >>>>> Thank you for your time,<br>> > >>>>><br>> > >>>>> Quentin<br>> > >>>><br>> > >><br>> ></div>