Hi Nidish, I may not fully understand your problem, but it sounds like you could benefit from continuation methods. Have you looked into this? If it's helpful, I have some experience with this and I can discuss with you by email.
Cheers, Zak On Mon, Aug 17, 2020 at 2:31 PM Nidish <n...@rice.edu> wrote: > Thankfully for this step, the matrix is not dense. But thank you. > > Nidish > > On 8/17/20 8:52 AM, Jose E. Roman wrote: > > > >> El 17 ago 2020, a las 14:27, Barry Smith <bsm...@petsc.dev> escribió: > >> > >> > >> Nidish, > >> > >> Your matrix is dense, correct? MUMPS is for sparse matrices. > >> > >> Then I guess you could use Scalapack > http://netlib.org/scalapack/slug/node48.html#SECTION04323200000000000000 > to do the SVD. The work is order N^3 and parallel efficiency may not be > great but it might help you solve your problem. > >> > >> I don't know if SLEPc has an interface to Scalapack for SVD or > not. > > Yes, SLEPc (master) has interfaces for ScaLAPACK and Elemental for both > SVD and (symmetric) eigenvalues. > > > > Jose > > > >> Barry > >> > >> > >> > >> > >> > >> > >> > >>> On Aug 17, 2020, at 2:51 AM, Jose E. Roman <jro...@dsic.upv.es> wrote: > >>> > >>> You can use SLEPc's SVD to compute the nullspace, but it has pitfalls: > make sure you use an absolute convergence test (not relative); for the > particular case of zero singular vectors, accuracy may not be very good and > convergence may be slow (with the corresponding high computational cost). > >>> > >>> MUMPS has functionality to get a basis of the nullspace, once you have > computed the factorization. But I don't know if this is easily accessible > via PETSc. > >>> > >>> Jose > >>> > >>> > >>> > >>>> El 17 ago 2020, a las 3:10, Nidish <n...@rice.edu> escribió: > >>>> > >>>> Oh damn. Alright, I'll keep trying out the different options. > >>>> > >>>> Thank you, > >>>> Nidish > >>>> > >>>> On 8/16/20 8:05 PM, Barry Smith wrote: > >>>>> SVD is enormously expensive, needs to be done on a full dense matrix > so completely impractical. You need the best tuned iterative method, Jose > is the by far the most knowledgeable about that. > >>>>> > >>>>> Barry > >>>>> > >>>>> > >>>>>> On Aug 16, 2020, at 7:46 PM, Nidish <n...@rice.edu> wrote: > >>>>>> > >>>>>> Thank you for the suggestions. > >>>>>> > >>>>>> I'm getting a zero pivot error for the LU in slepc while > calculating the rest of the modes. > >>>>>> > >>>>>> Would conducting an SVD for just the stiffness matrix and then > using the singular vectors as bases for the nullspace work? I haven't tried > this out just yet, but I'm wondering if you could provide me insights into > whether this will. > >>>>>> > >>>>>> Thanks, > >>>>>> Nidish > >>>>>> > >>>>>> On 8/16/20 2:50 PM, Barry Smith wrote: > >>>>>>> If you know part of your null space explicitly (for example the > rigid body modes) I would recommend you always use that information > explicitly since it is extremely expensive numerically to obtain. Thus > rather than numerically computing the entire null space compute the part > orthogonal to the part you already know. Presumably SLEPc has tools to help > do this, naively I would just orthogonalized against the know subspace > during the computational process but there are probably better ways. > >>>>>>> > >>>>>>> Barry > >>>>>>> > >>>>>>> > >>>>>>> > >>>>>>> > >>>>>>>> On Aug 16, 2020, at 11:26 AM, Nidish <n...@rice.edu> wrote: > >>>>>>>> > >>>>>>>> Well some of the zero eigenvectors are rigid body modes, but > there are some more which are introduced by lagrange-multiplier based > constraint enforcement, which are non trivial. > >>>>>>>> > >>>>>>>> My final application is for a nonlinear simulation, so I don't > mind the extra computational effort initially. Could you have me the > suggested solver configurations to get this type of eigenvectors in slepc? > >>>>>>>> > >>>>>>>> Nidish > >>>>>>>> On Aug 16, 2020, at 00:17, Jed Brown <j...@jedbrown.org> wrote: > >>>>>>>> It's possible to use this or a similar algorithm in SLEPc, but > keep in mind that it's more expensive to compute these eigenvectors than to > solve a linear system. Do you have a sequence of systems with the same > null space? > >>>>>>>> > >>>>>>>> You referred to the null space as "rigid body modes". Why can't > those be written down? Note that PETSc has convenience routines for > computing rigid body modes from coordinates. > >>>>>>>> > >>>>>>>> Nidish < > >>>>>>>> n...@rice.edu > >>>>>>>>> writes: > >>>>>>>> > >>>>>>>> I just use the standard eigs function ( > https://www.mathworks.com/help/matlab/ref/eigs.html > >>>>>>>> ) as a black box. I think it uses a lanczos type method under the > hood. > >>>>>>>> > >>>>>>>> Nidish > >>>>>>>> > >>>>>>>> On Aug 15, 2020, 21:42, at 21:42, Barry Smith < > >>>>>>>> bsm...@petsc.dev > >>>>>>>>> wrote: > >>>>>>>> > >>>>>>>> Exactly what algorithm are you using in Matlab to get the 10 > smallest > >>>>>>>> eigenvalues and their corresponding eigenvectors? > >>>>>>>> > >>>>>>>> Barry > >>>>>>>> > >>>>>>>> > >>>>>>>> > >>>>>>>> On Aug 15, 2020, at 8:53 PM, Nidish <n...@rice.edu > >>>>>>>>> wrote: > >>>>>>>> The section on solving singular systems in the manual starts with > >>>>>>>> > >>>>>>>> assuming that the singular eigenvectors are already known. > >>>>>>>> > >>>>>>>> > >>>>>>>> I have a large system where finding the singular eigenvectors is > not > >>>>>>>> > >>>>>>>> trivially written down. How would you recommend I proceed with > making > >>>>>>>> initial estimates? In MATLAB (with MUCH smaller matrices), I > conduct an > >>>>>>>> eigensolve for the first 10 smallest eigenvalues and take the > >>>>>>>> eigenvectors corresponding to the zero eigenvalues from this. This > >>>>>>>> approach doesn't work here since I'm unable to use SLEPc for > solving > >>>>>>>> > >>>>>>>> > >>>>>>>> K.v = lam*M.v > >>>>>>>> > >>>>>>>> for cases where K is positive semi-definite (contains a few "rigid > >>>>>>>> > >>>>>>>> body modes") and M is strictly positive definite. > >>>>>>>> > >>>>>>>> > >>>>>>>> I'd appreciate any assistance you may provide with this. > >>>>>>>> > >>>>>>>> Thank you, > >>>>>>>> Nidish > >>>>>>>> > >>>>>> -- > >>>>>> Nidish > >>>> -- > >>>> Nidish > -- > Nidish >