[sdpb] 85/233: First attempt at stabilized cholesky decomposition

Tobias Hansen thansen at moszumanska.debian.org
Thu Mar 9 04:06:22 UTC 2017


This is an automated email from the git hooks/post-receive script.

thansen pushed a commit to branch master
in repository sdpb.

commit 3588eb052f22e53785c3711952bd40b6883abaf7
Author: David Simmons-Duffin <dsd at athena.sns.ias.edu>
Date:   Mon Sep 29 18:42:59 2014 -0400

    First attempt at stabilized cholesky decomposition
---
 src/Matrix.cpp                 |  18 +++++
 src/Matrix.h                   |   2 +
 src/main.cpp                   |   4 +-
 src/mpack/Rpotf2Stabilized.cpp | 162 +++++++++++++++++++++++++++++++++++++++++
 src/mpack/RpotrfStabilized.cpp | 147 +++++++++++++++++++++++++++++++++++++
 src/mpack/mlapack_gmp.h        |   6 ++
 src/tests.cpp                  |  26 +++++--
 src/tests.h                    |   6 ++
 8 files changed, 362 insertions(+), 9 deletions(-)

diff --git a/src/Matrix.cpp b/src/Matrix.cpp
index b1b32de..34d1179 100644
--- a/src/Matrix.cpp
+++ b/src/Matrix.cpp
@@ -238,6 +238,24 @@ void choleskyDecomposition(Matrix &A, Matrix &L) {
       L.elements[i + j*dim] = 0;
 }
 
+void choleskyDecompositionStabilized(Matrix &A, Matrix &L, vector<Integer> &stabilizeIndices, vector<Real> &stabilizeLambdas) {
+  int dim = A.rows;
+  assert(A.cols == dim);
+  assert(L.rows == dim);
+  assert(L.cols == dim);
+
+  // Set lower-triangular part of L to cholesky decomposition
+  L.copyFrom(A);
+  Integer info;
+  RpotrfStabilized("Lower", dim, &L.elements[0], dim, &info, stabilizeIndices, stabilizeLambdas);
+  assert(info == 0);
+
+  // Set the upper-triangular part of the L to zero
+  for (int j = 0; j < dim; j++)
+    for (int i = 0; i < j; i++)
+      L.elements[i + j*dim] = 0;
+}
+
 // L' (lower triangular) such that L' L'^T = L L^T + v v^T. i.e., if L
 // is a cholesky decomposition of A, then L' is a cholesky
 // decomposition of A + v v^T.  This function dominates the running
diff --git a/src/Matrix.h b/src/Matrix.h
index d9616d1..6a99466 100644
--- a/src/Matrix.h
+++ b/src/Matrix.h
@@ -174,6 +174,8 @@ void solveWithLUDecomposition(Matrix &LU, vector<Integer> &pivots, Vector &b);
 // - L : dim x dim lower-triangular matrix
 void choleskyDecomposition(Matrix &A, Matrix &L);
 
+void choleskyDecompositionStabilized(Matrix &A, Matrix &L, vector<Integer> &stabilizeIndices, vector<Real> &stabilizeLambdas);
+
 // L' (lower triangular) such that L' L'^T = L L^T + v v^T. i.e., if L
 // is a cholesky decomposition of A, then L' is a cholesky
 // decomposition of A + v v^T.  This function dominates the running
diff --git a/src/main.cpp b/src/main.cpp
index 295254d..28faac9 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -9,6 +9,7 @@
 #include "Timers.h"
 #include "SDP.h"
 #include "parse.h"
+#include "tests.h"
 #include "SDPSolver.h"
 
 using std::cout;
@@ -183,5 +184,6 @@ int main(int argc, char** argv) {
     return 1; 
   } 
 
-  return solveSDP(sdpFile, outFile, checkpointFile, parameters);
+  //return solveSDP(sdpFile, outFile, checkpointFile, parameters);
+  testCholeskyStabilize();
 }
