Fast matrix copying and array appending
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponses (2)
Jan
le 30 Jan 2013
Modifié(e) : Jan
le 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?
4 commentaires
Jan
le 30 Jan 2013
Modifié(e) : Jan
le 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.
Matt J
le 30 Jan 2013
Modifié(e) : Matt J
le 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.
1 commentaire
Matt J
le 30 Jan 2013
Modifié(e) : Matt J
le 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.
Voir également
Catégories
En savoir plus sur Logical dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!