3
$\begingroup$

I don't have much experience with SparseArray, and the search here did not give me an answer to this.

I'm trying to numerically compute a certain spectrum, and the method creates a tridiagonal matrix. Using SparseArray to create the matrix$^*$ makes the computation run many times faster, and the result is correct, however Mathematica gives out warnings:

Part::pkspec1: The expression -1+mr cannot be used as a part specification.
Part::pkspec1: The expression -1+mr cannot be used as a part specification.
Part::pkspec1: The expression -1+mr cannot be used as a part specification.
General::stop: Further output of Part::pkspec1 will be suppressed during this calculation.

$^*$ I say "to create the matrix", because I need a large number of Eigenvalues, so Mathematica converts the matrix to a dense one anyway.

    Nr = 1500;(*Number of grid points*)
    R = 150;(*Size of grid*)
    hr = N[R/Nr, 10];(*Grid step*)
    B = 0;(*Parameter*)
    Anr = Table[N[Sqrt[1 - 1/2/nr], 10], {nr, 1, Nr + 1}];
(*Precomputed values for the matrix start*)
    Fgs[nr_] := Exp[-2 nr hr];
    Bnr = Table[
       N[2 + 1/hr^2 - 
         1/2 (Anr[[nr]] Fgs[nr - 1] + Anr[[nr + 1]] Fgs[nr + 1])/Fgs[nr]/
           hr^2, 10], {nr, 1, Nr}];
    Bnr0 = N[2 + 1/hr^2 - 1/2 Sqrt[1/2]*Fgs[1]/hr^2, 10];
    Bnr[[Nr]] = N[2 + 1/hr^2 - 1/2 Anr[[Nr]]*Fgs[Nr - 1]/Fgs[Nr]/hr^2, 10];
    Hmm = Table[N[1/hr^2 - Bnr[[mr]] + 1/8 B^2 mr^2 hr^2], {mr, 1, Nr}];
    Hnm1 = Table[N[-1/2 Anr[[mr]]/hr^2], {mr, 1, Nr}];
    Hnm2 = Table[N[-1/2 Anr[[mr + 1]]/hr^2], {mr, 1, Nr}];
(*Precomputed values for the matrix end*)
    Hnm = SparseArray[{{1, 1} -> N[1/hr^2 - Bnr0], {2, 1} -> 
         N[-1/2 Sqrt[1/2]/hr^2], {mr_, mr_} -> 
         Hmm[[mr - 1]], {nr_, mr_} /; nr - mr == -1 -> 
         Hnm1[[mr - 1]], {nr_, mr_} /; nr - mr == 1 -> 
         Hnm2[[mr - 1]]}, {Nr + 1, Nr + 1}, 0];(*Creating the matrix*)
    T1 = AbsoluteTiming[Hnm = Normal[Hnm];
       {En, Cn} = Eigensystem[Hnm]];(*Eigensystem*)
    T1[[1]]
    g = N[1/10, 10];(*Broadening*)
    Cl = Table[Cn[[l, 1]]^2, {l, 1, Length[En]}];
    HiIm = Table[{w, g Total[Cl/(g^2 + (w - En)^2)]}, {w, -5, 15, 0.05}];
(*Computing the spectrum*)
    ListPlot[HiIm, ScalingFunctions -> "Log", Joined -> True, ImageSize -> Large]

As a result, I get a very good time (I measure only the time to solve for Eigensystem, because creating the matrix itself takes way less time):

0.535721

And a correct looking plot:

enter image description here

How do I correctly define the Sparse Array, if I need to use a list of pre-computed values? Why does Mathematica give out warnings in this case, but still gives a correct result?

This is just a model problem, I need to understand how to work here so I can solve the more complicated cases.

$\endgroup$
4
$\begingroup$

Use RuleDelayed instead of simply Rule. No error if you modify your code like this:

Hnm=SparseArray[{
    {1,1}->N[1/hr^2-Bnr0],{2,1}->N[-1/2 Sqrt[1/2]/hr^2],
    {mr_,mr_}:>Hmm[[mr-1]],
    {nr_,mr_}/;nr-mr==-1:>Hnm1[[mr-1]],
    {nr_,mr_}/;nr-mr==1:>Hnm2[[mr-1]]},
    {Nr+1,Nr+1},0];
$\endgroup$
  • $\begingroup$ Thank you for this! $\endgroup$ – Yuriy S Sep 17 '18 at 21:07
3
$\begingroup$

Vitaliy has already given complete answer to your question. I would just like to add as an extended comment that using patterns in SparseArray creation comes at a certain cost, in particular if the values to write are already precomputed. Whenever you want to fill diagonals of a matrix, it is a good idea to use Band instead.

Here is Vitaliy's solution:

A = SparseArray[{
       {1, 1} -> N[1/hr^2 - Bnr0], 
       {2, 1} -> N[-1/2 Sqrt[1/2]/hr^2], 
       {mr_, mr_} :> Hmm[[mr - 1]], 
       {nr_, mr_} /; nr - mr == -1 :> Hnm1[[mr - 1]], 
       {nr_, mr_} /; nr - mr == 1 :> Hnm2[[mr - 1]]}, 
       {Nr + 1, Nr + 1}, 0.]; // RepeatedTiming // First

0.012

And here is the same using Band:

B = SparseArray[
     {
      {1, 1} -> N[1/hr^2 - Bnr0], 
      {2, 1} -> N[-1/2 Sqrt[1/2]/hr^2],
      Band[{2, 2}] -> Hmm,
      Band[{1, 2}] -> Hnm1,
      Band[{3, 2}] -> Most@Hnm2
      },
     {Nr + 1, Nr + 1}, 0.]; // RepeatedTiming // First

0.0027

Checking the correctness:

A == B

True

$\endgroup$
  • $\begingroup$ Thank you very much! I have read about Band, but all the examples listed are constants, and I didn't know how to use a list $\endgroup$ – Yuriy S Sep 18 '18 at 7:59
  • $\begingroup$ You're welcome. $\endgroup$ – Henrik Schumacher Sep 18 '18 at 8:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.