Index: /trunk/psphot/Makefile
===================================================================
--- /trunk/psphot/Makefile	(revision 5652)
+++ /trunk/psphot/Makefile	(revision 5653)
@@ -69,5 +69,5 @@
 TEST = \
 $(SRC)/psphotTest.$(ARCH).o \
-$(SRC)/psphotSparseMatrix.$(ARCH).o
+$(SRC)/psSparse.$(ARCH).o
 
 psphot: $(BIN)/psphot.$(ARCH)
Index: /trunk/psphot/src/psSparse.c
===================================================================
--- /trunk/psphot/src/psSparse.c	(revision 5653)
+++ /trunk/psphot/src/psSparse.c	(revision 5653)
@@ -0,0 +1,217 @@
+# include "psphot.h"
+
+typedef struct {
+    psVector *Aij;
+    psVector *Bfj;
+    psVector *Qii;
+    psVector *Si;
+    psVector *Sj;
+    int Nelem;
+    int Nrows;
+} psSparse;
+
+void psSparseMatrixTest () {
+
+    // build a sparse matrix
+    psSparse *sparse = psSparseAlloc (3, 9);
+
+    psSparseMatrixElement (sparse, 0, 0, 1.0);
+    psSparseMatrixElement (sparse, 1, 1, 1.0);
+    psSparseMatrixElement (sparse, 2, 2, 1.0);
+
+    psSparseMatrixElement (sparse, 0, 1, 0.1);
+    psSparseMatrixElement (sparse, 0, 2, 0.1);
+
+    psSparseMatrixResort (sparse);
+    for (int i = 0; i < sparse->Nelem; i++) {
+	fprintf (stderr, "%d %d %f\n", 
+		 sparse->Si->data.S32[i],
+		 sparse->Sj->data.S32[i],
+		 sparse->Aij->data.F32[i]);
+    }
+
+    psVector *x = psVectorAlloc (3, PS_DATA_F32);
+    x->data.F32[0] = 3;
+    x->data.F32[0] = 5;
+    x->data.F32[0] = 7;
+
+    psVector B = psSparseMatrixTimesVector (NULL, sparse, x);
+    fprintf (stderr, "B: %f %f %f\n", B->data.F32[0], B->data.F32[1], B->data.F32[2]);
+
+    sparse->Bfj->data.F32[0] = B->data.F32[0];
+    sparse->Bfj->data.F32[1] = B->data.F32[1];
+    sparse->Bfj->data.F32[2] = B->data.F32[2];
+
+    x = psSparseSolve (x, sparse);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+    return;
+} 
+
+void psSparseResort (psSparse *sparse) {
+
+    Nelem = sparse->Nelem;
+
+    psVector *index = psVectorSortIndex (NULL, sparse->Sj);
+    psVector Aij = sparse->Aij;
+    psVector Si = sparse->Si;
+    psVector Sj = sparse->Sj;
+
+    // allocate new temporary vectors
+    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 = tAij;
+    sparse->Si = tSi;
+    sparse->Sj = tSj;
+    return;
+}
+
+// user should only add elements above the diagonal, but we don't check this
+void psSparseMatrixElement (psSparse *sparse, int i, int j, float value) {
+
+    if (i < j) {
+	fprintf (stderr, "*** error: subdiagonal element ***\n");
+	return;
+    }
+
+    if (i == j) {
+	// add to the diagonal
+	sparse->Qii->data.F32[i] = value;
+
+	// check vectors lengths and extend if needed
+	if (sparse->Nelem >= sparse->Aij->nalloc) {
+	    sparse->Aij->nalloc += 100;	    
+	    psVectorRealloc (sparse->Aij, sparse->Aij->nalloc);
+	    sparse->Si->nalloc += 100;	    
+	    psVectorRealloc (sparse->Si, sparse->Si->nalloc);
+	    sparse->Sj->nalloc += 100;	    
+	    psVectorRealloc (sparse->Sj, sparse->Sj->nalloc);
+	}
+
+	k = sparse->Nelem;
+	sparse->Aij->data.F32[k] = value;
+	sparse->Si->data.S32[k]  = i;
+	sparse->Sj->data.S32[k]  = j;
+	
+	sparse->Nelem ++;
+	sparse->Aij->n ++;
+	sparse->Si->n ++;
+	sparse->Sj->n ++;
+    } else {
+	// check vectors lengths and extend if needed
+	if (sparse->Nelem >= sparse->Aij->nalloc - 1) {
+	    sparse->Aij->nalloc += 100;	    
+	    psVectorRealloc (sparse->Aij, sparse->Aij->nalloc);
+	    sparse->Si->nalloc += 100;	    
+	    psVectorRealloc (sparse->Si, sparse->Si->nalloc);
+	    sparse->Sj->nalloc += 100;	    
+	    psVectorRealloc (sparse->Sj, sparse->Sj->nalloc);
+	}
+
+	k = sparse->Nelem;
+	sparse->Aij->data.F32[k] = value;
+	sparse->Si->data.S32[k]  = i;
+	sparse->Sj->data.S32[k]  = j;
+	k++;
+
+	sparse->Aij->data.F32[k] = value;
+	sparse->Si->data.S32[k]  = j;
+	sparse->Sj->data.S32[k]  = i;
+	
+	sparse->Nelem  += 2;
+	sparse->Aij->n += 2;
+	sparse->Si->n  += 2;
+	sparse->Sj->n  += 2;
+    }
+    return;
+}
+
+void psSparseVectorElement (psSparse *sparse, int i, float value) {
+
+    sparse->Bjf->data.F32[i] = value;
+    return;
+}
+
+// multiple A * x
+psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector) {
+
+    if (output == NULL) {
+	psVector *output = psVectorAlloc (vector->n, PS_DATA_F32);
+    }
+
+    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);
+}
+
+psVector *psSparseSolve (psVector *guess, psSparse *sparse) {
+
+    psF32 dG;
+
+    psVector *Qii = sparse->Qii;
+    psVector *Bfj = sparse->Bfj;
+
+    guess = psVectorCopy (guess, Bfj);
+
+    // temporary storage for intermediate results
+    psVector *dQ = psVectorAlloc (guess->n, PS_DATA_F32);
+
+    for (int j = 0; j < 2; j++) {
+	dQ = psSparseMatrixTimesVector (dQ, sparse, guess);
+	for (int i = 0; i < dG->n; i++) {
+	    dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i];
+	    guess->data.F32[i] -= dG;
+	}
+    }
+    psFree (dQ);
+    return (guess);
+}
+
+static void psSparseFree (psSparse *sparse) {
+    if (sparse == NULL) return;
+    psFree (sparse.Aij);
+    psFree (sparse.Bfj);
+    psFree (sparse.Qii);
+    psFree (sparse.Si);
+    psFree (sparse.Sj);
+    return;
+}
+
+// allocate a sparse matrix container for Nrows, with Nelem slots allocated
+psSparse *psSparseAlloc (int Nrows, int Nelem) {
+
+    psSparse *sparse = (psSparse *) psAlloc (sizeof(psSparse));
+    sparse->Aij = psVectorAlloc (Nelem, PS_DATA_F32);
+    sparse->Si  = psVectorAlloc (Nelem, PS_DATA_S32);
+    sparse->Sj  = psVectorAlloc (Nelem, PS_DATA_S32);
+
+    sparse->Aij->n = 0;
+    sparse->Si->n  = 0;
+    sparse->Sj->n  = 0;
+    sparse->Nelem = 0;
+
+    sparse->Bfj = psVectorAlloc (Nrows, PS_DATA_F32);
+    sparse->Qii = psVectorAlloc (Nrows, PS_DATA_F32);
+
+    sparse->Nrows = Nrows;
+    return (sparse);
+}
+
+
Index: /trunk/psphot/src/psphotSparseMatrix.c
===================================================================
--- /trunk/psphot/src/psphotSparseMatrix.c	(revision 5652)
+++ /trunk/psphot/src/psphotSparseMatrix.c	(revision 5653)
@@ -1,37 +1,3 @@
 # include "psphot.h"
-
-typedef struct {
-    psVector *Aij;
-    psVector *Bfj;
-    psVector *Qii;
-    psVector *Si;
-    psVector *Sj;
-    int Nelem;
-    int Nrow;
-} psSparse;
-
-void psphotSparseMatrixTest () {
-
-    // build a sparse matrix
-   
-
-} 
-
-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) {
@@ -136,24 +102,2 @@
     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);
-}