//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [C] = %hm_m_hm(A,B)
// Multiplication of 2 hypermatrices. 
// A is considered as a set of matrices : A = [A1, A2, ...Ak] 
// size(A) = (N,P,K) 
// B is considered as a set of matrices : B = [B1, A2, ...Bk]
// size(B) =(P,Q,K) 
// Then C is defined as : 
// C = [A1*B1, A2*B2, ... Ak*Bk]
// size(C) =(N,Q,K) 
// NOTA: %hm_m_hm has no default definition in Scilab


// Declarations:

// Code:

if (size(A.dims,'*') ~= 3 | size(B.dims,'*') ~= 3)
  CL__error("Only hypermatrices with 3 dimensions are handled");
end

[N,P,K] = size(A);
[P2,Q,K2] = size(B);

if (K ~= K2 | P ~= P2)
  CL__error("Hypermatrices have incompatible sizes");
end


C = hypermat([N, Q, K]); 

// Note : 
// we note Aik = A(i,:,k)(:), and Bjk = B(:,j,k)(:)
// (Aik and Bjk = column vectors of P elements)
// => A(i,:,:) contains IN THIS ORDER: Ai1; Ai2;... Aik
// which can be reformatted COLUMN BY COLUMN: 
// matrix(A(i,:,:),P,K) = [Ai1, Ai2... ]  -> size = PxK
// => B(:,j,:) contains IN THIS ORDER: Bj1; Bj2;... Bjk
// which can be reformatted COLUMN BY COLUMN:
// matrix(B(:,j,:),P,K) = [Bj1, Bj2,... ] -> size = PxK
// => easy to compute sum(Ai1 .* Bj1) = C(i,j,1), etc... 

for (i=1:N)
  for (j=1:Q)
    C(i,j,:) = sum(matrix(A(i,:,:),P,K) .* matrix(B(:,j,:),P,K), 'r');
  end
end


endfunction