Home > manopt > tools > multiprod.m

multiprod

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

 Multiplying 1-D or 2-D subarrays contained in two N-D arrays.
 
   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 DOT2, OUTER, CROSS2, BAXFUN, MULTITRANSP, MULTITRACE, MULTISCALE.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005