Asked by Florian
on 3 Feb 2013

Is there a way to shift every column c by (c-1)*n rows down, where n is size(A,1) in the following example, without a loop and without excessive memory usage? So to get

A = [ a11 a12 a13 a14 ... a21 a22 a23 a24 ... ... ... ... ... ... ]

to read

B = [ a11 0 0 0 ... a21 0 0 0 ... ... ... ... ... ... 0 a12 0 0 ... 0 a22 0 0 ... ... ... ... ... ... 0 0 a13 0 ... 0 0 a23 0 ... ... ... ... ... ... 0 0 0 a14 ... 0 0 0 a24 ... ... ... ... ... ... ]

I tried diag and sparse with no luck so far. I need Matrix A in an operation X=J\A afterwards. It is crucial that the zero elements do not consume memory as A can have up to 20e6 columns and probably 10 rows.

Thanks in advance.

*No products are associated with this question.*

Answer by Cedric Wannaz
on 3 Feb 2013

Edited by Cedric Wannaz
on 3 Feb 2013

Accepted answer

You could use `sparse()`, building appropriate vectors of indices and values, e.g.:

nCol = size(A, 2) ; nRow = nCol * size(A, 1) ; I = (1:nRow).' ; J = reshape(repmat(1:nCol, size(A, 1), 1), [], 1) ; B = sparse(I, J, A(:), nRow, nCol) ;

(that you can write in two lines actually, but I wanted to keep it clear)

Cedric Wannaz
on 3 Feb 2013

Well, one way to use sparse() is to provide 5 args: a vector or row indices, a vector of column indices, a vector of values, the number of rows and the number of columns. Here we build these "things", using one "trick" that simplifies a bit the problem: we create the vector of values using linear indexing on `A` (in this case, it means lining up columns of `A` and reading that as a vector). Example:

>> A = [1, 2, 3; 4, 5, 6] A = 1 2 3 4 5 6 >> A(:) ans = 1 4 2 5 3 6

Knowing that we will pass `A(:)` to sparse() as a vector of values, we have to build appropriate vectors of row and column indices..

Looking at what you depicted, row indices will just increase linearly, so we build `I` as `1:nRow`, transposed so it is a column vector (same dim. as `A(:)`).

For column indices, it is a little more tricky as it will be something like

`1 1 1 ... 2 2 2 ...` etc, with the number of each integer being the number of rows of `A`. We can obtain this by creating a matrix

1 2 3 4 ... 1 2 3 4 ... . . . . . . . .

that we ultimately index linearly. This matrix is actually the vertical repetition of `(1:nCol)` a number of times that matches the number of rows of `A`. You obtain that with `repmat()`. Example:

>> repmat(1:4, 3, 1) ans = 1 2 3 4 1 2 3 4 1 2 3 4

Now I could have written

... J = repmat(1:nCol, size(A, 1), 1)) ; B = sparse(I, J(:), A(:), nRow, nCol) ;

and avoid the reshape, but you would have had `I` given as a vector, and `J` as a matrix indexed linearly, which I thought would have brought confusion. The reshape transforms the output of `repmat()` into a column vector (last arg = 1). As `reshape()` is flexible, it allows to specify only one size (the other is left to an empty array), and computes the other. The reshape does therefore essentially what `J(:)` would have done.

Answer by bym
on 3 Feb 2013

s = spdiags(A',[0,-1],c*(n-1),c);

Florian
on 3 Feb 2013

Thank you for your answer.

Unfortunately it didn't do the trick. Actually no specific value is assign to `"c"`. In my question it was meant to be a variable to describe that the first column `(c=1)` should be moved `(c-1)*n=0` rows down, the second column `(c=2)` should be moved by `(c-2)*n=n` rows down, the third column `(c=3)` should be moved `(c-2)*n=2n` rows down, etc.

From what I see, `diags` or `spdiags` places its elements only on true diagonals, shifted or not. In my case, however, the elements are not placed on a straight diagonal line but in a rather zig-zag way.

Opportunities for recent engineering grads.

## 0 Comments