Matrix transposition: Difference between revisions

Playing more to c's strengths, the following implementation

transposes a matrix of any type and dimensions

in place with only O(1) space. See the [http://www.en.wikipedia.org/wiki/In-place_matrix_transposition Wikipedia article]

for more information.

<lang C>

void

*transpose_matrix(matrix,rows,cols,item_size)

void *matrix;

int rows;

int cols;

size_t item_size;

{

#define ALIGNMENT 16 /* power of 2 >= minimum array boundary alignment; maybe unnecessary but machine dependent */

char *cursor;

char carry[ALIGNMENT];

size_t block_size,remaining_size;

int nadir,lag,orbit,ents;

if ((rows == 1) ? 1 : (cols == 1))

return matrix;

ents = rows * cols;

cursor = (char *) matrix;

remaining_size = item_size;

block_size = ((ALIGNMENT < remaining_size) ? ALIGNMENT : remaining_size);

while (block_size)

{

nadir = 1; /* first and last entries are always fixed points so aren't visited */

while (nadir + 1 < ents)

{

memcpy(carry,&(cursor[(lag = nadir) * item_size]),block_size);

while ((orbit = (lag / rows) + cols * (lag % rows)) > nadir) /* follow a complete cycle */

{

memcpy(&(cursor[lag * item_size]),&(cursor[orbit * item_size]),block_size);

lag = orbit;

}

memcpy(&(cursor[lag * item_size]),carry,block_size);

orbit = nadir++;

while ((orbit < nadir) ? (nadir + 1 < ents) : 0) /* find the next unvisited index by an exhaustive search */

{

orbit = nadir;

while ((orbit = (orbit / rows) + cols * (orbit % rows)) > nadir);

nadir = ((orbit < nadir) ? nadir + 1 : nadir);

}

}

cursor = cursor + block_size;

remaining_size = remaining_size - block_size;

block_size = ((ALIGNMENT < remaining_size) ? ALIGNMENT : remaining_size);

}

return matrix;

}

</lang>

No extra storage allocation is required by the caller. Here are

usage examples for an array of doubles and an array of complex numbers.

<lang c>

a = (double *) transpose_matrix((void *) a, n, m, sizeof(double));

b = (complex *) transpose_matrix((void *) b, n, m, sizeof(complex));

</lang>

After execution, the memory maps of a and b will be those of m by n arrays instead

of n by m.