25

I am working on my data in a C/C++ program, which is 2 dimensional. Here my value is calculated for pair-wise and here values would be same for foo[i][j] and foo[j][i].

Thus if I implement it by using a simple 2 dimensional array, half of my space would be wasted. So what would be best data structure to represent this lower/upper triangular matrix.

Regards,

1

7 Answers 7

19

If you have N items then a lower triangular array without the main diagonal will have (N - 1) * N / 2 elements, or (N + 1) * N / 2 elements with the main diagonal. Without the main diagonal, (I, J) (I,J ∈ 0..N-1, I > J) ⇒ (I * (I - 1) / 2 + J). With the main diagonal, (I,J ∈ 0..N-1, I ≥ J) ⇒ ((I + 1) * I / 2 + J).

(And yes, when you're allocating 4 gigabytes on a 2.5 gigabyte machine, cutting it half does make a huge difference.)

3
  • This is actually Szudzik's "elegant pairing function", as can be found in this paper: szudzik.com/ElegantPairing.pdf
    – awllower
    Feb 3, 2022 at 15:35
  • @awllower I can't reach that site from here, but I came up with that off the top of my head, and if that's Matthew Peter Szudzik, born 1973, then I'm sure the first time it was discovered was before he was born.
    – prosfilaes
    Feb 3, 2022 at 16:30
  • Thanks for the information. I just found that paper recently, so it is the first thing that popped to my mind when I saw your answer.
    – awllower
    Feb 4, 2022 at 7:40
12

Really, you're best off just using a regular two dimensional matrix. RAM is pretty cheap. If you really don't want to do that, then you can build a one-dimensional array with the right number of elements and then figure out how to access each element. For example, if the array is structured like this:

    j
    1234
i 1 A
  2 BC
  3 DEF
  4 GHIJ

and you have it stored as a one dimensional array, left to right, you'd access element C (2, 2) with array[3]. You can work out a function to go from [i][j] to [n] but I won't spoil your fun. But you don't have to do this unless the triangular array in question is really huge or you're very concerned about space.

5

As Dan and Praxeolitic proposed for lower triangular matrix with diagonal but with corrected transition rule.

For matrix n by n you need array (n+1)*n/2 length and transition rule is Matrix[i][j] = Array[i*(i+1)/2+j].

#include<iostream>
#include<cstring>

struct lowerMatrix {
  double* matArray;
  int sizeArray;
  int matDim;

  lowerMatrix(int matDim) {
    this->matDim = matDim;
    sizeArray = (matDim + 1)*matDim/2;
    matArray = new double[sizeArray];
    memset(matArray, .0, sizeArray*sizeof(double));
  };

  double &operator()(int i, int j) {
    int position = i*(i+1)/2+j;
    return matArray[position];
  };
};

I did it with double but you can make it as template. This is just basic skeleton so don't forget to implement destructor.

3

Use a jagged array:

int N;
// populate N with size

int **Array = new Array[N];
for(int i = 0; i < N; i++)
{
    Array[i] = new Array[N - i];
}

it will create array like

   0 1 2 3 4 5
0 [           ]
1 [         ]
2 [       ]
3 [     ]
4 [   ]
5 [ ]
2
  • 12
    This will allocate the individual arrays independently, which may be bad for cache behaviour and memory fragmentation. This may be OK if you don't care too much about performance, but in that case you should probably just use a single NxN array. If you do decide that you want to use an array of pointers anyway, then allocate the N*(N+1)/2 elements in a single array and create the row pointers as offsets into that array.
    – Erik P.
    Oct 31, 2011 at 1:49
  • @ErikP. : I know making a class with continuous array and access methods that calculate the offset is better, but this is a much simpler way.
    – Daniel
    Oct 31, 2011 at 2:03
1

The number of unique elements, m, needed to be represented in an n by n symmetric matrix:

With the main diagonal

m = (n*(n + 1))/2

Without the diagonal (for symmetric matrix as the OP describes, main diagonal is needed, but just for good measure...)

m = (n*(n - 1))/2.

Not dividing by 2 until the last operation is important if integer arithmetic with truncation is used.

You also need to do some arithmetic to find the index, i, in the allocated memory corresponding to row x and column y in the diagonal matrix.

Index in allocated memory, i, of row x and column y in upper diagonal matrix:

With the diagonal

i = (y*(2*n - y + 1))/2 + (x - y - 1)

Without the diagonal

i = (y*(2*n - y - 1))/2 + (x - y -1)

For a lower diagonal matrix flip x and y in the equations. For a symmetric matrix just choose either x>=y or y>=x internally and have member functions flip as needed.

2
  • 3
    This seems not quite right - plugging (0,0) into your "with the diagonal" yields -1. Mar 14, 2016 at 17:02
  • i = (y* (2 * n - y + 1)) / 2 + x;
    – Catriel
    Aug 18, 2020 at 18:00
0

Riffing on Dani's answer...

Instead of allocating many arrays of various sizes, which could lead to memory fragmentation or weird cache access patterns, you could allocate one array to hold the data and one small array to hold pointers to the rows within the first allocation.

const int side = ...;
T *backing_data = new T[side * (side + 1) / 2];  // watch for overflow
T **table = new T*[side];
auto p = backing_data;
for (int row = 0; row < side; ++row) {
   table[row] = p;
   p += side - row;
}

Now you can use table as though it was a jagged array as shown in Dani's answer:

table[row][col] = foo;

But all the data is in a single block, which it might not otherwise be depending on your allocator's strategy.

Using the table of row pointers may or may not be faster than computing the offset using Praxeolitic's formula.

-1
#include <stdio.h>

// Large math problems running on massively parallel systems sometimes use a lot
// of upper triangular matrices.  Implemented naively, these waste 50% of the memory
// in the machine, which is not recoverable by virtual memory techniques because it
// is interspersed with data on each row.  By mapping the array elements into an
// array that is half the size and not actually storing the zeroes, we can do twice
// the computation in the same machine or use half as many machines in total.

// To implement a safety feature of making the zero-part of an upper triangular matrix
// read-only, we place all the zeroes in write-protected memory and cause a memory violation
// if the programmer attempts to write to them.  System dependent but relatively portable.
// Requires that you compile with the -Wno-discarded-qualifiers option.

// for the awkward case (an even-sized array bound):

//                         +--------/
//     row  0, 40 items -> |0      /
//     row  1, 39 items -> |      /
//     row 19, 21 items -> |     /
//     row 20, 20 items -> |----/ <------  cut and rotate here to form a rectangle.
//     row 21, 19 items -> |   /
//                         |  /
//     row 39,  1 item  -> | /
//     row 40,  0 items -> |/
//                         /


//    x y   x  y
//    0,0   39,0
//     +----/           |                     +--------/
//     |   /            | row  0, 40 items -> |0      /| <-- row 40, 0 items
//     |  / - 20,18     | row  1, 39 items -> |      /0| <-- row 39, 1 item
//     | /\             | row 19, 21 items -> |     /  | <-- row 21, 19 items
//     |/  19,19        | row 20, 20 items -> |    /???| <-- row 20, 20 items  half of row 20 is wasted...
//0,39 v                |                     ~~~~~~~~~~
//     |                |

// for odd-sized array bounds, there is no need for the wasted half-row marked '???' above...

// And once the mapping above is done, mirror the indexes in x to get a proper Upper Triangular Matrix which
// looks like this...
//     ____
//     \  |
//      \ |
//       \|
//

// Rather than store the underlying data in a 2D array, if we use a 1-D array,
// and map the indexes ourselves, it is possible to recover that final half-row...

// The implementation allows for the matrix elements to be any scalar type.

#define DECLARE_TRIANGULAR_MATRIX(type, name, bound, zero)                      \
  type _##name[bound * (bound+1) / 2 + 1]; /* +1 is for a debugging tombstone */ \
  type *__##name(int x, int y) { \
    static const type Zero = zero; /* writing to the lower half of the matrix will segfault */ \
    x = (bound-1)-x; /* mirror */ \
    if (x+y >= bound) return &Zero; /* requires cc -Wno-discarded-qualifiers */ \
    if (y > bound/2) {x = (bound-1)-x; y = bound-y;} \
    return &_##name[y*bound+x]; /* standard mapping of x,y -> X */ \
  }
