qp4.c

Go to the documentation of this file.

 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "coiheader.h"
#include "comdecl.h"
 
/*  QP Model: Min (x-target)*Q*(x-target)/2 s.t.
              sum( x ) = 1;
    where x is a vector of non-negative variables of length NN
          target is a vector of +10, and
          Q is a matrix with +1 on the diagonal and 0.1 on the
            first upper and lower bi-diagonal                         */
 
#define NN 1000      /* Size of the QP model */
#define NQ (NN*2-1)
double Q[NQ];
double Target[NN];
int    Qrow[NQ];
int    Qcol[NQ];
 

int COI_CALLCONV QP_ReadMatrix( double LOWER[], double CURR[], double UPPER[], int VSTA[], int TYPE[], double RHS[],
                            int ESTA[], int COLSTA[], int ROWNO[], double VALUE[],
                            int NLFLAG[], int NUMVAR, int NUMCON, int NUMNZ, void* USRMEM )
{
  int i, j;
 /*                                                               */
 /*   Information about Variables:                                */
 /*   Default: Lower = -Inf, Curr = 0, and Upper = +inf.          */
 /*   Default: the status information in Vsta is not used.        */
 /*                                                               */
 /*   Lower bound = 0 on all variables:                           */
 /*                                                               */
   for ( i=0; i<NN; i++ )
      LOWER[i] = 0.0 ;
 /*                                                               */
 /*   Information about Constraints:                              */
 /*   Default: Rhs = 0                                            */
 /*   Default: the status information in Esta and the function    */
 /*            value in FV are not used.                          */
 /*   Default: Type: There is no default.                         */
 /*            0 = Equality,                                      */
 /*            1 = Greater than or equal,                         */
 /*            2 = Less than or equal,                            */
 /*            3 = Non binding.                                   */
 /*                                                               */
 /*   Constraint 0 (Objective)                                    */
 /*   Rhs = 0.0 and type Non binding                              */
 /*                                                               */
   TYPE[0] = 3    ;
 /*                                                               */
 /*   Constraint 1 (Sum to 1)                                     */
 /*   Rhs = 1 and type Equality                                   */
 /*                                                               */
   RHS[1]  = 1.0 ;
   TYPE[1] = 0   ;
 /*                                                               */
 /*   Information about the Jacobian. We use the standard method with */
 /*   Rowno, Value, Nlflag and Colsta and we do not use Colno.    */
 /*                                                               */
 /*   Colsta = Start of column indices (No Defaults):             */
 /*   Rowno  = Row indices                                        */
 /*   Value  = Value of derivative (by default only linear        */
 /*                                 derivatives are used)         */
 /*   Nlflag = 0 for linear and 1 for nonlinear derivative        */
 /*            (not needed for completely linear models)          */
 /*                                                               */
   j = 0;  /* counter for current nonzero   */
   for ( i=0; i<NN; i++ ) {
      COLSTA[i] = j;
      ROWNO[j]  = 0;
      NLFLAG[j] = 1;
      j++;
      ROWNO[j]  = 1;
      VALUE[j]  = 1.0;
      NLFLAG[j] = 0;
      j++;
      }
   COLSTA[NN] = j;
   return 0;
}
 

int COI_CALLCONV QP_FDEval( const double X[], double* G, double JAC[], int ROWNO, const int JACNUM[], int MODE,
                        int IGNERR, int* ERRCNT, int NUMVAR, int NUMJAC, int THREAD, void* USRMEM )
{
   int i, j, k;
   double sum;
 /*                                                                   */
 /*   Row 0: the objective function is nonlinear                      */
 /*                                                                   */
 
   if ( ROWNO == 0 ) {
 /*                                                                   */
 /*   Mode = 1 or 3. Function value: G = P * Out                      */
 /*                                                                   */
      if ( MODE == 1 || MODE == 3 ) {
         sum = 0.0;
         for ( k=0; k<NQ; k++ ) {
            i = Qrow[k]; j = Qcol[k];
            if ( i == j )
               sum += (X[i]-Target[i])*Q[k]*(X[i]-Target[i]);
            else
               sum += 2*(X[i]-Target[i])*Q[k]*(X[j]-Target[j]);
            }
         *G    = sum / 2;
         }
 /*                                                                   */
 /*   Mode = 2 or 3: Derivative values:                               */
 /*                                                                   */
      if ( MODE == 2 || MODE == 3 ) {
         for ( i=0; i<NN; i++ )
            JAC[i] = 0.0;
         for ( k=0; k<NQ; k++ ) {
            i = Qrow[k]; j = Qcol[k];
            if ( i == j )
               JAC[i] += Q[k]*(X[i]-Target[i]);
            else {
               JAC[i] += Q[k]*(X[j]-Target[j]);
               JAC[j] += Q[k]*(X[i]-Target[i]);
               }
            }
         }
      }
 /*                                                                   */
 /*   Row = 2: The row is linear and will not be called.              */
 /*                                                                   */
 
   return 0;
}
 

