Improving large array manipulation

Question

I am doing an iterative program, which uses large arrays/vectors and where the output of one round is used as an input for the next steps but works slowly when dimensions are bigger. In brief, we create firstly the complex matrix s (dimension dos x dos) using SparseArray

    dos=101;
    s = SparseArray[{{i_, i_} :> 1. + I, Band[{2, 1}] -> 2., 
    Band[{1, 2}] -> 2., Band[{dos + 1, 1}] -> 2.0, 
    Band[{1, dos + 1}] -> 2.0}, {dos*dos, dos*dos}, 0.];

And initial condition

    input=ConstantArray[1,{dos,dos}];
    input=Flatten[input];
    u1=input;

And finally, we do the iterative steps following some conditions/rules related to which position in the array is used

    lr = Table[i, {i, 1 + dos, dos*dos - dos, dos}];
    ll = Table[i, {i, dos + dos, dos*dos - dos, dos}];
    lt = Table[i, {i, 1, dos, 1}];
    lb = Table[i, {i, dos*dos - dos + 1, dos*dos, 1}];
    lll = Sort[Join[lr, ll, lt, lb]]; (*rules*)
    
    steps=50;
    Do[f = SparseArray[{{i_} :> 
       If[MemberQ[lll, i], 
        0, -2.0*u1[[i - 1]] - 2.0*u1[[i + 1]] - 2.0*u1[[i + dos]] - 
         2.0*u1[[i - dos]] + (1.+I)*u1[[i]]]}, {dos*dos}];
    u1 = LinearSolve[s, f];, {j, 1, steps, 1}]

This takes around 15 sec on my computer but with increasing dimensions, the time needed grows quickly. I noticed that applying LinearSolve is not slow, but creating the array f is where almost all the time is consumed. Is there a way to create f more efficiently? I tried to use Compile but not good results. Thanks in advance

Answer 1

First OP's implementation with timing on my machine:

First@AbsoluteTiming[
 u1 = input;
 Do[f = SparseArray[{{i_} :> 
      If[MemberQ[lll, i], 
       0, -2.0*u1[[i - 1]] - 2.0*u1[[i + 1]] - 2.0*u1[[i + dos]] - 
        2.0*u1[[i - dos]] + (1. + I)*u1[[i]]]}, {dos*dos}];
  u1 = LinearSolve[s, f];, {j, 1, steps, 1}];
 result0 = u1;
 ]

17.1247

The computation of f from u1 is linear and as such can be represented by the following sparse matrix:

A = Times[
    SparseArray[Partition[Complement[Range[dos dos], lll], 1] -> 1. + 0. I, {dos dos}, 0. I],
    SparseArray[{
      Band[{1, 1}] -> (1. + 1. I)
      , Band[{1, 1 + 1}] -> -2.0 + 0. I
      , Band[{1, 1 + dos}] -> -2.0 + 0. I
      , Band[{1 + 1, 1}] -> -2.0 + 0. I
      , Band[{1 + dos, 1}] -> -2.0 + 0. I
      },
     {dos dos, dos dos}, 0. I
     ]
    ]; // AbsoluteTiming // First

0.053559

So, using also one-time factorization of s, OP's implementation can be rewritten into

First@AbsoluteTiming[
 sinv = LinearSolve[s];
 u1 = input;
 Do[u1 = sinv[A.u1];, {j, 1, steps, 1}];
 result1 = u1;
 ]

0.141947

That is almost a 100-fold speedup.

This can also be formulated a bit shorter with Nest (which is not faster than Do):

First@AbsoluteTiming[
  sinv = LinearSolve[s];
  result2 = Nest[sinv[A.#] &, input, steps];
  ]

0.143743

Both lead to exactly the same result:

Max[Abs[result1 - result2]]

Comparison with OP's result:

Max[Abs[result1 - result0]]

541.03

Looks huge, but the relative error is much smaller:

Max[Abs[result1 - result0]]/Max[Abs[result0]]

1.64591*10^-14

What to learn from this: Pattern matching is all nice, but i maximizes rather programming speed, not execution speed.

Answer 2

The following (and my comment) overlooks the fact that u1 and f are dependent. Copy code, don't rewrite it -- also don't over edit once copied!

---

When I look at Normal@f for various values of dos, I find the "edge" elements of f (considered as a dos by dos 2d grid) are 0 and the "inner" elements of f are 1.0I-7.0. Here, then, is a quicker way to construct such arrays

f=Join[
  0&/@Range@dos,
  Join@@(Append[Prepend[1.0I-7.0&/@Range[dos-2],0],0]&/@Range[dos-2]),
  0&/@Range@dos]
]

There are definitely faster/more efficient ways to construct f, but I think this may be a case of over-early-optimization. If you're interested in situations for arbitrary rules lll, it may simply not be possible to speed things up.

A good rule of thumb is to present the computationally hardest example you can think of in the question (i.e. a randomized set of rules, although this is probably not what you're interested in). Of course, simple cases are useful for understanding the problem at hand. And then there's the necessity of splitting up your problem into small, single-question-size parts.

I think there's more that can be elucidated with further code exchange. The crossed out text detracts from this post.