#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 (ARRAYSIZE-x+y <= ARRAYSIZE) 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 (ARRAYSIZE-x+y <= ARRAYSIZE) { // 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;
}
#ifdef IMP
begin
integer name Force an error { being uninitialised, Imp will raise a signal on access through the pointer... }
const integer ARRAY SIZE = 11
! for compatibility with the C version we will base our array at 0
! even though in Imp we would probably normally base it at 1...
integer array triangular fred(0 : ARRAY SIZE * (ARRAY SIZE+1)//2 + 1)
! Actually, in imp, we would more likely declare the array as a two dimensional array
! similar to: %integer %array triangular fred(0 : ARRAY SIZE, 0 : (ARRAY SIZE+1)//2)
integer map fred(integer x, y)
x = (ARRAY SIZE-1) - x
result == Force an error if x+y >= bound
if y > ARRAY SIZE//2 start
x = (ARRAY SIZE-1) - x; y = ARRAY SIZE - y
finish
result == triangular fred(y*ARRAY SIZE + x)
! or if we had used a 2d array, ... %result == triangular fred(x, y)
end
integer x, y, element
select output(0)
! treat as fully populated 2D array...
for y = 0, 1, ARRAY SIZE-1 cycle
for x = 0, 1, ARRAY SIZE-1 cycle
if ARRAY SIZE-x+y <= ARRAY SIZE start
fred(x,y) = (x+1) * 100 + (y+1) { upper triangle test }
finish
repeat
repeat
print string("Upper Triangular matrix:".snl.snl)
spaces(4)
write(x, 5) for x = 0, 1, ARRAY SIZE-1
newline; spaces(4);
print string("_____") for x = 0, 1, ARRAY SIZE-1
newline
for y = 0, 1, ARRAY SIZE-1 cycle
write(y, 2); print string(" |")
for x = 0, 1, ARRAY SIZE-1 cycle
element = fred(x,y)
write(element, 5);
if ARRAY SIZE-x+y <= ARRAY SIZE start { upper triangle test }
if element # (x+1) * 100 + (y+1) start
print string("Mismatch! at ")
write(x, 0); print symbol(','); write(y, 0)
print string(" ("); write(element, 0); print string(" # ")
write(x * 100 + y, 0); print symbol(')'); newline
finish
else if element # 0
print string("Mismatch! at ")
write(x, 0); print symbol(','); write(y, 0)
print string(" ("); write(element, 0); print string(" # 0)".snl)
finish
repeat
newline
repeat
end of program
#endif