Enables single precision use with configure option '--enable-single-p…

…recision'. Pardiso is only interfaced solver

configure

@@ -648,6 +648,8 @@ IPOPTAMPLINTERFACELIB_CFLAGS_NOPC
 IPOPTAMPLINTERFACELIB_LFLAGS_NOPC
 IPOPTLIB_CFLAGS_NOPC
 IPOPTLIB_LFLAGS_NOPC
+IPOPT_DOUBLE_FALSE
+IPOPT_DOUBLE_TRUE
 BUILD_SIPOPT_FALSE
 BUILD_SIPOPT_TRUE
 BUILD_LINEARSOLVERLOADER_FALSE
@@ -886,6 +888,7 @@ enable_inexact_solver
 enable_java
 enable_linear_solver_loader
 enable_sipopt
+enable_single_precision
 '
       ac_precious_vars='build_alias
 host_alias
@@ -1554,6 +1557,8 @@ Optional Features:
   --disable-linear-solver-loader
                           disable build of linear solver loader
   --disable-sipopt        disable build of sIpopt
+  --enable-single-precision
+                          use single precision in this build of IPOPT
 
 Optional Packages:
   --with-PACKAGE[=ARG]    use PACKAGE [ARG=yes]
@@ -24898,7 +24903,7 @@ else
 JAVA_TEST=Test.java
 CLASS_TEST=Test.class
 cat << \EOF > $JAVA_TEST
-/* #line 24901 "configure" */
+/* #line 24906 "configure" */
 public class Test {
 }
 EOF
@@ -25384,7 +25389,7 @@ else
 JAVA_TEST=Test.java
 CLASS_TEST=Test.class
 cat << \EOF > $JAVA_TEST
-/* #line 25387 "configure" */
+/* #line 25392 "configure" */
 public class Test {
 }
 EOF
@@ -25418,7 +25423,7 @@ JAVA_TEST=Test.java
 CLASS_TEST=Test.class
 TEST=Test
 cat << \EOF > $JAVA_TEST
