Asked by nedo nodo
on 30 Jan 2013

Hi,

I would like to know if there exists a fast way to copy a subpart of a matrix into a new matrix and to append an array to the same matrix. For instance, something like this: if true A=[A(2:end,2:end),B;B',b]; end where A is nxn matrix, B is nx1 array and b is a scalar.

Thank you Best Regards Salvo

Answer by Jan Simon
on 30 Jan 2013

Edited by Jan Simon
on 30 Jan 2013

I'm convinced that this is a good idea already:

A = [A(2:end, 2:end), B; B', b];

You save some micro-seconds by avoiding `end`, because the explicit definition if the size is faster:

[s1, s2] = size(A); A = [A(2:s1, 2:s2), B; B', b];

In my measurements this is remarkably slower:

A(1:s1-1, 1:s2-1) = A(2:s1, 2:s2); A(s1, 1:s2-1) = B'; A(1:s1-1, s2) = B; A(s1, s2) = b;

This might save some temporary memory, but the indexing needs more time.

**[EDITED]**

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double *A, *B, *Bp, b, *R; mwSize i, n, n_m1;

% Get inputs: A = mxGetPr(prhs[0]); B = mxGetPr(prhs[1]); b = mxGetScalar(prhs[2]); n = mxGetM(prhs[0]);

% Create outputs: plhs[0] = mxCreateDoubleMatrix(n, n, mxReal); R = mxGetPr(plhs[0]);

Bp = B; n_m1 = n - 1; for (i = 0; i < n_m1; i++) { memcpy(R, A+1, n_m1 * sizeof(double)); R[n_m1] = *Bp++; R += n; A += n; }

memcpy(R, B, n_m1 * sizeof(double)); R[n_m1] = b; }

**Use this with great care and after exhaustive testing only!** This is written in the forum's editor and neither compiled nor tested. Do you know how to compile a C-Mex file?

José-Luis
on 30 Jan 2013

*your* computer.

Jan Simon
on 30 Jan 2013

@nedo nodo: I do not believe, that this assignment needs the most time in your program. Optimizing a not relevant part of the code does not lead to a significant acceleration of the whole program. So please use the profiler to find the bottlenecks, if you did not do this already.

You could create a C-MEX function for the creation of the matrix, but it is not sure if this is faster at all, because most of Matlab's indexing functions are very efficient already. If you find out, that this operation needs 90% of the overall processing time, I could offer to create a small C-code.

nedo nodo
on 30 Jan 2013

Log in to comment.

Answer by Matt J
on 30 Jan 2013

Edited by Matt J
on 30 Jan 2013

If this really is the performance critical part of the code, you should consider redesigning so that A(2:end,2:end) gets prepended instead of appended. It leads to orders of magnitude speed difference, as shown below. The redesign should be trivial if you are doing likely operations like A*x. Just do A*[x(end);x(1:end-1)] instead.

N=3000; A=rand(N); B=rand(N-1,1); b=1;

#########

tic %Appending A=[A(2:end,2:end),B;B',b]; toc %Elapsed time is 0.086694 seconds.

############

tic%Prepending

A(1)=b; A(2:end,1)=B; A(1,2:end)=B'; toc %Elapsed time is 0.000218 seconds.

Matt J
on 30 Jan 2013

You could also use my ProdCascade class so that even though you're prepending A, it will automatically pre- and post-permute things you want to multiply with A to give the same I/O behavior as if you had appended to A. The class introduces some overhead, but it's still close to three times as fast as your original approach when doing A*x operations:

tic %Appending A=[A(2:end,2:end),B;B',b]; z1=A*x; toc

Elapsed time is 0.079278 seconds.

P=circshift(speye(N),[1,0]); AA=ProdCascade({P.',A,P});

tic;

A(1)=b; A(2:end,1)=B; A(1,2:end)=B'; AA{2}=A; z2=AA*x; toc;

Elapsed time is 0.027279 seconds.

Log in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Log in to comment.