Comment by mlochbaum

I prefer "middle of the road" languages that are high-level AND readable AND have decent performance optimisation for bulk operations. Python with C libraries suffices for a lot of people, Julia similarly is getting popular.

Even the older Mathematica language blows both K and C out of the water for readability and performance:

    m = Table[RandomReal[], 100, 100];
    t = RepeatedTiming[MatrixPower[m, 2]];
    First[t]*1000*1000

    18.0245 

You would have to know literally zero about Mathematica's Wolfram Language to be able to read that clearly! From a standing start you could understand what another person has created. For K, you'd have to have memorized the K-specific syntax. For C, if you hadn't seen the standard patterns for matrix multiplication you'd have to read the code carefully. A lot of it is just noise, like the verbose for loop syntax.

Oh, and Mathematica's matrix power function is:

- Parallel! If I use a 10K x 10K matrix as an input, it uses about 75% CPU on my 8-core laptop. It can complete a single multiplication in 5.3 seconds. For laughs, try that with either K or C and see what you get...

- Extends to negative or fraction powers.

- Has optimisations for applying the matrix power directly to a vector.

- Is extensively documented, unlike the terse K snippet or the hand-rolled C code: https://reference.wolfram.com/language/ref/MatrixPower.html....

Essentially what I'm trying to say is that terseness is an anti-pattern, and doesn't even begin to approach the utility of a well designed high-level language intended for teams of collaborating humans.

When you compile the C example with the compiler option "-Wall", you'll get a warning that the variable 'r' is set but not used, so the compiler would be free to simply skip the for loops. In fact, if you compile with clang instead of gcc, the compiler will do just that and you'll get almost zero computation time.

It would be better to do something with the computed result so the compiler does not remove the computation, e.g. print a randomly selected value.

I also benchmarked the fixed C code against Python and "Python" (which uses OpenBLAS under the hood in my case) was 10 times faster:

    import time
    import numpy as np

    a = np.random.rand(100, 100)
    t = time.perf_counter()

    for _ in range(1000): a @ a

    print(time.perf_counter() - t, "seconds")

Implementation matters a lot for matrix multiplication.

Naive implementations of stock matrix math can't get anywhere close to numpy or julia, which both use BLAS and automatically parallelize across cores.

  % python matrix.py
  Timing 10 squares of a random 10000 x 10000 matrix
  97.3976636590669 seconds
  python matrix.py  364.41s user 8.10s system 379% cpu 1:38.25 total

julia has more overhead, and the first multiply triggers code compilation so there's an additional warm-up square outside of the timing loop, but its "warm" performance is equivalent to numpy. Turning on extra optimizations (-O3) can even make it a couple seconds faster than numpy once warmed up.

  % julia matrix.jl
  Timing 10 squares of a random 10000 x 10000 matrix
   97.787679 seconds (31 allocations: 7.451 GiB, 0.33% gc time)
  julia matrix.jl  405.34s user 8.13s system 375% cpu 1:50.09 total

If you're going to wait for that C implementation, or the other comment's K implementation, to finish that loop, you'll want a book.

Nice benchmarking, thank you for your effort.

Also the scientist was leading a team, so my program would've been used by at least 20 people, 10 times per day. That was why Python one was a no go for them from beginning.