Case a is the simplest one to use and enforces B-orthogonality of eigenvectors. Case b can also be used if you have a good reason to do so, and in that case you can recover symmetry (and B-orthogonality) by explicitly setting the inner product matrix as illustrated in ex47.c https://slepc.upv.es/documentation/current/src/eps/tutorials/ex47.c.html
Jose > El 11 dic 2021, a las 1:10, Wang, Kuang-chung <[email protected]> > escribió: > > • I was able to use MatIsHermitian function successfully since my > matrix type is seqaij. Just hoped that it can return MatNorm (H -H^+) with > this function. But it wasn’t hard for me to implement that with the function > listed in the previous email. > • I have a related question. In user manual > https://slepc.upv.es/documentation/slepc.pdf 2.1, it says that Ax=\lambda B > x problem is usually reformulated to B^-1 Ax =\lambda x. If A matrix is > Hermitian, B is diagonal but Bii and Bjj can be different. > a) Will solving “Ax=\lambda B x” directly with slepc guarantees users > receiving orthogonal eigenvectors? Namely, does xi^T*B*xj=delta_ij hold true? > > b) if we reformulate, (B^-1A) x = \lambda x will yield (B^-1 A) to be > non-hermitian and therefore doesn’t give orthogonal eigenvectors ( pointed > out by Jose). Xi^T*xj != delta_ij. What about xi^T*B*xj=delta_ij, will this > be guaranteed(since this is the same problem as “a” )? Currently, my test is > that xi^T*B*xj=delta_ij is no longer true. > To help with visibility: > <image002.jpg> > > Although a and b are the same problem formulated differently, but the > orthogonality isn’t ensured in the case b while in case a is ensured(?) . > If so, does it mean that we should be encouraged to use case a ? > > Best, > Kuang > > > From: Zhang, Hong <[email protected]> > Sent: Thursday, December 2, 2021 2:18 PM > To: Wang, Kuang-chung <[email protected]>; Jose E. Roman > <[email protected]> > Cc: [email protected]; Obradovic, Borna <[email protected]>; > Cea, Stephen M <[email protected]> > Subject: Re: [petsc-users] Orthogonality of eigenvectors in SLEPC > > Kuang, > PETSc supports MatIsHermitian() for SeqAIJ, IS and SeqSBAIJ matrix types. > What is your matrix type? > We should be able to add this support to other mat types. > Hong > From: petsc-users <[email protected]> on behalf of Wang, > Kuang-chung <[email protected]> > Sent: Thursday, December 2, 2021 2:06 PM > To: Jose E. Roman <[email protected]> > Cc: [email protected] <[email protected]>; Obradovic, Borna > <[email protected]>; Cea, Stephen M <[email protected]> > Subject: Re: [petsc-users] Orthogonality of eigenvectors in SLEPC > > Thanks Jose for your prompt reply. > I did find my matrix highly non-hermitian. By forcing the solver to be > hermtian, the orthogonality was restored. > But I do need to root cause why my matrix is non-hermitian in the first > place. > Along the way, I highly recommend MatIsHermitian() function or combining > functions like MatHermitianTranspose () MatAXPY MatNorm to determine the > hermiticity to safeguard our program. > > Best, > Kuang > > -----Original Message----- > From: Jose E. Roman <[email protected]> > Sent: Wednesday, November 24, 2021 6:20 AM > To: Wang, Kuang-chung <[email protected]> > Cc: [email protected]; Obradovic, Borna <[email protected]>; > Cea, Stephen M <[email protected]> > Subject: Re: [petsc-users] Orthogonality of eigenvectors in SLEPC > > In Hermitian eigenproblems orthogonality of eigenvectors is > guaranteed/enforced. But you are solving the problem as non-Hermitian. > > If your matrix is Hermitian, make sure you solve it as a HEP, and make sure > that your matrix is numerically Hermitian. > > If your matrix is non-Hermitian, then you cannot expect the eigenvectors to > be orthogonal. What you can do in this case is get an orthogonal basis of the > computed eigenspace, > seehttps://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSGetInvariantSubspace.html > > > By the way, version 3.7 is more than 5 years old, it is better if you can > upgrade to a more recent version. > > Jose > > > > > El 24 nov 2021, a las 7:15, Wang, Kuang-chung <[email protected]> > > escribió: > > > > Dear Jose : > > I came across this thread describing issue using krylovschur and finding > > eigenvectors non-orthogonal. > > https://lists.mcs.anl.gov/pipermail/petsc-users/2014-October/023360.ht > > ml > > > > I furthermore have tested by reducing the tolerance as highlighted below > > from 1e-12 to 1e-16 with no luck. > > Could you please suggest options/sources to try out ? > > Thanks a lot for sharing your knowledge! > > > > Sincere, > > Kuang-Chung Wang > > > > ======================================================= > > Kuang-Chung Wang > > Computational and Modeling Technology > > Intel Corporation > > Hillsboro OR 97124 > > ======================================================= > > > > Here are more info: > > • slepc/3.7.4 > > • output message from by doing EPSView(eps,PETSC_NULL): > > EPS Object: 1 MPI processes > > type: krylovschur > > Krylov-Schur: 50% of basis vectors kept after restart > > Krylov-Schur: using the locking variant > > problem type: non-hermitian eigenvalue problem > > selected portion of the spectrum: closest to target: 20.1161 (in > > magnitude) > > number of eigenvalues (nev): 40 > > number of column vectors (ncv): 81 > > maximum dimension of projected problem (mpd): 81 > > maximum number of iterations: 1000 > > tolerance: 1e-12 > > convergence test: relative to the eigenvalue BV Object: 1 MPI > > processes > > type: svec > > 82 columns of global length 2988 > > vector orthogonalization method: classical Gram-Schmidt > > orthogonalization refinement: always > > block orthogonalization method: Gram-Schmidt > > doing matmult as a single matrix-matrix product DS Object: 1 MPI > > processes > > type: nhep > > ST Object: 1 MPI processes > > type: sinvert > > shift: 20.1161 > > number of matrices: 1 > > KSP Object: (st_) 1 MPI processes > > type: preonly > > maximum iterations=1000, initial guess is zero > > tolerances: relative=1.12005e-09, absolute=1e-50, divergence=10000. > > left preconditioning > > using NONE norm type for convergence test > > PC Object: (st_) 1 MPI processes > > type: lu > > LU: out-of-place factorization > > tolerance for zero pivot 2.22045e-14 > > matrix ordering: nd > > factor fill ratio given 0., needed 0. > > Factored matrix follows: > > Mat Object: 1 MPI processes > > type: seqaij > > rows=2988, cols=2988 > > package used to perform factorization: mumps > > total: nonzeros=614160, allocated nonzeros=614160 > > total number of mallocs used during MatSetValues calls =0 > > MUMPS run parameters: > > SYM (matrix type): 0 > > PAR (host participation): 1 > > ICNTL(1) (output for error): 6 > > ICNTL(2) (output of diagnostic msg): 0 > > ICNTL(3) (output for global info): 0 > > ICNTL(4) (level of printing): 0 > > ICNTL(5) (input mat struct): 0 > > ICNTL(6) (matrix prescaling): 7 > > ICNTL(7) (sequential matrix ordering):7 > > ICNTL(8) (scaling strategy): 77 > > ICNTL(10) (max num of refinements): 0 > > ICNTL(11) (error analysis): 0 > > ICNTL(12) (efficiency control): 1 > > ICNTL(13) (efficiency control): 0 > > ICNTL(14) (percentage of estimated workspace increase): 20 > > ICNTL(18) (input mat struct): 0 > > ICNTL(19) (Schur complement info): 0 > > ICNTL(20) (rhs sparse pattern): 0 > > ICNTL(21) (solution struct): 0 > > ICNTL(22) (in-core/out-of-core facility): 0 > > ICNTL(23) (max size of memory can be allocated locally):0 > > ICNTL(24) (detection of null pivot rows): 0 > > ICNTL(25) (computation of a null space basis): 0 > > ICNTL(26) (Schur options for rhs or solution): 0 > > ICNTL(27) (experimental parameter): -24 > > ICNTL(28) (use parallel or sequential ordering): 1 > > ICNTL(29) (parallel ordering): 0 > > ICNTL(30) (user-specified set of entries in inv(A)): 0 > > ICNTL(31) (factors is discarded in the solve phase): 0 > > ICNTL(33) (compute determinant): 0 > > CNTL(1) (relative pivoting threshold): 0.01 > > CNTL(2) (stopping criterion of refinement): 1.49012e-08 > > CNTL(3) (absolute pivoting threshold): 0. > > CNTL(4) (value of static pivoting): -1. > > CNTL(5) (fixation for null pivots): 0. > > RINFO(1) (local estimated flops for the elimination after > > analysis): > > [0] 8.15668e+07 > > RINFO(2) (local estimated flops for the assembly after > > factorization): > > [0] 892584. > > RINFO(3) (local estimated flops for the elimination after > > factorization): > > [0] 8.15668e+07 > > INFO(15) (estimated size of (in MB) MUMPS internal data for > > running numerical factorization): > > [0] 16 > > INFO(16) (size of (in MB) MUMPS internal data used during > > numerical factorization): > > [0] 16 > > INFO(23) (num of pivots eliminated on this processor after > > factorization): > > [0] 2988 > > RINFOG(1) (global estimated flops for the elimination after > > analysis): 8.15668e+07 > > RINFOG(2) (global estimated flops for the assembly after > > factorization): 892584. > > RINFOG(3) (global estimated flops for the elimination after > > factorization): 8.15668e+07 > > (RINFOG(12) RINFOG(13))*2^INFOG(34) (determinant): > > (0.,0.)*(2^0) > > INFOG(3) (estimated real workspace for factors on all > > processors after analysis): 614160 > > INFOG(4) (estimated integer workspace for factors on all > > processors after analysis): 31971 > > INFOG(5) (estimated maximum front size in the complete > > tree): 246 > > INFOG(6) (number of nodes in the complete tree): 197 > > INFOG(7) (ordering option effectively use after analysis): 2 > > INFOG(8) (structural symmetry in percent of the permuted > > matrix after analysis): 100 > > INFOG(9) (total real/complex workspace to store the matrix > > factors after factorization): 614160 > > INFOG(10) (total integer space store the matrix factors > > after factorization): 31971 > > INFOG(11) (order of largest frontal matrix after > > factorization): 246 > > INFOG(12) (number of off-diagonal pivots): 0 > > INFOG(13) (number of delayed pivots after factorization): 0 > > INFOG(14) (number of memory compress after factorization): 0 > > INFOG(15) (number of steps of iterative refinement after > > solution): 0 > > INFOG(16) (estimated size (in MB) of all MUMPS internal > > data for factorization after analysis: value on the most memory consuming > > processor): 16 > > INFOG(17) (estimated size of all MUMPS internal data for > > factorization after analysis: sum over all processors): 16 > > INFOG(18) (size of all MUMPS internal data allocated during > > factorization: value on the most memory consuming processor): 16 > > INFOG(19) (size of all MUMPS internal data allocated during > > factorization: sum over all processors): 16 > > INFOG(20) (estimated number of entries in the factors): > > 614160 > > INFOG(21) (size in MB of memory effectively used during > > factorization - value on the most memory consuming processor): 14 > > INFOG(22) (size in MB of memory effectively used during > > factorization - sum over all processors): 14 > > INFOG(23) (after analysis: value of ICNTL(6) effectively > > used): 0 > > INFOG(24) (after analysis: value of ICNTL(12) effectively > > used): 1 > > INFOG(25) (after factorization: number of pivots modified > > by static pivoting): 0 > > INFOG(28) (after factorization: number of null pivots > > encountered): 0 > > INFOG(29) (after factorization: effective number of entries > > in the factors (sum over all processors)): 614160 > > INFOG(30, 31) (after solution: size in Mbytes of memory > > used during solution phase): 13, 13 > > INFOG(32) (after analysis: type of analysis done): 1 > > INFOG(33) (value used for ICNTL(8)): 7 > > INFOG(34) (exponent of the determinant if determinant is > > requested): 0 > > linear system matrix = precond matrix: > > Mat Object: 1 MPI processes > > type: seqaij > > rows=2988, cols=2988 > > total: nonzeros=151488, allocated nonzeros=151488 > > total number of mallocs used during MatSetValues calls =0 > > using I-node routines: found 996 nodes, limit used is 5 >
