diff --git a/Changes.md b/Changes.md
--- a/Changes.md
+++ b/Changes.md
@@ -1,5 +1,11 @@
 # Change log for the `lapack` package
 
+## 0.3.1
+
+ * `Matrix.Symmetric`:
+   You can now import many functions for symmetric matrices from this module.
+   This is more natural than importing them from `Triangular`.
+
 ## 0.3
 
  * Matrix data family
diff --git a/lapack.cabal b/lapack.cabal
--- a/lapack.cabal
+++ b/lapack.cabal
@@ -1,5 +1,5 @@
 Name:             lapack
-Version:          0.3.0.1
+Version:          0.3.1
 License:          BSD3
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -48,7 +48,7 @@
   Makefile
 
 Source-Repository this
-  Tag:         0.3.0.1
+  Tag:         0.3.1
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/lapack/
 
@@ -89,6 +89,7 @@
     Numeric.LAPACK.Matrix.Square
     Numeric.LAPACK.Matrix.Hermitian
     Numeric.LAPACK.Matrix.HermitianPositiveDefinite
+    Numeric.LAPACK.Matrix.Symmetric
     Numeric.LAPACK.Matrix.Triangular
     Numeric.LAPACK.Matrix.Banded
     Numeric.LAPACK.Matrix.BandedHermitian
@@ -126,6 +127,7 @@
     Numeric.LAPACK.Matrix.Hermitian.Linear
     Numeric.LAPACK.Matrix.Hermitian.Eigen
     Numeric.LAPACK.Matrix.HermitianPositiveDefinite.Linear
+    Numeric.LAPACK.Matrix.Symmetric.Basic
     Numeric.LAPACK.Matrix.Symmetric.Private
     Numeric.LAPACK.Matrix.Banded.Basic
     Numeric.LAPACK.Matrix.Banded.Linear
@@ -194,6 +196,7 @@
     Test.Matrix
     Test.Square
     Test.Triangular
+    Test.Symmetric
     Test.Hermitian
     Test.Orthogonal
     Test.LowerUpper
diff --git a/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs b/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
@@ -545,7 +545,8 @@
 congruence b =
    skipCheckCongruence Basic.mapWidth $ \a ->
       scaledAnticommutator one a $
