Cannot get the <\> operator to work

[denjoh]

I have just started to learn Haskell and I am trying to build a simple tool for analysing a beam on multiple spring supports. I have been able to get the lusolve to work but not the <>.

I keep on getting this error

*** Exception: givens: no compatible rows for entry (3,2) CallStack (from HasCallStack): error, called at src\Numeric\LinearAlgebra\Sparse.hs:168:28 in sparse-linear-algebra-0.2.2.0-CbAalqT73THJ7zUSs35jqX:Numeric.LinearAlgebra.Sparse

My code

import System.IO
import Numeric.LinearAlgebra.Sparse
import Numeric.LinearAlgebra.Class
import Data.Sparse.SpMatrix
import Data.Sparse.Common

eMod = 210000.0

k 0 0 i l = 12* eMod * i /(l**3)
k 0 1 i l = 6 * eMod * i /(l**2)
k 0 2 i l = (- (k 0 0 i l))
k 0 3 i l = k 0 1 i l 
k 1 0 i l = k 0 1 i l 
k 1 1 i l = 4 * eMod *i / l
k 1 2 i l = - k 0 1 i l 
k 1 3 i l = 2* eMod * i /l
k 2 0 i l = - k 0 0 i l
k 2 1 i l = - k 0 1 i l
k 2 2 i l = k 0 0 i l
k 2 3 i l = - k 0 1 i l
k 3 0 i l = k 0 1 i l 
k 3 1 i l = k 1 3 i l
k 3 2 i l = - k 0 1 i l
k 3 3 i l = k 1 1 i l
k _ _ _ _ = error "Index out of bounds: Local Stiffness Matrix is only 4x4"

localStiffness (i1:[]) (l1:[]) 0 = [ (a , b, (k a b i1 l1) ) | a <- [0..3], b <- [0..3], not(a>1 && b> 1)]
localStiffness (i1:i2:[]) (l1:l2:[]) rn = 
    firstQ ++ rest
    where 
        firstQ = [ (rn + a, rn + b, ((k (a+2) (b+2) i1 l1) + (k a b i2 l2) )) | a <- [0..1], b <- [0..1]]
        rest = [ (rn + a, rn + b, ( k a b i2 l2 )) | a <- [0..3], b <- [0..3], not(a<2 && b<2)]

localStiffness (i1:i2:it) (l1:l2:lt) rn = 
    firstQ ++ sndTrdQ
    where 
        firstQ = [ (rn + a, rn + b, ((k (a+2) (b+2) i1 l1) + (k a b i2 l2) )) | a <- [0..1], b <- [0..1]]
        sndTrdQ = [ (rn + a, rn + b, ( k a b i2 l2 )) | a <- [0..3], b <- [0..3], not(a<2 && b<2), not(b>1 && a >1)]

createStiffnessMatrix (i1:i2:iTail) (l1:l2:lTail) [] rn = 
   createStiffnessMatrix (i1:i2:iTail) (l1:l2:lTail) (localStiffness [i1] [l1] 0) 2
createStiffnessMatrix (i1:i2:[]) (l1:l2:[]) matrix rn =
    fromListSM (rn+4,rn+4) ( matrix ++ localStiffness [i1,i2] [l1,l2] (rn))
createStiffnessMatrix (i1:iTail) (l1:lTail) matrix rn = 
   createStiffnessMatrix (iTail) (lTail) ((localStiffness (i1:iTail) (l1:lTail) rn)++ matrix) (rn+2)

insertSprings :: SpMatrix Double -> [(Int, Int, Double)] -> SpMatrix Double
insertSprings kMat springs = 
        foldl (addSpring) kMat springs
        where
                addSpring k (a, b, v) = insertSpMatrix a b ( ( k @@! (a,b) ) + v) k 

main = do
        let l=replicate 500 100
        let i=replicate 500 400000000
        let s = [(x,x,2000)|x<- [200,400..1000]]
        let k=createStiffnessMatrix i l [] 0
        let knew = insertSprings  k s 
        let f= [ (x , 0) | x<-[0..9]] ++ [(1000,5000.0),(1001,819000000.0)]
        let fVec = (fromListSV 1002 f)
        --let res= (\(lo, up)-> luSolve lo up (dropSV 2 fVec)) $ lu (extractSubmatrixSM (\x -> x-2) (\x -> x-2) knew (2,1001) (2,1001))
        let res = (extractSubmatrixSM (\x -> x-2) (\x -> x-2) knew (2,1001) (2,1001)) <\> (dropSV 2 fVec)
        let myres = map (\x -> lookupSV x res) [198,398..998]
        print myres

[ocramz]

Hi, thank you for reporting this.

That error means that the QR factorization (which is used internally in <\>) fails due to some bad matrix shape.

I am not familiar with this particular problem domain so it would take me some time to refactor your code; below you can find a slightly more idiomatic version.

Please upgrade to the latest version of the library and notice that methods such as lu and <\> now return in a monad in order to raise exceptions (this is captured in the code below).

