Abstract : Hyperspectral data unmixing aims at identifying the components (endmembers) of an observed surface and at determining their fractional abundances inside each pixel area. Assuming that the spectral signatures of the surface components have been previously determined by an endmember extraction algorithm, or to be part of an available spectral library, the main problem is reduced to the estimation of the fractional abundances. For large hyperspectral image data sets, the estimation of the abundance maps requires the resolution of a large-scale optimization problem subject to linear constraints such as non-negativity and sum less or equal to one. This paper proposes a primal-dual interior-point optimization algorithm allowing a constrained least squares estimation approach. In comparison with existing methods, the proposed algorithm is more flexible since it can handle any linear equality and/or inequality constraint and has the advantage of a reduced computational cost. It also presents an algorithmic structure suitable for a parallel implementation on modern intensive computing devices such as Graphics Processing Units (GPU). The implementation issues are discussed and the applicability of the proposed approach is illustrated with the help of examples on synthetic and real hyperspectral data.