Index: /trunk/psphot/src/psphotSparseMatrix.c
===================================================================
--- /trunk/psphot/src/psphotSparseMatrix.c	(revision 5650)
+++ /trunk/psphot/src/psphotSparseMatrix.c	(revision 5650)
@@ -0,0 +1,152 @@
+# include "psphot.h"
+
+typedef struct {
+    psVector *Aij;
+    psVector *Bfj;
+    psVector *Qii;
+    psVector *Si;
+    psVector *Sj;
+    int Nelem;
+    int Nrow;
+} psSparse;
+
+psVector *psSparseMatrixTimesVector (psSparse *matrix, psVector *vector) {
+
+    psVector *output = psVectorAlloc (vector->n, PS_DATA_F32);
+    psVectorInit (output, 0);
+
+    Nelem = 0;
+    for (int j = 0; j < vector->n; j++) {
+	F = 0;
+	while (matrix->Sj->data.S32[Nelem] == j) {
+	    i = matrix->Sj->data.S32[Nelem];
+	    F += vector->data.F32[i] * matrix->Aij->data.F32[i];
+	}
+	output->data.F32[j] = F;
+    }
+    return (output);
+}
+
+psSparse *psphotStarOverlaps (psArray *sources) {
+
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
+    psVector *params = model->params;
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, psphotSortByY);
+
+    // fill out the models for each object
+    // XXX need to assign the model parameters using the PSF function, with unit peak
+    psArray  *pixlist = psArrayAlloc (sources->n);
+    for (int k = 0; k < sources->n; k++) {
+	src = sources->data[k];
+	nx = src->pixels->numCols;
+	ny = src->pixels->numRows;
+	im = psImageAlloc (nx, ny, PS_DATA_F32);
+	for (int i = 0; i < nx; i++) {
+	    for (int j = 0; j < ny; j++) {
+		x->data.F32[0] = (float) i + pix->col0;
+		x->data.F32[1] = (float) j + pix->row0;
+		im->data.F32[j][i] = modelFunc (NULL, params, x);
+	    }
+	}
+	pixlist->data[k] = im;
+    }
+    
+    psVector *Aij = psVectorAlloc (100, PS_DATA_F32);
+    psVector *Bfj = psVectorAlloc (100, PS_DATA_F32);
+    psVector *Si  = psVectorAlloc (100, PS_DATA_S32);
+    psVector *Sj  = psVectorAlloc (100, PS_DATA_S32);
+
+    int Ri, Rj;
+
+    Nelem = 0;
+    for (i = 0; i < pixlist->n; i++) {
+	pf = ((psSource *)sources->data[k])->pixels;
+	pi = pixlist->data[i];
+
+	// diagonal element (auto-cross-product)
+	f = psphotCrossProduct (pi, pi);
+	Aij->data.F32[Nelem] = f;
+	Si->data.S32[Nelem] = i;
+	Sj->data.S32[Nelem] = j;
+	Nelem ++;
+
+	Qii->data.F32[i] = f;
+
+	// find the image x model value
+	f = psphotCrossProduct (pi, pf);
+	Bfj->data.F32[i] = f;
+
+	// loop over all other stars following this one
+	for (j = i + 1; j < sources->n; j++) {
+	    pj = ((pmSource *)sources->data[j])->pixels;
+
+	    if (pi->col0 + pi->numCols < pj->col0) break;
+	    if (pj->col0 + pj->numCols < pi->col0) continue;
+	    if (pi->row0 + pi->numRows < pj->row0) continue;
+	    if (pj->row0 + pj->numRows < pi->row0) continue;
+	    
+	    // got an overlap; calculate cross-product and add to output array
+	    f = psphotCrossProduct (pi, pj);
+	    
+	    Aij->data.F32[Nelem] = f;
+	    Si->data.S32[Nelem] = i;
+	    Sj->data.S32[Nelem] = j;
+	    Nelem ++;
+
+	    Aij->data.F32[Nelem] = f;
+	    Si->data.S32[Nelem] = j;
+	    Sj->data.S32[Nelem] = i;
+	    Nelem ++;
+	}
+    }
+
+    index = psVectorSortIndex (NULL, Sj);
+
+    psVector *tAij = psVectorAlloc (Nelem, PS_DATA_F32);
+    psVector *tSi  = psVectorAlloc (Nelem, PS_DATA_S32);
+    psVector *tSj  = psVectorAlloc (Nelem, PS_DATA_S32);
+    for (i = 0; i < Nelem; i++) {
+	j = index->data.U32[i];
+	tAij->data.F32[i] = Aij->data.F32[j];
+	tSi->data.S32[i]  = Si->data.S32[j];
+	tSj->data.S32[i]  = Sj->data.S32[j];
+    }
+    psFree (Aij);
+    psFree (Si);
+    psFree (Sj);
+
+    sparse->Aij = Aij;
+    sparse->Bfj = Bfj;
+    sparse->Qii = Qii;
+    sparse->Si = Si;
+    sparse->Sj = Sj;
+    sparse->Nelem = Nelem;
+    sparse->Nrows = sources->n;
+
+    return (match);
+}
+
+psVector *psSparseSolve (sparse) {
+
+    psVector *Qii = sparse->Qii;
+    psVector *Bfj = sparse->Bfj;
+
+    psVector *Gii = psVectorCopy (NULL, Bfj);
+    psVector *dQ = psVectorAlloc (Gii->n, PS_DATA_F32);
+    psVector *dG = psVectorAlloc (Gii->n, PS_DATA_F32);
+
+    for (int j = 0; j < 2; j++) {
+	dQ = psSparseMatrixTimesVector (dQ, sparse, Gii);
+	for (int i = 0; i < dG->n; i++) {
+	    dG->data.F32[i] = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i];
+	    Gii->data.F32[i] -= dG->data.F32[i];
+	}
+    }
+    psFree (dQ);
+    psFree (dG);
+    
+    return (Gii);
+}