One very useful thing you could do right away is getting rid of all the hardcoded constants and make this code parametric (so that we can more quickly explore smaller problem sizes).

import Numeric.LinearAlgebra.Sparse
import Data.Sparse.Common

eMod = 210000.0

k :: Int -> Int -> Double -> Double -> Double
k 0 0 i l = 12* eMod * i /(l**3)
k 0 1 i l = 6 * eMod * i /(l**2)
k 0 2 i l = (- (k 0 0 i l))
k 0 3 i l = k 0 1 i l 
k 1 0 i l = k 0 1 i l 
k 1 1 i l = 4 * eMod *i / l
k 1 2 i l = - k 0 1 i l 
k 1 3 i l = 2* eMod * i /l
k 2 0 i l = - k 0 0 i l
k 2 1 i l = - k 0 1 i l
k 2 2 i l = k 0 0 i l
k 2 3 i l = - k 0 1 i l
k 3 0 i l = k 0 1 i l 
k 3 1 i l = k 1 3 i l
k 3 2 i l = - k 0 1 i l
k 3 3 i l = k 1 1 i l
k _ _ _ _ = error "Index out of bounds: Local Stiffness Matrix is only 4x4" 

localStiffness :: 
     [Double] -> [Double] -> Int -> [(Int, Int, Double)]
localStiffness [i1] [l1] 0 = [ (a , b, (k a b i1 l1) ) | a <- [0..3], b <- [0..3], not(a>1 && b> 1)]
localStiffness [i1,i2] [l1,l2] rn = 
    firstQ ++ rest
    where 
        firstQ = [ (rn + a, rn + b, ((k (a+2) (b+2) i1 l1) + (k a b i2 l2) )) | a <- [0..1], b <- [0..1]]
        rest = [ (rn + a, rn + b, ( k a b i2 l2 )) | a <- [0..3], b <- [0..3], not(a<2 && b<2)]
localStiffness (i1:i2:it) (l1:l2:lt) rn = 
    firstQ ++ sndTrdQ
    where 
        firstQ = [ (rn + a, rn + b, ((k (a+2) (b+2) i1 l1) + (k a b i2 l2) )) | a <- [0..1], b <- [0..1]]
        sndTrdQ = [ (rn + a, rn + b, ( k a b i2 l2 )) | a <- [0..3], b <- [0..3], not(a<2 && b<2), not(b>1 && a >1)]

createStiffnessMatrix :: [Double]
     -> [Double] -> [(Int, Int, Double)] -> Int -> SpMatrix Double
createStiffnessMatrix (i1:i2:iTail) (l1:l2:lTail) [] _ = 
   createStiffnessMatrix (i1:i2:iTail) (l1:l2:lTail) (localStiffness [i1] [l1] 0) 2
createStiffnessMatrix [i1, i2] [l1, l2] matrix rn =
    fromListSM (rn+4,rn+4) ( matrix ++ localStiffness [i1,i2] [l1,l2] (rn))
createStiffnessMatrix (i1:iTail) (l1:lTail) matrix rn = 
   createStiffnessMatrix iTail lTail ((localStiffness (i1:iTail) (l1:lTail) rn)++ matrix) (rn+2)

insertSprings :: (Foldable t, Num a) =>
     SpMatrix a -> t (IxRow, IxCol, a) -> SpMatrix a
insertSprings kMat springs = 
        foldl (addSpring) kMat springs
        where
                addSpring k (a, b, v) = insertSpMatrix a b ( ( k @@! (a,b) ) + v) k 

main = do
        let l=replicate 500 100
        let i=replicate 500 400000000
        let s = [(x,x,2000)|x<- [200,400..1000]]
        let k=createStiffnessMatrix i l [] 0
        let knew = insertSprings  k s
        let f= [ (x , 0) | x<-[0..9]] ++ [(1000,5000.0),(1001,819000000.0)] :: [(Int, Double)]
        let fVec = (fromListSV 1002 f) 

        res <- extractSubmatrixSM (subtract 2) (subtract 2) knew (2,1001) (2,1001) <\> dropSV 2 fVec
        let myres = map (`lookupSV` res) [198,398..998]
        print myres

[denjoh]

Thank you for your feedback.

Is there something that I can do to help you with this library? I have just started learning Haskell so it should not be any complicated.

[ocramz]

Hi @denjoh , thank you for the offer; I am currently working on a new release that would make it more robust and easier to use.

If you could provide a parametric version of your code above (e.g. in which the number of finite elements can be changed), that would already be a very useful first contribution, since it would enable to characterize the linear solver performance.

[denjoh]

Sounds good, I will do so and get back to you when I have done that.

[ocramz]

Hi again @denjoh , in case you haven't noticed you can also use linSolve0, which is the interface to the various iterative methods. If the system matrix is positive semidefinite or close to it, the CGS method works quite well. The syntax would be xhat <- linSolve0 CGS_ aa b x0, where x0 is a starting vector (you can create one with constv n 0.1 where n is the size).

[ocramz]

Closing this for now, still waiting for @denjoh 's input