A fast method is described for evaluating multiple exchange matrices in a Gaussian atomic orbital basis. For insulators, it is asymptotically linear scaling, and is a generalization of the linear scaling exchange (LinK) method, which was formulated for a single exchange matrix [J. Chem. Phys. 109 (1998) 1663]. It is employed to evaluate exchange-type contractions of all derivative density matrices with two-electron integrals for a series of linear alkanes, linear polyacenes, and water clusters using STO-3G, 3-21G, and 6-31G* basis sets. Significant computational savings are obtained for molecules with as few as 10 non-hydrogen atoms. (C) 2002 Published by Elsevier Science B.V.