# Error calculating eigs of pentadiagonal matrix.

8 messages
Open this post in threaded view
|

## Error calculating eigs of pentadiagonal matrix.

 Hello.I am working with a pentadiagonal sparse matrix that represents a 2D Schrodinger's time-independent equation. I first work the laplacian expressed in Finite Differences form and then I apply the potential on the same matrix.So far I've been able to validate my results for both an electron in a box as well as a harmonic oscillator, but when I change to the following potential of a dipole, Julia pretty much quits on me when I try to obtain the eigenvalues and eigenvectors:  O         = [round(L/2)-hx round(L/2)-hy]# --el ORIGEN (centro -- x,y) del potencial.  Eps_o  = 8.854187817e10-12# --F*m^-1  C         = 1/(4*pi*Eps_o)  D         = 1e-21#C*m^2/s# --Debyes)  pe        = 1.8*D   P(X,Y)  = -(C)*pe*(Y/X)*(1/( (X)^2 + (Y)^2    )   )# --How the potential gets described. #--I'm aware there's singularities in the potential.#--and here's how I apply the potential to my sparse matrix.  Vi      = Float64[]# --container for the potential.  for j=Y  for i=X  push!(Vi,P(i,j))   end  end@ --applying the potential.I use this command: l, v = eigs(M,nev=15,which = :SM ,ritzvec=true)My problem seems to be that there's an error that I can't get past:ERROR: LoadError: ArgumentError: matrix has one or more zero pivots in #ldltfact!#10(::Float64, ::Function, ::Base.SparseArrays.CHOLMOD.Factor{Float64}, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1350 in (::Base.LinAlg.#kw##ldltfact!)(::Array{Any,1}, ::Base.LinAlg.#ldltfact!, ::Base.SparseArrays.CHOLMOD.Factor{Float64}, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./:0 in #ldltfact#12(::Float64, ::Array{Int64,1}, ::Function, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1386 in #ldltfact#13(::Array{Any,1}, ::Function, ::Hermitian{Float64,SparseMatrixCSC{Float64,Int64}}) at ./sparse/cholmod.jl:1426 in factorize(::SparseMatrixCSC{Float64,Int64}) at ./sparse/linalg.jl:897 in #_eigs#62(::Int64, ::Int64, ::Symbol, ::Float64, ::Int64, ::Void, ::Array{Float64,1}, ::Bool, ::Base.LinAlg.#_eigs, ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at ./linalg/arnoldi.jl:251 in (::Base.LinAlg.#kw##_eigs)(::Array{Any,1}, ::Base.LinAlg.#_eigs, ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at ./:0 in #eigs#55(::Array{Any,1}, ::Function, ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at ./linalg/arnoldi.jl:78 in (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs, ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at ./:0 in #eigs#54(::Array{Any,1}, ::Function, ::SparseMatrixCSC{Float64,Int64}) at ./linalg/arnoldi.jl:77 in (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs, ::SparseMatrixCSC{Float64,Int64}) at ./:0 in include_from_node1(::String) at ./loading.jl:488while loading /home/alejandro/Desktop/ACAD/PROG/ACADEMIC_PROGRAMMING/FDM (Finite_Difference_Method)/2D/SCHROD_DIP_2D/SCRATCH-2D-SI-DIPOLO2.jl, in expression starting on line 106My question is, is there a way to work around it, or, am I completely screwed?Thanks so much in advance.
Open this post in threaded view
|

