[julia-users] Re: Error calculating eigs of pentadiagonal matrix.
As of now, I'm leaning towards maybe my feeding the matrix some Inf values
due to the grid spacing I'm using but I'm not sure.
Would it help if I uploaded the whole of my code, in here? Mind you, the
comments are in Spanish...
Cheers.
El miércoles, 2 de noviembre de 2016, 3:43:48 (UTC-7), Alejandro
Castellanos escribió:
>
> 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:488
>>>
>>> while 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 106
>>>
>>>
>>> My question is, is there a way to work around it, or, am I completely
> screwed?
>
> Thanks so much in advance.
>
[julia-users] 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.
[julia-users] Re: Error calculating eigs of pentadiagonal matrix.
Eigs uses shift-and-invert for the :sm variant, so it tries to solve M x = b. Try adding a small sigma parameter.
[julia-users] Re: Error calculating eigs of pentadiagonal matrix.
Hi, Steven. This is how I'm building the matrix:
#==#
# OPERATING THE MATRIX #
#==#
# --Creating HAMILTONIAN container.
M = spzeros(n*n, n*n)#--Creating the SPARSE matrix (n is
the number of points per row).
# --Working kinetic energy into matrix M (using 2D Laplacian scheme)...
# --Central part:
trmC = -4#--el termino central de la matriz.
M[1,1] = trmC
for i = 2:n*n M[i,i] = trmC; M[i-1,i] = 1 end#--diags. central and 1
above.
for i = 1:(n*n)-1 M[i+1,i] = 1 end#--diag. below.
#--Exterior diagonals (2 in 1 loop).
for i = 1:n*(n-1) M[i,i+n] = 1; M[i+n,i] = 1 end
#--"poking" zeroes into the matrix to make "holes."
for i = 1:n*(n-1) if i%n==0 M[i+1,i]=0;M[i,i+1] = 0 endend
After that I multiply M by the appropriate constants to get the full
kinetic energy. Then I use the routine I posted previously so as to add the
potential, all in the same matrix.
I assume my matrix creation scheme may not be the most efficient, but it
did get the job done (up til now).
Cheers.
El miércoles, 2 de noviembre de 2016, 3:43:48 (UTC-7), Alejandro
Castellanos escribió:
>
> 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:488
>>>
>>> while 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 106
>>>
>>>
>>> My question is, is there a way to work around it, or, am I completely
> screwed?
>
> Thanks so much in advance.
>
[julia-users] 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.
[julia-users] Re: Error calculating eigs of pentadiagonal matrix.
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.)
[julia-users] 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:488
>>>
>>> while 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 106
>>>
>>>
>>> My question is, is there a way to work around it, or, am I completely
> screwed?
>
> Thanks so much in advance.
>
