diff --git a/docs/man/futhark-run.rst b/docs/man/futhark-run.rst
--- a/docs/man/futhark-run.rst
+++ b/docs/man/futhark-run.rst
@@ -9,14 +9,14 @@
 SYNOPSIS
 ========
 
-futhark run <program.fut>
+futhark run [options...] <program.fut>
 
 DESCRIPTION
 ===========
 
-Execute the given program by evaluating the ``main`` function with
-arguments read from standard input, and write the results on standard
-output.
+Execute the given program by evaluating an entry point (``main`` by
+default) with arguments read from standard input, and write the
+results on standard output.
 
 ``futhark run`` is very slow, and in practice only useful for testing,
 teaching, and experimenting with the language.  The ``#[trace]`` and
diff --git a/futhark.cabal b/futhark.cabal
--- a/futhark.cabal
+++ b/futhark.cabal
@@ -1,6 +1,6 @@
 cabal-version: 2.4
 name:           futhark
-version:        0.21.2
+version:        0.21.3
 synopsis:       An optimising compiler for a functional, array-oriented language.
 
 description:    Futhark is a small programming language designed to be compiled to
@@ -308,7 +308,7 @@
       Paths_futhark
   hs-source-dirs:
       src
-  ghc-options: -Wall -Wcompat -Wredundant-constraints -Wincomplete-record-updates -Wmissing-export-lists
+  ghc-options: -Wall -Wcompat -Wno-incomplete-uni-patterns -Wredundant-constraints -Wincomplete-record-updates -Wmissing-export-lists
   build-tool-depends:
       alex:alex
     , happy:happy
@@ -366,7 +366,7 @@
   main-is: src/futhark.hs
   other-modules:
       Paths_futhark
-  ghc-options: -Wall -Wcompat -Wredundant-constraints -Wincomplete-record-updates -Wmissing-export-lists -threaded -rtsopts "-with-rtsopts=-N -qg1 -A16M"
+  ghc-options: -Wall -Wcompat -Wno-incomplete-uni-patterns -Wredundant-constraints -Wincomplete-record-updates -Wmissing-export-lists -threaded -rtsopts "-with-rtsopts=-N -qg1 -A16M"
   build-depends:
       base
     , futhark
diff --git a/src/Futhark/Analysis/UsageTable.hs b/src/Futhark/Analysis/UsageTable.hs
--- a/src/Futhark/Analysis/UsageTable.hs
+++ b/src/Futhark/Analysis/UsageTable.hs
@@ -18,8 +18,10 @@
     inResultUsage,
     sizeUsage,
     sizeUsages,
+    withoutU,
     Usages,
     consumedU,
+    presentU,
     usageInStm,
   )
 where
@@ -138,7 +140,7 @@
 matches :: Usages -> Usages -> Bool
 matches (Usages x) (Usages y) = x == (x .&. y)
 
--- | x - y, but for Usages.
+-- | x - y, but for 'Usages'.
 withoutU :: Usages -> Usages -> Usages
 withoutU (Usages x) (Usages y) = Usages $ x .&. complement y
 
diff --git a/src/Futhark/CLI/Autotune.hs b/src/Futhark/CLI/Autotune.hs
--- a/src/Futhark/CLI/Autotune.hs
+++ b/src/Futhark/CLI/Autotune.hs
@@ -8,6 +8,7 @@
 import qualified Data.ByteString.Char8 as SBS
 import Data.Function (on)
 import Data.List (elemIndex, intersect, isPrefixOf, minimumBy, sort, sortOn)
+import qualified Data.Map as M
 import Data.Maybe
 import qualified Data.Set as S
 import qualified Data.Text as T