## Re: Error calculating eigs of pentadiagonal matrix.

 I have one suggestion: see if Andreas' LinearAlgebra.jl does any better.`Pkg.clone(""https://github.com/andreasnoack/LinearAlgebra.jl")using LinearAlgebra# ...`If there is still difficulty, you may look at eigenGeneral.jl and eigenSelfAdjoint.jl to find which more specific versions of eigvals/eigvals! is appropriate for your scenario.On Wednesday, November 2, 2016 at 6:43:48 AM UTC-4, Alejandro Castellanos wrote:Hello.I am working with a pentadiagonal sparse matrix that represents a 2D Schrodinger's time-independent equation. I first work the laplacian expressed in Finite Differences form and then I apply the potential on the same matrix.So far I've been able to validate my results for both an electron in a box as well as a harmonic oscillator, but when I change to the following potential of a dipole, Julia pretty much quits on me when I try to obtain the eigenvalues and eigenvectors:  O         = [round(L/2)-hx round(L/2)-hy]# --el ORIGEN (centro -- x,y) del potencial.  Eps_o  = 8.854187817e10-12# --F*m^-1  C         = 1/(4*pi*Eps_o)  D         = 1e-21#C*m^2/s# --Debyes)  pe        = 1.8*D   P(X,Y)  = -(C)*pe*(Y/X)*(1/( (X)^2 + (Y)^2    )   )# --How the potential gets described. #--I'm aware there's singularities in the potential.#--and here's how I apply the potential to my sparse matrix.  Vi      = Float64[]# --container for the potential.  for j=Y  for i=X  push!(Vi,P(i,j))   end  end@ --applying the potential.I use this command: l, v = eigs(M,nev=15,which = :SM ,ritzvec=true)My problem seems to be that there's an error that I can't get past:ERROR: LoadError: ArgumentError: matrix has one or more zero pivots in #ldltfact!#10(::Float64, ::Function, ::Base.SparseArrays.CHOLMOD.Factor{Float64}, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1350 in (::Base.LinAlg.#kw##ldltfact!)(::Array{Any,1}, ::Base.LinAlg.#ldltfact!, ::Base.SparseArrays.CHOLMOD.Factor{Float64}, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./:0 in #ldltfact#12(::Float64, ::Array{Int64,1}, ::Function, ::Base.SparseArrays.CHOLMOD.Sparse{Float64}) at ./sparse/cholmod.jl:1386 in #ldltfact#13(::Array{Any,1}, ::Function, ::Hermitian{Float64,SparseMatrixCSC{Float64,Int64}}) at ./sparse/cholmod.jl:1426 in factorize(::SparseMatrixCSC{Float64,Int64}) at ./sparse/linalg.jl:897 in #_eigs#62(::Int64, ::Int64, ::Symbol, ::Float64, ::Int64, ::Void, ::Array{Float64,1}, ::Bool, ::Base.LinAlg.#_eigs, ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at ./linalg/arnoldi.jl:251 in (::Base.LinAlg.#kw##_eigs)(::Array{Any,1}, ::Base.LinAlg.#_eigs, ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at ./:0 in #eigs#55(::Array{Any,1}, ::Function, ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at ./linalg/arnoldi.jl:78 in (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs, ::SparseMatrixCSC{Float64,Int64}, ::UniformScaling{Int64}) at ./:0 in #eigs#54(::Array{Any,1}, ::Function, ::SparseMatrixCSC{Float64,Int64}) at ./linalg/arnoldi.jl:77 in (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs, ::SparseMatrixCSC{Float64,Int64}) at ./:0 in include_from_node1(::String) at ./loading.jl:488while loading /home/alejandro/Desktop/ACAD/PROG/ACADEMIC_PROGRAMMING/FDM (Finite_Difference_Method)/2D/SCHROD_DIP_2D/SCRATCH-2D-SI-DIPOLO2.jl, in expression starting on line 106My question is, is there a way to work around it, or, am I completely screwed?Thanks so much in advance.
Open this post in threaded view
|

## Re: Error calculating eigs of pentadiagonal matrix.

 In reply to this post by Alejandro Castellanos On Wednesday, November 2, 2016 at 6:43:48 AM UTC-4, Alejandro Castellanos wrote:I use this command: l, v = eigs(M,nev=15,which = :SM ,ritzvec=true)How are you constructing the matrix M? Apparently eigs thinks that your matrix is symmetric positive-definite (SPD), because it is trying to use a Cholesky factorization.    Schrodinger operators are indefinite unless the potential V is >= 0, so this won't work.  (If you're calling cholfact, don't.)
Open this post in threaded view
|

## Re: Error calculating eigs of pentadiagonal matrix.

 On Wednesday, November 2, 2016 at 9:48:38 AM UTC-4, Steven G. Johnson wrote:On Wednesday, November 2, 2016 at 6:43:48 AM UTC-4, Alejandro Castellanos wrote:I use this command: l, v = eigs(M,nev=15,which = :SM ,ritzvec=true)How are you constructing the matrix M? Apparently eigs thinks that your matrix is symmetric positive-definite (SPD), because it is trying to use a Cholesky factorization.    Schrodinger operators are indefinite unless the potential V is >= 0, so this won't work.  (If you're calling cholfact, don't.)Oh, nevermind, it is using an LDTL factorization, which is appropriate for symmetric-indefinite operators.Seems like it might be a bug in the eigs function — it shouldn't have a problem with singular matrices.
Open this post in threaded view
|

## Re: Error calculating eigs of pentadiagonal matrix.

Open this post in threaded view
|

## Re: Error calculating eigs of pentadiagonal matrix.

 In reply to this post by Steven G. Johnson Eigs uses shift-and-invert for the :sm variant, so it tries to solve M x = b. Try adding a small sigma parameter.
Open this post in threaded view
|

## Re: Error calculating eigs of pentadiagonal matrix.

 On Wednesday, November 2, 2016 at 1:40:17 PM UTC-4, Ralph Smith wrote:Eigs uses shift-and-invert for the :sm variant, so it tries to solve M x = b. Try adding a small sigma parameter.Yes, but it really should be able to handle matrices with a nontrivial nullspace ... the nullspace is (part of) exactly what you want, and in principle you can get it efficiently and then project it out.
Open this post in threaded view
|