#define TRIANGULAR_MATRIX(name, x, y)  *__##name(x,y)


// ----------------------------------------------------------------------------------------


// Simulate 'int fred[11][11];' as an upper triangular matrix:
#define ARRAYSIZE 11
DECLARE_TRIANGULAR_MATRIX(int, fred, ARRAYSIZE, 0)
#define fred(x, y) TRIANGULAR_MATRIX(fred, x, y)
/* unfortunately we can't #define fred[x][y] here ... In the Imp language which used () both
   for array indexes and procedure parameters, we could write a mapping function fred(x,y)
   which made the indirected function call indistinguishable from a normal array access.
   We attempt to do something similar here using macros, but C is not as cooperative. */



int main(int argc, char **argv) {
  int x,y, element;

  // treat as fully populated 2D array...
  for (y = 0; y < ARRAYSIZE; y++) {
    for (x = 0; x < ARRAYSIZE; x++) {
      if (y <= x) fred(x,y) = (x+1) * 100 + (y+1); // upper triangle test
    }
  }

  fprintf(stdout, "Upper Triangular matrix:\n\n");
  fprintf(stdout, "    ");
  for (x = 0; x < ARRAYSIZE; x++) fprintf(stdout, "%5d", x);
  fprintf(stdout, "\n    ");
  for (x = 0; x < ARRAYSIZE; x++) fprintf(stdout, "_____");
  fprintf(stdout, "\n");
  for (y = 0; y < ARRAYSIZE; y++) {
    fprintf(stdout, "%2d |", y);
    for (x = 0; x < ARRAYSIZE; x++) {
      element = fred(x,y);
      fprintf(stdout, "%5d", element);
      if (y <= x) { // upper triangle test
        if (element != (x+1) * 100 + (y+1)) {
          fflush(stdout); fprintf(stderr, "Mismatch! at %d,%d (%d != %d)\n", x, y, element, x * 100 + y);
        }
      } else if (element != 0) {
        fflush(stdout); fprintf(stderr, "Mismatch! at %d,%d (%d != 0)\n", x, y, element);
      }
    }
    fprintf(stdout, "\n");
  }

  return 0;
}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.