* Bug 15428 fixed: repmat() was slow: rewritten => > 7x faster
[scilab.git] / scilab / modules / elementary_functions / macros / repmat.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) - 2011 -  INRIA, Serge Steer
3 // Copyright (C) 2012 - 2016 - Scilab Enterprises
4 // Copyright (C) - 2018 - Samuel GOUGEON
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12
13 function B = repmat(A,varargin)
14
15     // CHECKING INPUT ARGUMENTS
16     // ------------------------
17     rhs = argn(2);
18     if rhs < 2 then
19         msg = _("%s: Wrong number of input arguments: at least %d expected.\n")
20         error(msprintf(msg, "repmat", 2))
21     end
22
23     if A == [] then
24         B = A
25         return
26     end
27
28     narg = size(varargin);
29
30     // Reading & setting replication numbers
31     if narg == 1 then
32         // Scalar of vector needed
33         v = varargin(1);
34         if type(v) <> 1 | ~isreal(v, 0) | ~and(v == int(v)) | v == [] | v <0 then
35             msg = _("%s: Argument #%d: Non-negative integers expected.\n")
36             error(msprintf(msg, "repmat", 2))
37         else
38             v = real(v);    // in case of complex encoding with null imag
39         end
40         if ~isscalar(v) & ~isvector(v) then
41             msg = _("%s: Wrong size for input argument #%d: A scalar or vector expected.\n")
42             error(msprintf(msg, "repmat", 2))
43         end
44         sizes = v;
45         if length(v) == 1 then
46             sizes = [v,v];
47         end
48     else
49         sizes = [];
50         for i = 1:narg
51             v = varargin(i);
52             if type(v) <> 1 | ~isreal(v, 0) | ~and(v == int(v)) | v == [] | v < 0 then
53                 msg = _("%s: Argument #%d: Non-negative integers expected.\n")
54                 error(msprintf(msg, "repmat", i+1))
55             end
56             if length(v)<>1 then
57                 msg = _("%s: Argument #%d: Scalar (1 element) expected.\n")
58                 error(msprintf(msg, "repmat", i+1))
59             end
60             sizes = [sizes v];
61         end
62     end
63
64     // Preprocessing sizes = replication factors
65     // -----------------------------------------
66     // Trimming ending ones in sizes
67     for s = length(sizes):-1:3
68         if sizes(s) == 1 then
69             sizes(s) = []
70         else
71             break
72         end
73     end
74     // B expected sizes
75     sA = size(A)
76     ndimsA = length(sA)
77     ndimsR = length(sizes)
78     if ndimsA >= ndimsR then
79         Bsizes = size(A).*[sizes ones(1,(ndimsA-ndimsR))]
80     else
81         Bsizes = [sA ones(1,(ndimsR-ndimsA))] .* sizes
82     end
83     // Sparse matrices are only 2D (limited output dimensions):
84     if or(type(A) == [5 6]) & length(Bsize) > 2 then
85         msg = _("%s: Wrong dimensions for input arguments: The sparse result must be 2D, not %dD.\n");
86         error(msprintf(msg, "repmat", length(Bsizes)));
87     end
88     // Replication factors
89     Rfactors = list();
90     for s = sizes
91         Rfactors($+1) = s
92     end
93
94     // PROCESSING
95     // ----------
96     if find(sizes == 0) then
97         B = []
98         return
99     end
100
101     // Special case where the input A is scalar: special efficient processing
102     if size(A,"*") == 1 then
103         i = uint8(0);
104         i(Rfactors(:)) = 0;
105         i = i + 1;
106         if typeof(A) <> "rational" then
107             B = matrix(A(i), sizes);
108         else
109             B = matrix(A(i,0), sizes);
110         end
111         return
112     end
113
114     // Numerical data: double, integers, polynomials, sparse : kron() is used
115     if or(type(A) == [1 2 5 7 8]) then
116         B = ones(Rfactors(:)) .*. A;
117         return
118     end
119
120     // Other regular types : booleans, textes..:
121     // => we replicate and use the matrix of linearized indices of A components
122     if or(type(A) == [4 6 9 10]) then
123         // Processing (without intermediate variable, to save memory)
124         B = matrix(A(ones(Rfactors(:)) .*. matrix(1:size(A,"*"), size(A))), Bsizes);
125         return
126     end
127
128     // Rationals, other mlists, overloads
129     if typeof(A) == "rational" then
130         B = rlist(repmat(A.num,sizes), repmat(A.den,sizes), A.dt)
131     else
132         execstr("B=%" + typeof(A) + "_repmat(A, sizes)")
133     end
134 endfunction