How can you quickly map the indices of a banded matrix to a 1-dimensional array?

Question

This is closely related to a the question: How to map the indexes of a matrix to a 1-dimensional array (C++)?

I need to assign a reversible index to each non-zero element in a banded matrix.
In the normal, full matrix it is easy to do:

           |-------- 5 ---------|
  Row      ______________________   _ _
   0      |0    1    2    3    4 |   |
   1      |5    6    7    8    9 |   4
   2      |10   11   12   13   14|   |
   3      |15   16   17   18   19|  _|_
          |______________________|
Column     0    1    2    3    4

To find the array index we just use the following bijective formula:

matrix[ i ][ j ] = array[ i*m + j ]

In my case, we have a symmetrically banded matrix with some constraint on distance from the diagonal. For example, the following uses an upper and lower bound of 1:

           |-------- 5 ---------|
  Row      ______________________   _ _
   0      |0    1    X    X    X |   |
   1      |2    3    4    X    X |   4
   2      |X    5    6    7    X |   |
   3      |X    X    8    9    10|  _|_
          |______________________|
Column     0    1    2    3    4

In this case, I want to assign an index position to each element within the bandwidth, and ignore everything outside. There are a couple of ways to do this, one of which is to create a list of all the acceptable indices ix's, and then use map lookups to quickly go back and forth between a (row,col) pair and a singular index:

ix's :: [(Int,Int)] -- List of all valid indices 

lkup :: Map (Int,Int) Int
lkup = M.fromList $ zip ix's [0..]

rlkup :: Map Int (Int, Int)
rlkup = M.fromList $ zip [0..] ix's

fromTup :: (Int, Int) -> Int
fromTup tup = fromMaybe 0 $ M.lookup tup lkup

toTup :: Int -> (Int, Int)
toTup i = fromMaybe (0,0) $ M.lookup i rlkup

For large matrices, this leads to a huge number of map lookups, which causes a bottleneck. Is there a more efficient formula to translate between the valid addresses, k, and (row,col) pairs?

K. A. Buhr · Answer 1 · 2017-09-29T02:44:34.167

You might find it more straightforward to "waste" a few indexes at the beginning and end of the matrix, and so assign:

  Row         ______________________   _ _
   0    (0)  |1    2    X    X    X |   |
   1         |3    4    5    X    X |   4
   2         |X    6    7    8    X |   |
   3         |X    X    9   10   11 |  _|_
             |______________________|
Column        0    1    2    3    4

where (0) is an ignored index.

This is similar to the band matrix representation used by the highly respected LAPACK library.

You just need to take care that the unused elements are properly ignored when performing operations where they might affect used elements. (For example, a fast fill routine can be written without regard to which elements are used or unused; but a matrix multiplication would need to take a little more more care.)

If you take this approach, then the bijections are pretty simple:

import Data.Char
import Data.Maybe

type Index = Int

-- |(row,col) coordinate: (0,0) is top level
type Coord = (Int, Int)

-- |Matrix dimensions: (rows, cols, edges) where edges gives
-- the count of auxiliary diagonals to *each side* of the main
-- diagonal (i.e., what you call the maximum distance), so the
-- total band width is 1+2*edges
type Dims = (Int, Int, Int)

-- |Get index for (row,col)
idx :: Dims -> Coord -> Index
idx (m, n, e) (i, j) = let w = 1+2*e in w*i+(j-i+e)

-- |Get (row,col) for index
ij :: Dims -> Index -> Coord
ij  (m, n, e) idx    = let w = 1+2*e
                           (i, j') = idx `quotRem` w
                       in (i, j'+i-e)

--
-- test code
--

showCoords :: Dims -> [(Coord, Char)] -> String
showCoords (m, n, _) cs =
  unlines $
  for [0..m-1] $ \i ->
    for [0..n-1] $ \j ->
      fromMaybe '.' $ lookup (i,j) cs
  where for = flip map

test :: Dims -> IO ()
test dm@(m,n,_) = do
  putStrLn $ "Testing " ++ show dm
  let idxs = [0..]
  -- get valid index/coordinates for this matrix
  let cs = takeWhile (\(_, (i,j)) -> i<m || j<n)
           $ filter (\(_, (i,j)) -> i>=0 && j>=0)
           $ map (\ix -> (ix, ij dm ix)) idxs
  -- prove the coordinates are right
  putStr $ showCoords dm (map (\(ix, (i,j)) -> ((i,j), chr (ord 'A' + ix))) cs)
  -- prove getIndex inverts getCoord
  print $ all (\(ix, (i,j)) -> idx dm (i,j) == ix) cs
  putStrLn ""

main = do test (4, 5, 1)   -- your example
          test (3, 8, 2)   -- another example

I actually did try going this route originally, but the band storage matrix, while simpler, really just shifts the problems from one set of corners to the other. It's harder to see that in this toy example, but you also have to throw away certain indices in the middle of the array as well, not just those at the beginning and edge. — jkeuhlen, Sep 29 '17 at 02:36
I'm not sure I understand. The point of the LAPACK-like storage is that you don't worry about those things. There can be unused indexes in either corner. — K. A. Buhr, Sep 29 '17 at 02:40
My problem has been efficiently marching through the acceptable indices. You still have to worry about the unused positions when deciding a formula to maybe between `(i,j)` pairs and your `ix` right? I'll take a closer look at your code tomorrow morning and see if there's something I'm missing, but I did go down this path before (even looking at the linked page). — jkeuhlen, Sep 29 '17 at 02:47

How can you quickly map the indices of a banded matrix to a 1-dimensional array?

1 Answers1