We investigate the thermodynamic properties of a polyelectrolyte solution in
a presence of {\it multivalent} salts. The polyions are modeled as rigid
cylinders with the charge distributed uniformly along the major axis. The
solution, besides the polyions, contain monovalent and divalent counterions as
well as monovalent coions. The strong electrostatic attraction existing between
the polyions and the counterions results in formation of clusters consisting of
one polyion and a number of associated monovalent and divalent counterions. The
theory presented in the paper allows us to explicitly construct the Helmholtz
free energy of a polyelectrolyte solution. The characteristic cluster size, as
well as any other thermodynamic property can then be determined by an
appropriate operation on the free energy