Index: trunk/psphot/src/psSparse.c
===================================================================
--- trunk/psphot/src/psSparse.c	(revision 5653)
+++ trunk/psphot/src/psSparse.c	(revision 5654)
@@ -1,13 +1,3 @@
 # include "psphot.h"
-
-typedef struct {
-    psVector *Aij;
-    psVector *Bfj;
-    psVector *Qii;
-    psVector *Si;
-    psVector *Sj;
-    int Nelem;
-    int Nrows;
-} psSparse;
 
 void psSparseMatrixTest () {
@@ -16,12 +6,12 @@
     psSparse *sparse = psSparseAlloc (3, 9);
 
-    psSparseMatrixElement (sparse, 0, 0, 1.0);
-    psSparseMatrixElement (sparse, 1, 1, 1.0);
+    psSparseMatrixElement (sparse, 0, 0, 3.0);
+    psSparseMatrixElement (sparse, 1, 1, 2.0);
     psSparseMatrixElement (sparse, 2, 2, 1.0);
 
-    psSparseMatrixElement (sparse, 0, 1, 0.1);
-    psSparseMatrixElement (sparse, 0, 2, 0.1);
-
-    psSparseMatrixResort (sparse);
+    psSparseMatrixElement (sparse, 1, 0, 0.1);
+    psSparseMatrixElement (sparse, 2, 0, -0.1);
+
+    psSparseResort (sparse);
     for (int i = 0; i < sparse->Nelem; i++) {
 	fprintf (stderr, "%d %d %f\n", 
@@ -33,8 +23,9 @@
     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);
+    x->data.F32[1] = 5;
+    x->data.F32[2] = 7;
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    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]);
 
@@ -43,5 +34,17 @@
     sparse->Bfj->data.F32[2] = B->data.F32[2];
 
-    x = psSparseSolve (x, sparse);
+    x = psSparseSolve (x, sparse, 0);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, sparse, 1);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, sparse, 2);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, sparse, 3);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, sparse, 4);
     fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     return;
@@ -50,10 +53,10 @@
 void psSparseResort (psSparse *sparse) {
 
-    Nelem = sparse->Nelem;
+    int Nelem = sparse->Nelem;
 
     psVector *index = psVectorSortIndex (NULL, sparse->Sj);
-    psVector Aij = sparse->Aij;
-    psVector Si = sparse->Si;
-    psVector Sj = sparse->Sj;
+    psVector *Aij = sparse->Aij;
+    psVector *Si = sparse->Si;
+    psVector *Sj = sparse->Sj;
 
     // allocate new temporary vectors
@@ -61,6 +64,6 @@
     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];
+    for (int i = 0; i < Nelem; i++) {
+	int j = index->data.U32[i];
 	tAij->data.F32[i] = Aij->data.F32[j];
 	tSi->data.S32[i]  = Si->data.S32[j];
@@ -80,4 +83,6 @@
 void psSparseMatrixElement (psSparse *sparse, int i, int j, float value) {
 
+    int k;
+
     if (i < j) {
 	fprintf (stderr, "*** error: subdiagonal element ***\n");
@@ -91,10 +96,7 @@
 	// 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);
+	    psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100);
+	    psVectorRealloc (sparse->Si,  sparse->Si->nalloc + 100);
+	    psVectorRealloc (sparse->Sj,  sparse->Sj->nalloc + 100);
 	}
 
@@ -111,10 +113,7 @@
 	// 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);
+	    psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100);
+	    psVectorRealloc (sparse->Si,  sparse->Si->nalloc + 100);
+	    psVectorRealloc (sparse->Sj,  sparse->Sj->nalloc + 100);
 	}
 
@@ -139,5 +138,5 @@
 void psSparseVectorElement (psSparse *sparse, int i, float value) {
 
-    sparse->Bjf->data.F32[i] = value;
+    sparse->Bfj->data.F32[i] = value;
     return;
 }
@@ -146,6 +145,9 @@
 psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector) {
 
+    int i, Nelem;
+    float F;
+
     if (output == NULL) {
-	psVector *output = psVectorAlloc (vector->n, PS_DATA_F32);
+	output = psVectorAlloc (vector->n, PS_DATA_F32);
     }
 
@@ -154,6 +156,7 @@
 	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];
+	    i = matrix->Si->data.S32[Nelem];
+	    F += vector->data.F32[i] * matrix->Aij->data.F32[Nelem];
+	    Nelem++;
 	}
 	output->data.F32[j] = F;
@@ -162,5 +165,5 @@
 }
 
-psVector *psSparseSolve (psVector *guess, psSparse *sparse) {
+psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter) {
 
     psF32 dG;
@@ -169,12 +172,12 @@
     psVector *Bfj = sparse->Bfj;
 
-    guess = psVectorCopy (guess, Bfj);
+    guess = psVectorCopy (guess, Bfj, PS_DATA_F32);
 
     // temporary storage for intermediate results
     psVector *dQ = psVectorAlloc (guess->n, PS_DATA_F32);
 
-    for (int j = 0; j < 2; j++) {
+    for (int j = 0; j < Niter; j++) {
 	dQ = psSparseMatrixTimesVector (dQ, sparse, guess);
-	for (int i = 0; i < dG->n; i++) {
+	for (int i = 0; i < dQ->n; i++) {
 	    dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i];
 	    guess->data.F32[i] -= dG;
@@ -187,9 +190,9 @@
 static void psSparseFree (psSparse *sparse) {
     if (sparse == NULL) return;
-    psFree (sparse.Aij);
-    psFree (sparse.Bfj);
-    psFree (sparse.Qii);
-    psFree (sparse.Si);
-    psFree (sparse.Sj);
+    psFree (sparse->Aij);
+    psFree (sparse->Bfj);
+    psFree (sparse->Qii);
+    psFree (sparse->Si);
+    psFree (sparse->Sj);
     return;
 }
@@ -212,4 +215,5 @@
 
     sparse->Nrows = Nrows;
+    psMemSetDeallocator(sparse, (psFreeFunc) psSparseFree);
     return (sparse);
 }