Index: trunk/psLib/src/math/psSparse.c
===================================================================
--- trunk/psLib/src/math/psSparse.c	(revision 10017)
+++ trunk/psLib/src/math/psSparse.c	(revision 10121)
@@ -141,4 +141,5 @@
 }
 
+// solve Ax = B
 psVector *psSparseSolve(psVector *output, psSparseConstraint constraint, const psSparse *sparse, int Niter)
 {
@@ -151,5 +152,9 @@
     psVector *Bfj = sparse->Bfj;
 
+    // initial guess is B / Qii
     output = psVectorCopy(output, Bfj, PS_DATA_F32);
+    for (int i = 0; i < output->n; i++) {
+        output->data.F32[i] /= Qii->data.F32[i];
+    }
 
     // temporary storage for intermediate results
@@ -429,5 +434,5 @@
 
         // solve Sx = f
-        xVec = psSparseSolve(xVec, constraint, border->sparse, 4);
+        xVec = psSparseSolve(xVec, constraint, border->sparse, Niter);
 
         dG = psSparseBorderLowerProduct (dG, border, xVec);