% new vech
% this vectorization is constructed on the following fact
% define z=vec(M); then if need to compute y=vech(M), notice that among the
% elements of z, y selects the elements 1,
% m+1:m+2,2*m+1:2*m+3,...,(j-1)*m+1:(j-1)*m+j,...,(n-1)*m+1:(n-1)*m+n
% usage: [vv,Dn]=vech3(size of M);
% Giovanni.Lombardo@ecb.int 2005
% vech3 returns the indeces of the vec that are taken into the vech.
% Optionally it returns the Dn matrix that transforms vech back into vec

function [elvh,Dn]=vech3(m);

% if nargin<2
%     yesym=0;
% end
elvh=0;
Dn=0;


%m=size(M,1);
%s=size(M,2);
n=m*(m+1)/2;

      
  
%if m~=s
 %  disp(' matrix must be square')
 %  return
%end
nvec=m*m;
elvh=[];
% do the following separately from D so if D not requested can return
for jj=1:m
    elvh=[elvh;[1+(jj-1)*m:1+(jj-1)*(m+1)]'];
end
%conversion matrix
if nargout>1
 Dn=(zeros(m*m,n));   
z=1;
u=1;
for jj=1:m
%    xxv=[xxv;[(jj-1)*m+1:(jj-1)*m+jj]'];
    for ii=1:jj
        Dn((jj-1)*m+ii,z)=1;
        Dn((ii-1)*m+jj,z)=1;
        z=z+1;
    end
end

% if yesym==1
% Dn=sym(Dna);
% else
%    Dn=Dna;
%end
end %ifnargout
return