12 files changed, 217 insertions, 64 deletions

diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG
index 93b2594..c3e118c 100644
--- a/packages/base/CHANGELOG
+++ b/packages/base/CHANGELOG
@@ -5,7 +5,7 @@
 
     * basic support of Int32 and Int64 elements
     
-    * remap, ccompare, sortIndex
+    * remap, more general cond
 
     * experimental support of type safe modular arithmetic, including linear
       systems and lu factorization
diff --git a/packages/base/src/Internal/CG.hs b/packages/base/src/Internal/CG.hs
index fd14212..758d130 100644
--- a/packages/base/src/Internal/CG.hs
+++ b/packages/base/src/Internal/CG.hs
@@ -45,13 +45,13 @@ cg sym at a (CGState p r r2 x _) = CGState p' r' r'2 x' rdx
     ap1 = a p
     ap  | sym       = ap1
         | otherwise = at ap1
-    pap | sym       = p <·> ap1
+    pap | sym       = p <.> ap1
         | otherwise = norm2 ap1 ** 2
     alpha = r2 / pap
     dx = scale alpha p
     x' = x + dx
     r' = r - scale alpha ap
-    r'2 = r' <·> r'
+    r'2 = r' <.> r'
     beta = r'2 / r2
     p' = r' + scale beta p
 
@@ -75,9 +75,9 @@ solveG mat ma meth rawb x0' ϵb ϵx
     b = mat rawb
     x0  = if x0' == 0 then konst 0 (dim b) else x0'
     r0  = b - a x0
-    r20 = r0 <·> r0
+    r20 = r0 <.> r0
     p0  = r0
-    nb2 = b <·> b
+    nb2 = b <.> b
     ok CGState {..}
         =  cgr2 <nb2*ϵb**2
         || cgdx < ϵx
diff --git a/packages/base/src/Internal/Container.hs b/packages/base/src/Internal/Container.hs
index 7fe5758..1c158ff 100644
--- a/packages/base/src/Internal/Container.hs
+++ b/packages/base/src/Internal/Container.hs
@@ -53,28 +53,23 @@ linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n
 
 --------------------------------------------------------------------------------
 
-infixl 7 <.>
--- | An infix synonym for 'dot'
-(<.>) :: Numeric t => Vector t -> Vector t -> t
-(<.>) = dot
+infixr 8 <.>
+{- | An infix synonym for 'dot'
 
+>>> vector [1,2,3,4] <.> vector [-2,0,1,1]
+5.0
 
-infixr 8 <·>, #>
+>>> let 𝑖 = 0:+1 :: C
+>>> fromList [1+𝑖,1] <.> fromList [1,1+𝑖]
+2.0 :+ 0.0
 
-{- | infix synonym for 'dot'
+-}
 
->>> vector [1,2,3,4] <·> vector [-2,0,1,1]
-5.0
+(<.>) :: Numeric t => Vector t -> Vector t -> t
+(<.>) = dot
 
->>> let 𝑖 = 0:+1 :: ℂ
->>> fromList [1+𝑖,1] <·> fromList [1,1+𝑖]
-2.0 :+ 0.0
 
-(the dot symbol "·" is obtained by Alt-Gr .)
 
--}
-(<·>) :: Numeric t => Vector t -> Vector t -> t
-(<·>) = dot
 
 
 {- | infix synonym for 'app'
@@ -91,6 +86,7 @@ infixr 8 <·>, #>
 fromList [140.0,320.0]
 
 -}