int COI_CALLCONV QP_2DDir( const double X[], const double* DX, double D2G[], int ROWNO,
                const int JACNUM[], int* NODRV, int NUMVAR, int NUMJAC, int THREAD, void* USRMEM )
{
int i, j, k;
/*                                                                    */
/*   Is only called for Irow = 1, the nonlinear objective             */
/*   Return H*Dx = Q*Dx in D2G                                        */
/*                                                                    */
   if ( ROWNO == 0 ) {
      for ( i=0; i<NN; i++ )
         D2G[i] = 0.0;
      for ( k=0; k<NQ; k++ ) {
         i = Qrow[k]; j = Qcol[k];
         if ( i == j )
            D2G[i] += Q[k]*DX[i];
         else {
            D2G[i] += Q[k]*DX[j];
            D2G[j] += Q[k]*DX[i];
            }
         }
      }
   return 0;
}
 

int COI_CALLCONV QP_2DLagrStr( int HSRW[], int HSCL[], int* NODRV,
                               int NUMVAR, int NUMCON, int NHESS, void* USRMEM )
{
   int k;
   for ( k=0; k<NQ; k++ ) {
      HSRW[k] = Qrow[k];
      HSCL[k] = Qcol[k];
      }
   return 0;
}
 

int COI_CALLCONV QP_2DLagrVal( const double X[], const double U[], const int HSRW[], const int HSCL[],
                double HSVL[], int* NODRV, int NUMVAR, int NUMCON, int NHESS, void* USRMEM )
{
   int k;
 
   for ( k=0; k<NQ; k++ )
      HSVL[k] = Q[k]*U[0];
   return 0;
}
 
#include "std.c"

int main(int argc, char** argv)
{
   int i,j;
 
   c_log( "Starting to execute", START );
/*
   Initialize the Q matrix
 */
   j = 0;
   for ( i=0; i<NN; i++ ) {
      Target[i] = 10.0;
      Q[j]    = 1.0;
      Qrow[j] = i;
      Qcol[j] = i;
      j++;
      if ( i < NN-1 ) {
         Q[j]    = 0.1;
         Qrow[j] = i+1;
         Qcol[j] = i;
         j++;
         }
      }
 /*
   Tell CONOPT about the sizes in the model
 */
   COI_Error +=  COIDEF_NumVar    ( CntVect, NN   );       /* NN variables                      */
   COI_Error +=  COIDEF_NumCon    ( CntVect, 2    );       /* 2 constraints                     */
   COI_Error +=  COIDEF_NumNz     ( CntVect, 2*NN );       /* 2*NN nonzeros in the Jacobian     */
   COI_Error +=  COIDEF_NumNlNz   ( CntVect, NN   );       /* NN of which are nonlinear         */
   COI_Error +=  COIDEF_NumHess   ( CntVect, NQ   );       /* NQ nonzero Hessian elements       */
   COI_Error +=  COIDEF_OptDir    ( CntVect, -1   );       /* Minimize                          */
   COI_Error +=  COIDEF_ObjCon    ( CntVect, 0    );       /* Objective is constraint 0         */
   COI_Error +=  COIDEF_DebugFV   ( CntVect, 0    );       /* 0 means no debugging, 1 each iter */
   COI_Error +=  COIDEF_StdOut    ( CntVect, 0    );       /* 1 means Allow output to StdOut    */
 /*
   Register the necessary callback routines with CONOPT
 */
   COI_Error +=  COIDEF_Message   ( CntVect, &Std_Message  );  /* Register the callback Message     */
   COI_Error +=  COIDEF_ErrMsg    ( CntVect, &Std_ErrMsg   );  /* Register the callback ErrMsg      */
   COI_Error +=  COIDEF_Status    ( CntVect, &Std_Status   );  /* Register the callback Status      */
   COI_Error +=  COIDEF_Solution  ( CntVect, &Std_Solution );  /* Register the callback Solution    */
   COI_Error +=  COIDEF_ReadMatrix( CntVect, &QP_ReadMatrix);  /* Register the callback ReadMatrix  */
   COI_Error +=  COIDEF_FDEval    ( CntVect, &QP_FDEval);      /* Register the callback FDEval      */
   COI_Error +=  COIDEF_2DDir     ( CntVect, &QP_2DDir );      /* Register the callback 2DDir       */
   COI_Error +=  COIDEF_2DLagrStr ( CntVect, &QP_2DLagrStr);   /* Register the callback 2DLagrStr   */
   COI_Error +=  COIDEF_2DLagrVal ( CntVect, &QP_2DLagrVal);   /* Register the callback 2DLagrVal   */
 
#if defined(LICENSE_INT_1) && defined(LICENSE_INT_2) && defined(LICENSE_INT_3) && defined(LICENSE_TEXT)
   COI_Error +=  COIDEF_License   ( CntVect, LICENSE_INT_1, LICENSE_INT_2, LICENSE_INT_3, LICENSE_TEXT);
#endif
 
   if ( COI_Error ) {
      printf("Skipping COI_Solve due to setup errors. COI_Error = %d\n",COI_Error);
      c_log( "Skipping Solve due to setup errors", COI_Error); }
   COI_Error =  COI_Solve        ( CntVect );           /* Optimize                      */
   printf("After solving. COI_Error =%d\n",COI_Error);
   if ( COI_Error ) {
      c_log( "Errors encountered during first solve", COI_Error); }
   else if ( stacalls == 0 || solcalls == 0 ) {
      c_log( "Status or Solution routine was not called", -1); }
   else if ( mstat != 2 || sstat != 1 ) {
      c_log( "Incorrect Model or Solver Status", -1); }
   else if ( fabs( OBJ-59978.000 ) > 0.001 ) {
      c_log( "Incorrect objective returned", -1); };
   c_log( "Successful Solve", OK );
}