Imperative Sieve of Eratosthenes using destructive updates in Haskell

I've implemented the prime sieve of Eratosthenes in Haskell using a mutable array of unboxed Bools instead of the usual pure and immutable datastructures.

My main concern with this code is speed. For large upper bounds, it runs approximately 70% slower than C code implementing the same algorithm.
Memory is already optimal (1 bit per number plus some small overhead).

I would also like to know if the code is reasonably idiomatic, considering the fact that is impure, and basically non-functional (meaning imperative).
Any nitpicks are appreciated.

import Prelude
import System.Environment (getArgs)
import Data.Array.Unboxed (UArray, (!))
import Data.Array.Storable (newArray, readArray, writeArray, getBounds)
import Data.Array.ST (STUArray, runSTUArray)
import Control.Monad (when)
import Control.Monad.ST (ST)

markMultiplesAsNotPrime :: STUArray s Int Bool -> Int -> ST s ()
markMultiplesAsNotPrime isPrime k = do
  (_, n)  writeArray isPrime i False) multiples

sieve :: Int -> UArray Int Bool
sieve n | n  Int
isqrt' = ceiling . (sqrt :: Double -> Double) . fromIntegral

primesSmallerThanOrEqualTo :: Int -> [Int]
primesSmallerThanOrEqualTo n | n < 2 = []
primesSmallerThanOrEqualTo n = let 
    isPrime = sieve n
  in
    [i | i <- [2 .. n], isPrime ! i]

main :: IO ()
main = do
  args <- getArgs
  mapM_ (print . last . primesSmallerThanOrEqualTo . read) args

In case it is relevant:
I compiled and executed the code on an amd64 CPU using "The Glorious Glasgow Haskell Compilation System, version 7.10.3" with the optimization switch "-O2".

Some timings (as measured with the runtime system switch "-s", best of five runs): $$ \begin{array}{4c} n & π(n) & t [s] & \log\left(\frac{t_2}{t_1}\right) \bigg/ \log\left(\frac{π(n_2)}{π(n_1)}\right) \\ 10^4 & 1.23 \cdot 10^3 & 2 \cdot 10^{-3} & - \\ 10^5 & 9.59 \cdot 10^3 & 4 \cdot 10^{-3} & 0.3 \\ 10^6 & 7.85 \cdot 10^4 & 1.0 \cdot 10^{-2} & 0.4 \\ 10^7 & 6.65 \cdot 10^5 & 9.3 \cdot 10^{-2} & 1.04 \\ 10^8 & 5.76 \cdot 10^6 & 1

takeWhile (<=n) $ iterate (+k) (k+k) can be written as [2k, 3k..n] which, here, can be replaced with [kk, kk+k..n].

Similarly, takeWhile (<= kMax) . iterate (+1) $ 2 can be written as [2..kMax].

Algorithmically, we can work with odds only, and switch to [kk, kk+2*k..n] and [3,5..kMax] (special-casing the 2). Taking it further to pre-ignoring the multiples of 3 (5,7...) also, this is known as wheel factorization optimization.

mapM and mapM_, though very idiomatic, can be faster if coded manually as loops (recursive functions). This would also enumerate the numbers manually, instead of creating the interim lists and hoping the compiler won't actually create them.

Using unsafeRead and unsafeWrite should be faster than the safe variants.

Lastly, last . primesSmallerThanOrEqualTo can be replaced with lastPrime,

lastPrime n = let arr = sieve n in
              take 1 [ i | i <- [n, n-1..2], arr ! i]

Even faster is to move the last prime-finding code into the sieve function itself, and return that prime from it, instead of the array.

StackExchange Code Review Q#144124, answer score: 3