/***************************************************************************/
/* */
/* ACM Algorithm 302 - TRANSPOSE VECTOR STORED ARRAY. */
/* */
/* J. Boothroyd, University of Tasmania, Hobart, Tasmania, Australia. */
/* */
/***************************************************************************/
/* Translated to PL/I by R. A. Vowels, 11 February 2006. */
/* This procedure performs an in-situ transposition of an m x n array
A(1:m, 1:n) stored by rows in the vector a(1:m*n). The method is
essentially that of Windley [P. F. Windley, "Transposing matrices
in a digital computer", Comp. J. 2 (1959), pp. 47-48], modified
for use with vectors having unit lower bounds.
The algorithm processes only elements A(1,2) through
A(m,n-1) since A(1,1) and A(m,n) retain their original positions.
Elements A(q,p) of the transposed matrix are placed in
a(i), in the order i = 2, 3, ..., mn-2, by an exchanging process.
At the last step two elements are correctly placed which
accounts for the value mn-2 as the upper bound on i. Valid
subscripts of the vectora(1:m*n) are elements in the 1-origin
index of set [1, 2, ..., mn]. Computationally, however, it is more
convenient to use the zero-origin set [0, 1, ..., mn-1]. Denoting
by i(0) (i(0) = i-1) the corresponding zero-origin index of
a(i), to be occupied by A(q,p), we have i = m(q-1) + (p-1).
The corresponding zero-origin index j(0) of the A(p,q) element
now in a(j), which must be transferred to a(i), is:
j(0) = j-1 = n(p-1) + (q-1) = n*i(0) mod(mn-1).
For each value of i = 2, 3, ..., mn-2 (or i(0) =
1, 2, ..., mn-3) we compute the index j of a(j) and exchange
a(i) and a(j) provided j >= i (i.e., j(o) >= i(0)). The case j < i indicates
that an element originally in a(j) is now elsewhere following
previous exchanges. Its present position is given by the first
j(r) > i(0) in the series of zero-origin indices:
j(0), j(r+1) = n * j(r) mod(mn-1).
The two sequences modulo(mn-1) are generated by different
methods. An additive process generates the first, using k to
dulplicate the function of j, in case this is adjusted in the second
recurrence-generated sequence if j < i.
{one paragraph omitted]
*/
/* The algorithm was tested by correctly transposing a 3 by 4 matrix */
transpose:
procedure (a, m, n) options (reorder);
declare A(*) float, (m, n) fixed binary;
declare (i, j, k, iless1, mnless1, jn, modlessn, it) fixed binary,
t float;;
mnless1 = m*n-1; modlessn = mnless1 - n;
k = 0; iless1 = 1;
do i = 2 to mnless1-1;
/* computes j, k = n*i(0) mod(mn-1) */
if k <= modlessn then k = k + n; else k = k - modlessn;
j= k; /* [can be shortened to a multiple ass. stmt] */
do while (j < iless1);
/* computes j(r+1) = n * j(r)mod(mn-1) */
jn = j*n;
it = jn/mnless1;
j = jn - it * mnless1;
end;
/* avoid unnecessary exchanges */
if j ^= iless1 then
do;
j = j + 1;
t = a(i); a(i) = a(j); a(j) = t;
end;
iless1 = i;
end;
end transpose;