@@ -240,31 +241,36 @@
             return ((v, cmp), ancestors')
        in ((parent, parent_cmp), mapMaybe isChild thresholds)
 
+-- | The performance difference in percentage that triggers a non-monotonicity
+-- warning. This is to account for slight variantions in run-time.
+epsilon :: Double
+epsilon = 1.02
+
 --- Doing the atual tuning
 
 tuneThreshold ::
   AutotuneOptions ->
   Server ->
   [(DatasetName, RunDataset, T.Text)] ->
-  Path ->
+  (Path, M.Map DatasetName Int) ->
   (String, Path) ->
-  IO Path
-tuneThreshold opts server datasets already_tuned (v, _v_path) = do
-  tune_result <-
-    foldM tuneDataset Nothing datasets
+  IO (Path, M.Map DatasetName Int)
+tuneThreshold opts server datasets (already_tuned, best_runtimes0) (v, _v_path) = do
+  (tune_result, best_runtimes) <-
+    foldM tuneDataset (Nothing, best_runtimes0) datasets
   case tune_result of
     Nothing ->
-      return $ (v, thresholdMin) : already_tuned
+      return ((v, thresholdMin) : already_tuned, best_runtimes)
     Just (_, threshold) ->
-      return $ (v, threshold) : already_tuned
+      return ((v, threshold) : already_tuned, best_runtimes)
   where
-    tuneDataset :: Maybe (Int, Int) -> (DatasetName, RunDataset, T.Text) -> IO (Maybe (Int, Int))
-    tuneDataset thresholds (dataset_name, run, entry_point) =
+    tuneDataset :: (Maybe (Int, Int), M.Map DatasetName Int) -> (DatasetName, RunDataset, T.Text) -> IO (Maybe (Int, Int), M.Map DatasetName Int)
+    tuneDataset (thresholds, best_runtimes) (dataset_name, run, entry_point) =
       if not $ isPrefixOf (T.unpack entry_point ++ ".") v
         then do
           when (optVerbose opts > 0) $
             putStrLn $ unwords [v, "is irrelevant for", T.unpack entry_point]
-          return thresholds
+          return (thresholds, best_runtimes)
         else do
           putStrLn $
             unwords
@@ -290,7 +296,7 @@
               when (optVerbose opts > 0) $
                 putStrLn $
                   "Sampling run failed:\n" ++ err
-              return thresholds
+              return (thresholds, best_runtimes)
             Right (cmps, t) -> do
               let (tMin, tMax) = fromMaybe (thresholdMin, thresholdMax) thresholds
               let ePars =
@@ -311,10 +317,30 @@
               when (optVerbose opts > 1) $
                 putStrLn $ unwords ("Got ePars: " : map show ePars)
 
-              newMax <- binarySearch runner (t, tMax) ePars
-              let newMinIdx = pred <$> elemIndex newMax ePars
-              let newMin = maxinum $ catMaybes [Just tMin, newMinIdx]
-              return $ Just (newMin, newMax)
+              (best_t, newMax) <- binarySearch runner (t, tMax) ePars
+              let newMinIdx = do
+                    i <- pred <$> elemIndex newMax ePars
+                    if i < 0 then fail "Invalid lower index" else return i
+              let newMin = maxinum $ catMaybes [Just tMin, fmap (ePars !!) newMinIdx]
+              best_runtimes' <-
+                case dataset_name `M.lookup` best_runtimes of
+                  Just rt
+                    | fromIntegral rt * epsilon < fromIntegral best_t -> do
+                      putStrLn $
+                        unwords
+                          [ "WARNING! Possible non-monotonicity detected. Previous best run-time for dataset",
+                            dataset_name,
+                            " was",
+                            show rt,
+                            "but after tuning threshold",
+                            v,
+                            "it is",
+                            show best_t
+                          ]
+                      return best_runtimes
+                  _ ->
+                    return $ M.insertWith min dataset_name best_t best_runtimes
+              return (Just (newMin, newMax), best_runtimes')
 
     bestPair :: [(Int, Int)] -> (Int, Int)
     bestPair = minimumBy (compare `on` fst)
@@ -327,7 +353,7 @@
     candidateEPar (tMin, tMax) (threshold, ePar) =
       ePar > tMin && ePar < tMax && threshold == v
 
-    binarySearch :: (Int -> Int -> IO (Maybe Int)) -> (Int, Int) -> [Int] -> IO Int
+    binarySearch :: (Int -> Int -> IO (Maybe Int)) -> (Int, Int) -> [Int] -> IO (Int, Int)
     binarySearch runner best@(best_t, best_e_par) xs =
       case splitAt (length xs `div` 2) xs of
         (lower, middle : middle' : upper) -> do
@@ -363,14 +389,14 @@
                       "and",
                       show middle'
                     ]
-              return best_e_par
+              return (best_t, best_e_par)
         (_, _) -> do
           when (optVerbose opts > 0) $
             putStrLn $ unwords ["Trying e_pars", show xs]
           candidates <-
             catMaybes . zipWith (fmap . flip (,)) xs
               <$> mapM (runner $ timeout best_t) xs
-          return $ snd $ bestPair $ best : candidates
+          return $ bestPair $ best : candidates
 
 --- CLI
 
@@ -389,7 +415,8 @@
   putStrLn $ "Running with options: " ++ unwords (serverOptions opts)
 
   withServer (futharkServerCfg progbin (serverOptions opts)) $ \server ->
-    foldM (tuneThreshold opts server datasets) [] $ tuningPaths forest
+    fmap fst . foldM (tuneThreshold opts server datasets) ([], mempty) $
+      tuningPaths forest
 
 runAutotuner :: AutotuneOptions -> FilePath -> IO ()
 runAutotuner opts prog = do
diff --git a/src/Futhark/CLI/Pkg.hs b/src/Futhark/CLI/Pkg.hs
--- a/src/Futhark/CLI/Pkg.hs
+++ b/src/Futhark/CLI/Pkg.hs
@@ -185,7 +185,6 @@
 
 instance Monad PkgM where
   PkgM m >>= f = PkgM $ m >>= unPkgM . f
-  return = PkgM . return
 
 instance MonadFail PkgM where
   fail s = liftIO $ do
diff --git a/src/Futhark/CodeGen/ImpGen/GPU/Base.hs b/src/Futhark/CodeGen/ImpGen/GPU/Base.hs
--- a/src/Futhark/CodeGen/ImpGen/GPU/Base.hs
+++ b/src/Futhark/CodeGen/ImpGen/GPU/Base.hs
@@ -1189,7 +1189,7 @@
   barrier
 
   sComment "restore correct values for first block" $
-    sWhen is_first_block $
+    sWhen (is_first_block .&&. ltid_in_bounds) $
       forM_ (zip3 x_params y_params arrs) $ \(x, y, arr) ->
         if primType (paramType y)
           then copyDWIM arr [DimFix ltid] (Var $ paramName y) []
diff --git a/src/Futhark/CodeGen/ImpGen/GPU/SegRed.hs b/src/Futhark/CodeGen/ImpGen/GPU/SegRed.hs
--- a/src/Futhark/CodeGen/ImpGen/GPU/SegRed.hs
+++ b/src/Futhark/CodeGen/ImpGen/GPU/SegRed.hs
@@ -162,7 +162,7 @@
     forM (lambdaReturnType lam) $ \t -> do
       let pt = elemType t
           extra_dim
-            | primType t = intConst Int64 1
+            | primType t, shapeRank shape == 0 = intConst Int64 1
             | otherwise = group_size
           full_shape = Shape [extra_dim, virt_num_groups] <> shape <> arrayShape t
           -- Move the groupsize dimension last to ensure coalesced
@@ -446,7 +446,7 @@
             `rem` (sExt64 (unCount group_size') * groups_per_segment)
 
       let first_group_for_segment = sExt64 flat_segment_id * groups_per_segment
-      dIndexSpace (zip segment_gtids (init dims')) $sExt64 flat_segment_id
+      dIndexSpace (zip segment_gtids (init dims')) $ sExt64 flat_segment_id
       dPrim_ (last gtids) int64
       let num_elements = Imp.elements $ toInt64Exp w
 
@@ -695,11 +695,8 @@
     reductionStageZero constants ispace num_elements global_tid elems_per_thread threads_per_segment slugs body
 
   case slugsComm slugs of
-    Noncommutative ->
-      forM_ slugs $ \slug ->
-        forM_ (zip (accParams slug) (slugAccs slug)) $ \(p, (acc, acc_is)) ->
-          copyDWIMFix (paramName p) [] (Var acc) acc_is
-    _ -> doTheReduction
+    Noncommutative -> pure ()
+    Commutative -> doTheReduction
 
   return slugs_op_renamed
 
diff --git a/src/Futhark/IR/Mem/Simplify.hs b/src/Futhark/IR/Mem/Simplify.hs
--- a/src/Futhark/IR/Mem/Simplify.hs
+++ b/src/Futhark/IR/Mem/Simplify.hs
@@ -170,6 +170,7 @@
         Just tse <- maybeNth j $ bodyResult tbranch,
         Just fse <- maybeNth j $ bodyResult fbranch,
         mem `onlyUsedIn` patElemName pat_elem,
+        length (IxFun.base ixfun) == shapeRank shape, -- See #1325
         all knownSize (shapeDims shape),
         not $ freeIn ixfun `namesIntersect` namesFromList (patNames pat),
         fse /= tse =
diff --git a/src/Futhark/IR/Primitive.hs b/src/Futhark/IR/Primitive.hs
--- a/src/Futhark/IR/Primitive.hs
+++ b/src/Futhark/IR/Primitive.hs
@@ -113,6 +113,16 @@
 import qualified Data.Binary.Get as G
 import qualified Data.Binary.Put as P
 import Data.Bits
+  ( complement,
+    countLeadingZeros,
+    countTrailingZeros,
+    popCount,
+    shift,
+    shiftR,
+    xor,
+    (.&.),
+    (.|.),
+  )
 import Data.Fixed (mod') -- Weird location.
 import Data.Int (Int16, Int32, Int64, Int8)
 import qualified Data.Map as M
diff --git a/src/Futhark/Optimise/Simplify/Engine.hs b/src/Futhark/Optimise/Simplify/Engine.hs
--- a/src/Futhark/Optimise/Simplify/Engine.hs
+++ b/src/Futhark/Optimise/Simplify/Engine.hs
@@ -535,7 +535,8 @@
         zip (patNames pat) (patAliases pat)
     usageThroughBindeeAliases (name, aliases) = do
       uses <- UT.lookup name utable
-      return $ mconcat $ map (`UT.usage` uses) $ namesToList aliases
+      pure . mconcat $
+        map (`UT.usage` (uses `UT.withoutU` UT.presentU)) $ namesToList aliases
 
 type BlockPred rep = ST.SymbolTable rep -> UT.UsageTable -> Stm rep -> Bool
 
@@ -716,8 +717,7 @@
 usageFromDiet Consume = UT.consumedU
 usageFromDiet _ = mempty
 
--- | Simplify a single 'Result'.  The @[Diet]@ only covers the value
--- elements, because the context cannot be consumed.
+-- | Simplify a single 'Result'.
 simplifyResult ::
   SimplifiableRep rep => [UT.Usages] -> Result -> SimpleM rep (Result, UT.UsageTable)
 simplifyResult usages res = do
@@ -725,7 +725,11 @@
   vtable <- askVtable
   let more_usages = mconcat $ do
         (u, Var v) <- zip usages $ map resSubExp res
-        map (`UT.usage` u) $ v : namesToList (ST.lookupAliases v vtable)
+        let als_usages =
+              map
+                (`UT.usage` (u `UT.withoutU` UT.presentU))
+                (namesToList (ST.lookupAliases v vtable))
+        UT.usage v u : als_usages
   return (res', UT.usages (freeIn res') <> more_usages)
 
 isDoLoopResult :: Result -> UT.UsageTable
diff --git a/src/Futhark/Pass/ExpandAllocations.hs b/src/Futhark/Pass/ExpandAllocations.hs
--- a/src/Futhark/Pass/ExpandAllocations.hs
+++ b/src/Futhark/Pass/ExpandAllocations.hs
@@ -84,15 +84,15 @@
 -- code versions.  If so, we can remove the offending branch.  Only if
 -- both versions fail do we propagate the error.
 transformStm (Let pat aux (If cond tbranch fbranch (IfDec ts IfEquiv))) = do
-  tbranch' <- (Right <$> transformBody tbranch) `catchError` (return . Left)
-  fbranch' <- (Right <$> transformBody fbranch) `catchError` (return . Left)
+  tbranch' <- (Right <$> transformBody tbranch) `catchError` (pure . Left)
+  fbranch' <- (Right <$> transformBody fbranch) `catchError` (pure . Left)
   case (tbranch', fbranch') of
     (Left _, Right fbranch'') ->
-      return $ useBranch fbranch''
+      pure $ useBranch fbranch''
     (Right tbranch'', Left _) ->
-      return $ useBranch tbranch''
+      pure $ useBranch tbranch''
     (Right tbranch'', Right fbranch'') ->
-      return $ oneStm $ Let pat aux $ If cond tbranch'' fbranch'' (IfDec ts IfEquiv)
+      pure $ oneStm $ Let pat aux $ If cond tbranch'' fbranch'' (IfDec ts IfEquiv)
     (Left e, _) ->
       throwError e
   where
@@ -104,7 +104,7 @@
         <> stmsFromList (zipWith bindRes (patElems pat) (bodyResult b))
 transformStm (Let pat aux e) = do
   (stms, e') <- transformExp =<< mapExpM transform e
-  return $ stms <> oneStm (Let pat aux e')
+  pure $ stms <> oneStm (Let pat aux e')
   where
     transform =
       identityMapper
@@ -114,7 +114,7 @@
 transformExp :: Exp GPUMem -> ExpandM (Stms GPUMem, Exp GPUMem)
 transformExp (Op (Inner (SegOp (SegMap lvl space ts kbody)))) = do
   (alloc_stms, (_, kbody')) <- transformScanRed lvl space [] kbody
-  return
+  pure
     ( alloc_stms,
       Op $ Inner $ SegOp $ SegMap lvl space ts kbody'
     )
@@ -122,7 +122,7 @@
   (alloc_stms, (lams, kbody')) <-
     transformScanRed lvl space (map segBinOpLambda reds) kbody
   let reds' = zipWith (\red lam -> red {segBinOpLambda = lam}) reds lams
-  return
+  pure
     ( alloc_stms,
       Op $ Inner $ SegOp $ SegRed lvl space reds' ts kbody'
     )
@@ -130,14 +130,14 @@
   (alloc_stms, (lams, kbody')) <-
     transformScanRed lvl space (map segBinOpLambda scans) kbody
   let scans' = zipWith (\red lam -> red {segBinOpLambda = lam}) scans lams
-  return
+  pure
     ( alloc_stms,
       Op $ Inner $ SegOp $ SegScan lvl space scans' ts kbody'
     )
 transformExp (Op (Inner (SegOp (SegHist lvl space ops ts kbody)))) = do
   (alloc_stms, (lams', kbody')) <- transformScanRed lvl space lams kbody
   let ops' = zipWith onOp ops lams'
-  return
+  pure
     ( alloc_stms,
       Op $ Inner $ SegOp $ SegHist lvl space ops' ts kbody'
     )
@@ -172,7 +172,7 @@
             "Cannot handle un-sliceable allocation size: " ++ pretty v
               ++ "\nLikely cause: irregular nested operations inside accumulator update operator."
         [] ->
-          return ()
+          pure ()
 
       let num_is = shapeRank shape
           is = map paramName $ take num_is $ lambdaParams op_lam
@@ -184,9 +184,9 @@
       either throwError pure $
         runOffsetM scope' alloc_offsets $ do
           op_lam'' <- offsetMemoryInLambda op_lam'
-          return (alloc_stms, (shape, arrs, Just (op_lam'', nes)))
+          pure (alloc_stms, (shape, arrs, Just (op_lam'', nes)))
 transformExp e =
-  return (mempty, e)
+  pure (mempty, e)
 
 transformScanRed ::
   SegLevel ->
@@ -213,19 +213,19 @@
         "Cannot handle un-sliceable allocation size: " ++ pretty v
           ++ "\nLikely cause: irregular nested operations inside parallel constructs."
     Nothing ->
-      return ()
+      pure ()
 
   case lvl of
     SegGroup {}
       | not $ null variant_allocs ->
         throwError "Cannot handle invariant allocations in SegGroup."
     _ ->
-      return ()
+      pure ()
 
   allocsForBody variant_allocs invariant_allocs lvl space kbody' $ \alloc_stms kbody'' -> do
     ops'' <- forM ops' $ \op' ->
       localScope (scopeOf op') $ offsetMemoryInLambda op'
-    return (alloc_stms, (ops'', kbody''))
+    pure (alloc_stms, (ops'', kbody''))
   where
     bound_in_kernel =
       namesFromList (M.keys $ scopeOfSegSpace space)
@@ -289,7 +289,7 @@
         kstms
         variant_allocs
 
-  return
+  pure
     ( invariant_alloc_offsets <> variant_alloc_offsets,
       num_threads_stms <> invariant_alloc_stms <> variant_alloc_stms
     )
@@ -374,7 +374,7 @@
       -- scalar allocations.
       || (boundInKernel size && notScalar space) = do
     tell $ M.singleton (patElemName patElem) (user, size, space)
-    return Nothing
+    pure Nothing
   where
     expandableSize (Var v) = v `nameIn` bound_outside || v `nameIn` bound_kernel
     expandableSize Constant {} = True
@@ -382,7 +382,7 @@
     boundInKernel Constant {} = False
 extractStmAllocations user bound_outside bound_kernel stm = do
   e <- mapExpM (expMapper user) $ stmExp stm
-  return $ Just $ stm {stmExp = e}
+  pure $ Just $ stm {stmExp = e}
   where
     expMapper user' =
       identityMapper
@@ -393,14 +393,14 @@
     onBody user' body = do
       let (body', allocs) = extractBodyAllocations user' bound_outside bound_kernel body
       tell allocs
-      return body'
+      pure body'
 
     onOp (_, user_ids) (Inner (SegOp op)) =
       Inner . SegOp <$> mapSegOpM (opMapper user'') op
       where
         user'' =
           (segLevel op, user_ids ++ [le64 (segFlat (segSpace op))])
-    onOp _ op = return op
+    onOp _ op = pure op
 
     opMapper user' =
       identitySegOpMapper
@@ -411,11 +411,11 @@
     onKernelBody user' body = do
       let (body', allocs) = extractKernelBodyAllocations user' bound_outside bound_kernel body
       tell allocs
-      return body'
+      pure body'
 
     onLambda user' lam = do
       body <- onBody user' $ lambdaBody lam
-      return lam {lambdaBody = body}
+      pure lam {lambdaBody = body}
 
 genericExpandedInvariantAllocations ::
   (User -> (Shape, [TPrimExp Int64 VName])) -> Extraction -> ExpandM (Stms GPUMem, RebaseMap)
@@ -424,7 +424,7 @@
   -- equal to the number of kernel threads.
   (rebases, alloc_stms) <- runBuilder $ mapM expand $ M.toList invariant_allocs
 
-  return (alloc_stms, mconcat rebases)
+  pure (alloc_stms, mconcat rebases)
   where
     expand (mem, (user, per_thread_size, space)) = do
       let num_users = fst $ getNumUsers user
@@ -476,7 +476,7 @@
   Extraction ->
   ExpandM (Stms GPUMem, RebaseMap)
 expandedVariantAllocations _ _ _ variant_allocs
-  | null variant_allocs = return (mempty, mempty)
+  | null variant_allocs = pure (mempty, mempty)
 expandedVariantAllocations num_threads kspace kstms variant_allocs = do
   let sizes_to_blocks = removeCommonSizes variant_allocs
       variant_sizes = map fst sizes_to_blocks
@@ -502,11 +502,11 @@
   -- equal to the sum of the sizes required by different threads.
   (alloc_stms, rebases) <- unzip <$> mapM expand variant_allocs'
 
-  return (slice_stms' <> stmsFromList alloc_stms, mconcat rebases)
+  pure (slice_stms' <> stmsFromList alloc_stms, mconcat rebases)
   where
     expand (mem, (offset, total_size, space)) = do
       let allocpat = Pat [PatElem mem $ MemMem space]
-      return
+      pure
         ( Let allocpat (defAux ()) $ Op $ Alloc total_size space,
           M.singleton mem $ newBase offset
         )
@@ -563,7 +563,7 @@
 lookupNewBase :: VName -> ([TPrimExp Int64 VName], PrimType) -> OffsetM (Maybe IxFun)
 lookupNewBase name x = do
   offsets <- askRebaseMap
-  return $ ($ x) <$> M.lookup name offsets
+  pure $ ($ x) <$> M.lookup name offsets
 
 offsetMemoryInKernelBody :: KernelBody GPUMem -> OffsetM (KernelBody GPUMem)
 offsetMemoryInKernelBody kbody = do
@@ -574,7 +574,7 @@
         (\scope' -> localScope scope' . offsetMemoryInStm)
         scope
         (stmsToList $ kernelBodyStms kbody)
-  return kbody {kernelBodyStms = stms'}
+  pure kbody {kernelBodyStms = stms'}
 
 offsetMemoryInBody :: Body GPUMem -> OffsetM (Body GPUMem)
 offsetMemoryInBody (Body dec stms res) = do
@@ -585,7 +585,7 @@
         (\scope' -> localScope scope' . offsetMemoryInStm)
         scope
         (stmsToList stms)
-  return $ Body dec stms' res
+  pure $ Body dec stms' res
 
 offsetMemoryInStm :: Stm GPUMem -> OffsetM (Scope GPUMem, Stm GPUMem)
 offsetMemoryInStm (Let pat dec e) = do
@@ -598,7 +598,7 @@
   let pat'' = Pat $ zipWith pick (patElems pat') rts
       stm = Let pat'' dec e'
   let scope' = scopeOf stm <> scope
-  return (scope', stm)
+  pure (scope', stm)
   where
     pick ::
       PatElemT (MemInfo SubExp NoUniqueness MemBind) ->
@@ -615,7 +615,7 @@
     instantiateIxFun = traverse (traverse inst)
       where
         inst Ext {} = Nothing
-        inst (Free x) = return x
+        inst (Free x) = pure x
 
 offsetMemoryInPat :: Pat GPUMem -> [ExpReturns] -> OffsetM (Pat GPUMem)
 offsetMemoryInPat (Pat pes) rets = do
@@ -635,30 +635,30 @@
 offsetMemoryInParam :: Param (MemBound u) -> OffsetM (Param (MemBound u))
 offsetMemoryInParam fparam = do
   fparam' <- offsetMemoryInMemBound $ paramDec fparam
-  return fparam {paramDec = fparam'}
+  pure fparam {paramDec = fparam'}
 
 offsetMemoryInMemBound :: MemBound u -> OffsetM (MemBound u)
 offsetMemoryInMemBound summary@(MemArray pt shape u (ArrayIn mem ixfun)) = do
   new_base <- lookupNewBase mem (IxFun.base ixfun, pt)
-  return . fromMaybe summary $ do
+  pure . fromMaybe summary $ do
     new_base' <- new_base
-    return $ MemArray pt shape u $ ArrayIn mem $ IxFun.rebase new_base' ixfun
-offsetMemoryInMemBound summary = return summary
+    pure $ MemArray pt shape u $ ArrayIn mem $ IxFun.rebase new_base' ixfun
+offsetMemoryInMemBound summary = pure summary
 
 offsetMemoryInBodyReturns :: BodyReturns -> OffsetM BodyReturns
 offsetMemoryInBodyReturns br@(MemArray pt shape u (ReturnsInBlock mem ixfun))
   | Just ixfun' <- isStaticIxFun ixfun = do
     new_base <- lookupNewBase mem (IxFun.base ixfun', pt)
-    return . fromMaybe br $ do
+    pure . fromMaybe br $ do
       new_base' <- new_base
-      return . MemArray pt shape u . ReturnsInBlock mem $
+      pure . MemArray pt shape u . ReturnsInBlock mem $
         IxFun.rebase (fmap (fmap Free) new_base') ixfun
-offsetMemoryInBodyReturns br = return br
+offsetMemoryInBodyReturns br = pure br
 
 offsetMemoryInLambda :: Lambda GPUMem -> OffsetM (Lambda GPUMem)
 offsetMemoryInLambda lam = inScopeOf lam $ do
   body <- offsetMemoryInBody $ lambdaBody lam
-  return $ lam {lambdaBody = body}
+  pure $ lam {lambdaBody = body}
 
 -- A loop may have memory parameters, and those memory blocks may
 -- be expanded.  We assume (but do not check - FIXME) that if the
@@ -687,7 +687,7 @@
       localScope
         (scopeOfFParams (map fst merge') <> scopeOf form)
         (offsetMemoryInBody body)
-    return $ DoLoop merge' form body'
+    pure $ DoLoop merge' form body'
 offsetMemoryInExp e = mapExpM recurse e
   where
     recurse =
@@ -705,7 +705,7 @@
             { mapOnSegOpBody = offsetMemoryInKernelBody,
               mapOnSegOpLambda = offsetMemoryInLambda
             }
-    onOp op = return op
+    onOp op = pure op
 
 ---- Slicing allocation sizes out of a kernel.
 
@@ -723,24 +723,19 @@
 
     unAllocStm nested stm@(Let _ _ (Op Alloc {}))
       | nested = throwError $ "Cannot handle nested allocation: " ++ pretty stm
-      | otherwise = return Nothing
+      | otherwise = pure Nothing
     unAllocStm _ (Let pat dec e) =
       Just <$> (Let <$> unAllocPat pat <*> pure dec <*> mapExpM unAlloc' e)
 
     unAllocLambda (Lambda params body ret) =
-      Lambda (unParams params) <$> unAllocBody body <*> pure ret
-
-    unParams = mapMaybe $ traverse unMem
+      Lambda (map unParam params) <$> unAllocBody body <*> pure ret
 
-    unAllocPat pat@(Pat merge) =
-      Pat <$> maybe bad return (mapM (rephrasePatElem unMem) merge)
-      where
-        bad = Left $ "Cannot handle memory in pattern " ++ pretty pat
+    unAllocPat (Pat pes) =
+      Pat <$> mapM (rephrasePatElem (Right . unMem)) pes
 
     unAllocOp Alloc {} = Left "unAllocOp: unhandled Alloc"
     unAllocOp (Inner OtherOp {}) = Left "unAllocOp: unhandled OtherOp"
-    unAllocOp (Inner (SizeOp op)) =
-      return $ SizeOp op
+    unAllocOp (Inner (SizeOp op)) = pure $ SizeOp op
     unAllocOp (Inner (SegOp op)) = SegOp <$> mapSegOpM mapper op
       where
         mapper =
@@ -749,43 +744,37 @@
               mapOnSegOpBody = unAllocKernelBody
             }
 
-    unParam p = maybe bad return $ traverse unMem p
-      where
-        bad = Left $ "Cannot handle memory-typed parameter '" ++ pretty p ++ "'"
+    unParam = fmap unMem
 
-    unT t = maybe bad return $ unMem t
-      where
-        bad = Left $ "Cannot handle memory type '" ++ pretty t ++ "'"
+    unT = Right . unMem
 
     unAlloc' =
       Mapper
         { mapOnBody = const unAllocBody,
           mapOnRetType = unT,
           mapOnBranchType = unT,
-          mapOnFParam = unParam,
-          mapOnLParam = unParam,
+          mapOnFParam = Right . unParam,
+          mapOnLParam = Right . unParam,
           mapOnOp = unAllocOp,
           mapOnSubExp = Right,
           mapOnVName = Right
         }
 
-unMem :: MemInfo d u ret -> Maybe (TypeBase (ShapeBase d) u)
-unMem (MemPrim pt) = Just $ Prim pt
-unMem (MemArray pt shape u _) = Just $ Array pt shape u
-unMem (MemAcc acc ispace ts u) = Just $ Acc acc ispace ts u
-unMem MemMem {} = Nothing
+unMem :: MemInfo d u ret -> TypeBase (ShapeBase d) u
+unMem (MemPrim pt) = Prim pt
+unMem (MemArray pt shape u _) = Array pt shape u
+unMem (MemAcc acc ispace ts u) = Acc acc ispace ts u
+unMem MemMem {} = Prim Unit
 
 unAllocScope :: Scope GPUMem -> Scope GPU.GPU
-unAllocScope = M.mapMaybe unInfo
+unAllocScope = M.map unInfo
   where
-    unInfo (LetName dec) = LetName <$> unMem dec
-    unInfo (FParamName dec) = FParamName <$> unMem dec
-    unInfo (LParamName dec) = LParamName <$> unMem dec
-    unInfo (IndexName it) = Just $ IndexName it
+    unInfo (LetName dec) = LetName $ unMem dec
+    unInfo (FParamName dec) = FParamName $ unMem dec
+    unInfo (LParamName dec) = LParamName $ unMem dec
+    unInfo (IndexName it) = IndexName it
 
-removeCommonSizes ::
-  Extraction ->
-  [(SubExp, [(VName, Space)])]
+removeCommonSizes :: Extraction -> [(SubExp, [(VName, Space)])]
 removeCommonSizes = M.toList . foldl' comb mempty . M.toList
   where
     comb m (mem, (_, size, space)) = M.insertWith (++) size [(mem, space)] m
@@ -797,7 +786,7 @@
   Stms GPUMem ->
   ExpandM (Stms GPU.GPU, [VName], [VName])
 sliceKernelSizes num_threads sizes space kstms = do
-  kstms' <- either throwError return $ unAllocGPUStms kstms
+  kstms' <- either throwError pure $ unAllocGPUStms kstms
   let num_sizes = length sizes
       i64s = replicate num_sizes $ Prim int64
 
@@ -811,7 +800,7 @@
         forM (zip xs ys) $ \(x, y) ->
           fmap subExpRes . letSubExp "z" . BasicOp $
             BinOp (SMax Int64) (Var $ paramName x) (Var $ paramName y)
-    return $ Lambda (xs ++ ys) (mkBody stms zs) i64s
+    pure $ Lambda (xs ++ ys) (mkBody stms zs) i64s
 
   flat_gtid_lparam <- newParam "flat_gtid" (Prim (IntType Int64))
 
@@ -831,7 +820,7 @@
         zipWithM_ letBindNames (map pure kspace_gtids) =<< mapM toExp new_inds
 
         mapM_ addStm kstms'
-        return $ subExpsRes sizes
+        pure $ subExpsRes sizes
 
     localScope (scopeOfSegSpace space) $
       GPU.simplifyLambda (Lambda [flat_gtid_lparam] (Body () stms zs) i64s)
@@ -846,8 +835,7 @@
 
     thread_space_iota <-
       letExp "thread_space_iota" $
-        BasicOp $
-          Iota w (intConst Int64 0) (intConst Int64 1) Int64
+        BasicOp $ Iota w (intConst Int64 0) (intConst Int64 1) Int64
     let red_op =
           SegBinOp
             Commutative
@@ -863,6 +851,6 @@
       letExp "size_sum" $
         BasicOp $ BinOp (Mul Int64 OverflowUndef) (Var threads_max) num_threads
 
-    return (patNames pat, size_sums)
+    pure (patNames pat, size_sums)
 
-  return (slice_stms, maxes_per_thread, size_sums)
+  pure (slice_stms, maxes_per_thread, size_sums)
diff --git a/src/Futhark/Test.hs b/src/Futhark/Test.hs
--- a/src/Futhark/Test.hs
+++ b/src/Futhark/Test.hs
@@ -78,20 +78,18 @@
 -- executable, and the second is the directory which file paths are
 -- read relative to.
 getValues :: (MonadFail m, MonadIO m) => FutharkExe -> FilePath -> Values -> m [V.Value]
-getValues _ _ (Values vs) =
-  pure vs
+getValues _ _ (Values vs) = pure vs
 getValues futhark dir v = do
   s <- getValuesBS futhark dir v
-  case valuesFromByteString file s of
+  case valuesFromByteString (fileName v) s of
     Left e -> fail e
     Right vs -> pure vs
   where
-    file = case v of
-      Values {} -> "<values>"
-      InFile f -> f
-      GenValues {} -> "<randomly generated>"
-      ScriptValues {} -> "<FutharkScript expression>"
-      ScriptFile f -> f
+    fileName Values {} = "<values>"
+    fileName GenValues {} = "<randomly generated>"
+    fileName ScriptValues {} = "<FutharkScript expression>"
+    fileName (InFile f) = f
+    fileName (ScriptFile f) = f
 
 readAndDecompress :: FilePath -> IO (Either DecompressError BS.ByteString)
 readAndDecompress file = E.try $ do
diff --git a/src/Language/Futhark/Parser/Lexer.x b/src/Language/Futhark/Parser/Lexer.x
--- a/src/Language/Futhark/Parser/Lexer.x
+++ b/src/Language/Futhark/Parser/Lexer.x
@@ -382,6 +382,22 @@
                     , alex_scd = 0}) of Left msg -> Left msg
                                         Right ( _, a ) -> Right a
 
+
+getToken :: FilePath -> Alex (Lexeme Token)
+getToken file = do
+  inp__@(_,_,_,n) <- alexGetInput
+  sc <- alexGetStartCode
+  case alexScan inp__ sc of
+    AlexEOF -> alexEOF
+    AlexError ((AlexPn _ line column),_,_,_) ->
+      alexError $ "Error at " ++ file ++ ":" ++ show line ++ ":" ++ show column ++ ": lexical error."
+    AlexSkip  inp__' _len -> do
+      alexSetInput inp__'
+      getToken file
+    AlexToken inp__'@(_,_,_,n') _ action -> let len = n'-n in do
+      alexSetInput inp__'
+      action (ignorePendingBytes inp__) len
+
 -- | Given a starting position, produce tokens from the given text (or
 -- a lexer error).  Returns the final position.
 scanTokensText :: Pos -> T.Text -> Either String ([L Token], Pos)
@@ -391,7 +407,7 @@
 scanTokens (Pos file start_line start_col start_off) str =
   runAlex' (AlexPn start_off start_line start_col) str $ do
   fix $ \loop -> do
-    tok <- alexMonadScan
+    tok <- getToken file
     case tok of
       (start, end, EOF) ->
         return ([], posnToPos end)
diff --git a/src/Language/Futhark/Parser/Parser.y b/src/Language/Futhark/Parser/Parser.y
--- a/src/Language/Futhark/Parser/Parser.y
+++ b/src/Language/Futhark/Parser/Parser.y
@@ -53,6 +53,7 @@
 
 %tokentype { L Token }
 %error { parseError }
+%errorhandlertype explist
 %monad { ParserMonad }
 %lexer { lexer } { L _ EOF }
 
@@ -1195,13 +1196,18 @@
       putTokens (xs, pos)
       cont x
 
-parseError :: L Token -> ParserMonad a
-parseError (L loc EOF) =
-  parseErrorAt (srclocOf loc) $ Just "unexpected end of file."
-parseError (L loc DOC{}) =
+parseError :: (L Token, [String]) -> ParserMonad a
+parseError (L loc EOF, expected) =
+  parseErrorAt (srclocOf loc) $ Just $
+  unlines ["unexpected end of file.",
+           "Expected one of the following: " ++ unwords expected]
+parseError (L loc DOC{}, _) =
   parseErrorAt (srclocOf loc) $
   Just "documentation comments ('-- |') are only permitted when preceding declarations."
-parseError tok = parseErrorAt (srclocOf tok) Nothing
+parseError (L loc tok, expected) =
+  parseErrorAt loc $ Just $
+  unlines ["unexpected " ++ show tok,
+          "Expected one of the following: " ++ unwords expected]
 
 parseErrorAt :: SrcLoc -> Maybe String -> ParserMonad a
 parseErrorAt loc Nothing = throwError $ "Error at " ++ locStr loc ++ ": Parse error."
diff --git a/src/Language/Futhark/TypeChecker/Terms/Monad.hs b/src/Language/Futhark/TypeChecker/Terms/Monad.hs
--- a/src/Language/Futhark/TypeChecker/Terms/Monad.hs
+++ b/src/Language/Futhark/TypeChecker/Terms/Monad.hs
@@ -356,7 +356,7 @@
     where
       expected' = commasep (map ppr expected)
   ppr (CheckingBranches t1 t2) =
-    "Conditional branches differ in type."
+    "Branches differ in type."
       </> "Former:" <+> ppr t1
       </> "Latter:" <+> ppr t2
 
diff --git a/src/Language/Futhark/TypeChecker/Unify.hs b/src/Language/Futhark/TypeChecker/Unify.hs
--- a/src/Language/Futhark/TypeChecker/Unify.hs
+++ b/src/Language/Futhark/TypeChecker/Unify.hs
@@ -266,6 +266,23 @@
     . S.toList
     . typeDimNames
 
+typeVarNotes :: MonadUnify m => VName -> m Notes
+typeVarNotes v = maybe mempty (aNote . note . snd) . M.lookup v <$> getConstraints
+  where
+    note (HasConstrs cs _) =
+      pprName v <+> "="
+        <+> mconcat (map ppConstr (M.toList cs))
+        <+> "..."
+    note (Overloaded ts _) =
+      pprName v <+> "must be one of" <+> mconcat (punctuate ", " (map ppr ts))
+    note (HasFields fs _) =
+      pprName v <+> "="
+        <+> braces (mconcat (punctuate ", " (map ppField (M.toList fs))))
+    note _ = mempty
+
+    ppConstr (c, _) = "#" <> ppr c <+> "..." <+> "|"
+    ppField (f, _) = pprName f <> ":" <+> "..."
+
 -- | Monads that which to perform unification must implement this type
 -- class.
 class Monad m => MonadUnify m where
@@ -725,29 +742,32 @@
             modifyConstraints $
               M.insert vn (lvl, Constraint (RetType ext tp') usage)
             unifySharedConstructors onDims usage bound bcs required_cs ts
-        Scalar (TypeVar _ _ (TypeName [] v) [])
-          | not $ isRigid v constraints -> do
-            link
-            case M.lookup v constraints of
-              Just (_, HasConstrs v_cs _) ->
-                unifySharedConstructors onDims usage bound bcs required_cs v_cs
-              _ -> pure ()
-            modifyConstraints $
-              M.insertWith
-                combineConstrs
-                v
-                (lvl, HasConstrs required_cs old_usage)
+        Scalar (TypeVar _ _ (TypeName [] v) []) -> do
+          case M.lookup v constraints of
+            Just (_, HasConstrs v_cs _) -> do
+              unifySharedConstructors onDims usage bound bcs required_cs v_cs
+            Just (_, NoConstraint {}) -> pure ()
+            Just (_, Equality {}) -> pure ()
+            _ -> do
+              notes <- (<>) <$> typeVarNotes vn <*> typeVarNotes v
+              noSumType notes
+          link
+          modifyConstraints $
+            M.insertWith
+              combineConstrs
+              v
+              (lvl, HasConstrs required_cs old_usage)
           where
             combineConstrs (_, HasConstrs cs1 usage1) (_, HasConstrs cs2 _) =
               (lvl, HasConstrs (M.union cs1 cs2) usage1)
             combineConstrs hasCs _ = hasCs
-        _ -> noSumType
+        _ -> noSumType mempty
     _ -> link
   where
-    noSumType =
+    noSumType notes =
       unifyError
         usage
-        mempty
+        notes
         bcs
         "Cannot unify a sum type with a non-sum type"
 
@@ -998,9 +1018,8 @@
           | otherwise ->
             unifyError usage mempty noBreadCrumbs $
               "Different arity for constructor" <+> pquote (ppr c) <+> "."
-    _ -> do
+    _ ->
       unify usage t $ Scalar $ Sum $ M.singleton c fs
-      return ()
 
 mustHaveFieldWith ::
   MonadUnify m =>