-      Split.tallMultiplyR NonTransposed (takeHalf b) a
+      Split.tallMultiplyR NonTransposed
+         (Split.takeHalf MatrixShape.hermitianOrder b) a
 
 congruenceAdjoint ::
    (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
@@ -553,31 +554,8 @@
 congruenceAdjoint =
    flip $ \b -> skipCheckCongruence Basic.mapHeight $ \a ->
       scaledAnticommutatorAdjoint one a $
-      Basic.swapMultiply (Split.tallMultiplyR Transposed) a (takeHalf b)
-
-
-data Corrupt = Corrupt
-   deriving (Eq)
-
-{- |
-> let b = takeHalf a
-> ==>
-> isTriangular b && a == addAdjoint b
--}
-takeHalf ::
-   (Shape.C sh, Class.Floating a) =>
-   Hermitian sh a -> Split.Square Corrupt sh a
-takeHalf (Array (MatrixShape.Hermitian order sh) a) =
-   Array.unsafeCreate (MatrixShape.Split Corrupt order (Extent.square sh)) $
-      \bPtr -> evalContT $ do
-   let n = Shape.size sh
-   aPtr <- ContT $ withForeignPtr a
-   nPtr <- Call.cint n
-   alphaPtr <- Call.number 0.5
-   incxPtr <- Call.cint (n+1)
-   liftIO $ do
-      unpack order n aPtr bPtr
-      BlasGen.scal nPtr alphaPtr bPtr incxPtr
+      Basic.swapMultiply (Split.tallMultiplyR Transposed)
+         a (Split.takeHalf MatrixShape.hermitianOrder b)
 
 
 scaledAnticommutator ::
diff --git a/src/Numeric/LAPACK/Matrix/Symmetric.hs b/src/Numeric/LAPACK/Matrix/Symmetric.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Symmetric.hs
@@ -0,0 +1,178 @@
+{-# LANGUAGE TypeOperators #-}
+module Numeric.LAPACK.Matrix.Symmetric (
+   Symmetric,
+   size,
+   fromList, autoFromList,
+   identity,
+   diagonal,
+   takeDiagonal,
+   transpose,
+   adjoint,
+
+   stack, (#%%%#),
+   split,
+
+   toSquare,
+
+   gramian,            gramianTransposed,
+   congruenceDiagonal, congruenceDiagonalTransposed,
+   congruence,         congruenceTransposed,
+   anticommutator,     anticommutatorTransposed,
+   ) where
+
+import qualified Numeric.LAPACK.Matrix.Symmetric.Basic as Basic
+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
+
+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Extent as Extent
+import Numeric.LAPACK.Matrix.Array.Triangular (Symmetric)
+import Numeric.LAPACK.Matrix.Array (Full, General, Square)
+import Numeric.LAPACK.Matrix.Shape.Private (Order)
+import Numeric.LAPACK.Matrix.Private (ShapeInt)
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Scalar (one)
+
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Shape ((:+:))
+
+import Foreign.Storable (Storable)
+
+
+size :: Symmetric sh a -> sh
+size = MatrixShape.triangularSize . ArrMatrix.shape
+
+transpose :: Symmetric sh a -> Symmetric sh a
+transpose = Triangular.transpose
+
+adjoint :: (Shape.C sh, Class.Floating a) => Symmetric sh a -> Symmetric sh a
+adjoint = Triangular.adjoint
+
+
+fromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Symmetric sh a
+fromList = Triangular.symmetricFromList
+
+autoFromList :: (Storable a) => Order -> [a] -> Symmetric ShapeInt a
+autoFromList = Triangular.autoSymmetricFromList
+
+
+toSquare :: (Shape.C sh, Class.Floating a) => Symmetric sh a -> Square sh a
+toSquare = Triangular.toSquare
+
+
+identity :: (Shape.C sh, Class.Floating a) => Order -> sh -> Symmetric sh a
+identity order = Triangular.relaxUnitDiagonal . Triangular.identity order
+
+diagonal ::
+   (Shape.C sh, Class.Floating a) => Order -> Vector sh a -> Symmetric sh a
+diagonal = Triangular.diagonal
+
+takeDiagonal :: (Shape.C sh, Class.Floating a) => Symmetric sh a -> Vector sh a
+takeDiagonal = Triangular.takeDiagonal
+
+
+stack ::
+   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>
+   Symmetric sh0 a ->
+   General sh0 sh1 a ->
+   Symmetric sh1 a ->
+   Symmetric (sh0:+:sh1) a
+stack = Triangular.stackSymmetric
+
+infixr 2 #%%%#
+
+(#%%%#) ::
+   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>
+   (Symmetric sh0 a, General sh0 sh1 a) ->
+   Symmetric sh1 a ->
+   Symmetric (sh0:+:sh1) a
+(#%%%#) = uncurry stack
+
+
+split ::
+   (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1, Class.Floating a) =>
+   Symmetric (sh0:+:sh1) a ->
+   (Symmetric sh0 a, General sh0 sh1 a, Symmetric sh1 a)
+split = Triangular.splitSymmetric
+
+
+
+{- |
+gramian A = A^T * A
+-}
+gramian ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   General height width a -> Symmetric width a
+gramian = ArrMatrix.lift1 Basic.gramian
+
+{- |
+gramianTransposed A = A * A^T = gramian (A^T)
+-}
+gramianTransposed ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   General height width a -> Symmetric height a
+gramianTransposed = ArrMatrix.lift1 Basic.gramianTransposed
+
+{- |
+congruenceDiagonal D A = A^T * D * A
+-}
+congruenceDiagonal ::
+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
+   Vector height a -> General height width a -> Symmetric width a
+congruenceDiagonal = ArrMatrix.lift1 . Basic.congruenceDiagonal
+
+{- |
+congruenceDiagonalTransposed A D = A * D * A^T
+-}
+congruenceDiagonalTransposed ::
+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
+   General height width a -> Vector width a -> Symmetric height a
+congruenceDiagonalTransposed a =
+   ArrMatrix.lift0 . Basic.congruenceDiagonalTransposed (ArrMatrix.toVector a)
+
+{- |
+congruence B A = A^T * B * A
+-}
+congruence ::
+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
+   Symmetric height a -> General height width a -> Symmetric width a
+congruence = ArrMatrix.lift2 Basic.congruence
+
+{- |
+congruenceTransposed B A = A * B * A^T
+-}
+congruenceTransposed ::
+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
+   General height width a -> Symmetric width a -> Symmetric height a
+congruenceTransposed = ArrMatrix.lift2 Basic.congruenceTransposed
+
+
+{- |
+anticommutator A B  =  A^T * B + B^T * A
+
+Not exactly a matrix anticommutator,
+thus I like to call it Symmetric anticommutator.
+-}
+anticommutator ::
+   (Extent.C vert, Extent.C horiz,
+    Shape.C height, Eq height, Shape.C width, Eq width,
+    Class.Floating a) =>
+   Full vert horiz height width a ->
+   Full vert horiz height width a -> Symmetric width a
+anticommutator = ArrMatrix.lift2 $ Basic.scaledAnticommutator one
+
+{- |
+anticommutatorTransposed A B
+   = A * B^T + B * A^T
+   = anticommutator (transpose A) (transpose B)
+-}
+anticommutatorTransposed ::
+   (Extent.C vert, Extent.C horiz,
+    Shape.C height, Eq height, Shape.C width, Eq width,
+    Class.Floating a) =>
+   Full vert horiz height width a ->
+   Full vert horiz height width a -> Symmetric height a
+anticommutatorTransposed =
+   ArrMatrix.lift2 $ Basic.scaledAnticommutatorTransposed one
diff --git a/src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs b/src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs
@@ -0,0 +1,209 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Numeric.LAPACK.Matrix.Symmetric.Basic (
+   Symmetric,
+
+   gramian,              gramianTransposed,
+   congruenceDiagonal,   congruenceDiagonalTransposed,
+   congruence,           congruenceTransposed,
+   scaledAnticommutator, scaledAnticommutatorTransposed,
+   ) where
+
+import qualified Numeric.LAPACK.Matrix.Private as Matrix
+import qualified Numeric.LAPACK.Matrix.Basic as Basic
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
+import qualified Numeric.LAPACK.Split as Split
+import Numeric.LAPACK.Matrix.Triangular.Private (pack, recheck)
+import Numeric.LAPACK.Matrix.Shape.Private
+         (Order(RowMajor,ColumnMajor), NonUnit(NonUnit))
+import Numeric.LAPACK.Matrix.Modifier
+         (Transposition(NonTransposed, Transposed))
+import Numeric.LAPACK.Matrix.Private (Full)
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Scalar (zero, one)
+import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))
+
+import qualified Numeric.BLAS.FFI.Generic as BlasGen
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+
+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
+import Foreign.Ptr (Ptr)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+
+
+type Symmetric sh = Array (MatrixShape.FlexSymmetric NonUnit sh)
+
+
+-- cf. Hermitian.Basic
+gramian ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   Matrix.General height width a -> Symmetric width a
+gramian (Array (MatrixShape.Full order extent) a) =
+   Array.unsafeCreate (symmetricShape order $ Extent.width extent) $
+   \bPtr -> gramianIO order a bPtr $ gramianParameters order extent
+
+gramianParameters ::
+   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>
+   Order ->
+   Extent.Extent vert horiz height width ->
+   ((Int, Int), (Char, Char, Int))
+gramianParameters order extent =
+   let (height, width) = Extent.dimensions extent
+       n = Shape.size width
+       k = Shape.size height
+    in ((n,k),
+         case order of
+            ColumnMajor -> ('U', 'T', k)
+            RowMajor -> ('L', 'N', n))
+
+
+gramianTransposed ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   Matrix.General height width a -> Symmetric height a
+gramianTransposed (Array (MatrixShape.Full order extent) a) =
+   Array.unsafeCreate (symmetricShape order $ Extent.height extent) $
+   \bPtr -> gramianIO order a bPtr $ gramianTransposedParameters order extent
+
+gramianTransposedParameters ::
+   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>
+   Order ->
+   Extent.Extent vert horiz height width ->
+   ((Int, Int), (Char, Char, Int))
+gramianTransposedParameters order extent =
+   let (height, width) = Extent.dimensions extent
+       n = Shape.size height
+       k = Shape.size width
+   in ((n,k),
+         case order of
+            ColumnMajor -> ('U', 'N', n)
+            RowMajor -> ('L', 'T', k))
+
+gramianIO ::
+   (Class.Floating a) =>
+   Order ->
+   ForeignPtr a -> Ptr a ->
+   ((Int, Int), (Char, Char, Int)) -> IO ()
+gramianIO order a bPtr ((n,k), (uplo,trans,lda)) =
+   evalContT $ do
+      uploPtr <- Call.char uplo
+      transPtr <- Call.char trans
+      nPtr <- Call.cint n
+      kPtr <- Call.cint k
+      alphaPtr <- Call.number one
+      aPtr <- ContT $ withForeignPtr a
+      ldaPtr <- Call.leadingDim lda
+      betaPtr <- Call.number zero
+      cPtr <- Call.allocaArray (n*n)
+      ldcPtr <- Call.leadingDim n
+      liftIO $ do
+         BlasGen.syrk uploPtr transPtr
+            nPtr kPtr alphaPtr aPtr ldaPtr betaPtr cPtr ldcPtr
+         pack order n cPtr bPtr
+
+skipCheckCongruence ::
+   ((sh -> Unchecked sh) -> matrix0 -> matrix1) ->
+   (matrix1 -> Symmetric (Unchecked sh) a) -> matrix0 -> Symmetric sh a
+skipCheckCongruence mapSize f a =
+   recheck $ f $ mapSize Unchecked a
+
+
+congruenceDiagonal ::
+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
+   Vector height a -> Matrix.General height width a -> Symmetric width a
+congruenceDiagonal d =
+   skipCheckCongruence Basic.mapWidth $ \a ->
+      scaledAnticommutator 0.5 a $ Basic.scaleRows d a
+
+congruenceDiagonalTransposed ::
+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
+   Matrix.General height width a -> Vector width a -> Symmetric height a
+congruenceDiagonalTransposed =
+   flip $ \d -> skipCheckCongruence Basic.mapHeight $ \a ->
+      scaledAnticommutatorTransposed 0.5 a $ Basic.scaleColumns d a
+
+
+congruence ::
+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
+   Symmetric height a -> Matrix.General height width a -> Symmetric width a
+congruence b =
+   skipCheckCongruence Basic.mapWidth $ \a ->
+      scaledAnticommutator one a $
+      Split.tallMultiplyR NonTransposed
+         (Split.takeHalf MatrixShape.triangularOrder b) a
+
+congruenceTransposed ::
+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
+   Matrix.General height width a -> Symmetric width a -> Symmetric height a
+congruenceTransposed =
+   flip $ \b -> skipCheckCongruence Basic.mapHeight $ \a ->
+      scaledAnticommutatorTransposed one a $
+      Basic.swapMultiply (Split.tallMultiplyR Transposed)
+         a (Split.takeHalf MatrixShape.triangularOrder b)
+
+
+scaledAnticommutator ::
+   (Extent.C vert, Extent.C horiz,
+    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
+   a ->
+   Full vert horiz height width a ->
+   Full vert horiz height width a -> Symmetric width a
+scaledAnticommutator alpha arr (Array (MatrixShape.Full order extentB) b) = do
+   let (Array (MatrixShape.Full _ extentA) a) = Basic.forceOrder order arr
+   Array.unsafeCreate (symmetricShape order $ Extent.width extentB) $
+         \cpPtr -> do
+      Call.assert "Symmetric.anticommutator: extents mismatch"
+         (extentA==extentB)
+      scaledAnticommutatorIO alpha order a b cpPtr $
+         gramianParameters order extentB
+
+scaledAnticommutatorTransposed ::
+   (Extent.C vert, Extent.C horiz,
+    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
+   a ->
+   Full vert horiz height width a ->
+   Full vert horiz height width a -> Symmetric height a
+scaledAnticommutatorTransposed
+      alpha arr (Array (MatrixShape.Full order extentB) b) = do
+   let (Array (MatrixShape.Full _ extentA) a) = Basic.forceOrder order arr
+   Array.unsafeCreate (symmetricShape order $ Extent.height extentB) $
+         \cpPtr -> do
+      Call.assert "Symmetric.anticommutatorTransposed: extents mismatch"
+         (extentA==extentB)
+      scaledAnticommutatorIO alpha order a b cpPtr $
+         gramianTransposedParameters order extentB
+
+scaledAnticommutatorIO ::
+   (Class.Floating a) =>
+   a ->
+   Order -> ForeignPtr a -> ForeignPtr a -> Ptr a ->
+   ((Int, Int), (Char, Char, Int)) -> IO ()
+scaledAnticommutatorIO alpha order a b cpPtr ((n,k), (uplo,trans,lda)) =
+   evalContT $ do
+      uploPtr <- Call.char uplo
+      transPtr <- Call.char trans
+      nPtr <- Call.cint n
+      kPtr <- Call.cint k
+      alphaPtr <- Call.number alpha
+      aPtr <- ContT $ withForeignPtr a
+      ldaPtr <- Call.leadingDim lda
+      bPtr <- ContT $ withForeignPtr b
+      let ldbPtr = ldaPtr
+      betaPtr <- Call.number zero
+      cPtr <- Call.allocaArray (n*n)
+      ldcPtr <- Call.leadingDim n
+      liftIO $ do
+         BlasGen.syr2k uploPtr transPtr nPtr kPtr alphaPtr
+            aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr
+         pack order n cPtr cpPtr
+
+
+symmetricShape :: Order -> size -> MatrixShape.Symmetric size
+symmetricShape = MatrixShape.Triangular NonUnit MatrixShape.autoUplo
diff --git a/src/Numeric/LAPACK/Matrix/Triangular/Basic.hs b/src/Numeric/LAPACK/Matrix/Triangular/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Triangular/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Triangular/Basic.hs
@@ -52,7 +52,7 @@
 import qualified Numeric.LAPACK.Vector as Vector
 import Numeric.LAPACK.Matrix.Triangular.Private
          (Triangular, FlexDiagonal, diagonalPointers, diagonalPointerPairs,
-          pack, packRect, unpack, unpackZero, unpackToTemp)
+          pack, packRect, unpack, unpackZero, unpackToTemp, uncheck, recheck)
 import Numeric.LAPACK.Matrix.Shape.Private
          (Order(RowMajor,ColumnMajor),
           flipOrder, transposeFromOrder, uploFromOrder, uploOrder,
@@ -587,18 +587,6 @@
    Array.unsafeCreate shape $
       Symmetric.square NonConjugated order (Shape.size sh) a
 
-
-uncheck :: Triangular lo diag up sh a -> Triangular lo diag up (Unchecked sh) a
-uncheck =
-   Array.mapShape $
-      \(MatrixShape.Triangular diag uplo order sh) ->
-         MatrixShape.Triangular diag uplo order (Unchecked sh)
-
-recheck :: Triangular lo diag up (Unchecked sh) a -> Triangular lo diag up sh a
-recheck =
-   Array.mapShape $
-      \(MatrixShape.Triangular diag uplo order (Unchecked sh)) ->
-         MatrixShape.Triangular diag uplo order sh
 
 {-
 Requires frequent unpacking and packing of triangles.
diff --git a/src/Numeric/LAPACK/Matrix/Triangular/Private.hs b/src/Numeric/LAPACK/Matrix/Triangular/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Triangular/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Triangular/Private.hs
@@ -13,6 +13,7 @@
 import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated))
 import Numeric.LAPACK.Matrix.Private (Full)
 import Numeric.LAPACK.Scalar (zero)
+import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))
 import Numeric.LAPACK.Private
          (pointerSeq, copyBlock, copyCondConjugateToTemp,
           pokeCInt, fill, withInfo, errorCodeMsg)
@@ -138,6 +139,19 @@
    zPtr <- Call.number z
    liftIO $ LapackGen.laset uploPtr nPtr nPtr zPtr zPtr aPtr nPtr
 
+
+
+uncheck :: Triangular lo diag up sh a -> Triangular lo diag up (Unchecked sh) a
+uncheck =
+   Array.mapShape $
+      \(MatrixShape.Triangular diag uplo order sh) ->
+         MatrixShape.Triangular diag uplo order (Unchecked sh)
+
+recheck :: Triangular lo diag up (Unchecked sh) a -> Triangular lo diag up sh a
+recheck =
+   Array.mapShape $
+      \(MatrixShape.Triangular diag uplo order (Unchecked sh)) ->
+         MatrixShape.Triangular diag uplo order sh
 
 
 stack ::
diff --git a/src/Numeric/LAPACK/Split.hs b/src/Numeric/LAPACK/Split.hs
--- a/src/Numeric/LAPACK/Split.hs
+++ b/src/Numeric/LAPACK/Split.hs
@@ -1,12 +1,14 @@
+{-# LANGUAGE TypeFamilies #-}
 module Numeric.LAPACK.Split where
 
 import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Shape.Box as Box
 import qualified Numeric.LAPACK.Matrix.Triangular.Private as TriPriv
 import qualified Numeric.LAPACK.Matrix.Triangular.Basic as Tri
 import qualified Numeric.LAPACK.Matrix.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import qualified Numeric.LAPACK.Private as Private
-import Numeric.LAPACK.Matrix.Triangular.Private (diagonalPointers)
+import Numeric.LAPACK.Matrix.Triangular.Private (diagonalPointers, unpack)
 import Numeric.LAPACK.Matrix.Triangular.Basic (UnitLower, Upper)
 import Numeric.LAPACK.Matrix.Shape.Private
          (Order(RowMajor, ColumnMajor), transposeFromOrder,
@@ -238,3 +240,34 @@
          withInfo "trtrs" $
             LapackGen.trtrs uploPtr transPtr diagPtr
                nPtr nrhsPtr aPtr ldaPtr xPtr ldxPtr
+
+
+data Corrupt = Corrupt
+   deriving (Eq)
+
+{-
+We could use Plain.Class.ShapeOrder
+but this would currently cause an import cycle.
+-}
+{- |
+> let b = takeHalf a
+> ==>
+> isTriangular b && a == addTransposed b
+-}
+takeHalf ::
+   (Box.Box symShape, Box.HeightOf symShape ~ sh, Shape.C sh,
+    Class.Floating a) =>
+   (symShape -> Order) -> Array symShape a -> Square Corrupt sh a
+takeHalf shapeOrder (Array symShape a) =
+   let sh = Box.height symShape
+       order = shapeOrder symShape
+   in Array.unsafeCreate (MatrixShape.Split Corrupt order (Extent.square sh)) $
+         \bPtr -> evalContT $ do
+      let n = Shape.size sh
+      aPtr <- ContT $ withForeignPtr a
+      nPtr <- Call.cint n
+      alphaPtr <- Call.number 0.5
+      incxPtr <- Call.cint (n+1)
+      liftIO $ do
+         unpack order n aPtr bPtr
+         BlasGen.scal nPtr alphaPtr bPtr incxPtr
diff --git a/test/Main.hs b/test/Main.hs
--- a/test/Main.hs
+++ b/test/Main.hs
@@ -5,6 +5,7 @@
 import qualified Test.Matrix as Matrix
 import qualified Test.Square as Square
 import qualified Test.Triangular as Triangular
+import qualified Test.Symmetric as Symmetric
 import qualified Test.Hermitian as Hermitian
 import qualified Test.Banded as Banded
 import qualified Test.BandedHermitian as BandedHermitian
@@ -39,6 +40,7 @@
    prefix "Matrix" Matrix.testsVar ++
    prefix "Square" Square.testsVar ++
    prefix "Triangular" Triangular.testsVar ++
+   prefix "Symmetric" Symmetric.testsVar ++
    prefix "Hermitian" Hermitian.testsVar ++
    prefix "Banded" Banded.testsVar ++
    prefix "BandedHermitian" BandedHermitian.testsVar ++
diff --git a/test/Test/Symmetric.hs b/test/Test/Symmetric.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Symmetric.hs
@@ -0,0 +1,205 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Test.Symmetric (testsVar) where
+
+import qualified Test.Generator as Gen
+import qualified Test.Utility as Util
+import Test.Generator ((<-*#>), (<#*|>), (<#*#>), (<#=#>))
+import Test.Utility
+         (approxArray, approxMatrix, equalArray, Tagged, genOrder, (!===))
+
+import qualified Numeric.LAPACK.Matrix.Symmetric as Symmetric
+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
+import qualified Numeric.LAPACK.Matrix as Matrix
+import qualified Numeric.LAPACK.Vector as Vector
+import Numeric.LAPACK.Matrix.Symmetric (Symmetric)
+import Numeric.LAPACK.Matrix.Shape (Order)
+import Numeric.LAPACK.Matrix (General, ShapeInt, (#+#), (|||))
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Scalar (RealOf)
+
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Shape ((:+:))
+
+import Control.Applicative ((<$>))
+
+import Data.Semigroup ((<>))
+import Data.Tuple.HT (uncurry3)
+
+import qualified Test.QuickCheck as QC
+
+
+
+generalFromSymmetric ::
+   (Shape.C sh, Class.Floating a) => Symmetric sh a -> General sh sh a
+generalFromSymmetric = Matrix.fromFull . Triangular.toSquare
+
+stack ::
+   (Class.Floating a) =>
+   (Symmetric ShapeInt a, General ShapeInt ShapeInt a, Symmetric ShapeInt a) ->
+   Bool
+stack (a,b,c) =
+   let abc = generalFromSymmetric $ Symmetric.stack a b c
+   in equalArray abc $
+         (Matrix.fromFull (Triangular.toSquare a) ||| b
+          !===
+          Matrix.transpose b ||| Matrix.fromFull (Triangular.toSquare c))
+
+split :: (Class.Floating a) => Symmetric (ShapeInt:+:ShapeInt) a -> Bool
+split abc = equalArray abc $ uncurry3 Symmetric.stack $ Symmetric.split abc
+
+
+gramian ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   General ShapeInt ShapeInt a -> Bool
+gramian x =
+   approxArray
+      (generalFromSymmetric $ Symmetric.gramian x)
+      (Matrix.transpose x <> x)
+
+gramianTransposed ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   General ShapeInt ShapeInt a -> Bool
+gramianTransposed x =
+   approxArray
+      (generalFromSymmetric $ Symmetric.gramianTransposed x)
+      (Matrix.adaptOrder x $ x <> Matrix.transpose x)
+
+gramianNonTransposed ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   General ShapeInt ShapeInt a -> Bool
+gramianNonTransposed x =
+   approxArray
+      (Matrix.forceOrder (ArrMatrix.shapeOrder $ ArrMatrix.shape x) $
+       Symmetric.gramian $ Matrix.transpose x)
+      (Symmetric.gramianTransposed x)
+
+congruenceDiagonal ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Vector ShapeInt a, General ShapeInt ShapeInt a) -> Bool
+congruenceDiagonal (d,a) =
+   approxArray
+      (generalFromSymmetric $ Symmetric.congruenceDiagonal d a)
+      (Matrix.transpose a <> Matrix.scaleRows d a)
+
+congruenceDiagonalTransposed ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (General ShapeInt ShapeInt a, Vector ShapeInt a) -> Bool
+congruenceDiagonalTransposed (a,d) =
+   approxMatrix 1e-5
+      (generalFromSymmetric $ Symmetric.congruenceDiagonalTransposed a d)
+      (Matrix.scaleColumns d a <> Matrix.transpose a)
+
+congruenceDiagonalGramian ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   General ShapeInt ShapeInt a -> Bool
+congruenceDiagonalGramian a =
+   approxArray
+      (Symmetric.congruenceDiagonal (Vector.one $ Matrix.height a) a)
+      (Symmetric.gramian a)
+
+congruence ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Symmetric ShapeInt a, General ShapeInt ShapeInt a) -> Bool
+congruence (b,a) =
+   approxArray
+      (generalFromSymmetric $ Symmetric.congruence b a)
+      (Matrix.transpose a <> generalFromSymmetric b <> a)
+
+congruenceTransposed ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (General ShapeInt ShapeInt a, Symmetric ShapeInt a) -> Bool
+congruenceTransposed (a,b) =
+   approxMatrix 1e-5
+      (generalFromSymmetric $ Symmetric.congruenceTransposed a b)
+      (a <> generalFromSymmetric b <> Matrix.transpose a)
+
+congruenceCongruenceDiagonal ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   Order -> (Vector ShapeInt a, General ShapeInt ShapeInt a) -> Bool
+congruenceCongruenceDiagonal order (d,a) =
+   approxArray
+      (Symmetric.congruenceDiagonal d a)
+      (Symmetric.congruence (Triangular.diagonal order d) a)
+
+anticommutator ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (General ShapeInt ShapeInt a, General ShapeInt ShapeInt a) -> Bool
+anticommutator (a,b) =
+   approxArray
+      (generalFromSymmetric $ Symmetric.anticommutator a b)
+      ((Matrix.transpose b <> a) #+# (Matrix.transpose a <> b))
+
+anticommutatorCommutative ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (General ShapeInt ShapeInt a, General ShapeInt ShapeInt a) -> Bool
+anticommutatorCommutative (a,b) =
+   approxMatrix 1e-5
+      (Symmetric.anticommutator a b)
+      (Symmetric.anticommutator b a)
+
+anticommutatorTransposed ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (General ShapeInt ShapeInt a, General ShapeInt ShapeInt a) -> Bool
+anticommutatorTransposed (a,b) =
+   approxArray
+      (Matrix.forceOrder (ArrMatrix.shapeOrder $ ArrMatrix.shape b) $
+       Symmetric.anticommutator (Matrix.transpose a) (Matrix.transpose b))
+      (Symmetric.anticommutatorTransposed a b)
+
+
+
+checkForAll ::
+   (Show a, QC.Testable test) =>
+   Gen.T dim tag a -> (a -> test) -> Tagged tag QC.Property
+checkForAll gen = Util.checkForAll (Gen.run gen 3 5)
+
+checkForAllExtra ::
+   (Show a, Show b, QC.Testable test) =>
+   QC.Gen a -> Gen.T dim tag b ->
+   (a -> b -> test) -> Tagged tag QC.Property
+checkForAllExtra = Gen.withExtra checkForAll
+
+
+testsVar ::
+   (Show a, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>
+   [(String, Tagged a QC.Property)]
+testsVar =
+   ("stack",
+      checkForAll (Gen.stack3 Gen.triangular Gen.matrix Gen.triangular) stack) :
+   ("split",
+      checkForAll Gen.triangular split) :
+
+   ("gramian",
+      checkForAll Gen.matrix gramian) :
+   ("gramianTransposed",
+      checkForAll Gen.matrix gramianTransposed) :
+   ("gramianNonTransposed",
+      checkForAll Gen.matrix gramianNonTransposed) :
+   ("congruenceDiagonal",
+      checkForAll ((,) <$> Gen.vector <-*#> Gen.matrix) congruenceDiagonal) :
+   ("congruence",
+      checkForAll ((,) <$> Gen.triangular <#*#> Gen.matrix) congruence) :
+   ("congruenceDiagonalTransposed",
+      checkForAll ((,) <$> Gen.matrix <#*|> Gen.vector)
+         congruenceDiagonalTransposed) :
+   ("congruenceDiagonalGramian",
+      checkForAll Gen.matrix congruenceDiagonalGramian) :
+   ("congruenceTransposed",
+      checkForAll ((,) <$> Gen.matrix <#*#> Gen.triangular)
+         congruenceTransposed) :
+   ("congruenceCongruenceDiagonal",
+      checkForAllExtra genOrder
+         ((,) <$> Gen.vector <-*#> Gen.matrix) congruenceCongruenceDiagonal) :
+   ("anticommutator",
+      checkForAll ((,) <$> Gen.matrix <#=#> Gen.matrix) anticommutator) :
+   ("anticommutatorCommutative",
+      checkForAll ((,) <$> Gen.matrix <#=#> Gen.matrix)
+         anticommutatorCommutative) :
+   ("anticommutatorTransposed",
+      checkForAll ((,) <$> Gen.matrix <#=#> Gen.matrix)
+         anticommutatorTransposed) :
+   []
diff --git a/test/Test/Triangular.hs b/test/Test/Triangular.hs
--- a/test/Test/Triangular.hs
+++ b/test/Test/Triangular.hs
@@ -116,8 +116,7 @@
     Triangular.FlexUpper diag ShapeInt a) ->
    Bool
 stackUpper (a,b,c) =
-   let abc = Matrix.fromFull $
-         Triangular.toSquare $ Triangular.stackUpper a b c
+   let abc = expandTriangle $ Triangular.stackUpper a b c
    in equalArray abc $
          (expandTriangle a ||| b
           !===