0
$\begingroup$

The following is a question in regards to the matrix $A$ (originally the zero $4n \times 4n$ matrix) below:

 In[288]:= For[j = 1, j <= 2 n, j++,
   For[k = 1, k <= 2 n, k++,
     A[[2 j - 1]][[2 k - 1]] = -2 I*H[[j]][[k]] - M[[k]][[j]] + 
       M[[j]][[k]];
     A[[2 j - 1]][[2 k]] = 2*I*M[[j]][[k]];
     A[[2 j]][[2 k - 1]] = -2*I*M[[k]][[j]];
     A[[2 j]][[2 k]] = -2 I*H[[j]][[k]] + M[[k]][[j]] - M[[j]][[k]];
     ];
   ]; // AbsoluteTiming

Out[288]= {25.767, Null}

$H$ and $M$ are both matrices of size $2n \times 2n$ and $n = 640$.

I have tried using Table (code below) but it seems to run much slower than the above code.

In[291]:= A = Table[
    p = Ceiling[j/2]; q = Ceiling[k/2];

    If[OddQ[j] && OddQ[k], -2 I*H[[p]][[q]] - M[[q]][[p]] + 
       M[[p]][[q]], 0]
     If[OddQ[j] && EvenQ[k], 2*I*M[[p]][[q]], 0]
     If[EvenQ[j] && OddQ[k], -2*I*M[[q]][[p]], 0]
     If[EvenQ[j] && EvenQ[k], -2 I*H[[p]][[q]] + M[[q]][[p]] - 
       M[[p]][[q]], 0]

    , {j, 1, 4 n}, {k, 1, 4 n}]; // AbsoluteTiming

Out[291]= {52.4175, Null}

I was wondering what possible speed ups there are for the generation of $A$?

$\endgroup$

1 Answer 1

3
$\begingroup$
n = 640;
A1 = A2 = A3 = ConstantArray[0. + 0. I, {4 n, 4 n}];
{H, M} = RandomComplex[{-1 - I, 1 + I}, {2, 2 n, 2 n}];

The original code and its timing from my machine:

For[j = 1, j <= 2 n, j++,
    For[k = 1, k <= 2 n, k++,
      A1[[2 j - 1]][[2 k - 1]] = -2 I*H[[j]][[k]] - M[[k]][[j]] + 
        M[[j]][[k]];
      A1[[2 j - 1]][[2 k]] = 2*I*M[[j]][[k]];
      A1[[2 j]][[2 k - 1]] = -2*I*M[[k]][[j]];
      A1[[2 j]][[2 k]] = -2 I*H[[j]][[k]] + M[[k]][[j]] - 
        M[[j]][[k]];
      ];
    ]; // AbsoluteTiming // First

39.3159

The same with a Do loop and with simplified Part-syntax (main speedup is indeed achieved by the change of the Part-syntax!):

Do[
    A2[[2 j - 1, 2 k - 1]] = -2 I H[[j, k]] - M[[k, j]] + M[[j, k]];
    A2[[2 j - 1, 2 k]] = 2 I M[[j, k]];
    A2[[2 j, 2 k - 1]] = -2 I M[[k, j]];
    A2[[2 j, 2 k]] = -2 I H[[j, k]] + M[[k, j]] - M[[j, k]];
    , {j, 1, 2 n}, {k, 1, 2 n}]; // AbsoluteTiming // First

18.3622

And a vectorized approach:

First@RepeatedTiming[
  A3[[1 ;; ;; 2, 1 ;; ;; 2]] = -2 I H - Transpose[M] + M;
  A3[[1 ;; ;; 2, 2 ;; ;; 2]] = 2 I M;
  A3[[2 ;; ;; 2, 1 ;; ;; 2]] = -2 I Transpose[M];
  A3[[2 ;; ;; 2, 2 ;; ;; 2]] = -2 I H + Transpose[M] - M;
  ]
A1 == A2 == A3

0.150

For the meaning of ;;, lookup Span in the documentation.

Also notice that it is important to initialize A3 as matrix of complex zeroes 0. + 0. I in order to keep A3 a packed array.

$\endgroup$
1
  • 1
    $\begingroup$ +1 for the vectorized part, this is how it should be done IMHO. $\endgroup$ Aug 5, 2018 at 22:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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