<div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr">Okay i think i found my error. I now used this code snippet:</div><div dir="ltr"><div dir="ltr">if nconv > 0:</div><div dir="ltr"><span style="white-space:pre"> </span># Create the results vectors</div><div dir="ltr"><span style="white-space:pre"> </span>vr, wr = DynMatrix_nn.getVecs()</div><div dir="ltr"><span style="white-space:pre"> </span>vi, wi = DynMatrix_nn.getVecs()</div><div dir="ltr"><span style="white-space:pre"> </span></div><div dir="ltr"><span style="white-space:pre"> </span># Array for the eigenvectos</div><div dir="ltr"><span style="white-space:pre"> </span>vecs = np.zeros((69999,nconv))<span style="white-space:pre"> </span></div><div dir="ltr"><span style="white-space:pre"> </span>print("Shape of output array: ", vecs.shape)<span style="white-space:pre"> </span></div><div dir="ltr"><br></div><div dir="ltr"><span style="white-space:pre"> </span>Print()</div><div dir="ltr"><span style="white-space:pre"> </span>Print(" k ||Ax-kx||/||kx|| ")</div><div dir="ltr"><span style="white-space:pre"> </span>Print("-----------------------------------------")</div><div dir="ltr"><span style="white-space:pre"> </span>for i in range(nconv):</div><div dir="ltr"><span style="white-space:pre"> </span>k = E.getEigenpair(i, vr, vi)</div><div dir="ltr"><span style="white-space:pre"> </span>error = E.computeError(i)</div><div dir="ltr"><span style="white-space:pre"> </span>if k.imag != 0.0:</div><div dir="ltr"><span style="white-space:pre"> </span>Print("%9f%+9f j %12g" % (k.real, k.imag, error))</div><div dir="ltr"><span style="white-space:pre"> </span>else:</div><div dir="ltr"><span style="white-space:pre"> </span>Print("%12f %12g" %(k.real,error))</div><div dir="ltr"><span style="white-space:pre"> </span># Save all eigenvalues to one array<span style="white-space:pre"> </span></div><div dir="ltr"><span style="white-space:pre"> </span>eigenvals_n[i,0] = E.getEigenvalue(i).real</div><div dir="ltr"><span style="white-space:pre"> </span>eigenvals_n[i,1] = error </div><div dir="ltr"><span style="white-space:pre"> </span></div><div dir="ltr"><span style="white-space:pre"> </span># Insert eigenvectors <span style="white-space:pre"> </span></div><div dir="ltr"><span style="white-space:pre"> </span>scatter, vec_full = PETSc.Scatter().toZero(vr)</div><div dir="ltr"><span style="white-space:pre"> </span>scatter.scatter(vr, vec_full, PETSc.InsertMode.INSERT, PETSc.ScatterMode.FORWARD)</div><div dir="ltr"><span style="white-space:pre"> </span>if rank == 0:<span style="white-space:pre"> </span></div><div dir="ltr"><span style="white-space:pre"> </span>vecs[:,i] = vec_full.getArray()<span style="white-space:pre"> </span></div><div dir="ltr"><div dir="ltr">E.view()</div><div dir="ltr"># Save eigenvalues to file </div><div dir="ltr">if rank == 0:</div><div dir="ltr"><span style="white-space:pre"> </span>np.savetxt("Eigenvalues.txt", eigenvals_n, header="Eigenvalues, ascending ordering due to real part")</div><div dir="ltr"><span style="white-space:pre"> </span>np.savetxt("Eigenvec2.txt", vecs, header="Eigenvalues, ascending ordering due to real part")</div><div>The problem was that all MPI processes create the array vecs. The scattering worked fine but in the end every MPI process wrote the eigenvec2.txt file and not only rank0. So i assume that 19 processes writing an empty array and 1 process writing the correct one, was leading to this strange kind of .txt file. Now only rank0 with the correct eigenvectors write the array to file. I still have to check if the rest is correct but i appears to me that it now works and it makes sense.</div><div>Thank you very much for your help Matt!!</div></div><div dir="ltr"><br></div></div></div></div></div></div></div><br><div class="gmail_quote"><div dir="ltr">Am Di., 18. Sep. 2018 um 18:19 Uhr schrieb Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>>:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_quote"><div dir="ltr">On Tue, Sep 18, 2018 at 11:53 AM Jan Grießer <<a href="mailto:griesser.jan@googlemail.com" target="_blank">griesser.jan@googlemail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div dir="ltr"><div dir="ltr">But i'am still struggling a little bit with copying the vector entries from the processes to one sequentiel proces. I tried to use the comment with Scatter().toZero() in Python and tried the following code:<div><div>vecs = np.zeros((69999,nconv))</div><div><span style="white-space:pre-wrap"> </span>for i in range(nconv):</div><div><span style="white-space:pre-wrap"> </span>val = E.getEigenpair(i, vr, vi)</div><div><span style="white-space:pre-wrap"> </span>scatter, vec_full = PETSc.Scatter().toZero(vr)</div><div><span style="white-space:pre-wrap"> </span>scatter.scatter(vr, vec_full, False, PETSc.ScatterMode.FORWARD)</div><div><span style="white-space:pre-wrap"> </span>if rank == 0:</div><div><span style="white-space:pre-wrap"> </span>vecs[:,i] = vec_full.getArray()<span style="white-space:pre-wrap"> </span></div></div><div>but when i take a look at the eigenvectors there are some entries which are strange, because in this exampled i printed 4 eigenvalues and vectors and get nconv=4, but when i print the vectors i get lines </div><div>with 4 entries but they are all zero and also get lines with strange entries e.g.</div><div><div>6.386079897185445800e-05 3.675310683411162830e-04 9.068495707804999221e-06 7.43627600.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00</div></div><div>where i have more than 4 entries. Can somebody explain to me why this happens or what is wrong with the code. (Sorry for asking such beginner questions, but i just started to work with PETSc and SLEPc)</div></div></div></div></blockquote><div><br></div><div>We can definitely figure that out. However, can you try</div><div><br></div><div> vr.view(PETSC_VIEWER_STDOUT_WORLD)</div><div><br></div><div>just to check your vectors initially?</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="gmail_quote"><div dir="ltr">Am Di., 18. Sep. 2018 um 09:26 Uhr schrieb Jan Grießer <<a href="mailto:griesser.jan@googlemail.com" target="_blank">griesser.jan@googlemail.com</a>>:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Now it works. Thank you very much!</div><br><div class="gmail_quote"><div dir="ltr">Am Mo., 17. Sep. 2018 um 10:54 Uhr schrieb Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>>:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">For this you have to set a number of partitions equal to the number of MPI processes. If you get this error it is because these values are different.<br>
Jose<br>
<br>
<br>
> El 17 sept 2018, a las 10:50, Jan Grießer <<a href="mailto:griesser.jan@googlemail.com" target="_blank">griesser.jan@googlemail.com</a>> escribió:<br>
> <br>
> Is this relly necessary, because in the last sentences of the chapter it states that:<br>
> An additional benefit of multi-communicator support is that it enables parallel spectrum slicing runs without the need to install a parallel direct solver (MUMPS). The following commandline example uses sequential linear solves in 4 partitions, one process each:<br>
> Therefore i assumed that it is not necessary to compile PETsc4py with an external solver e.g. MUMPS<br>
> <br>
> Am Mo., 17. Sep. 2018 um 10:47 Uhr schrieb Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>>:<br>
> You need a parallel direct solver such as MUMPS. This is explained in section 3.4.5.<br>
> Jose<br>
> <br>
> <br>
> > El 17 sept 2018, a las 10:41, Jan Grießer <<a href="mailto:griesser.jan@googlemail.com" target="_blank">griesser.jan@googlemail.com</a>> escribió:<br>
> > <br>
> > def solve_eigensystem(DynMatrix_nn, Unity_nn, Dimension, LowerLimit, UpperLimit):<br>
> > # Create the EPS solver <br>
> > E = SLEPc.EPS().create()<br>
> > <br>
> > # Create the preconditioner and set it to Cholesky<br>
> > pc = PETSc.PC().create()<br>
> > pc.setType(pc.Type.CHOLESKY)<br>
> > <br>
> > # Create the KSP object<br>
> > ksp = PETSc.KSP().create()<br>
> > ksp.setType(ksp.Type.PREONLY)<br>
> > ksp.setPC(pc)<br>
> > <br>
> > # Set up the spectral transformations <br>
> > st = SLEPc.ST().create()<br>
> > st.setType("sinvert")<br>
> > st.setKSP(ksp)<br>
> > # Setup spectral transformation <br>
> > E.setST(st)<br>
> > <br>
> > # Eigenvalues should be real, therefore we start to order them from the smallest real value |l.real|<br>
> > E.setWhichEigenpairs(E.Which.ALL)<br>
> > # Set the interval of spectrum slicing<br>
> > E.setInterval(LowerLimit, UpperLimit)<br>
> > # Since the dynamical matrix is symmetric and real it is hermitian. Use GHEP for the spectrum slicing. Operatormatrix B is just a unit matrix <br>
> > E.setProblemType(SLEPc.EPS.ProblemType.GHEP)<br>
> > # Use the Krylov Schur method to solve the eigenvalue problem<br>
> > E.setType(E.Type.KRYLOVSCHUR)<br>
> > # Partition the Krylov schnur problem in npart procceses<br>
> > E.setKrylovSchurPartitions(10)<br>
> > # Set the convergence criterion to relative to the eigenvalue and the maximal number of iterations<br>
> > E.setConvergenceTest(E.Conv.REL)<br>
> > E.setTolerances(tol = 1e-7, max_it = 1000)<br>
> > # Set the matrix in order to solve <br>
> > E.setOperators(DynMatrix_nn, Unity_nn)<br>
> > # Sets EPS options from the options database. <br>
> > E.setFromOptions()<br>
> > # Sets up all the internal data structures necessary for the execution of the eigensolver.<br>
> > E.setUp()<br>
> > <br>
> > # Solve eigenvalue problem <br>
> > startClock = time.clock()<br>
> > startTime = time.time()<br>
> > E.solve()<br>
> > <br>
> > Has maybe one of you any idea why this happens and where the problem is ? <br>
> > <br>
> > Am Mo., 17. Sep. 2018 um 10:40 Uhr schrieb Jan Grießer <<a href="mailto:griesser.jan@googlemail.com" target="_blank">griesser.jan@googlemail.com</a>>:<br>
> > I am aware that SLEPc is not supposed to calculate all eigenvalues and eigenvectors, my problem is simply that i want for a physical large enough system all of them before i can make the transition to go to the smallest ones. <br>
> > Competitiveness is of secondary importance at the moment. <br>
> > But ihave a problem connected with spectrum slicing. I followed the instructions in the manual of Chap. 3.4.5 Spectrum Slicing and converted them to the python package. <br>
> > But now i get the following error. It appears to me that it is not able to find the ksp object, but i actually do not know why this is the case. <br>
> > aceback (most recent call last):<br>
> > File "Eigensolver_spectrum_slicing.py", line 216, in <module><br>
> > solve_eigensystem(DynMatrix_nn, Unity_nn, D_nn.shape, opt_dict.LowLimit, opt_dict.UpperLimit)<br>
> > File "Eigensolver_spectrum_slicing.py", line 121, in solve_eigensystem<br>
> > E.setUp()<br>
> > File "SLEPc/EPS.pyx", line 1099, in slepc4py.SLEPc.EPS.setUp<br>
> > petsc4py.PETSc.Error: error code 92<br>
> > [14] EPSSetUp() line 165 in /tmp/pip-install-golhudw7/slepc/src/eps/interface/epssetup.c<br>
> > [14] EPSSetUp_KrylovSchur() line 146 in /tmp/pip-install-golhudw7/slepc/src/eps/impls/krylov/krylovschur/krylovschur.c<br>
> > [14] EPSSetUp_KrylovSchur_Slice() line 410 in /tmp/pip-install-golhudw7/slepc/src/eps/impls/krylov/krylovschur/ks-slice.c<br>
> > [14] EPSSliceGetEPS() line 300 in /tmp/pip-install-golhudw7/slepc/src/eps/impls/krylov/krylovschur/ks-slice.c<br>
> > [14] EPSSetUp() line 165 in /tmp/pip-install-golhudw7/slepc/src/eps/interface/epssetup.c<br>
> > [14] EPSSetUp_KrylovSchur() line 146 in /tmp/pip-install-golhudw7/slepc/src/eps/impls/krylov/krylovschur/krylovschur.c<br>
> > [14] EPSSetUp_KrylovSchur_Slice() line 461 in /tmp/pip-install-golhudw7/slepc/src/eps/impls/krylov/krylovschur/ks-slice.c<br>
> > [14] EPSSliceGetInertia() line 331 in /tmp/pip-install-golhudw7/slepc/src/eps/impls/krylov/krylovschur/ks-slice.c<br>
> > [14] STSetUp() line 271 in /tmp/pip-install-golhudw7/slepc/src/sys/classes/st/interface/stsolve.c<br>
> > [14] STSetUp_Sinvert() line 132 in /tmp/pip-install-golhudw7/slepc/src/sys/classes/st/impls/sinvert/sinvert.c<br>
> > [14] KSPSetUp() line 381 in /tmp/pip-install-xmiaat2t/petsc/src/ksp/ksp/interface/itfunc.c<br>
> > [14] PCSetUp() line 923 in /tmp/pip-install-xmiaat2t/petsc/src/ksp/pc/interface/precon.c<br>
> > [14] PCSetUp_Cholesky() line 86 in /tmp/pip-install-xmiaat2t/petsc/src/ksp/pc/impls/factor/cholesky/cholesky.c<br>
> > [14] MatGetFactor() line 4318 in /tmp/pip-install-xmiaat2t/petsc/src/mat/interface/matrix.c<br>
> > [14] See <a href="http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html</a> for possible LU and Cholesky solvers<br>
> > [14] Could not locate a solver package. Perhaps you must ./configure with --download-<package><br>
> > <br>
> > The code i used to solve the problem is <br>
> > <br>
> > Am Fr., 14. Sep. 2018 um 18:34 Uhr schrieb Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>>:<br>
> > On Fri, Sep 14, 2018 at 12:19 PM Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>> wrote:<br>
> > El 14 sept 2018, a las 17:45, Jan Grießer <<a href="mailto:griesser.jan@googlemail.com" target="_blank">griesser.jan@googlemail.com</a>> escribió:<br>
> > <br>
> >> Hey there,<br>
> >> first i want to say thanks to Satish and Matt for helping with with my last problem with the mpi compilation. I have two questions related to solving a big, hermitian, standard eigenvalue problem using SLEPc4py., compiled with Intel MKL and Intel MPI. <br>
> >> <br>
> >> - I am using slepc4py with <br>
> >> mpi and run it with around -n 20 cores at the moment and how i wanted to ask if there is an easy way to retrieve the eigenvectors? When i run my code and print for i in range(nconv):<br>
> >> for i in range(nconv):<br>
> >> <br>
> >> val = E.<br>
> >> getEigenpair(i, vr<br>
> >> , vi)<br>
> >> Print(<br>
> >> vr.getArray())<br>
> >> i get the parts of the eigenvectors according to the partition of the matrix. Is there any easy way to put them together in an array and write them to file ? (I am struggling a little bit with the building them in the correct order)<br>
> > <br>
> > You need VecScatterCreateToZero. There must be an equivalent in python.<br>
> > <br>
> > An alternative to this which you should consider, because it is simpler, is to write the vector to a file<br>
> > using some format that PETSc understands, Then you just need vr.view(viewer) for a viewer like<br>
> > the binary viewer or some ASCII format you like.<br>
> > <br>
> > Thanks,<br>
> > <br>
> > Matt<br>
> >> - I need to solve eigenvalue problems up to a dimension of 100000 degrees of freedom and i need all eigenvalues and eigenvectors. I think solving all eigenvalues in one process is far too much and i thought about if it is possible to apply the spectrum slicing described in Chap. 3.4.5. Due to the nature of my problem, i am able to simulate smaller systems of 10000 DOF and extract the biggest eigenvalue, which will be the same for larger systems sizes. Is this in general possible since i have a standard HEP problem or is there a better and faster possibility to do this?<br>
> > <br>
> > In general, SLEPc is not intended for computing the whole spectrum. You can try with spectrum slicing but this will be competitive if computing just a percentage of eigenvalues, 50% say. <br>
> > <br>
> > Jose<br>
> > <br>
> >> <br>
> >> Thank you very much!<br>
> > <br>
> > <br>
> > -- <br>
> > 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<br>
> > <br>
> > <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
> <br>
<br>
</blockquote></div>
</blockquote></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="m_7688103653174234089gmail_signature" data-smartmail="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>