function [Kr1,Mr1,Rcb1,phiN,psiC,K1,M1]=riduzioneCB(K1,M1,masters1,slaves1,dofcontrolled,Nmodes)
% partiziono le matrici
% kp1=K1([masters1;slaves1],[masters1;slaves1]);
% mp1=M1([masters1;slaves1],[masters1;slaves1]);
disp(['[' datestr(now,'HH:MM:SS') '] Riduzione CB'])

disp(['[' datestr(now,'HH:MM:SS') '] calcolo constraint modes psiC'])
tic
% Constraint modes - passaggio lungo
psiC=-(K1(slaves1,slaves1)\K1(slaves1,masters1));
% Kslave=K1(slaves1,slaves1);
% perm=amd(Kslave);
% Kslave=Kslave(perm,perm);
% L=chol(Kslave,'lower');
% slaves1=slaves1(perm);
% psiC=-(L.'\(L\K1(slaves1,masters1)));
% interm=(L\K1(slaves1,masters1));
% psiC=L.'\interm;
toc


disp(['[' datestr(now,'HH:MM:SS') '] slave modes phiN ', num2str(Nmodes)])
% Slave modes
[av,ev]=eigs(K1(slaves1,slaves1),M1(slaves1,slaves1),Nmodes,'sm');
evd=diag(ev);
[evds,ind]=sort(evd);
evds=abs(evds);

% phiN=av(:,ind);

phiNn=av(:,ind);
% normalizzo autovettori
mpr=phiNn'*M1(slaves1,slaves1)*phiNn;
mnorm=abs(diag(mpr)).^0.5;

phiN=zeros(size(phiNn));
for mn=1:Nmodes
    phiN(:,mn)=phiNn(:,mn)/mnorm(mn);
end
phiN=sparse(phiN);
%%%% Guyan
% phiN=[];
% Nmodes=0;
%%%%

% % % CB_controlled
c=length(dofcontrolled);
psiC(1:c,:)=[];
phiN(1:c,:)=[];
% psiC(:,[end-c+1:end]')=[];


% % Matrice di trasformazione
Rcb1=[speye(length(masters1)+c), sparse(length(masters1)+c,Nmodes);
    psiC , zeros(length(slaves1)-c,c), phiN];
% 
% Rcb1=[speye(length(masters1)), sparse(length(masters1),Nmodes);
%     psiC , zeros(length(slaves1),c), phiN];
% % % CB_controlled


%     % Queste due righe sono la riduzione vera e propria delle matrici di
%     % massa e di rigidezza, però per risparmiare tempo si può fare in
%     % maniera diversa che richiede meno moltiplicazioni e quindi meno tempo
%     Kr1=Rcb1'*kp1*Rcb1; % Rigidezza ridotta
%     Mr1=Rcb1'*kp1*Rcb1; % Massa ridotta

% Il risultato è sempre Kr1 e Mr1 di sopra, però questo dovrebbe essere più
% veloce
% Partizioni della matrice di massa e rigidezza
Kmm=K1(masters1,masters1);
Kss=K1(slaves1,slaves1);
Kms=K1(masters1,slaves1);
Ksm=K1(slaves1,masters1);
Mmm=M1(masters1,masters1); 
Mss=M1(slaves1,slaves1);
Mms=M1(masters1,slaves1);
Msm=M1(slaves1,masters1);

K1=[Kmm Kms; Ksm Kss];
M1=[Mmm Mms; Msm Mss];
% 
    Kr1=Rcb1'*K1*Rcb1; 
    Mr1=Rcb1'*M1*Rcb1;
% % CB_base
% disp(['[' datestr(now,'HH:MM:SS') '] calcoli per M'])
% % Calcoli per M
% % prealloc
% Mbb=zeros(size(masters1,1));
% Mbm=zeros(size(masters1,1),Nmodes);
% tic
% % Mbb=Mbb+Mmm+Mms*psiC+psiC'*Msm+psiC'*Mss*psiC;
% 
% Mbb=Mbb+Mmm+psiC.'*Mss*psiC;
% toc
% tic
% % Mbm=Mbm+(Mms+psiC'*Mss)*phiN;
% 
% %%%%% Guyan
% Mbm=Mbm+psiC.'*Mss*phiN;
% %%%%%
% toc
% %     Mmb=phiN'*(Msm+Mss*psiC); % trasposta di Mbm
% % Assemblaggio Mrid
% Mr1=[Mbb,Mbm;
%     Mbm.',speye(Nmodes)];
% 
% if ~isreal(Mr1)
%     warning('Complex Mrid')
%     keyboard
% end
% 
% % Calcoli per K
% disp(['[' datestr(now,'HH:MM:SS') '] calcoli per K'])
% % Kbb=zeros(size(masters1,1));
% tic
% Kbb=Kmm+Kms*psiC;
% toc
% % Assemblaggio Krid
% %%%% Guyan
% Kr1=[Kbb, sparse(length(Kbb),Nmodes);
%     sparse(Nmodes,length(Kbb)), diag(evds)];
% % Kr1=Kbb;
% %%%% Guyan
% 
% 
% if ~isreal(Kr1)
%     warning('Complex Krid')
%     keyboard
% end
% % CB_base
% end % function end