+infixr 8 #>
 (#>) :: Numeric t => Matrix t -> Vector t -> Vector t
 (#>) = mXv
 
@@ -264,6 +260,24 @@ instance Numeric Z
 sortVector :: (Ord t, Element t) => Vector t -> Vector t
 sortVector = sortV
 
+{- |
+
+>>> m <- randn 4 10
+>>> disp 2 m
+4x10
+-0.31   0.41   0.43  -0.19  -0.17  -0.23  -0.17  -1.04  -0.07  -1.24
+ 0.26   0.19   0.14   0.83  -1.54  -0.09   0.37  -0.63   0.71  -0.50
+-0.11  -0.10  -1.29  -1.40  -1.04  -0.89  -0.68   0.35  -1.46   1.86
+ 1.04  -0.29   0.19  -0.75  -2.20  -0.01   1.06   0.11  -2.09  -1.58
+
+>>> disp 2 $ m ?? (All, Pos $ sortIndex (m!1))
+4x10
+-0.17  -1.04  -1.24  -0.23   0.43   0.41  -0.31  -0.17  -0.07  -0.19
+-1.54  -0.63  -0.50  -0.09   0.14   0.19   0.26   0.37   0.71   0.83
+-1.04   0.35   1.86  -0.89  -1.29  -0.10  -0.11  -0.68  -1.46  -1.40
+-2.20   0.11  -1.58  -0.01   0.19  -0.29   1.04   1.06  -2.09  -0.75
+
+-}
 sortIndex :: (Ord t, Element t) => Vector t -> Vector I
 sortIndex = sortI
 
@@ -273,11 +287,50 @@ ccompare = ccompare'
 cselect :: (Container c t) => c I -> c t -> c t -> c t -> c t
 cselect = cselect'
 
+{- | Extract elements from positions given in matrices of rows and columns.
+
+>>> r
+(3><3)
+ [ 1, 1, 1
+ , 1, 2, 2
+ , 1, 2, 3 ]
+>>> c
+(3><3)
+ [ 0, 1, 5
+ , 2, 2, 1
+ , 4, 4, 1 ]
+>>> m
+(4><6)
+ [  0,  1,  2,  3,  4,  5
+ ,  6,  7,  8,  9, 10, 11
+ , 12, 13, 14, 15, 16, 17
+ , 18, 19, 20, 21, 22, 23 ]
+
+>>> remap r c m
+(3><3)
+ [  6,  7, 11
+ ,  8, 14, 13
+ , 10, 16, 19 ]
+
+The indexes are autoconformable.
+
+>>> c'
+(3><1)
+ [ 1
+ , 2
+ , 4 ]
+>>> remap r c' m
+(3><3)
+ [  7,  7,  7
+ ,  8, 14, 14
+ , 10, 16, 22 ]
+
+-}
 remap :: Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t
 remap i j m
     | minElement i >= 0 && maxElement i < fromIntegral (rows m) &&
       minElement j >= 0 && maxElement j < fromIntegral (cols m) = remapM i' j' m
-    | otherwise = error $ "out of range index in rmap"
+    | otherwise = error $ "out of range index in remap"
   where
     [i',j'] = conformMs [i,j]
     
diff --git a/packages/base/src/Internal/Element.hs b/packages/base/src/Internal/Element.hs
index 4007491..51d5686 100644
--- a/packages/base/src/Internal/Element.hs
+++ b/packages/base/src/Internal/Element.hs
@@ -80,7 +80,7 @@ breakAt c l = (a++[c],tail b) where
     (a,b) = break (==c) l
 
 --------------------------------------------------------------------------------
-
+-- | Specification of indexes for the operator '??'.
 data Extractor
     = All
     | Range Int Int Int
@@ -102,7 +102,32 @@ ppext (Drop n) = printf "Drop %d" n
 ppext (TakeLast n) = printf "TakeLast %d" n
 ppext (DropLast n) = printf "DropLast %d" n
 
+{- | General matrix slicing.
+
+>>> m
+(4><5)
+ [  0,  1,  2,  3,  4
+ ,  5,  6,  7,  8,  9
+ , 10, 11, 12, 13, 14
+ , 15, 16, 17, 18, 19 ]
+
+>>> m ?? (Take 3, DropLast 2)
+(3><3)
+ [  0,  1,  2
+ ,  5,  6,  7
+ , 10, 11, 12 ]
+
+>>> m ?? (Pos (idxs[2,1]), All)
+(2><5)
+ [ 10, 11, 12, 13, 14
+ ,  5,  6,  7,  8,  9 ]
+
+>>> m ?? (PosCyc (idxs[-7,80]), Range 4 (-2) 0)
+(2><3)
+ [ 9, 7, 5
+ , 4, 2, 0 ]
 
+-}
 infixl 9 ??
 (??)  :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t
 
@@ -328,27 +353,30 @@ r >< c = f where
 
 ----------------------------------------------------------------
 
--- | Creates a matrix with the first n rows of another matrix
 takeRows :: Element t => Int -> Matrix t -> Matrix t
 takeRows n mt = subMatrix (0,0) (n, cols mt) mt
+
 -- | Creates a matrix with the last n rows of another matrix
 takeLastRows :: Element t => Int -> Matrix t -> Matrix t
 takeLastRows n mt = subMatrix (rows mt - n, 0) (n, cols mt) mt
--- | Creates a copy of a matrix without the first n rows
+
 dropRows :: Element t => Int -> Matrix t -> Matrix t
 dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt
+
 -- | Creates a copy of a matrix without the last n rows
 dropLastRows :: Element t => Int -> Matrix t -> Matrix t
 dropLastRows n mt = subMatrix (0,0) (rows mt - n, cols mt) mt
--- |Creates a matrix with the first n columns of another matrix
+
 takeColumns :: Element t => Int -> Matrix t -> Matrix t
 takeColumns n mt = subMatrix (0,0) (rows mt, n) mt
+
 -- |Creates a matrix with the last n columns of another matrix
 takeLastColumns :: Element t => Int -> Matrix t -> Matrix t
 takeLastColumns n mt = subMatrix (0, cols mt - n) (rows mt, n) mt
--- | Creates a copy of a matrix without the first n columns
+
 dropColumns :: Element t => Int -> Matrix t -> Matrix t
 dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt
+
 -- | Creates a copy of a matrix without the last n columns
 dropLastColumns :: Element t => Int -> Matrix t -> Matrix t
 dropLastColumns n mt = subMatrix (0,0) (rows mt, cols mt - n) mt
diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs
index e4b1226..75e92a5 100644
--- a/packages/base/src/Internal/Matrix.hs
+++ b/packages/base/src/Internal/Matrix.hs
@@ -378,7 +378,6 @@ foreign import ccall unsafe "transL" ctransL :: TMM Z
 
 ----------------------------------------------------------------------
 
--- | Extracts a submatrix from a matrix.
 subMatrix :: Element a
           => (Int,Int) -- ^ (r0,c0) starting position
           -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs
index 3b27310..0d906bb 100644
--- a/packages/base/src/Internal/Modular.hs
+++ b/packages/base/src/Internal/Modular.hs
@@ -9,8 +9,8 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE GADTs #-}
-{-# LANGUAGE TypeFamilies #-}
-
+{-# LANGUAGE TypeFamilies  #-}
+{-# LANGUAGE TypeOperators #-}
 
 {- |
 Module      :  Internal.Modular
@@ -23,7 +23,7 @@ Proof of concept of statically checked modular arithmetic.
 -}
 
 module Internal.Modular(
-    Mod
+    Mod, type (./.)
 ) where
 
 import Internal.Vector
@@ -49,6 +49,9 @@ import Data.Complex
 newtype Mod (n :: Nat) t = Mod {unMod:: t}
   deriving (Storable)
 
+infixr 5 ./.
+type (./.) x n = Mod n x
+
 instance (Integral t, Enum t, KnownNat m) => Enum (Mod m t)
   where
     toEnum = l0 (\m x -> fromIntegral $ x `mod` (fromIntegral m))
diff --git a/packages/base/src/Internal/Numeric.hs b/packages/base/src/Internal/Numeric.hs
index ca17c23..efcde2c 100644
--- a/packages/base/src/Internal/Numeric.hs
+++ b/packages/base/src/Internal/Numeric.hs
@@ -481,14 +481,16 @@ step = step'
 -- ,   0.0, 100.0,   7.0,  8.0
 -- ,   0.0,   0.0, 100.0, 12.0 ]
 --
+-- >>> let chop x = cond (abs x) 1E-6 0 0 x
+--
 cond
-    :: (Ord e, Container c e)
+    :: (Ord e, Container c e, Container c x)
     => c e -- ^ a
     -> c e -- ^ b
-    -> c e -- ^ l
-    -> c e -- ^ e
-    -> c e -- ^ g
-    -> c e -- ^ result
+    -> c x -- ^ l
+    -> c x -- ^ e
+    -> c x -- ^ g
+    -> c x -- ^ result
 cond a b l e g = cselect' (ccompare' a b) l e g
 
 
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index 09ba21c..f08f710 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -85,7 +85,7 @@ type ℤ = Int
 type ℂ = Complex Double
 
 -- | imaginary unit
-iC :: ℂ
+iC :: C
 iC = 0:+1
 
 {- | Create a real vector.
@@ -94,7 +94,7 @@ iC = 0:+1
 fromList [1.0,2.0,3.0,4.0,5.0]
 
 -}
-vector :: [ℝ] -> Vector ℝ
+vector :: [R] -> Vector R
 vector = fromList
 
 {- | Create a real matrix.
@@ -108,8 +108,8 @@ vector = fromList
 -}
 matrix
   :: Int -- ^ number of columns
-  -> [ℝ] -- ^ elements in row order
-  -> Matrix ℝ
+  -> [R] -- ^ elements in row order
+  -> Matrix R
 matrix c = reshape c . fromList
 
 
@@ -260,34 +260,34 @@ norm = pnorm PNorm2
 
 class Normed a
   where
-    norm_0   :: a -> ℝ
-    norm_1   :: a -> ℝ
-    norm_2   :: a -> ℝ
-    norm_Inf :: a -> ℝ
+    norm_0   :: a -> R
+    norm_1   :: a -> R
+    norm_2   :: a -> R
+    norm_Inf :: a -> R
 
 
-instance Normed (Vector ℝ)
+instance Normed (Vector R)
   where
     norm_0 v = sumElements (step (abs v - scalar (eps*normInf v)))
     norm_1 = pnorm PNorm1
     norm_2 = pnorm PNorm2
     norm_Inf = pnorm Infinity
 
-instance Normed (Vector ℂ)
+instance Normed (Vector C)
   where
     norm_0 v = sumElements (step (fst (fromComplex (abs v)) - scalar (eps*normInf v)))
     norm_1 = pnorm PNorm1
     norm_2 = pnorm PNorm2
     norm_Inf = pnorm Infinity
 
-instance Normed (Matrix ℝ)
+instance Normed (Matrix R)
   where
     norm_0 = norm_0 . flatten
     norm_1 = pnorm PNorm1
     norm_2 = pnorm PNorm2
     norm_Inf = pnorm Infinity
 
-instance Normed (Matrix ℂ)
+instance Normed (Matrix C)
   where
     norm_0 = norm_0 . flatten
     norm_1 = pnorm PNorm1
@@ -323,12 +323,22 @@ instance Normed (Vector (Complex Float))
     norm_Inf = norm_Inf . double
 
 
-norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> ℝ
+norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R
 norm_Frob = norm_2 . flatten
 
-norm_nuclear :: Field t => Matrix t -> ℝ
+norm_nuclear :: Field t => Matrix t -> R
 norm_nuclear = sumElements . singularValues
 
+{- | Check if the absolute value or complex magnitude is greater than a given threshold
+
+>>> magnit 1E-6 (1E-12 :: R)
+False
+>>> magnit 1E-6 (3+iC :: C)
+True
+>>> magnit 0 (3 :: I ./. 5)
+True
+
+-}
 magnit :: (Element t, Normed (Vector t)) => R -> t -> Bool
 magnit e x = norm_1 (fromList [x]) > e
 
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index fe524cc..2e6b8ca 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -47,7 +47,7 @@ module Numeric.LinearAlgebra (
 
     -- * Products
     -- ** dot
-    dot, (<·>),
+    dot, (<.>),
     -- ** matrix-vector
     app, (#>), (<#), (!#>),
     -- ** matrix-matrix
@@ -240,7 +240,53 @@ nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m)
 -- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'.
 orth m = orthSVD (Left (1*eps)) m (leftSV m)
 
+{- | Experimental implementation of LU factorization
+     working on any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.
+
+>>> let m = ident 5 + (5><5) [0..] :: Matrix (Z ./. 17)
+(5><5)
+ [  1,  1,  2,  3,  4
+ ,  5,  7,  7,  8,  9
+ , 10, 11, 13, 13, 14
+ , 15, 16,  0,  2,  2
+ ,  3,  4,  5,  6,  8 ]
+
+>>> let (l,u,p,s) = luFact $ luPacked' m
+>>> l
+(5><5)
+ [  1,  0, 0,  0, 0
+ ,  6,  1, 0,  0, 0
+ , 12,  7, 1,  0, 0
+ ,  7, 10, 7,  1, 0
+ ,  8,  2, 6, 11, 1 ]
+>>> u
+(5><5)
+ [ 15, 16,  0,  2,  2
+ ,  0, 13,  7, 13, 14
+ ,  0,  0, 15,  0, 11
+ ,  0,  0,  0, 15, 15
+ ,  0,  0,  0,  0,  1 ]
+
+-}
 luPacked' x = mutable (luST (magnit 0)) x
 
+{- | Experimental implementation of gaussian elimination
+     working on any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.
+
+>>> let a = (2><2) [1,2,3,5] :: Matrix (Z ./. 13)
+(2><2)
+ [ 1, 2
+ , 3, 5 ]
+>>> b
+(2><3)
+ [ 5, 1, 3
+ , 8, 6, 3 ]
+
+>>> let x = linearSolve' a b
+(2><3)
+ [ 4,  7, 4
+ , 7, 10, 6 ]
+
+-}
 linearSolve' x y = gaussElim x y
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index fffc2bd..2a45771 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE TypeOperators #-}
+
 --------------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Data
@@ -17,11 +19,11 @@ module Numeric.LinearAlgebra.Data(
     -- | 1D arrays are storable vectors from the vector package. There is no distinction
     --   between row and column vectors.
 
-    fromList, toList, vector, (|>),
+    fromList, toList, (|>), vector, range, idxs,
 
     -- * Matrix
 
-    matrix, (><), tr, tr',
+    (><), matrix, tr, tr',
 
     -- * Dimensions
 
@@ -43,7 +45,7 @@ module Numeric.LinearAlgebra.Data(
     Indexable(..),
 
     -- * Construction
-    scalar, Konst(..), Build(..), assoc, accum, linspace, range, idxs, -- ones, zeros,
+    scalar, Konst(..), Build(..), assoc, accum, linspace,  -- ones, zeros,
 
     -- * Diagonal
     ident, diag, diagl, diagRect, takeDiag,
@@ -53,16 +55,19 @@ module Numeric.LinearAlgebra.Data(
 
     -- * Matrix extraction
     Extractor(..), (??),
+
     takeRows, dropRows, takeColumns, dropColumns,
-    subMatrix, (?), (¿), fliprl, flipud, remap,
+    subMatrix, (?), (¿), fliprl, flipud,
+
+    remap,
 
     -- * Block matrix
     fromBlocks, (|||), (===), diagBlock, repmat, toBlocks, toBlocksEvery,
 
     -- * Mapping functions
     conj, cmap, cmod,
-    
-    step, cond, ccompare, cselect,
+
+    step, cond,
 
     -- * Find elements
     find, maxIndex, minIndex, maxElement, minElement,
@@ -87,7 +92,7 @@ module Numeric.LinearAlgebra.Data(
     separable,
     fromArray2D,
     module Data.Complex,
-    R,C,I,Z,Mod,
+    R,C,I,Z,Mod, type(./.),
     Vector, Matrix, GMatrix, nRows, nCols
 
 ) where
diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
index 11c2487..8adaaaf 100644
--- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
@@ -13,10 +13,13 @@ compatibility with previous version, to be removed
 
 module Numeric.LinearAlgebra.HMatrix (
     module Numeric.LinearAlgebra,
-    (¦),(——),ℝ,ℂ,
+    (¦),(——),ℝ,ℂ,(<·>)
 ) where
 
 import Numeric.LinearAlgebra
 import Internal.Util
 
+infixr 8 <·>
+(<·>) :: Numeric t => Vector t -> Vector t -> t
+(<·>) = dot
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
index a657bd0..d0a790d 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs
@@ -44,7 +44,7 @@ module Numeric.LinearAlgebra.Static(
     -- * Complex
     C, M, Her, her, 𝑖,
     -- * Products
-    (<>),(#>),(<·>),
+    (<>),(#>),(<.>),
     -- * Linear Systems
     linSolve, (<\>),
     -- * Factorizations
@@ -55,13 +55,13 @@ module Numeric.LinearAlgebra.Static(
     Disp(..), Domain(..),
     withVector, withMatrix,
     toRows, toColumns,
-    Sized(..), Diag(..), Sym, sym, mTm, unSym
+    Sized(..), Diag(..), Sym, sym, mTm, unSym, (<·>)
 ) where
 
 
 import GHC.TypeLits
 import Numeric.LinearAlgebra hiding (
-    (<>),(#>),(<·>),Konst(..),diag, disp,(===),(|||),
+    (<>),(#>),(<.>),Konst(..),diag, disp,(===),(|||),
     row,col,vector,matrix,linspace,toRows,toColumns,
     (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',
     eigenvalues,eigenvaluesSH,eigenvaluesSH',build,
@@ -207,6 +207,10 @@ infixr 8 <·>
 (<·>) :: R n -> R n -> ℝ
 (<·>) = dotR
 
+infixr 8 <.>
+(<.>) :: R n -> R n -> ℝ
+(<.>) = dotR
+
 --------------------------------------------------------------------------------
 
 class Diag m d | m -> d
@@ -496,7 +500,7 @@ appC m v = mkC (extract m LA.#> extract v)
 dotC :: KnownNat n => C n -> C n -> ℂ
 dotC (unwrap -> u) (unwrap -> v)
     | singleV u || singleV v = sumElements (conj u * v)
-    | otherwise = u LA.<·> v
+    | otherwise = u LA.<.> v
 
 
 crossC :: C 3 -> C 3 -> C 3
@@ -584,12 +588,12 @@ test = (ok,info)
       where
         q = tm :: L 10 3
 
-    thingD = vjoin [ud1 u, 1] LA.<·> tr m LA.#> m LA.#> ud1 v
+    thingD = vjoin [ud1 u, 1] LA.<.> tr m LA.#> m LA.#> ud1 v
       where
         m = LA.matrix 3 [1..30]
 
     precS = (1::Double) + (2::Double) * ((1 :: R 3) * (u & 6)) <·> konst 2 #> v
-    precD = 1 + 2 * vjoin[ud1 u, 6] LA.<·> LA.konst 2 (LA.size (ud1 u) +1, LA.size (ud1 v)) LA.#> ud1 v
+    precD = 1 + 2 * vjoin[ud1 u, 6] LA.<.> LA.konst 2 (LA.size (ud1 u) +1, LA.size (ud1 v)) LA.#> ud1 v
 
 
 splittest