Quantcast

Error calculating eigs of pentadiagonal matrix.

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Error calculating eigs of pentadiagonal matrix.

Alejandro Castellanos
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 ./<missing>: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 ./<missing>: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 ./<missing>: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 ./<missing>: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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Error calculating eigs of pentadiagonal matrix.

Jeffrey Sarnoff
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 ./<missing>: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 ./<missing>: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 ./<missing>: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 ./<missing>: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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Error calculating eigs of pentadiagonal matrix.

Steven G. Johnson
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.)
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Error calculating eigs of pentadiagonal matrix.

Steven G. Johnson


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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Error calculating eigs of pentadiagonal matrix.

Alejandro Castellanos
In reply to this post by Alejandro Castellanos
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  end    end

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 ./<missing>: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 ./<missing>: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 ./<missing>: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 ./<missing>: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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Error calculating eigs of pentadiagonal matrix.

Ralph Smith
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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Error calculating eigs of pentadiagonal matrix.

Steven G. Johnson


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.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Error calculating eigs of pentadiagonal matrix.

Alejandro Castellanos
In reply to this post by Alejandro Castellanos
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 ./<missing>: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 ./<missing>: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 ./<missing>: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 ./<missing>: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.
Loading...