function [indsrows,pcols,pcolsinv,Mfinal,Mtemp,pinvMfinal,pinvMtemp,...
           Mtempnull,finalblock]=matrixrep(m,n,A)
%function [indsrows,pcols,pcolsinv,Mfinal,Mtemp,pinvMfinal,pinvMtemp,...
%           Mtempnull,finalblock]=matrixrep(m,n,A)
%% find the matrix representation for completely positive T constraints
%%  T maps Hn ---> Hm  with constraints T(Ai)=Bi for all i
%% This always includes the trace preserving constraint, TP.
%% The unital constraint is assume to be one of the T(Ai)=Bi
%INPUT:
%     m,n     order of Hermitian inputs
%     A cell of size k for input matrices (size is k+1 as unital
%             constraint is added in run.m file.)
%OUTPUT 
%   indsrows  row indices for the permuted rows of L2 and RHS
%   colindices column indices/permutation for block representation
%   colindices inverse column indices/permutation for block representation
%   Lsparse       from original matrix representation, sparse, with permuted
%             rows without redundant rows, but original columns!
%             permuted columns are needed in altproj.m
%  pinvL    gen inverse with both permuted rows and permuted cols
%         

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Mar 22/14 ... code can be sped up by vectorizing the for kk= ...  loop
%%%      using a kronecker product ... some details are in previous versions
%i=sqrt(-1);     % better to use   1i
%%Apr 5/14 ... completed the sparse representation of the input
%%             i.e. still looping for each column but just saving
%               nonzeros as a sparse representation
%% Apr10/14: start implementing better block structure
%% Apr22/14: proper block structure now works ... the explict expression
%       for the entire Moore-Penrose inverse is given below - but the block
%       structure is exploited in the evaluation and so neither Lsparse, pinvL
%        are constructed or passed to other files.
% R is the RHS of the constraints used in altproj.m
% The other global variables are evaluated/described below.
% jun21/14   global variables are now passed

%global indsrows 
%global pcols pcolsinv Mfinal Mtemp pinvMfinal
%global pinvMtemp Mtempnull finalblock
sqrt2=sqrt(2);
k=size(A,1);  % number of Ai,Bi constraints
mn=m*n;
tmnm=mn*(mn-1)/2;
tmm=m*(m-1)/2;  % m=n so tmm=tnm throught
tnm=n*(n-1)/2;
M=zeros(k,m^2);
for ii=1:k,
	M(ii,:)=sHvectwo(A{ii})';    % do NOT use sqrt(2)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%testing
%%main diagonal block matrix rep.
%% rows are ordered by elements - not by B_i
M1=M(:,end-n+1:end);     % diagonal parts of A_i
M2=M(:,1:tnm);     % real parts of A_i
M3=M(:,tnm+1:2*tnm);     % imag parts of A_i
%%% for off diagonal parts of L:
%%%%%%%%%%%%%%%%%%%%%%%Following for Mfinal
Mtest=[M2  M2 -M3 -M3             % blocks for real parts
      -M3  M3 -M2  M2 ];           % blocks for imaginary parts
Mtest=...
[Mtest [M1 zeros(k,n); zeros(k,n) M1]];
%%% permute the rows first; for interlacing of real/imag parts
tempinds=1:k;
tempinds=[tempinds; tempinds+k];
tempinds=tempinds(:);
Mtest=Mtest(tempinds,:);   % interlacing rows real/imag
%%%now for column type interlacing based on triangular matrix structure
tempinds=[];
for ii=1:n-1,
	newinds=(1+ii*(ii-1)/2):(ii*(ii+1)/2);
	tempinds=[tempinds newinds newinds+tnm];
