Index: /branches/pap/psLib/src/math/psMatrix.c
===================================================================
--- /branches/pap/psLib/src/math/psMatrix.c	(revision 25864)
+++ /branches/pap/psLib/src/math/psMatrix.c	(revision 25865)
@@ -1033,10 +1033,11 @@
 }
 
-psVector *psMatrixSolveSVD(const psImage *matrix, const psVector *vector)
-{
-    PS_ASSERT_IMAGE_NON_NULL(matrix, NULL);
-    PS_CHECK_DIMEN_AND_TYPE(matrix, PS_DIMEN_IMAGE, {return NULL;});
-    PS_ASSERT_VECTOR_NON_NULL(vector, NULL);
-    PS_CHECK_DIMEN_AND_TYPE(vector, PS_DIMEN_VECTOR, {return NULL;});
+psVector *psMatrixSolveSVD(psVector *out, const psImage *matrix, const psVector *vector)
+{
+    #define psMatrixSolveSVD_EXIT {psFree(out); return NULL; }
+    PS_ASSERT_GENERAL_IMAGE_NON_NULL(matrix, psMatrixSolveSVD_EXIT);
+    PS_CHECK_DIMEN_AND_TYPE(matrix, PS_DIMEN_IMAGE, psMatrixSolveSVD_EXIT);
+    PS_ASSERT_GENERAL_VECTOR_NON_NULL(RHS, psMatrixSolveSVD_EXIT);
+    PS_CHECK_DIMEN_AND_TYPE(vector, PS_DIMEN_VECTOR, psMatrixSolveSVD_EXIT);
 
     int numCols = matrix->numCols, numRows = matrix->numRows; // Size of matrix
@@ -1084,5 +1085,5 @@
     gsl_vector_free(b);
 
-    psVector *solution = psVectorAlloc(numCols, PS_TYPE_F64);
+    out = psVectorRecycle(out, numCols, PS_TYPE_F64);
 
     vectorGSLtoPS(solution, x);
Index: /branches/pap/psLib/src/math/psMatrix.h
===================================================================
--- /branches/pap/psLib/src/math/psMatrix.h	(revision 25864)
+++ /branches/pap/psLib/src/math/psMatrix.h	(revision 25865)
@@ -190,4 +190,5 @@
 /// Solves Ax = b for x
 psVector *psMatrixSolveSVD(
+    psVector *solution,                 ///< Solution to output, or NULL
     const psImage *matrix,              ///< Matrix to be solved
     const psVector *vector              ///< Vector of values