Index: /trunk/psLib/test/math/tap_psSparse.c
===================================================================
--- /trunk/psLib/test/math/tap_psSparse.c	(revision 10169)
+++ /trunk/psLib/test/math/tap_psSparse.c	(revision 10170)
@@ -8,13 +8,13 @@
 int main (void)
 {
-    plan_tests(23);
+    plan_tests(26);
 
     diag("psSparse() tests");
 
-    // test psSparseSolve for a simple example matrix
+    // test psSparseSolve for a simple normal example matrix
     {
         psMemId id = psMemGetId();
 
-        diag ("solve a simple, small matrix equation ");
+        diag ("solve a normalized matrix equation with psSparseSolve");
 
         // the basic equation is Ax = b
@@ -31,4 +31,77 @@
         {
             psSparseMatrixElement (matrix, i, i, 1.0);
+            if (i + 1 < 100) {
+                psSparseMatrixElement (matrix, i + 1, i, 0.1);
+            }
+        }
+
+        // incoming matrix elements do not need to be in order; sort before
+        // applying sparse matrix
+        psSparseResort (matrix);
+
+        psVector *xRef = psVectorAlloc (100, PS_TYPE_F32);
+        for (int i = 0; i < 100; i++)
+        {
+            xRef->data.F32[i] = 1.0;
+        }
+
+        psVector *bVec = psSparseMatrixTimesVector (NULL, matrix, xRef);
+
+        for (int i = 0; i < 100; i++)
+        {
+            psSparseVectorElement (matrix, i, bVec->data.F32[i]);
+        }
+
+        psSparseConstraint constraint;
+        constraint.paramMin   = 0.0;
+        constraint.paramMax   = 1e8;
+        constraint.paramDelta = 1e8;
+
+        // solve for normalization terms (need include local sky?)
+        psVector *xFit = psSparseSolve (NULL, constraint, matrix, 4);
+
+        // measure stdev between xFit and xRes
+        float dS = 0;
+        float dS2 = 0;
+        for (int i = 0; i < 100; i++)
+        {
+            float dX = xRef->data.F32[i] - xFit->data.F32[i];
+            // fprintf (stderr, "%5.3f %5.3f %5.3f %5.3f\n", bVec->data.F32[i], xRef->data.F32[i], xFit->data.F32[i], dX);
+            dS2 += PS_SQR(dX);
+            dS += dX;
+        }
+        dS /= 100.0;
+        dS2 = sqrt (dS2/100.0 - dS*dS);
+
+        ok_float_tol (dS2, 0.0, 1e-4, "scatter: %.20f", dS2);
+
+        psFree (matrix);
+        psFree (xRef);
+        psFree (bVec);
+        psFree (xFit);
+
+        skip_end();
+        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
+    }
+
+    // test psSparseSolve for a simple non-normal example matrix
+    {
+        psMemId id = psMemGetId();
+
+        diag ("solve a non-normalized matrix equation with psSparseSolve");
+
+        // the basic equation is Ax = b
+
+        // create a matrix A with diagonals of 1 and a small number of off diagonal elements.
+        // construct a vector x, construct the corresponding vector b by multiplication. solve
+        // Ax = b for x using psSparseSolve.  compare with the input values for x.
+
+        psSparse *matrix = psSparseAlloc (100, 100);
+        ok(matrix != NULL, "psSparse successfully allocated");
+        skip_start(matrix == NULL, 5, "Skipping tests because psSparseAlloc() failed");
+
+        for (int i = 0; i < 100; i++)
+        {
+            psSparseMatrixElement (matrix, i, i, 5.0);
             if (i + 1 < 100) {
                 psSparseMatrixElement (matrix, i + 1, i, 0.1);