# Specify positions in a SparseArray using patterns with functions

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:

Nx=2;
Ny=1;
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}];


Array:

delta =
Table[
Total[Flatten[V[[μ, ν, ;; , ;; , ;; , ;;]]*
Table[
Table[
Which[
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}, {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:

 A={{0,1},{1,0}}


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

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

Thanks

• Try Condition: A = Normal[SparseArray[{{i_, j_} /; j == N1[i][]} :> 1, {nn, nn}]] 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},
Compile[{},
Table[
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[
SparseArray[
Table[Band[{i, Mod[i + 1, Nx, 1]}, Automatic, {Nx, Nx}] -> 2, {i, Nx}],
{nn, nn}],
{Nband, Nband, Nband, Nband}]; // AbsoluteTiming

(* {0.006551, Null} *)


Check:

 Vc == spV
(* True *)


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

• 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!!! Dec 16 '13 at 15:39