Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I spent some time trying to figure out what was what was making so slow the calculation of the following array(actually this is only a SUPERsimple toy model):

Input to the array:

nn =Nx*Ny; 
Nband = 2;
Nstates = 2*nn*Nband;

eigenvectors = Table[Range[Nstates] + i, {i, 1, Nstates}];

InverseFlatten[l_, dimensions_] := Fold[Partition, Flatten@l,Most[Reverse[dimensions]]];
N1[i_] := Module[{ix, iy, n1x},
ix = Mod[i, Nx, 1];
iy = Quotient[i, Nx, 1] + 1;
If[ix + 1 > Nx, n1x = 1, n1x = ix + 1];
{n1x + (iy - 1)*Nx}];
uv = InverseFlatten[eigenvectors, {Nstates, Nband, 2, nn}];
u = uv[[1 ;; Nstates, 1 ;; Nband, 1]];
v = uv[[1 ;; Nstates, 1 ;; Nband, 2]];
f = Range[Nstates];
V = Table[Which[MemberQ[N1[i], j] == True, 2], 
{l, 1, Nband}, {m, 1, Nband}, {s, 1, Nband},{q, 1, Nband}, {i, 1, nn}, {j, 1, nn}];


delta = 
Total[Flatten[V[[μ, ν, ;; , ;; , ;; , ;;]]*
            MemberQ[N1[i], j] == True,
            f.(u[[;; , q, i]]*Conjugate[v[[;; , s, j]]])    ], 
          {i, 1, nn}, {j, 1, nn}],
      {q, 1, Nband}, {s, 1, Nband}], 1]], 
 {μ, 1, Nband}, {ν, 1, Nband}];

After having tried several things like using Packed arrays, using Compile...and seen that nothing really helped, I thing that what it actually slows down the thing is Which. I think, that the process should be a lot faster if instead of using Which, which gives a condition for which matrix elements to calculate, I directly specify the positions for which I would like to calculate the matrix element; I think so because the number of elements to be calculated is really small compared to the number of elements(for big nn).

Then looking at the SparseArray list manipulating tutorial I saw that one can do this using patterns, for example:

Normal[SparseArray[{{i_, i_} :> 1}, {nn, nn}]]

Now, what specifies the second index of my matrix is the function N1 defined above, therefore, I would like to use something like:

 A=Normal[SparseArray[{{i_, N1[i_][[1]]} :> 1}, {nn, nn}]]  
 SparseArray::posd: The left-hand side of {i_,n1x$671+2 Quotient[i_,2,1]}:>1  in {{i_,n1x$671+2 Quotient[i_,2,1]}:>1} is not a position or a pattern that will match the position of an element in an array with dimensions {2,2}. >>

but it complaints. I expect to get:


since N1[1][[1]]=2 and N1[2][[1]]=1.

Does anybody now how could I make a "pattern SparseArray" having this function N1 in the pattern?


share|improve this question
Try Condition: A = Normal[SparseArray[{{i_, j_} /; j == N1[i][[1]]} :> 1, {nn, nn}]] – Michael E2 Dec 16 '13 at 15:17

If you want a fast way, then you don't really want to pursue SparseArrays with patterns and condition. Built-in patterns like Band are optimized and work relatively fast.

Slow ways

With the parameters,

Nx = 10;
Ny = 20;
nn = Nx*Ny;
Nband = 5;
(* etc. *)

the following took nearly 90 sec., while I went to get some coffee.

SparseArray[{l_, m_, s_, q_, i_, j_} /; 
      Mod[i, Nx] + Quotient[i, Nx, 1]*Nx + 1 == j -> 2,
   {Nband, Nband, Nband, Nband, nn, nn}]; // AbsoluteTiming

(* {87.753965, Null} *)

The OP's Table method took too long for me to wait. A compiled version took almost 6 seconds:

Vc = With[{Nband = Nband, nn = nn, Nx = Nx, Ny = Ny},
       If[Mod[i, Nx] + Quotient[i, Nx, 1]*Nx + 1 == j, 2, 0],
       {l, 1, Nband}, {m, 1, Nband}, {s, 1, Nband}, {q, 1, Nband}, {i, 1, nn}, {j, 1, nn}]]
      ][]; // AbsoluteTiming

(* {5.772137, Null} *)

Much faster

This SparseArray approach was several hundred times faster:

spV = SparseArray@ConstantArray[
      Table[Band[{i, Mod[i + 1, Nx, 1]}, Automatic, {Nx, Nx}] -> 2, {i, Nx}],
      {nn, nn}],
     {Nband, Nband, Nband, Nband}]; // AbsoluteTiming

(* {0.006551, Null} *)


 Vc == spV
 (* True *)

P.S. I replaced Null with 0. That seemed to be acceptable in one of the OP's earlier posts.

share|improve this answer
thanks so much for answering(always there lately :)), I will see if I can understand what you did and apply it to calculate delta in a fast way, which is my problem because, although the computation of V is also slow, V is not in an iterating convergence process, whereas delta is, so I really need delta to be fast!!! – Mencia Dec 16 '13 at 15:39

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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