end
tempinds=[tempinds tempinds+length(tempinds)];
tempinds=[tempinds length(tempinds)+(1:(2*n))];
Mfinal=Mtest(:,tempinds);  % final block in Kronecker for upper part of Lsparse
%Mfinal works for kronecker for offdiag part
% Now for lower/diagonal related part of Lsparse
%   sqrt2 is needed as not used in original M2,M3
Mtemp=[sqrt2*M2 -sqrt2*M3 M1];      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The following evaluations are not needed but are kept here
% for documentation purposes.
%Loffdiagpart=kron(speye(tnm),Mfinal);  % tnm*(2*k) by tnm*2*n^2
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%diag part 
 % Ldiagpart=...
 %   [kron(speye(n-1),Mtemp) spalloc((n-1)*k,n^2,0)];
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%TP part 
 %LTPpart=(kron(ones(n,1),speye(n^2)))';
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%bottom part (diag-TP together)
 %Lbottom=[Ldiagpart;LTPpart];
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sparse block diag structure matrix rep
 %%Lsparse is 'large' block-diag structured form of matrix representation
 %Lsparse=blkdiag(Loffdiagpart,[Ldiagpart;LTPpart]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%now find the corresponding column ordering for P basis corresp. to Lsparse
%% upper triangular parts of blocks first
Lindsrealoff=[];
Lindsrealdiag=[];
for jj=2:m,  % column with element ii,jj for block 1,2
	for ii=1:jj-1,  % row with starting element ii,jj
		indsdiag=ii+(jj-2)*(jj-1)/2;
		newinds=[];
		for blockjj=2:n,
			 %% upper position is (ii,jj)
			 %% start with reflected position (jj,ii)
			 pcol=m*(blockjj-1)+ii;  % column position
			 elempos=(pcol-2)*(pcol-1)/2;
			 newinds=[newinds elempos+(jj:m:(blockjj-1)*m)];
			 %% continue with position (ii,jj)
			 pcol=m*(blockjj-1)+jj;  % column position
			 elempos=(pcol-2)*(pcol-1)/2;
			 elements=elempos+(ii:m:blockjj*m);
			 indsdiag=[indsdiag elements(end)];
			 newinds=[newinds elements(1:end-1)];
		end
		Lindsrealdiag=[Lindsrealdiag;indsdiag]; %diag els at end
	        Lindsrealoff=[Lindsrealoff;newinds];
	end
end
 % real part cols first; imag part cols next; diag parts last
  Lindsreal=[Lindsrealoff Lindsrealoff+tmnm Lindsrealdiag Lindsrealdiag+tmnm];
  Lindsreal=Lindsreal';
  tempLindstest=(Lindsreal(:))'; % all indices for uppertriangular parts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% find the indices in matrix rep L for the basis change
%% now find indices for diagonals of all blocks
  Lindsreal=[];
  for ii=1:n,   % diag element in the block
  for jj=2:n,    % starting in block 1,2, jj is the block column
	  tt=n*(jj-1)+ii; % current column number in matrix
	  %startelem: els from previous cols in previous blocks + 
	  startelem=((tt-1)*(tt-2)/2)+ii;
	  Lindsreal=[Lindsreal ; (startelem:n:startelem+(jj-2)*n)'];
  end
  end
  Lindsimag=Lindsreal+tmnm;
  Lindsreal=reshape(Lindsreal,tnm,m); 
  Lindsimag=reshape(Lindsimag,tnm,m);
  Lindsdiag=reshape(1:mn,n,m)+2*tmnm;
  Lindstest=[Lindsreal;Lindsimag;Lindsdiag'];
  Lindstest=(Lindstest(:))';   % indices in L for diagonal elements of
  % blocks: ordered by real ii els., emag ii els., diag ii els.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%change row order for Lsparse
  indsrows2=1:m^2:k*m^2;  %  first element from each of k blocks
  indsrows2=[indsrows2;indsrows2+tmm]; % two rows with first elements
       % i.e. first row has first element for reals for each of k blocks; 
        %    second row has first element for complexes for each of k blocks
  indsrows2=indsrows2(:)'; % puts real/complex in pairs (length 2*tmm)
  indsrows2=ones(tmm,1)*indsrows2; % tmm copies of the real/complex in pairs
  onesm=diag(0:tmm-1)*ones(tmm,2*k); % to allow for proper increase in vals
  indsrows2=indsrows2+onesm;  % increase the vals to proper pairs
  indsrows2=indsrows2';
  indsrows2=(indsrows2(:))'; % vector with real/complex pairs increasing order
 diagrows2=2*tmm+1:m^2:k*m^2;
 % remove last row since redundant constraints
 diagrows2=ones(m-1,1)*diagrows2;
 onesm=diag(0:m-2)*ones(m-1,k);
 diagrows2=diagrows2+onesm;
 diagrows2=diagrows2';
 diagrows2=(diagrows2(:))';
 indsrows=[indsrows2 diagrows2 k*m^2+1:k*m^2+n^2]; %%% all row indices
  colindices=[tempLindstest Lindstest]; % all col indices
[~,colindicesinv]=sort(colindices); % inverse permutation

%  %L2 is sparse form for original L with rows permuted/without deg rows
  %L2=Lsparse(:,colindicesinv);  % leave rows permuted 
  %% L2 is needed to evaluate residual in altproj.m
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%%gen inverse
  %pinvMfinal=(inv(Mfinal*Mfinal')*Mfinal)'; % is slower
  pinvMfinal=(Mfinal')*inv(Mfinal*Mfinal');  % is faster
%pinvLoffdiagpart=kron(speye(tnm),pinvMfinal);  % not used; block structure
pinvMtemp=pinv(Mtemp);
%% to find explicit expression for gen inv, use fact that range of gen
%    inv is same as range of transpose and use constraints of gen inv
%  to show that we need    Mtemp*Mtempnull=0
%  The last RHS blocks form a resolution of I
Mtempnull=null(Mtemp);  % need to use the orthog basis of nullspace
Mtempnull=(Mtempnull*Mtempnull')/n;
finalblock=speye(n^2)-(n-1)*Mtempnull;  % to satisfy spectral resol. of I
%pinvLbottom=...
%    [kron(speye(n-1),pinvMtemp) kron(ones(n-1,1),Mtempnull)
%     kron(ones(1,n-1),-pinvMtemp)  finalblock ];
%pinvL=blkdiag(pinvLoffdiagpart,pinvLbottom);
pcols=colindices;pcolsinv=colindicesinv;  % rename for global variables
