Home > manopt > tools > multiprod_legacy.m

# multiprod_legacy

## PURPOSE

Multiplying 1-D or 2-D subarrays contained in two N-D arrays.

## SYNOPSIS

function c = multiprod_legacy(a, b, idA, idB)

## DESCRIPTION

``` 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

## CROSS-REFERENCE INFORMATION

This function calls:
This function is called by:
• multiprod Matrix multiply 2-D slices of N-D arrays

## SOURCE CODE

```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);
0379     extsizeA = adjsizeA([1:idA1-1, idA2+1:end]); % Size of EDs
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
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
0437         elseif NidB==2        % 3) 1-D BLOCKS BY 2-D BLOCKS
0439         else                  % 4) 1-D BLOCKS BY 1-D BLOCKS
0440             if declared_outer
0442             else
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
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
0475             elseif scalarsinA % scal-by-vec
0476                 delA=idA1;
0477             elseif scalarsinB % mat-by-scal
0479             elseif p==1 % vec-by-vec ("inner")
0480                 delA=idA1; sumOK=idC1; checktimes=true;
0481             else % mat-by-vec
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
0491             elseif scalarsinB % vec-by-scal
0492                 delB=idB2;
0493             elseif scalarsinA % scal-by-mat
0495             elseif s==1 % vec-by-vec ("inner")
0496                 delB=idB2; sumOK=idC1; checktimes=true;
0497             else % vec-by-mat
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")
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
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
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.)
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;```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005