Sometimes in matrix computations you need to obtain the values of certain submatrices such as the diagonal elements or the super- or subdiagonal elements. About a year ago, I showed one way to do that: convert subscripts to indices and vice-versa.

However, a tip from @RLangTip on Twitter got me thinking about a different way to solve the same problem. The tweet was the following:

A single index can be used to select from a matrix. For example, X[row(X)==col(X)] is the same as diag(X) http://t.co/EY1uMfny

proc iml;
/* return matrix m such that m[i,j] = j */
start col(x);
   return( repeat(1:ncol(x), nrow(x)) );
finish;
/* return matrix m such that m[i,j] = i */
start row(x);
   return( repeat( T(1:nrow(x)), 1, ncol(x) ));
finish;
 
x=shape(1:16, 4); /* define 4x4 matrix */
r = row(x);
c = col(x);
print x, r, c;

Update: The ROW function and COL function are distributed with SAS/IML 12.3, which was released with SAS 9.4.

When I started thinking about these matrices, I realized that they can be used to assign, extract, and copy super- and subdiagonal elements, as well as upper and lower triangular portions of a matrix. For example, the following statements extract the diagonal, the superdiagonal, and the subdiagonal elements of the matrix x:

/* find indices of diagonal, superdiagonal, and subdiagonal */
diagIdx = loc(r=c);
superIdx = loc(r=c-1);
subIdx = loc(r=c+1);
 
/* extract elements */
diag = x[diagIdx];
superDiag = x[superIdx];
subDiag = x[subIdx];
print diag subDiag superDiag;

You can also use the diagIdx, superIdx, and subIdx vectors to assign values to the x matrix. For example,
x[diagIdx]=1:4;
assigns the values 1, 2, 3, and 4 to the diagonal elements of x.

You can also use the ROW and COL functions to extract or assign values for the upper or lower triangular portions of a matrix. For example, the following statements fill a matrix with uniform random numbers and copy values from the lower triangular portion into the upper triangular portion, which results in a random symmetric matrix:

/* create random symmetric matrix */
call randseed(1234);
x = j(4, 4);
call randgen(x, "Uniform");
 
upperTri = loc(r <= c); /* upper tri indices in row major order */
v = vech(x);            /* lower tri element in col major order */
x[uppertri] = v;        /* copy elements */
print x[format=4.3];

The advantage of using this approach is its convenience. This technique makes it simple to specify elements of a banded or triangular matrix. For dealing with small matrices, I like this method and I will add it to my programming toolkit.

The disadvantage of this method is its inefficiency, especially for large matrices. The ROW and COL functions each return matrices that are the same size as x. If want to extract the diagonal of a 1000x1000 matrix, I need to create two additional matrices of the same size, plus there are further (temporary) matrices created by expressions such as r=c or r > c, which create matrices of zeros and ones. Although I like the convenience of the ROW and COL functions, for the sake of efficiency I don't plan to discard my functions that convert subscripts to indices.