Index: trunk/psphot/src/psSparse.c
===================================================================
--- trunk/psphot/src/psSparse.c	(revision 6379)
+++ trunk/psphot/src/psSparse.c	(revision 6481)
@@ -34,17 +34,22 @@
     sparse->Bfj->data.F32[2] = B->data.F32[2];
 
-    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);
+    psSparseConstraint constraint;
+    constraint.paramMin   = -1e8;
+    constraint.paramMax   = +1e8;
+    constraint.paramDelta = +1e8;
+
+    x = psSparseSolve (x, constraint, sparse, 0);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 1);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 2);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 3);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 4);
     fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     return;
@@ -166,5 +171,5 @@
 }
 
-psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter) {
+psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter) {
 
     psF32 dG;
@@ -182,5 +187,14 @@
 	for (int i = 0; i < dQ->n; i++) {
 	    dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i];
+	    if (fabs (dG) > constraint.paramDelta) {
+		if (dG > 0) {
+		    dG = +constraint.paramDelta;
+		} else {
+		    dG = -constraint.paramDelta;
+		}
+	    }
 	    guess->data.F32[i] -= dG;
+	    guess->data.F32[i] = PS_MAX (guess->data.F32[i], constraint.paramMin);
+	    guess->data.F32[i] = PS_MIN (guess->data.F32[i], constraint.paramMax);
 	}
     }
@@ -220,3 +234,2 @@
 }
 
-