% This is called from the main program

% Form the P matrices
k=[1:ngen;(ngen+1):2*ngen];
kk=k(:);

PGmod1=sparse(zeros(nb,ngen));
PGmod1((gen(:,1)),[1:ngen])=diag(ones(1,ngen));
PGmod=[PGmod1 zeros(nb,ngen);zeros(nb,ngen) PGmod1];
PG=PGmod(ll,kk);

clear PGmod1;
clear PGmod;

m=[1:nload;(nload+1):2*nload];
mmm=m(:);

PLmod1=sparse(zeros(nb,nload));
PLmod1((ld(:,1)),[1:nload])=diag(ones(1,nload));
PLmod=[PLmod1 zeros(nb,nload);zeros(nb,nload) PLmod1];
PL=PLmod(ll,mmm);

clear PLmod1;
clear PLmod;

%pack;