Multiplying 1-D or 2-D subarrays contained in two N-D arrays.
THIS ORIGINAL MULTIPROD IS NOW CALLED MULTIPROD_LEGACY IN MANOPT
C = MULTIPROD(A,B) is equivalent to C = MULTIPROD(A,B,[1 2],[1 2])
C = MULTIPROD(A,B,[D1 D2]) is eq. to C = MULTIPROD(A,B,[D1 D2],[D1 D2])
C = MULTIPROD(A,B,D1) is equival. to C = MULTIPROD(A,B,D1,D1)
MULTIPROD performs multiple matrix products, with array expansion (AX)
enabled. Its first two arguments A and B are "block arrays" of any
size, containing one or more 1-D or 2-D subarrays, called "blocks" (*).
For instance, a 5x6x3 array may be viewed as an array containing five
6x3 blocks. In this case, its size is denoted by 5x(6x3). The 1 or 2
adjacent dimensions along which the blocks are contained are called the
"internal dimensions" (IDs) of the array (°).
1) 2-D by 2-D BLOCK(S) (*)
C = MULTIPROD(A, B, [DA1 DA2], [DB1 DB2]) contains the products
of the PxQ matrices in A by the RxS matrices in B. [DA1 DA2] are
the IDs of A; [DB1 DB2] are the IDs of B.
2) 2-D by 1-D BLOCK(S) (*)
C = MULTIPROD(A, B, [DA1 DA2], DB1) contains the products of the
PxQ matrices in A by the R-element vectors in B. The latter are
considered to be Rx1 matrices. [DA1 DA2] are the IDs of A; DB1 is
the ID of B.
3) 1-D by 2-D BLOCK(S) (*)
C = MULTIPROD(A, B, DA1, [DB1 DB2]) contains the products of the
Q-element vectors in A by the RxS matrices in B. The vectors in A
are considered to be 1xQ matrices. DA1 is the ID of A; [DB1 DB2]
are the IDs of B.
4) 1-D BY 1-D BLOCK(S) (*)
(a) If either SIZE(A, DA1) == 1 or SIZE(B, DB1) == 1, or both,
C = MULTIPROD(A, B, DA1, DB1) returns products of scalars by
vectors, or vectors by scalars or scalars by scalars.
(b) If SIZE(A, DA1) == SIZE(B, DB1),
C = MULTIPROD(A, B, [0 DA1], [DB1 0]) or
C = MULTIPROD(A, B, DA1, DB1) virtually turns the vectors
contained in A and B into 1xP and Px1 matrices, respectively,
then returns their products, similar to scalar products.
Namely, C = DOT2(A, B, DA1, DB1) is equivalent to
C = MULTIPROD(CONJ(A), B, [0 DA1], [DB1 0]).
(c) Without limitations on the length of the vectors in A and B,
C = MULTIPROD(A, B, [DA1 0], [0 DB1]) turns the vectors
contained in A and B into Px1 and 1xQ matrices, respectively,
then returns their products, similar to outer products.
Namely, C = OUTER(A, B, DA1, DB1) is equivalent to
C = MULTIPROD(CONJ(A), B, [DA1 0], [0 DB1]).
Common constraints for all syntaxes:
The external dimensions of A and B must either be identical or
compatible with AX rules. The internal dimensions of each block
array must be adjacent (DA2 == DA1 + 1 and DB2 == DB1 + 1 are
required). DA1 and DB1 are allowed to be larger than NDIMS(A) and
NDIMS(B). In syntaxes 1, 2, and 3, Q == R is required, unless the
blocks in A or B are scalars.
Array expansion (AX):
AX is a powerful generalization to N-D of the concept of scalar
expansion. Indeed, A and B may be scalars, vectors, matrices or
multi-dimensional arrays. Scalar expansion is the virtual
replication or annihilation of a scalar which allows you to combine
it, element by element, with an array X of any size (e.g. X+10,
X*10, or []-10). Similarly, in MULTIPROD, the purpose of AX is to
automatically match the size of the external dimensions (EDs) of A
and B, so that block-by-block products can be performed. ED matching
is achieved by means of a dimension shift followed by a singleton
expansion:
1) Dimension shift (see SHIFTDIM).
Whenever DA1 ~= DB1, a shift is applied to impose DA1 == DB1.
If DA1 > DB1, B is shifted to the right by DA1 - DB1 steps.
If DB1 > DA1, A is shifted to the right by DB1 - DA1 steps.
2) Singleton expansion (SX).
Whenever an ED of either A or B is singleton and the
corresponding ED of the other array is not, the mismatch is
fixed by virtually replicating the array (or diminishing it to
length 0) along that dimension.
MULTIPROD is a generalization for N-D arrays of the matrix
multiplication function MTIMES, with AX enabled. Vector inner, outer,
and cross products generalized for N-D arrays and with AX enabled are
performed by DOT2, OUTER, and CROSS2 (MATLAB Central, file #8782).
Elementwise multiplications (see TIMES) and other elementwise binary
operations with AX enabled are performed by BAXFUN (MATLAB Central,
file #23084). Together, these functions make up the "ARRAYLAB toolbox".
Input and output format:
The size of the EDs of C is determined by AX. Block size is
determined as follows, for each of the above-listed syntaxes:
1) C contains PxS matrices along IDs MAX([DA1 DA2], [DB1 DB2]).
2) Array Block size ID(s)
----------------------------------------------------
A PxQ (2-D) [DA1 DA2]
B R (1-D) DB1
C (a) P (1-D) MAX(DA1, DB1)
C (b) PxQ (2-D) MAX([DA1 DA2], [DB1 DB1+1])
----------------------------------------------------
(a) The 1-D blocks in B are not scalars (R > 1).
(b) The 1-D blocks in B are scalars (R = 1).
3) Array Block size ID(s)
----------------------------------------------------
A Q (1-D) DA1
B RxS (2-D) [DB1 DB2]
C (a) S (1-D) MAX(DA1, DB1)
C (b) RxS (2-D) MAX([DA1 DA1+1], [DB1 DB2])
----------------------------------------------------
(a) The 1-D blocks in A are not scalars (Q > 1).
(b) The 1-D blocks in A are scalars (Q = 1).
4) Array Block size ID(s)
--------------------------------------------------------------
(a) A P (1-D) DA1
B Q (1-D) DB1
C MAX(P,Q) (1-D) MAX(DA1, DB1)
--------------------------------------------------------------
(b) A P (1-D) DA1
B P (1-D) DB1
C 1 (1-D) MAX(DA1, DB1)
--------------------------------------------------------------
(c) A P (1-D) DA1
B Q (1-D) DB1
C PxQ (2-D) MAX([DA1 DA1+1], [DB1 DB1+1])
--------------------------------------------------------------
Terminological notes:
(*) 1-D and 2-D blocks are generically referred to as "vectors" and
"matrices", respectively. However, both may be also called
"scalars" if they have a single element. Moreover, matrices with a
single row or column (e.g. 1x3 or 3x1) may be also called "row
vectors" or "column vectors".
(°) Not to be confused with the "inner dimensions" of the two matrices
involved in a product X * Y, defined as the 2nd dimension of X and
the 1st of Y (DA2 and DB1 in syntaxes 1, 2, 3).
Examples:
1) If A is .................... a 5x(6x3)x2 array,
and B is .................... a 5x(3x4)x2 array,
C = MULTIPROD(A, B, [2 3]) is a 5x(6x4)x2 array.
A single matrix A pre-multiplies each matrix in B
If A is ........................... a (1x3) single matrix,
and B is ........................... a 10x(3x4) 3-D array,
C = MULTIPROD(A, B, [1 2], [3 4]) is a 10x(1x4) 3-D array.
Each matrix in A pre-multiplies each matrix in B (all possible
combinations)
If A is .................... a (6x3)x5 array,
and B is .................... a (3x4)x1x2 array,
C = MULTIPROD(A, B, [1 2]) is a (6x4)x5x2 array.
2a) If A is ........................... a 5x(6x3)x2 4-D array,
and B is ........................... a 5x(3)x2 3-D array,
C = MULTIPROD(A, B, [2 3], [2]) is a 5x(6)x2 3-D array.
2b) If A is ........................... a 5x(6x3)x2 4-D array,
and B is ........................... a 5x(1)x2 3-D array,
C = MULTIPROD(A, B, [2 3], [2]) is a 5x(6x3)x2 4-D array.
4a) If both A and B are .................. 5x(6)x2 3-D arrays,
C = MULTIPROD(A, B, 2) is .......... a 5x(1)x2 3-D array, while
4b) C = MULTIPROD(A, B, [2 0], [0 2]) is a 5x(6x6)x2 4-D array
See also MULTIPROD, DOT2, OUTER, CROSS2, BAXFUN, MULTITRANSP, MULTITRACE, MULTISCALE.0001 function c = multiprod_legacy(a, b, idA, idB) 0002 % Multiplying 1-D or 2-D subarrays contained in two N-D arrays. 0003 % 0004 % THIS ORIGINAL MULTIPROD IS NOW CALLED MULTIPROD_LEGACY IN MANOPT 0005 % 0006 % C = MULTIPROD(A,B) is equivalent to C = MULTIPROD(A,B,[1 2],[1 2]) 0007 % C = MULTIPROD(A,B,[D1 D2]) is eq. to C = MULTIPROD(A,B,[D1 D2],[D1 D2]) 0008 % C = MULTIPROD(A,B,D1) is equival. to C = MULTIPROD(A,B,D1,D1) 0009 % 0010 % MULTIPROD performs multiple matrix products, with array expansion (AX) 0011 % enabled. Its first two arguments A and B are "block arrays" of any 0012 % size, containing one or more 1-D or 2-D subarrays, called "blocks" (*). 0013 % For instance, a 5x6x3 array may be viewed as an array containing five 0014 % 6x3 blocks. In this case, its size is denoted by 5x(6x3). The 1 or 2 0015 % adjacent dimensions along which the blocks are contained are called the 0016 % "internal dimensions" (IDs) of the array (°). 0017 % 0018 % 1) 2-D by 2-D BLOCK(S) (*) 0019 % C = MULTIPROD(A, B, [DA1 DA2], [DB1 DB2]) contains the products 0020 % of the PxQ matrices in A by the RxS matrices in B. [DA1 DA2] are 0021 % the IDs of A; [DB1 DB2] are the IDs of B. 0022 % 0023 % 2) 2-D by 1-D BLOCK(S) (*) 0024 % C = MULTIPROD(A, B, [DA1 DA2], DB1) contains the products of the 0025 % PxQ matrices in A by the R-element vectors in B. The latter are 0026 % considered to be Rx1 matrices. [DA1 DA2] are the IDs of A; DB1 is 0027 % the ID of B. 0028 % 0029 % 3) 1-D by 2-D BLOCK(S) (*) 0030 % C = MULTIPROD(A, B, DA1, [DB1 DB2]) contains the products of the 0031 % Q-element vectors in A by the RxS matrices in B. The vectors in A 0032 % are considered to be 1xQ matrices. DA1 is the ID of A; [DB1 DB2] 0033 % are the IDs of B. 0034 % 0035 % 4) 1-D BY 1-D BLOCK(S) (*) 0036 % (a) If either SIZE(A, DA1) == 1 or SIZE(B, DB1) == 1, or both, 0037 % C = MULTIPROD(A, B, DA1, DB1) returns products of scalars by 0038 % vectors, or vectors by scalars or scalars by scalars. 0039 % (b) If SIZE(A, DA1) == SIZE(B, DB1), 0040 % C = MULTIPROD(A, B, [0 DA1], [DB1 0]) or 0041 % C = MULTIPROD(A, B, DA1, DB1) virtually turns the vectors 0042 % contained in A and B into 1xP and Px1 matrices, respectively, 0043 % then returns their products, similar to scalar products. 0044 % Namely, C = DOT2(A, B, DA1, DB1) is equivalent to 0045 % C = MULTIPROD(CONJ(A), B, [0 DA1], [DB1 0]). 0046 % (c) Without limitations on the length of the vectors in A and B, 0047 % C = MULTIPROD(A, B, [DA1 0], [0 DB1]) turns the vectors 0048 % contained in A and B into Px1 and 1xQ matrices, respectively, 0049 % then returns their products, similar to outer products. 0050 % Namely, C = OUTER(A, B, DA1, DB1) is equivalent to 0051 % C = MULTIPROD(CONJ(A), B, [DA1 0], [0 DB1]). 0052 % 0053 % Common constraints for all syntaxes: 0054 % The external dimensions of A and B must either be identical or 0055 % compatible with AX rules. The internal dimensions of each block 0056 % array must be adjacent (DA2 == DA1 + 1 and DB2 == DB1 + 1 are 0057 % required). DA1 and DB1 are allowed to be larger than NDIMS(A) and 0058 % NDIMS(B). In syntaxes 1, 2, and 3, Q == R is required, unless the 0059 % blocks in A or B are scalars. 0060 % 0061 % Array expansion (AX): 0062 % AX is a powerful generalization to N-D of the concept of scalar 0063 % expansion. Indeed, A and B may be scalars, vectors, matrices or 0064 % multi-dimensional arrays. Scalar expansion is the virtual 0065 % replication or annihilation of a scalar which allows you to combine 0066 % it, element by element, with an array X of any size (e.g. X+10, 0067 % X*10, or []-10). Similarly, in MULTIPROD, the purpose of AX is to 0068 % automatically match the size of the external dimensions (EDs) of A 0069 % and B, so that block-by-block products can be performed. ED matching 0070 % is achieved by means of a dimension shift followed by a singleton 0071 % expansion: 0072 % 1) Dimension shift (see SHIFTDIM). 0073 % Whenever DA1 ~= DB1, a shift is applied to impose DA1 == DB1. 0074 % If DA1 > DB1, B is shifted to the right by DA1 - DB1 steps. 0075 % If DB1 > DA1, A is shifted to the right by DB1 - DA1 steps. 0076 % 2) Singleton expansion (SX). 0077 % Whenever an ED of either A or B is singleton and the 0078 % corresponding ED of the other array is not, the mismatch is 0079 % fixed by virtually replicating the array (or diminishing it to 0080 % length 0) along that dimension. 0081 % 0082 % MULTIPROD is a generalization for N-D arrays of the matrix 0083 % multiplication function MTIMES, with AX enabled. Vector inner, outer, 0084 % and cross products generalized for N-D arrays and with AX enabled are 0085 % performed by DOT2, OUTER, and CROSS2 (MATLAB Central, file #8782). 0086 % Elementwise multiplications (see TIMES) and other elementwise binary 0087 % operations with AX enabled are performed by BAXFUN (MATLAB Central, 0088 % file #23084). Together, these functions make up the "ARRAYLAB toolbox". 0089 % 0090 % Input and output format: 0091 % The size of the EDs of C is determined by AX. Block size is 0092 % determined as follows, for each of the above-listed syntaxes: 0093 % 1) C contains PxS matrices along IDs MAX([DA1 DA2], [DB1 DB2]). 0094 % 2) Array Block size ID(s) 0095 % ---------------------------------------------------- 0096 % A PxQ (2-D) [DA1 DA2] 0097 % B R (1-D) DB1 0098 % C (a) P (1-D) MAX(DA1, DB1) 0099 % C (b) PxQ (2-D) MAX([DA1 DA2], [DB1 DB1+1]) 0100 % ---------------------------------------------------- 0101 % (a) The 1-D blocks in B are not scalars (R > 1). 0102 % (b) The 1-D blocks in B are scalars (R = 1). 0103 % 3) Array Block size ID(s) 0104 % ---------------------------------------------------- 0105 % A Q (1-D) DA1 0106 % B RxS (2-D) [DB1 DB2] 0107 % C (a) S (1-D) MAX(DA1, DB1) 0108 % C (b) RxS (2-D) MAX([DA1 DA1+1], [DB1 DB2]) 0109 % ---------------------------------------------------- 0110 % (a) The 1-D blocks in A are not scalars (Q > 1). 0111 % (b) The 1-D blocks in A are scalars (Q = 1). 0112 % 4) Array Block size ID(s) 0113 % -------------------------------------------------------------- 0114 % (a) A P (1-D) DA1 0115 % B Q (1-D) DB1 0116 % C MAX(P,Q) (1-D) MAX(DA1, DB1) 0117 % -------------------------------------------------------------- 0118 % (b) A P (1-D) DA1 0119 % B P (1-D) DB1 0120 % C 1 (1-D) MAX(DA1, DB1) 0121 % -------------------------------------------------------------- 0122 % (c) A P (1-D) DA1 0123 % B Q (1-D) DB1 0124 % C PxQ (2-D) MAX([DA1 DA1+1], [DB1 DB1+1]) 0125 % -------------------------------------------------------------- 0126 % 0127 % Terminological notes: 0128 % (*) 1-D and 2-D blocks are generically referred to as "vectors" and 0129 % "matrices", respectively. However, both may be also called 0130 % "scalars" if they have a single element. Moreover, matrices with a 0131 % single row or column (e.g. 1x3 or 3x1) may be also called "row 0132 % vectors" or "column vectors". 0133 % (°) Not to be confused with the "inner dimensions" of the two matrices 0134 % involved in a product X * Y, defined as the 2nd dimension of X and 0135 % the 1st of Y (DA2 and DB1 in syntaxes 1, 2, 3). 0136 % 0137 % Examples: 0138 % 1) If A is .................... a 5x(6x3)x2 array, 0139 % and B is .................... a 5x(3x4)x2 array, 0140 % C = MULTIPROD(A, B, [2 3]) is a 5x(6x4)x2 array. 0141 % 0142 % A single matrix A pre-multiplies each matrix in B 0143 % If A is ........................... a (1x3) single matrix, 0144 % and B is ........................... a 10x(3x4) 3-D array, 0145 % C = MULTIPROD(A, B, [1 2], [3 4]) is a 10x(1x4) 3-D array. 0146 % 0147 % Each matrix in A pre-multiplies each matrix in B (all possible 0148 % combinations) 0149 % If A is .................... a (6x3)x5 array, 0150 % and B is .................... a (3x4)x1x2 array, 0151 % C = MULTIPROD(A, B, [1 2]) is a (6x4)x5x2 array. 0152 % 0153 % 2a) If A is ........................... a 5x(6x3)x2 4-D array, 0154 % and B is ........................... a 5x(3)x2 3-D array, 0155 % C = MULTIPROD(A, B, [2 3], [2]) is a 5x(6)x2 3-D array. 0156 % 0157 % 2b) If A is ........................... a 5x(6x3)x2 4-D array, 0158 % and B is ........................... a 5x(1)x2 3-D array, 0159 % C = MULTIPROD(A, B, [2 3], [2]) is a 5x(6x3)x2 4-D array. 0160 % 0161 % 4a) If both A and B are .................. 5x(6)x2 3-D arrays, 0162 % C = MULTIPROD(A, B, 2) is .......... a 5x(1)x2 3-D array, while 0163 % 4b) C = MULTIPROD(A, B, [2 0], [0 2]) is a 5x(6x6)x2 4-D array 0164 % 0165 % See also MULTIPROD, DOT2, OUTER, CROSS2, BAXFUN, MULTITRANSP, MULTITRACE, MULTISCALE. 0166 0167 % $ Version: 2.1 $ 0168 % CODE by: Paolo de Leva 0169 % (Univ. of Rome, Foro Italico, IT) 2009 Jan 24 0170 % optimized by: Paolo de Leva 0171 % Jinhui Bai (Georgetown Univ., D.C.) 2009 Jan 24 0172 % COMMENTS by: Paolo de Leva 2009 Feb 24 0173 % OUTPUT tested by: Paolo de Leva 2009 Feb 24 0174 % ------------------------------------------------------------------------- 0175 0176 assert(nargin >= 2 && nargin <= 4, 'Takes from 2 to 4 inputs.'); 0177 0178 switch nargin % Setting IDA and/or IDB 0179 case 2, idA = [1 2]; idB = [1 2]; 0180 case 3, idB = idA; 0181 end 0182 0183 % ESC 1 - Special simple case (both A and B are 2D), solved using C = A * B 0184 0185 if ismatrix(a) && ismatrix(b) && ... 0186 isequal(idA,[1 2]) && isequal(idB,[1 2]) 0187 c = a * b; return 0188 end 0189 0190 % MAIN 0 - Checking and evaluating array size, block size, and IDs 0191 0192 sizeA0 = size(a); 0193 sizeB0 = size(b); 0194 [sizeA, sizeB, shiftC, delC, sizeisnew, idA, idB, ... 0195 squashOK, sxtimesOK, timesOK, mtimesOK, sumOK] = ... 0196 sizeval(idA,idB, sizeA0,sizeB0); 0197 0198 % MAIN 1 - Applying dimension shift (first step of AX) and 0199 % turning both A and B into arrays of either 1-D or 2-D blocks 0200 0201 if sizeisnew(1), a = reshape(a, sizeA); end 0202 if sizeisnew(2), b = reshape(b, sizeB); end 0203 0204 % MAIN 2 - Performing products with or without SX (second step of AX) 0205 0206 if squashOK % SQUASH + MTIMES (fastest engine) 0207 c = squash2D_mtimes(a,b, idA,idB, sizeA,sizeB, squashOK); 0208 elseif timesOK % TIMES (preferred w.r. to SX + TIMES) 0209 if sumOK, c = sum(a .* b, sumOK); 0210 else, c = a .* b; end 0211 elseif sxtimesOK % SX + TIMES 0212 if sumOK, c = sum(bsxfun(@times, a, b), sumOK); 0213 else, c = bsxfun(@times, a, b); end 0214 elseif mtimesOK % MTIMES (rarely used) 0215 c = a * b; 0216 end 0217 0218 % MAIN 3 - Reshaping C (by inserting or removing singleton dimensions) 0219 0220 [sizeC, sizeCisnew] = adjustsize(size(c), shiftC, false, delC, false); 0221 if sizeCisnew, c = reshape(c, sizeC); end 0222 0223 0224 function c = squash2D_mtimes(a, b, idA, idB, sizeA, sizeB, squashOK) 0225 % SQUASH2D_MTIMES Multiproduct with single-block expansion (SBX). 0226 % Actually, no expansion is performed. The multi-block array is 0227 % rearranged from N-D to 2-D, then MTIMES is applied, and eventually the 0228 % result is rearranged back to N-D. No additional memory is required. 0229 % One and only one of the two arrays must be single-block, and its IDs 0230 % must be [1 2] (MAIN 1 removes leading singletons). Both arrays 0231 % must contain 2-D blocks (MAIN 1 expands 1-D blocks to 2-D). 0232 0233 if squashOK == 1 % A is multi-block, B is single-block (squashing A) 0234 0235 % STEP 1 - Moving IDA(2) to last dimension 0236 nd = length(sizeA); 0237 d2 = idA(2); 0238 order = [1:(d2-1) (d2+1):nd d2]; % Partial shifting 0239 a = permute(a, order); % ...xQ 0240 0241 % STEP 2 - Squashing A from N-D to 2-D 0242 q = sizeB(1); 0243 s = sizeB(2); 0244 lengthorder = length(order); 0245 collapsedsize = sizeA(order(1:lengthorder-1)); 0246 n = prod(collapsedsize); 0247 a = reshape(a, [n, q]); % NxQ 0248 fullsize = [collapsedsize s]; % Size to reshape C back to N-D 0249 0250 else % B is multi-block, A is single-block (squashing B) 0251 0252 % STEP 1 - Moving IDB(1) to first dimension 0253 nd = length(sizeB); 0254 d1 = idB(1); 0255 order = [d1 1:(d1-1) (d1+1):nd]; % Partial shifting 0256 b = permute(b, order); % Qx... 0257 0258 % STEP 2 - Squashing B from N-D to 2-D 0259 p = sizeA(1); 0260 q = sizeA(2); 0261 lengthorder = length(order); 0262 collapsedsize = sizeB(order(2:lengthorder)); 0263 n = prod(collapsedsize); 0264 b = reshape(b, [q, n]); % QxN 0265 fullsize = [p collapsedsize]; % Size to reshape C back to N-D 0266 0267 end 0268 0269 % FINAL STEPS - Multiplication, reshape to N-D, inverse permutation 0270 invorder(order) = 1 : lengthorder; 0271 c = permute (reshape(a*b, fullsize), invorder); 0272 0273 0274 function [sizeA, sizeB, shiftC, delC, sizeisnew, idA, idB, ... 0275 squashOK, sxtimesOK, timesOK, mtimesOK, sumOK] = ... 0276 sizeval(idA0,idB0, sizeA0,sizeB0) 0277 %SIZEVAL Evaluation of array size, block size, and IDs 0278 % Possible values for IDA and IDB: 0279 % [DA1 DA2], [DB1 DB2] 0280 % [DA1 DA2], [DB1] 0281 % [DA1], [DB1 DB2] 0282 % [DA1], [DB1] 0283 % [DA1 0], [0 DB1] 0284 % [0 DA1], [DB1 0] 0285 % 0286 % sizeA/B Equal to sizeA0/B0 if RESHAPE is not needed in MAIN 1 0287 % shiftC, delC Variables controlling MAIN 3. 0288 % sizeisnew 1x2 logical array; activates reshaping of A and B. 0289 % idA/B May change only if squashOK ~= 0 0290 % squashOK If only A or B is a multi-block array (M-B) and the other 0291 % is single-block (1-B), it will be rearranged from N-D to 0292 % 2-D. If both A and B are 1-B or M-B arrays, squashOK = 0. 0293 % If only A (or B) is a M-B array, squashOK = 1 (or 2). 0294 % sxtimesOK, timesOK, mtimesOK Flags controlling MAIN 2 (TRUE/FALSE). 0295 % sumOK Dimension along which SUM is performed. If SUM is not 0296 % needed, sumOK = 0. 0297 0298 % Initializing output arguments 0299 0300 idA = idA0; 0301 idB = idB0; 0302 squashOK = 0; 0303 sxtimesOK = false; 0304 timesOK = false; 0305 mtimesOK = false; 0306 sumOK = 0; 0307 shiftC = 0; 0308 delC = 0; 0309 0310 % Checking for gross input errors 0311 0312 NidA = numel(idA); 0313 NidB = numel(idB); 0314 idA1 = idA(1); 0315 idB1 = idB(1); 0316 if NidA>2 || NidB>2 || NidA==0 || NidB==0 || ... 0317 ~isreal(idA1) || ~isreal(idB1) || ... 0318 ~isnumeric(idA1) || ~isnumeric(idB1) || ... 0319 0>idA1 || 0>idB1 || ... % negative 0320 idA1~=fix(idA1) || idB1~=fix(idB1) || ... % non-integer 0321 ~isfinite(idA1) || ~isfinite(idB1) % Inf or NaN 0322 error('MULTIPROD:InvalidDimensionArgument', ... 0323 ['Internal-dimension arguments (e.g., [IDA1 IDA2]) must\n', ... 0324 'contain only one or two non-negative finite integers']); 0325 end 0326 0327 % Checking Syntaxes containing zeros (4b/c) 0328 0329 declared_outer = false; 0330 idA2 = idA(NidA); % It may be IDA1 = IDA2 (1-D block) 0331 idB2 = idB(NidB); 0332 0333 if any(idA==0) || any(idB==0) 0334 0335 % "Inner products": C = MULTIPROD(A, B, [0 DA1], [DB1 0]) 0336 if idA1==0 && idA2>0 && idB1>0 && idB2==0 0337 idA1 = idA2; 0338 idB2 = idB1; 0339 % "Outer products": C = MULTIPROD(A, B, [DA1 0], [0 DB1]) 0340 elseif idA1>0 && idA2==0 && idB1==0 && idB2>0 0341 declared_outer = true; 0342 idA2 = idA1; 0343 idB1 = idB2; 0344 else 0345 error('MULTIPROD:InvalidDimensionArgument', ... 0346 ['Misused zeros in the internal-dimension arguments\n', ... 0347 '(see help heads 4b and 4c)']); 0348 end 0349 NidA = 1; 0350 NidB = 1; 0351 idA = idA1; 0352 idB = idB1; 0353 0354 elseif (NidA==2 && idA2~=idA1+1) || ... % Non-adjacent IDs 0355 (NidB==2 && idB2~=idB1+1) 0356 error('MULTIPROD:InvalidDimensionArgument', ... 0357 ['If an array contains 2-D blocks, its two internal dimensions', ... 0358 'must be adjacent (e.g. IDA2 == IDA1+1)']); 0359 end 0360 0361 % ESC - Case for which no reshaping is needed (both A and B are scalars) 0362 0363 scalarA = isequal(sizeA0, [1 1]); 0364 scalarB = isequal(sizeB0, [1 1]); 0365 if scalarA && scalarB 0366 sizeA = sizeA0; 0367 sizeB = sizeB0; 0368 sizeisnew = [false false]; 0369 timesOK = true; return 0370 end 0371 0372 % Computing and checking adjusted sizes 0373 % The lengths of ADJSIZEA and ADJSIZEB must be >= IDA(END) and IDB(END) 0374 0375 NsA = idA2 - length(sizeA0); % Number of added trailing singletons 0376 NsB = idB2 - length(sizeB0); 0377 adjsizeA = [sizeA0 ones(1,NsA)]; 0378 adjsizeB = [sizeB0 ones(1,NsB)]; 0379 extsizeA = adjsizeA([1:idA1-1, idA2+1:end]); % Size of EDs 0380 extsizeB = adjsizeB([1:idB1-1, idB2+1:end]); 0381 p = adjsizeA(idA1); 0382 q = adjsizeA(idA2); 0383 r = adjsizeB(idB1); 0384 s = adjsizeB(idB2); 0385 scalarsinA = (p==1 && q==1); 0386 scalarsinB = (r==1 && s==1); 0387 singleA = all(extsizeA==1); 0388 singleB = all(extsizeB==1); 0389 if q~=r && ~scalarsinA && ~scalarsinB && ~declared_outer 0390 error('MULTIPROD:InnerDimensionsMismatch', ... 0391 'Inner matrix dimensions must agree.'); 0392 end 0393 0394 % STEP 1/3 - DIMENSION SHIFTING (FIRST STEP OF AX) 0395 % Pipeline 1 (using TIMES) never needs left, and may need right shifting. 0396 % Pipeline 2 (using MTIMES) may need left shifting of A and right of B. 0397 0398 shiftA = 0; 0399 shiftB = 0; 0400 diffBA = idB1 - idA1; 0401 if scalarA % Do nothing 0402 elseif singleA && ~scalarsinB, shiftA = -idA1 + 1; % Left shifting A 0403 elseif idB1 > idA1, shiftA = diffBA; % Right shifting A 0404 end 0405 if scalarB % Do nothing 0406 elseif singleB && ~scalarsinA, shiftB = -idB1 + 1; % Left shifting B 0407 elseif idA1 > idB1, shiftB = -diffBA; % Right shifting B 0408 end 0409 0410 % STEP 2/3 - SELECTION OF PROPER ENGINE AND BLOCK SIZE ADJUSTMENTS 0411 0412 addA = 0; addB = 0; 0413 delA = 0; delB = 0; 0414 swapA = 0; swapB = 0; 0415 idC1 = max(idA1, idB1); 0416 idC2 = idC1 + 1; 0417 checktimes = false; 0418 0419 if (singleA||singleB) &&~scalarsinA &&~scalarsinB % Engine using MTIMES 0420 0421 if singleA && singleB 0422 mtimesOK = true; 0423 shiftC=idC1-1; % Right shifting C 0424 idC1=1; idC2=2; 0425 elseif singleA 0426 squashOK = 2; 0427 idB = [idB1, idB1+1] + shiftB; 0428 else % singleB 0429 squashOK = 1; 0430 idA = [idA1, idA1+1] + shiftA; 0431 end 0432 0433 if NidA==2 && NidB==2 % 1) 2-D BLOCKS BY 2-D BLOCKS 0434 % OK 0435 elseif NidA==2 % 2) 2-D BLOCKS BY 1-D BLOCKS 0436 addB=idB1+1; delC=idC2; 0437 elseif NidB==2 % 3) 1-D BLOCKS BY 2-D BLOCKS 0438 addA=idA1; delC=idC1; 0439 else % 4) 1-D BLOCKS BY 1-D BLOCKS 0440 if declared_outer 0441 addA=idA1+1; addB=idB1; 0442 else 0443 addA=idA1; addB=idB1+1; delC=idC2; 0444 end 0445 end 0446 0447 else % Engine using TIMES (also used if SCALARA || SCALARB) 0448 0449 sxtimesOK = true; 0450 0451 if NidA==2 && NidB==2 % 1) 2-D BLOCKS BY 2-D BLOCKS 0452 0453 if scalarA || scalarB 0454 timesOK=true; 0455 elseif scalarsinA && scalarsinB % scal-by-scal 0456 checktimes=true; 0457 elseif scalarsinA || scalarsinB || ... % scal-by-mat 0458 (q==1 && r==1) % vec-by-vec ("outer") 0459 elseif p==1 && s==1 % vec-by-vec ("inner") 0460 swapA=idA1; sumOK=idC1; checktimes=true; 0461 elseif s==1 % mat-by-vec 0462 swapB=idB1; sumOK=idC2; 0463 elseif p==1 % vec-by-mat 0464 swapA=idA1; sumOK=idC1; 0465 else % mat-by-mat 0466 addA=idA2+1; addB=idB1; sumOK=idC2; delC=idC2; 0467 end 0468 0469 elseif NidA==2 % 2) 2-D BLOCKS BY 1-D BLOCKS 0470 0471 if scalarA || scalarB 0472 timesOK=true; 0473 elseif scalarsinA && scalarsinB % scal-by-scal 0474 addB=idB1; checktimes=true; 0475 elseif scalarsinA % scal-by-vec 0476 delA=idA1; 0477 elseif scalarsinB % mat-by-scal 0478 addB=idB1; 0479 elseif p==1 % vec-by-vec ("inner") 0480 delA=idA1; sumOK=idC1; checktimes=true; 0481 else % mat-by-vec 0482 addB=idB1; sumOK=idC2; delC=idC2; 0483 end 0484 0485 elseif NidB==2 % 3) 1-D BLOCKS BY 2-D BLOCKS 0486 0487 if scalarA || scalarB 0488 timesOK=true; 0489 elseif scalarsinA && scalarsinB % scal-by-scal 0490 addA=idA1+1; checktimes=true; 0491 elseif scalarsinB % vec-by-scal 0492 delB=idB2; 0493 elseif scalarsinA % scal-by-mat 0494 addA=idA1+1; 0495 elseif s==1 % vec-by-vec ("inner") 0496 delB=idB2; sumOK=idC1; checktimes=true; 0497 else % vec-by-mat 0498 addA=idA1+1; sumOK=idC1; delC=idC1; 0499 end 0500 0501 else % 4) 1-D BLOCKS BY 1-D BLOCKS 0502 0503 if scalarA || scalarB 0504 timesOK=true; 0505 elseif declared_outer % vec-by-vec ("outer") 0506 addA=idA1+1; addB=idB1; 0507 elseif scalarsinA && scalarsinB % scal-by-scal 0508 checktimes=true; 0509 elseif scalarsinA || scalarsinB % vec-by-scal 0510 else % vec-by-vec 0511 sumOK=idC1; checktimes=true; 0512 end 0513 end 0514 end 0515 0516 % STEP 3/3 - Adjusting the size of A and B. The size of C is adjusted 0517 % later, because it is not known yet. 0518 0519 [sizeA, sizeisnew(1)] = adjustsize(sizeA0, shiftA, addA, delA, swapA); 0520 [sizeB, sizeisnew(2)] = adjustsize(sizeB0, shiftB, addB, delB, swapB); 0521 0522 if checktimes % Faster than calling BBXFUN 0523 diff = length(sizeB) - length(sizeA); 0524 if isequal([sizeA ones(1,diff)], [sizeB ones(1,-diff)]) 0525 timesOK = true; 0526 end 0527 end 0528 0529 0530 function [sizeA, sizeisnew] = adjustsize(sizeA0, shiftA, addA, delA, swapA) 0531 % ADJUSTSIZE Adjusting size of a block array. 0532 0533 % Dimension shifting (by adding or deleting trailing singleton dim.) 0534 if shiftA>0, [sizeA,newA1] = addsing(sizeA0, 1, shiftA); 0535 elseif shiftA<0, [sizeA,newA1] = delsing(sizeA0, 1,-shiftA); 0536 else, sizeA = sizeA0; newA1 = false; 0537 end 0538 % Modifying block size (by adding, deleting, or moving singleton dim.) 0539 if addA, [sizeA,newA2] = addsing(sizeA, addA+shiftA, 1); % 1D-->2D 0540 elseif delA, [sizeA,newA2] = delsing(sizeA, delA+shiftA, 1); % 2D-->1D 0541 elseif swapA, [sizeA,newA2] = swapdim(sizeA,swapA+shiftA); % ID Swapping 0542 else, newA2 = false; 0543 end 0544 sizeisnew = newA1 || newA2; 0545 0546 0547 function [newsize, flag] = addsing(size0, dim, ns) 0548 %ADDSING Adding NS singleton dimensions to the size of an array. 0549 % Warning: NS is assumed to be a positive integer. 0550 % Example: If the size of A is ..... SIZE0 = [5 9 3] 0551 % NEWSIZE = ADDSING(SIZE0, 3, 2) is [5 9 1 1 3] 0552 0553 if dim > length(size0) 0554 newsize = size0; 0555 flag = false; 0556 else 0557 newsize = [size0(1:dim-1), ones(1,ns), size0(dim:end)]; 0558 flag = true; 0559 end 0560 0561 0562 function [newsize, flag] = delsing(size0, dim, ns) 0563 %DELSING Removing NS singleton dimensions from the size of an array. 0564 % Warning: Trailing singletons are not removed 0565 % Example: If the size of A is SIZE0 = [1 1 1 5 9 3] 0566 % NEWSIZE = DELSING(SIZE, 1, 3) is [5 9 3] 0567 0568 if dim > length(size0)-ns % Trailing singletons are not removed 0569 newsize = size0; 0570 flag = false; 0571 else % Trailing singl. added, so NEWSIZE is guaranteed to be 2D or more 0572 newsize = size0([1:dim-1, dim+ns:end, dim]); 0573 flag = true; 0574 end 0575 0576 0577 function [newsize, flag] = swapdim(size0, dim) 0578 %SWAPDIM Swapping two adjacent dimensions of an array (DIM and DIM+1). 0579 % Used only when both A and B are multi-block arrays with 2-D blocks. 0580 % Example: If the size of A is .......... 5x(6x3) 0581 % NEWSIZE = SWAPIDS(SIZE0, 2) is 5x(3x6) 0582 0583 newsize = [size0 1]; % Guarantees that dimension DIM+1 exists. 0584 newsize = newsize([1:dim-1, dim+1, dim, dim+2:end]); 0585 flag = true;