diff --git a/src/mpack/Rpotf2Stabilized.cpp b/src/mpack/Rpotf2Stabilized.cpp
new file mode 100644
index 0000000..e4f5e47
--- /dev/null
+++ b/src/mpack/Rpotf2Stabilized.cpp
@@ -0,0 +1,162 @@
+/*
+ * Copyright (c) 2008-2010
+ *      Nakata, Maho
+ *      All rights reserved.
+ *
+ *  $Id: Rpotf2.cpp,v 1.10 2010/08/07 04:48:33 nakatamaho Exp $ 
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ */
+/*
+Copyright (c) 1992-2007 The University of Tennessee.  All rights reserved.
+
+$COPYRIGHT$
+
+Additional copyrights may follow
+
+$HEADER$
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+  notice, this list of conditions and the following disclaimer. 
+  
+- Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions and the following disclaimer listed
+  in this license in the documentation and/or other materials
+  provided with the distribution.
+  
+- Neither the name of the copyright holders nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+  
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT  
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT  
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
+*/
+
+#include <mblas.h>
+#include <mlapack.h>
+
+#include <vector>
+#include <assert.h>
+#include <iostream>
+
+using std::cout;
+using std::endl;
+
+const double CHOLESKY_STABILIZE_THRESHOLD = 1e-10;
+
+REAL lambdaGeometricMean(double totalLogLambda, INTEGER numLambdas) {
+  assert(numLambdas != 0);
+  return REAL(exp(totalLogLambda/numLambdas));
+}
+
+bool isSmallDiagonal(const REAL &ajj, double totalLogLambda, INTEGER numLambdas) {
+  double d = cast2double(ajj);
+  return (d <= 0) || (numLambdas != 0 && log(d)/2 < totalLogLambda/numLambdas + log(CHOLESKY_STABILIZE_THRESHOLD));
+}
+
+void correctDiagonal(REAL &ajj, INTEGER jIndex, double totalLogLambda, vector<INTEGER> &stabilizeIndices, vector<REAL> &stabilizeLambdas) {
+  REAL lambda = lambdaGeometricMean(totalLogLambda, jIndex);
+  ajj += lambda * lambda;
+  stabilizeIndices.push_back(jIndex);
+  stabilizeLambdas.push_back(lambda);
+}
+
+void Rpotf2Stabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGER * info, INTEGER indexStart, vector<INTEGER> &stabilizeIndices, vector<REAL> &stabilizeLambdas, double &totalLogLambda)
+{
+    INTEGER j, upper, success = 1;
+    REAL ajj;
+    REAL Zero = 0.0, One = 1.0;
+
+    *info = 0;
+    upper = Mlsame(uplo, "U");
+    if (!upper && !Mlsame(uplo, "L")) {
+	*info = -1;
+    } else if (n < 0) {
+	*info = -2;
+    } else if (lda < max((INTEGER) 1, n)) {
+	*info = -4;
+    }
+    if (*info != 0) {
+	Mxerbla("Rpotf2", -(*info));
+	return;
+    }
+//Quick return if possible
+    if (n == 0) {
+	return;
+    }
+    if (upper) {
+//Compute the Cholesky factorization A = U'*U.
+	for (j = 0; j < n; j++) {
+//Compute U(J,J) and test for non-positive-definiteness.
+	    ajj = A[j + j * lda] - Rdot(j, &A[j * lda], 1, &A[j * lda], 1);
+
+            if (isSmallDiagonal(ajj, totalLogLambda, indexStart + j))
+              correctDiagonal(ajj, indexStart + j, totalLogLambda, stabilizeIndices, stabilizeLambdas);
+
+	    ajj = sqrt(ajj);
+            totalLogLambda += log(cast2double(ajj));
+
+	    A[j + j * lda] = ajj;
+//Compute elements J+1:N of row J.
+	    if (j < n) {
+		Rgemv("Transpose", j, n - j - 1, -One, &A[(j + 1) * lda], lda, &A[j * lda], 1, One, &A[j + (j + 1) * lda], lda);
+		Rscal(n - j - 1, One / ajj, &A[j + (j + 1) * lda], lda);
+	    }
+	}
+    } else {
+//Compute the Cholesky factorization A = L*L'.
+	for (j = 0; j < n; j++) {
+//Compute L(J,J) and test for non-positive-definiteness.
+	    ajj = A[j + j * lda] - Rdot(j, &A[j], lda, &A[j], lda);
+
+            if (isSmallDiagonal(ajj, totalLogLambda, indexStart + j))
+              correctDiagonal(ajj, indexStart + j, totalLogLambda, stabilizeIndices, stabilizeLambdas);
+
+	    ajj = sqrt(ajj);
+            totalLogLambda += log(cast2double(ajj));
+
+	    A[j + j * lda] = ajj;
+//Compute elements J+1:N of column J.
+	    if (j < n) {
+		Rgemv("No transpose", n - j - 1, j, -One, &A[j + 1], lda, &A[j], lda, One, &A[j + 1 + j * lda], 1);
+		Rscal(n - j - 1, One / ajj, &A[j + 1 + j * lda], 1);
+	    }
+	}
+    }
+    if (!success)
+	*info = j + 1;
+    return;
+}
diff --git a/src/mpack/RpotrfStabilized.cpp b/src/mpack/RpotrfStabilized.cpp
new file mode 100644
index 0000000..cbf787b
--- /dev/null
+++ b/src/mpack/RpotrfStabilized.cpp
@@ -0,0 +1,147 @@
+/*
+ * Copyright (c) 2008-2010
+ *      Nakata, Maho
+ *      All rights reserved.
+ *
+ *  $Id: Rpotrf.cpp,v 1.9 2010/08/07 04:48:33 nakatamaho Exp $ 
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ */
+/*
+Copyright (c) 1992-2007 The University of Tennessee.  All rights reserved.
+
+$COPYRIGHT$
+
+Additional copyrights may follow
+
+$HEADER$
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+  notice, this list of conditions and the following disclaimer. 
+  
+- Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions and the following disclaimer listed
+  in this license in the documentation and/or other materials
+  provided with the distribution.
+  
+- Neither the name of the copyright holders nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+  
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT  
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT  
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
+*/
+
+#include <mblas.h>
+#include <mlapack.h>
+
+#include <vector>
+
+void RpotrfStabilized(const char *uplo, INTEGER n, REAL * A, INTEGER lda, INTEGER * info, vector<INTEGER> &stabilizeIndices, vector<REAL> &stabilizeLambdas)
+{
+    INTEGER upper;
+    INTEGER j, jb, nb;
+    REAL One = 1.0;
+    double totalLogLambda = 0;
+
+//Test the input parameters.
+    *info = 0;
+    upper = Mlsame(uplo, "U");
+    if (!upper && !Mlsame(uplo, "L")) {
+	*info = -1;
+    } else if (n < 0) {
+	*info = -2;
+    } else if (lda < max((INTEGER) 1, n)) {
+	*info = -4;
+    }
+    if (*info != 0) {
+	Mxerbla("Rpotrf", -(*info));
+	return;
+    }
+//Quick return if possible
+    if (n == 0) {
+	return;
+    }
+//Determine the block size for this environment.
+    nb = iMlaenv(1, "Rpotrf", uplo, n, -1, -1, -1);
+    if (nb <= 1 || nb >= n) {
+//Use unblocked code.
+      Rpotf2Stabilized(uplo, n, A, lda, info, 0, stabilizeIndices, stabilizeLambdas, totalLogLambda);
+    } else {
+//Use blocked code.
+	if (upper) {
+//Compute the Cholesky factorization A = U'*U.
+	    for (j = 1; j <= n; j = j + nb) {
+//Update and factorize the current diagonal block and test
+//for non-positive-definiteness.
+		jb = min(nb, n - j + 1);
+		Rsyrk("Upper", "Transpose", jb, j - 1, -One, &A[0 + (j - 1) * lda], lda, One, &A[(j - 1) + (j - 1) * lda], lda);
+		Rpotf2Stabilized("Upper", jb, &A[(j - 1) + (j - 1) * lda], lda, info, j-1, stabilizeIndices, stabilizeLambdas, totalLogLambda);
+		if (*info != 0) {
+		    goto L30;
+		}
+		if (j + jb <= n) {
+//Compute the current block row.
+		    Rgemm("Transpose", "No transpose", jb, n - j - jb + 1,
+			  j - 1, -One, &A[0 + (j - 1) * lda], lda, &A[0 + (j + jb - 1) * lda], lda, One, &A[(j - 1) + (j + jb - 1) * lda], lda);
+		    Rtrsm("Left", "Upper", "Transpose", "Non-unit", jb, n - j - jb + 1, One, &A[(j - 1) + (j - 1) * lda], lda, &A[(j - 1) + (j + jb - 1) * lda], lda);
+		}
+	    }
+	} else {
+//Compute the Cholesky factorization A = L*L'.
+	    for (j = 1; j <= n; j = j + nb) {
+//Update and factorize the current diagonal block and test
+//for non-positive-definiteness.
+		jb = min(nb, n - j + 1);
+		Rsyrk("Lower", "No transpose", jb, j - 1, -One, &A[(j - 1) + 0 * lda], lda, One, &A[(j - 1) + (j - 1) * lda], lda);
+		Rpotf2Stabilized("Lower", jb, &A[(j - 1) + (j - 1) * lda], lda, info, j-1, stabilizeIndices, stabilizeLambdas, totalLogLambda);
+		if (*info != 0) {
+		    goto L30;
+		}
+		if (j + jb <= n) {
+//Compute the current block column.
+		    Rgemm("No transpose", "Transpose", n - j - jb + 1, jb,
+			  j - 1, -One, &A[(j + jb - 1) + 0 * lda], lda, &A[(j - 1) + 0 * lda], lda, One, &A[(j + jb - 1) + (j - 1) * lda], lda);
+		    Rtrsm("Right", "Lower", "Transpose", "Non-unit", n - j - jb + 1, jb, One, &A[(j - 1) + (j - 1) * lda], lda, &A[(j + jb - 1) + (j - 1) * lda], lda);
+		}
+	    }
+	}
+    }
+    return;
+  L30:
+    *info = *info + j - 1;
+    return;
+}
diff --git a/src/mpack/mlapack_gmp.h b/src/mpack/mlapack_gmp.h
index d902d79..fdab543 100644
--- a/src/mpack/mlapack_gmp.h
+++ b/src/mpack/mlapack_gmp.h
@@ -30,6 +30,10 @@
 #ifndef _MLAPACK_GMP_H_
 #define _MLAPACK_GMP_H_
 