-/* [#]line 25421 "configure" */
+/* [#]line 25426 "configure" */
 public class Test {
 public static void main (String args[]) {
         System.exit (0);
@@ -25736,6 +25741,29 @@ else
 fi
 
 
+########################################################################
+##                          Single Precision                          ##
+########################################################################
+
+# Check whether --enable-single-precision was given.
+if test "${enable_single_precision+set}" = set; then :
+  enableval=$enable_single_precision;
+fi
+
+
+if test "x$enable_single_precision" = "xyes" ; then
+
+$as_echo "#define IPOPT_SINGLE 1" >>confdefs.h
+
+fi
+ if test "$enable_single_precision" != yes; then
+  IPOPT_DOUBLE_TRUE=
+  IPOPT_DOUBLE_FALSE='#'
+else
+  IPOPT_DOUBLE_TRUE='#'
+  IPOPT_DOUBLE_FALSE=
+fi
+
 ########################################################################
 ##             Create links for the test source files                 ##
 ########################################################################
@@ -26454,6 +26482,10 @@ if test -z "${BUILD_SIPOPT_TRUE}" && test -z "${BUILD_SIPOPT_FALSE}"; then
   as_fn_error $? "conditional \"BUILD_SIPOPT\" was never defined.
 Usually this means the macro was only invoked conditionally." "$LINENO" 5
 fi
+if test -z "${IPOPT_DOUBLE_TRUE}" && test -z "${IPOPT_DOUBLE_FALSE}"; then
+  as_fn_error $? "conditional \"IPOPT_DOUBLE\" was never defined.
+Usually this means the macro was only invoked conditionally." "$LINENO" 5
+fi
 
 : "${CONFIG_STATUS=./config.status}"
 ac_write_fail=0

configure.ac

@@ -409,6 +409,17 @@ AC_ARG_ENABLE([sipopt],
   [use_sipopt=yes])
 AM_CONDITIONAL(BUILD_SIPOPT, [test "$use_sipopt" = yes])
 
+########################################################################
+##                          Single Precision                          ##
+########################################################################
+
+AC_ARG_ENABLE([single-precision],
+			  AC_HELP_STRING([--enable-single-precision],[use single precision in this build of IPOPT]))
+
+if test "x$enable_single_precision" = "xyes" ; then
+   AC_DEFINE([IPOPT_SINGLE],[1],[Define to 1 if using single precision floating point])
+fi
+AM_CONDITIONAL([IPOPT_DOUBLE],[test "$enable_single_precision" != yes])
 ########################################################################
 ##             Create links for the test source files                 ##
 ########################################################################

contrib/sIPOPT/src/SensIndexPCalculator.cpp

@@ -90,7 +90,7 @@ bool IndexPCalculator::ComputeP()
          {
             comp_vec = dynamic_cast<const DenseVector*>(GetRawPtr(sol_vec->GetComp(j)));
             comp_values = comp_vec->Values();
-            IpBlasDcopy(comp_vec->Dim(), comp_values, 1, col_values + curr_dim, 1);
+            IpBlasCopy(comp_vec->Dim(), comp_values, 1, col_values + curr_dim, 1);
             curr_dim += comp_vec->Dim();
          }
          cols_[col] = new PColumn(col_values);

contrib/sIPOPT/src/SensIndexSchurData.cpp

@@ -298,7 +298,7 @@ void IndexSchurData::TransMultiply(
    {
       curr_dim = v.GetCompNonConst(i)->Dim();
       curr_val = dynamic_cast<DenseVector*>(GetRawPtr(v.GetCompNonConst(i)))->Values();
-      IpBlasDcopy(curr_dim, v_vals + v_idx, 1, curr_val, 1);
+      IpBlasCopy(curr_dim, v_vals + v_idx, 1, curr_val, 1);
 
       v_idx += curr_dim;
    }

contrib/sIPOPT/src/SensStdStepCalc.cpp

@@ -131,7 +131,7 @@ bool StdStepCalculator::Step(
          delta_u_space = new DenseVectorSpace(new_du_size); // create new delta_u space
          new_delta_u = new DenseVector(GetRawPtr(ConstPtr(delta_u_space)));
          new_du_values = new_delta_u->Values();
-         IpBlasDcopy(old_delta_u->Dim(), old_delta_u->Values(), 1, new_du_values, 1);
+         IpBlasCopy(old_delta_u->Dim(), old_delta_u->Values(), 1, new_du_values, 1);
          for( Index i = 0; i < (int) x_bound_violations_idx.size(); ++i )
          {
             //   printf("i=%d, delta_u_sort[i]=%d, x_bound_viol_du[i]=%f\n", i, delta_u_sort[i], x_bound_violations_du[i]);

examples/hs071_c/hs071_c.c

@@ -151,7 +151,11 @@ int main()
    /* Set some options.  Note the following ones are only examples,
     * they might not be suitable for your problem.
     */
+#ifdef IPOPT_SINGLE
+   AddIpoptNumOption(nlp, "tol", 3.82e-6); // reduce tolerance in single precision to avoid "EXIT: Search Direction is becoming Too Small."
+#else
    AddIpoptNumOption(nlp, "tol", 1e-7);
+#endif
    AddIpoptStrOption(nlp, "mu_strategy", "adaptive");
    AddIpoptStrOption(nlp, "output_file", "ipopt.out");
 

src/Algorithm/Inexact/IpInexactNormalTerminationTester.cpp

@@ -89,7 +89,7 @@ InexactNormalTerminationTester::ETerminationTest InexactNormalTerminationTester:
 
    ETerminationTest retval = CONTINUE;
 
-   double norm2_resid = IpBlasDnrm2(ndim, resid, 1);
+   Number norm2_resid = IpBlasNrm2(ndim, resid, 1);
    Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
                   "TTNormal: iter = %d ||resid|| = %23.16e ||rhs|| = %23.16e\n", iter,  norm2_resid, norm2_rhs);
 

src/Algorithm/Inexact/IpInexactPDTerminationTester.cpp

@@ -259,7 +259,7 @@ InexactPDTerminationTester::ETerminationTest InexactPDTerminationTester::TestTer
     return retval;
     }
     */
-   Number norm2_resid = IpBlasDnrm2(ndim, resid, 1);
+   Number norm2_resid = IpBlasNrm2(ndim, resid, 1);
    Number test_ratio = norm2_resid / norm2_rhs; // Min(norm2_resid/norm2_rhs, norm2_resid);
    Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
                   "TT: test ratio %e (norm2_rhs = %e norm2_resid = %e).\n",

src/Algorithm/Inexact/IpInexactTSymScalingMethod.cpp

@@ -28,8 +28,8 @@ bool InexactTSymScalingMethod::ComputeSymTScalingFactors(
    Index         nnz,
    const ipfint* airn,
    const ipfint* ajcn,
-   const double* a,
-   double*       scaling_factors
+   const Number* a,
+   Number*       scaling_factors
 )
 {
    DBG_START_METH("InexactTSymScalingMethod::ComputeTSymScalingFactors",

src/Algorithm/Inexact/IpInexactTSymScalingMethod.hpp

@@ -44,8 +44,8 @@ class InexactTSymScalingMethod: public TSymScalingMethod
       Index         nnz,
       const ipfint* airn,
       const ipfint* ajcn,
-      const double* a,
-      double*       scaling_factors
+      const Number* a,
+      Number*       scaling_factors
    );
 
 private:

src/Algorithm/Inexact/IpIterativePardisoSolverInterface.cpp

@@ -28,10 +28,10 @@ extern "C"
 {
    int IpoptTerminationTest(
       int     n,
-      double* sol,
-      double* resid,
+      Number* sol,
+      Number* resid,
       int     iter,
-      double  norm2_rhs
+      Number  norm2_rhs
    )
    {
       fflush(stdout);
@@ -55,10 +55,10 @@ extern "C"
    void SetIpoptCallbackFunction(
       int (*IpoptFunction)(
          int     n,
-         double* x,
-         double* r,
+         Number* x,
+         Number* r,
          int     k,
-         double  b)
+         Number  b)
    );
 }
 
@@ -70,7 +70,7 @@ extern "C"
       const ipfint* MTYPE,
       const ipfint* SOLVER,
       ipfint*       IPARM,
-      double*       DPARM,
+      Number*       DPARM,
       ipfint*       ERROR
    );
 
@@ -81,17 +81,17 @@ extern "C"
       const ipfint* MTYPE,
       const ipfint* PHASE,
       const ipfint* N,
-      const double* A,
+      const Number* A,
       const ipfint* IA,
       const ipfint* JA,
       const ipfint* PERM,
       const ipfint* NRHS,
       ipfint*       IPARM,
       const ipfint* MSGLVL,
-      double*       B,
-      double*       X,
+      Number*       B,
+      Number*       X,
       ipfint*       ERROR,
-      double*       DPARM
+      Number*       DPARM
    );
 }
 
@@ -121,7 +121,7 @@ IterativePardisoSolverInterface::IterativePardisoSolverInterface(
 
    PT_ = new void* [64];
    IPARM_ = new ipfint[64];
-   DPARM_ = new double[64];
+   DPARM_ = new Number[64];
 }
 
 IterativePardisoSolverInterface::~IterativePardisoSolverInterface()
@@ -137,7 +137,7 @@ IterativePardisoSolverInterface::~IterativePardisoSolverInterface()
       ipfint NRHS = 0;
       ipfint ERROR;
       ipfint idmy;
-      double ddmy;
+      Number ddmy;
       IPOPT_PARDISO_FUNC(pardiso, PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_, &PHASE, &N, &ddmy, &idmy, &idmy, &idmy, &NRHS, IPARM_,
                                      &MSGLVL_, &ddmy, &ddmy, &ERROR, DPARM_);
       DBG_ASSERT(ERROR == 0);
@@ -213,7 +213,7 @@ bool IterativePardisoSolverInterface::InitializeImpl(
       ipfint NRHS = 0;
       ipfint ERROR;
       ipfint idmy;
-      double ddmy;
+      Number ddmy;
       IPOPT_PARDISO_FUNC(pardiso, PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_, &PHASE, &N, &ddmy, &idmy, &idmy, &idmy, &NRHS, IPARM_,
                                      &MSGLVL_, &ddmy, &ddmy, &ERROR, DPARM_);
       DBG_ASSERT(ERROR == 0);
@@ -285,7 +285,11 @@ bool IterativePardisoSolverInterface::InitializeImpl(
    IPARM_[20] = 3; // Results in better accuracy
    IPARM_[23] = 1; // parallel fac
    IPARM_[24] = 1; // parallel solve
-   IPARM_[28] = 0; // 32-bit factorization
+#ifdef IPOPT_SINGLE
+   IPARM_[28] = 1; // 32-bit factorization (single precision)
+#else
+   IPARM_[28] = 0; // 64-bit factorization (double precision)
+#endif
    IPARM_[29] = 1; // we need this for IPOPT interface
 
    IPARM_[31] = 1;  // iterative solver
@@ -318,7 +322,7 @@ ESymSolverStatus IterativePardisoSolverInterface::MultiSolve(
    const Index* ia,
    const Index* ja,
    Index        nrhs,
-   double*      rhs_vals,
+   Number*      rhs_vals,
    bool         check_NegEVals,
    Index        numberOfNegEVals
 )
@@ -344,7 +348,7 @@ ESymSolverStatus IterativePardisoSolverInterface::MultiSolve(
    return Solve(ia, ja, nrhs, rhs_vals);
 }
 
-double* IterativePardisoSolverInterface::GetValuesArrayPtr()
+Number* IterativePardisoSolverInterface::GetValuesArrayPtr()
 {
    DBG_ASSERT(initialized_);
    DBG_ASSERT(a_);
@@ -367,7 +371,7 @@ ESymSolverStatus IterativePardisoSolverInterface::InitializeStructure(
    // Make space for storing the matrix elements
    delete[] a_;
    a_ = NULL;
-   a_ = new double[nonzeros_];
+   a_ = new Number[nonzeros_];
 
    // Do the symbolic facotrization
    ESymSolverStatus retval = SymbolicFactorization(ia, ja);
@@ -404,8 +408,8 @@ static void write_iajaa_matrix(
    int          N,
    const Index* ia,
    const Index* ja,
-   double*      a_,
-   double*      rhs_vals,
+   Number*      a_,
+   Number*      rhs_vals,
    int          iter_cnt,
    int          sol_cnt
 )
@@ -476,9 +480,9 @@ ESymSolverStatus IterativePardisoSolverInterface::Factorization(
    ipfint N = dim_;
    ipfint PERM;   // This should not be accessed by Pardiso
    ipfint NRHS = 0;
-   double B;  // This should not be accessed by Pardiso in factorization
+   Number B;  // This should not be accessed by Pardiso in factorization
    // phase
-   double X;  // This should not be accessed by Pardiso in factorization
+   Number X;  // This should not be accessed by Pardiso in factorization
    // phase
    ipfint ERROR;
 
@@ -669,7 +673,7 @@ ESymSolverStatus IterativePardisoSolverInterface::Solve(
    const Index* ia,
    const Index* ja,
    Index        nrhs,
-   double*      rhs_vals
+   Number*      rhs_vals
 )
 {
    DBG_START_METH("IterativePardisoSolverInterface::Solve", dbg_verbosity);
@@ -685,8 +689,8 @@ ESymSolverStatus IterativePardisoSolverInterface::Solve(
    ipfint N = dim_;
    ipfint PERM;   // This should not be accessed by Pardiso
    ipfint NRHS = nrhs;
-   double* X = new double[nrhs * dim_];
-   double* ORIG_RHS = new double[nrhs * dim_];
+   Number* X = new Number[nrhs * dim_];
+   Number* ORIG_RHS = new Number[nrhs * dim_];
    ipfint ERROR;
 
    // Initialize solution with zero and save right hand side
@@ -858,7 +862,7 @@ ESymSolverStatus IterativePardisoSolverInterface::Solve(
       // iterates to zero
       Index nvars = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
       const Number zero = 0.;
-      IpBlasDcopy(nvars, &zero, 0, rhs_vals, 1);
+      IpBlasCopy(nvars, &zero, 0, rhs_vals, 1);
    }
    return SYMSOLVER_SUCCESS;
 }