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]]} :> 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[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?
Thanks
Condition
:A = Normal[SparseArray[{{i_, j_} /; j == N1[i][[1]]} :> 1, {nn, nn}]]
– Michael E2 Dec 16 '13 at 15:17