+#include <vector>
+
+using std::vector;
+
 /* this is a subset of mpack for SDPA-GMP only */
 /* http://mplapack.sourceforge.net/ */
 
@@ -40,6 +44,7 @@ void
     Rsyev(const char *jobz, const char *uplo, mpackint n, mpf_class * A,
     mpackint lda, mpf_class * w, mpf_class * work, mpackint lwork, mpackint *info);
 void Rpotrf(const char *uplo, mpackint n, mpf_class * A, mpackint lda, mpackint *info);
+void RpotrfStabilized(const char *uplo, mpackint n, mpf_class * A, mpackint lda, mpackint * info, vector<mpackint> &stabilizeIndices, vector<mpf_class> &stabilizeLambdas);
 mpackint iMlaenv_gmp(mpackint ispec, const char *name, const char *opts, mpackint n1, mpackint n2,
     mpackint n3, mpackint n4);
 mpf_class Rlamch_gmp(const char *cmach);
@@ -77,6 +82,7 @@ void Rorg2r(mpackint m, mpackint n, mpackint k, mpf_class * A, mpackint lda, mpf
 void Rlarf(const char *side, mpackint m, mpackint n, mpf_class * v, mpackint incv,
     mpf_class tau, mpf_class * C, mpackint ldc, mpf_class * work);
 void Rpotf2(const char *uplo, mpackint n, mpf_class * A, mpackint lda, mpackint *info);
+void Rpotf2Stabilized(const char *uplo, mpackint n, mpf_class * A, mpackint lda, mpackint * info, mpackint indexStart, vector<mpackint> &stabilizeIndices, vector<mpf_class> &stabilizeLambdas, double &totalLogLambda);
 void Rlaset(const char *uplo, mpackint m, mpackint n, mpf_class alpha, mpf_class beta,
     mpf_class * A, mpackint lda);
 void Rlaev2(mpf_class a, mpf_class b, mpf_class c, mpf_class * rt1,
diff --git a/src/tests.cpp b/src/tests.cpp
index fc2dfa3..e99d96e 100644
--- a/src/tests.cpp
+++ b/src/tests.cpp
@@ -10,16 +10,17 @@ using std::cout;
 using std::endl;
 
 void testCholeskyStabilize() {
+  cout.precision(20);
   Matrix A(4,4);
   Matrix L(A);
-  A.elt(0,0) = 1e40;
-  A.elt(1,1) = 1e20;
-  A.elt(2,2) = 1;
-  A.elt(3,3) = 1e-20;
-  A.elt(1,0) = 1e15;
-  A.elt(0,1) = 1e15;
-  A.elt(1,2) = 2e8;
-  A.elt(2,1) = 2e8;
+  A.elt(0,0) = 1e20;
+  A.elt(1,1) = 2e18;
+  A.elt(2,2) = 1e-3;
+  A.elt(3,3) = 1e-5;
+  A.elt(1,0) = 1e2;
+  A.elt(0,1) = 1e2;
+  A.elt(1,2) = 2e2;
+  A.elt(2,1) = 2e2;
   vector<int> updateIndices;
   Vector updateVector(L.rows);
   Real lambdaGM;
@@ -30,6 +31,15 @@ void testCholeskyStabilize() {
   cout << "updateIndices = " << updateIndices << ";\n";
   cout << "L = " << L << "\n;";
   cout << "lambdaGM = " << lambdaGM << ";\n";
+
+
+  vector<Integer> stabilizeIndices;
+  vector<Real> stabilizeLambdas;
+  choleskyDecompositionStabilized(A, L, stabilizeIndices, stabilizeLambdas);
+  cout << "A = " << A << ";\n";
+  cout << "L = " << L << ";\n";
+  cout << "stabilizeIndices = " << stabilizeIndices << ";\n";
+  cout << "stabilizeLambdas = " << stabilizeLambdas << ";\n";
 }
 
 void testLinearlyIndependentRowIndices() {
diff --git a/src/tests.h b/src/tests.h
new file mode 100644
index 0000000..1e3ad2f
--- /dev/null
+++ b/src/tests.h
@@ -0,0 +1,6 @@
+#ifndef SDP_BOOTSTRAP_TESTS_H_
+#define SDP_BOOTSTRAP_TESTS_H_
+
+void testCholeskyStabilize();
+
+#endif  // SDP_BOOTSTRAP_TESTS_H_

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/sdpb.git



More information about the debian-science-commits mailing list