%%************************************************************
%% Prod3: compute the entries of Q = A*B*C specified in 
%%        nzlistQ. 
%% 
%% Q = Prod3(blk,A,B,C,sym,nzlistQ)
%% Important: (a) A is assumed to be symmetric if nzlistQ
%%                has 2 columns (since mexProd2nz computes A'*B). 
%%            (b) The 2nd column of nzlistQ must be sorted in 
%%                ascending order. 
%%
%% (optional) sym = 1, if Q is symmetric.
%%                = 0, otherwise. 
%% (optional) nzlistQ = list of non-zero elements of Q to be 
%%                      computed.
%%
%%************************************************************

  function Q = Prod3nz(blk,A,B,C,sym,nzlistQ)

  if (nargin<5); sym = 0; end;
  checkcell(1:3) = [iscell(A) iscell(B) iscell(C)]; 
  if (nargin==6);
     checkcell(4) = iscell(nzlistQ); 
  else 
     nzlistQ = inf; 
  end;
%%
  if any(checkcell-1)      
     if (size(blk,1) > 1); 
        error('Prod3nz: blk and A,B,C are not compatible'); 
     end;
     if strcmp(blk{1},'s')
        [len,len2] = size(nzlistQ);  
        if (len2 == 1) & (nzlistQ == inf) 
           tmp = Prod2(blk,A,B,0);   
           Q = Prod2(blk,tmp,C,sym); 
        else
           tmp = Prod2(blk,B,C,0);
           Q = mexProd2nz(blk,A,tmp,nzlistQ); 
           if sym; Q = 0.5*(Q+Q'); end; 
        end
     end
     if strcmp(blk{1},'l')
        Q = A.*B.*C;
     end
  else 
     error('Prod3nz: A,B,C,nzlistQ must all be matrices'); 
  end
%%
%%************************************************************
