17 List of examples

List of examples shipped in the distribution of Optimizer API for C:

Table 17.1 List of distributed examples

File

Description

acc1.c

A simple problem with one affine conic constraint (ACC)

acc2.c

A simple problem with two affine conic constraints (ACC)

blas_lapack.c

Demonstrates the MOSEK interface to BLAS/LAPACK linear algebra routines

callback.c

An example of data/progress callback

ceo1.c

A simple conic exponential problem

concurrent1.cc

Implementation of a concurrent optimizer for linear and mixed-integer problems

cqo1.c

A simple conic quadratic problem

djc1.c

A simple problem with disjunctive constraints (DJC)

errorreporting.c

Demonstrates how error reporting can be customized

feasrepairex1.c

A simple example of how to repair an infeasible problem

gp1.c

A simple geometric program (GP) in conic form

helloworld.c

A Hello World example

lo1.c

A simple linear problem

lo2.c

A simple linear problem

logistic.c

Implements logistic regression and simple log-sum-exp (CEO)

mico1.c

A simple mixed-integer conic problem

milo1.c

A simple mixed-integer linear problem

mioinitsol.c

A simple mixed-integer linear problem with an initial guess

opt_server_async.c

Uses MOSEK OptServer to solve an optimization problem asynchronously

opt_server_sync.c

Uses MOSEK OptServer to solve an optimization problem synchronously

parallel.c

Demonstrates parallel optimization using a batch method in MOSEK

parameters.c

Shows how to set optimizer parameters and read information items

pinfeas.c

Shows how to obtain and analyze a primal infeasibility certificate

portfolio_1_basic.c

Portfolio optimization - basic Markowitz model

portfolio_2_frontier.c

Portfolio optimization - efficient frontier

portfolio_3_impact.c

Portfolio optimization - market impact costs

portfolio_4_transcost.c

Portfolio optimization - transaction costs

portfolio_5_card.c

Portfolio optimization - cardinality constraints

portfolio_6_factor.c

Portfolio optimization - factor model

pow1.c

A simple power cone problem

qcqo1.c

A simple quadratically constrained quadratic problem

qo1.c

A simple quadratic problem

reoptimization.c

Demonstrate how to modify and re-optimize a linear problem

response.c

Demonstrates proper response handling

sdo1.c

A simple semidefinite problem with one matrix variable and a quadratic cone

sdo2.c

A simple semidefinite problem with two matrix variables

sdo_lmi.c

A simple semidefinite problem with an LMI using the SVEC domain.

sensitivity.c

Sensitivity analysis performed on a small linear problem

simple.c

A simple I/O example: read problem from a file, solve and write solutions

solutionquality.c

Demonstrates how to examine the quality of a solution

solvebasis.c

Demonstrates solving a linear system with the basis matrix

solvelinear.c

Demonstrates solving a general linear system

sparsecholesky.c

Shows how to find a Cholesky factorization of a sparse matrix

unicode.c

Demonstrates string conversion to Unicode

Additional examples can be found on the MOSEK website and in other MOSEK publications.

acc1.c

Listing 17.1 acc1.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      acc1.c

   Purpose :   Tutorial example for affine conic constraints.
               Models the problem:
 
               maximize c^T x
               subject to  sum(x) = 1
                           gamma >= |Gx+h|_2
 */
#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  MSKrescodee  r;
  MSKint32t i, j;

  MSKenv_t    env  = NULL;
  MSKtask_t   task = NULL;
  MSKint64t   quadDom;

  /* Input data dimensions */
  const MSKint32t n = 3,
                  k = 2;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, 0, 0, &task);

    if (r == MSK_RES_OK)
    {
      MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      /* Create n free variables */
      if (r == MSK_RES_OK)
        r = MSK_appendvars(task, n);
      if (r == MSK_RES_OK)
        r = MSK_putvarboundsliceconst(task, 0, n, MSK_BK_FR, -MSK_INFINITY, +MSK_INFINITY);

      /* Set up the objective */
      {
        double c[] = {2.0, 3.0, -1.0};

        if (r == MSK_RES_OK)
          r = MSK_putcslice(task, 0, n, c);
        if (r == MSK_RES_OK)
          r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE);
      }

      /* One linear constraint sum(x) == 1 */
      if (r == MSK_RES_OK)
        r = MSK_appendcons(task, 1);
      if (r == MSK_RES_OK)
        r = MSK_putconbound(task, 0, MSK_BK_FX, 1.0, 1.0);
      for(i = 0; i < n && r == MSK_RES_OK; i++)
        r = MSK_putaij(task, 0, i, 1.0);

      /* Append empty AFE rows for affine expression storage */
      if (r == MSK_RES_OK)
        r = MSK_appendafes(task, k + 1);

      {
        /* Fill in the affine expression storage with data */
        /* F matrix in sparse form */
        MSKint64t Fsubi[] = {1, 1, 2, 2};       /* G is placed from row 1 of F */
        MSKint32t Fsubj[] = {0, 1, 0, 2};
        double    Fval[]  = {1.5, 0.1, 0.3, 2.1};
        int       numEntries = 4;
        /* Other data */
        double h[]    = {0, 0.1};
        double gamma  = 0.03;

        /* Fill in F storage */
        if (r == MSK_RES_OK)
          r = MSK_putafefentrylist(task, numEntries, Fsubi, Fsubj, Fval);

        /* Fill in g storage */
        if (r == MSK_RES_OK)
          r = MSK_putafeg(task, 0, gamma);
        if (r == MSK_RES_OK)
          r = MSK_putafegslice(task, 1, k+1, h);
      }

      /* Define a conic quadratic domain */
      if (r == MSK_RES_OK)
        r = MSK_appendquadraticconedomain(task, k + 1, &quadDom);

      {
        /* Create the ACC */
        MSKint64t afeidx[] = {0, 1, 2};

        if (r == MSK_RES_OK)
          r = MSK_appendacc(task, 
                            quadDom,    /* Domain index */
                            k + 1,      /* Dimension */
                            afeidx,     /* Indices of AFE rows [0,...,k] */
                            NULL);      /* Ignored */
      }


      /* Begin optimization and fetching the solution */
      if (r == MSK_RES_OK)
      {
        MSKrescodee trmcode;

        /* Run optimizer */
        r = MSK_optimizetrm(task, &trmcode);

        /* Print a summary containing information
           about the solution for debugging purposes*/
        MSK_solutionsummary(task, MSK_STREAM_MSG);

        if (r == MSK_RES_OK)
        {
          MSKsolstae solsta;

          MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

          switch (solsta)
          {
            case MSK_SOL_STA_OPTIMAL:
              {
                double *xx, *doty, *activity = NULL;

                /* Fetch the solution */
                xx = calloc(n, sizeof(double));
                MSK_getxx(task,
                          MSK_SOL_ITR,    /* Request the interior solution. */
                          xx);

                printf("Optimal primal solution\n");
                for (j = 0; j < n; ++j)
                  printf("x[%d]: %e\n", j, xx[j]);

                free(xx);

                /* Fetch the doty dual of the ACC */
                doty = calloc(k + 1, sizeof(double));
                MSK_getaccdoty(task,
                               MSK_SOL_ITR,    /* Request the interior solution. */
                               0,              /* ACC index. */
                               doty);

                printf("Dual doty of the ACC\n");
                for (j = 0; j < k + 1; ++j)
                  printf("doty[%d]: %e\n", j, doty[j]);

                free(doty);

                /* Fetch the activity of the ACC */
                activity = calloc(k + 1, sizeof(double));
                MSK_evaluateacc(task,
                                MSK_SOL_ITR,    /* Request the interior solution. */
                                0,              /* ACC index. */
                                activity);

                printf("Activity of the ACC\n");
                for (j = 0; j < k + 1; ++j)
                  printf("activity[%d]: %e\n", j, activity[j]);

                free(activity);
              }
              break;
            case MSK_SOL_STA_DUAL_INFEAS_CER:
            case MSK_SOL_STA_PRIM_INFEAS_CER:
              printf("Primal or dual infeasibility certificate found.\n");
              break;
            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined. Termination code: %d.\n", trmcode);
              break;
            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          printf("Error while optimizing.\n");
        }
      }

      if (r != MSK_RES_OK)
      {
        /* In case of an error print error code and description. */
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];

        printf("An error occurred while optimizing.\n");
        MSK_getcodedesc(r,
                        symname,
                        desc);
        printf("Error %s - '%s'\n", symname, desc);
      }
    }
    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return (r);
} /* main */

acc2.c

Listing 17.2 acc2.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      acc2.c

   Purpose :   Tutorial example for affine conic constraints.
               Models the problem:
 
               maximize c^T x
               subject to  sum(x) = 1
                           gamma >= |Gx+h|_2

               This version inputs the linear constraint as an affine conic constraint.
 */
#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  MSKrescodee  r;
  MSKint32t i, j;

  MSKenv_t    env  = NULL;
  MSKtask_t   task = NULL;
  MSKint64t   zeroDom, quadDom;

  /* Input data dimensions */
  const MSKint32t n = 3,
                  k = 2;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, 0, 0, &task);

    if (r == MSK_RES_OK)
    {
      MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      /* Create n free variables */
      if (r == MSK_RES_OK)
        r = MSK_appendvars(task, n);
      if (r == MSK_RES_OK)
        r = MSK_putvarboundsliceconst(task, 0, n, MSK_BK_FR, -MSK_INFINITY, +MSK_INFINITY);

      /* Set up the objective */
      {
        double c[] = {2.0, 3.0, -1.0};

        if (r == MSK_RES_OK)
          r = MSK_putcslice(task, 0, n, c);
        if (r == MSK_RES_OK)
          r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE);
      }

      {
        /* Set AFE rows representing the linear constraint */
        if (r == MSK_RES_OK)
          r = MSK_appendafes(task, 1);
        for(i = 0; i < n && r == MSK_RES_OK; i++)
          r = MSK_putafefentry(task, 0, i, 1.0);
        if (r == MSK_RES_OK)
          r = MSK_putafeg(task, 0, -1.0);
      }

      {
        /* Set AFE rows representing the quadratic constraint */
        /* F matrix data in sparse form */
        MSKint64t Fsubi[] = {2, 2, 3, 3};       /* G is placed from row 2 of F */
        MSKint32t Fsubj[] = {0, 1, 0, 2};
        double    Fval[]  = {1.5, 0.1, 0.3, 2.1};
        int       numEntries = 4;
        /* Other data */
        double h[]    = {0, 0.1};
        double gamma  = 0.03;

        if (r == MSK_RES_OK)
          r = MSK_appendafes(task, k + 1);
        if (r == MSK_RES_OK)
          r = MSK_putafefentrylist(task, numEntries, Fsubi, Fsubj, Fval);
        if (r == MSK_RES_OK)
          r = MSK_putafeg(task, 1, gamma);
        if (r == MSK_RES_OK)
          r = MSK_putafegslice(task, 2, k+2, h);
      }

      /* Define a conic quadratic domain */
      if (r == MSK_RES_OK)
        r = MSK_appendrzerodomain(task, 1, &zeroDom);
      if (r == MSK_RES_OK)
        r = MSK_appendquadraticconedomain(task, k + 1, &quadDom);

      /* Append affine conic constraints */
      {
        /* Linear constraint */
        MSKint64t afeidx[] = {0};

        if (r == MSK_RES_OK)
          r = MSK_appendacc(task, 
                            zeroDom,    /* Domain index */
                            1,          /* Dimension */
                            afeidx,     /* Indices of AFE rows */
                            NULL);      /* Ignored */
      }

      {
        /* Quadratic constraint */
        MSKint64t afeidx[] = {1, 2, 3};

        if (r == MSK_RES_OK)
          r = MSK_appendacc(task, 
                            quadDom,    /* Domain index */
                            k + 1,      /* Dimension */
                            afeidx,     /* Indices of AFE rows */
                            NULL);      /* Ignored */
      }


      /* Begin optimization and fetching the solution */
      if (r == MSK_RES_OK)
      {
        MSKrescodee trmcode;

        /* Run optimizer */
        r = MSK_optimizetrm(task, &trmcode);

        /* Print a summary containing information
           about the solution for debugging purposes*/
        MSK_solutionsummary(task, MSK_STREAM_MSG);

        if (r == MSK_RES_OK)
        {
          MSKsolstae solsta;

          MSK_getsolsta(task, MSK_SOL_ITR, &solsta);
          MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

          switch (solsta)
          {
            case MSK_SOL_STA_OPTIMAL:
              {
                double *xx, *doty, *activity = NULL;

                /* Fetch the primal solution */
                xx = calloc(n, sizeof(double));
                MSK_getxx(task,
                          MSK_SOL_ITR,    /* Request the interior solution. */
                          xx);

                printf("Optimal primal solution\n");
                for (j = 0; j < n; ++j)
                  printf("x[%d]: %e\n", j, xx[j]);

                free(xx);

                /* Fetch the dual doty solution for the ACC */
                doty = calloc(k + 1, sizeof(double));
                MSK_getaccdoty(task,
                               MSK_SOL_ITR,    /* Request the interior solution. */
                               1,              /* ACC index of quadratic ACC. */
                               doty);

                printf("Dual doty of the quadratic ACC\n");
                for (j = 0; j < k + 1; ++j)
                  printf("doty[%d]: %e\n", j, doty[j]);

                free(doty);

                /* Fetch the activity of the ACC */
                activity = calloc(k + 1, sizeof(double));
                MSK_evaluateacc(task,
                                MSK_SOL_ITR,    /* Request the interior solution. */
                                0,              /* ACC index. */
                                activity);

                printf("Activity of the ACC\n");
                for (j = 0; j < k + 1; ++j)
                  printf("activity[%d]: %e\n", j, activity[j]);

                free(activity);                
              }
              break;
            case MSK_SOL_STA_DUAL_INFEAS_CER:
            case MSK_SOL_STA_PRIM_INFEAS_CER:
              printf("Primal or dual infeasibility certificate found.\n");
              break;
            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined. Termination code: %d.\n", trmcode);
              break;
            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          printf("Error while optimizing.\n");
        }
      }

      if (r != MSK_RES_OK)
      {
        /* In case of an error print error code and description. */
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];

        printf("An error occurred while optimizing.\n");
        MSK_getcodedesc(r,
                        symname,
                        desc);
        printf("Error %s - '%s'\n", symname, desc);
      }
    }
    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return (r);
} /* main */

blas_lapack.c

Listing 17.3 blas_lapack.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      blas_lapack.c

   Purpose: To demonstrate how to call BLAS/LAPACK routines for
            which MOSEK provides simplified interfaces.

 */
#include "mosek.h"
void print_matrix(MSKrealt* x, MSKint32t r, MSKint32t c)
{
  MSKint32t i, j;
  for (i = 0; i < r; i++)
  {
    for (j = 0; j < c; j++)
      printf("%f ", x[j * r + i]);

    printf("\n");
  }

}

int main(int argc, char*  argv[])
{

  MSKrescodee r   = MSK_RES_OK;
  MSKenv_t    env = NULL;

  const MSKint32t n = 3, m = 2, k = 3;

  MSKrealt alpha = 2.0, beta = 0.5;
  MSKrealt x[]    = {1.0, 1.0, 1.0};
  MSKrealt y[]    = {1.0, 2.0, 3.0};
  MSKrealt z[]    = {1.0, 1.0};
  MSKrealt A[]    = {1.0, 1.0, 2.0, 2.0, 3.0, 3.0};
  MSKrealt B[]    = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
  MSKrealt C[]    = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
  MSKrealt D[]    = {1.0, 1.0, 1.0, 1.0};
  MSKrealt Q[]    = {1.0, 0.0, 0.0, 2.0};
  MSKrealt v[]    = {0.0, 0.0, 0.0};

  MSKrealt xy;

  /* BLAS routines*/
  r = MSK_makeenv(&env, NULL);
  printf("n=%d m=%d k=%d\n", m, n, k);
  printf("alpha=%f\n", alpha);
  printf("beta=%f\n", beta);

  r = MSK_dot(env, n, x, y, &xy);
  printf("dot results= %f r=%d\n", xy, r);


  print_matrix(x, 1, n);
  print_matrix(y, 1, n);

  r = MSK_axpy(env, n, alpha, x, y);
  puts("axpy results is");
  print_matrix(y, 1, n);


  r = MSK_gemv(env, MSK_TRANSPOSE_NO, m, n, alpha, A, x, beta, z);
  printf("gemv results is (r=%d) \n", r);
  print_matrix(z, 1, m);

  r = MSK_gemm(env, MSK_TRANSPOSE_NO, MSK_TRANSPOSE_NO, m, n, k, alpha, A, B, beta, C);
  printf("gemm results is (r=%d) \n", r);
  print_matrix(C, m, n);

  r = MSK_syrk(env, MSK_UPLO_LO, MSK_TRANSPOSE_NO, m, k, 1., A, beta, D);
  printf("syrk results is (r=%d) \n", r);
  print_matrix(D, m, m);

  /* LAPACK routines*/

  r = MSK_potrf(env, MSK_UPLO_LO, m, Q);
  printf("potrf results is (r=%d) \n", r);
  print_matrix(Q, m, m);

  r = MSK_syeig(env, MSK_UPLO_LO, m, Q, v);
  printf("syeig results is (r=%d) \n", r);
  print_matrix(v, 1, m);

  r = MSK_syevd(env, MSK_UPLO_LO, m, Q, v);
  printf("syevd results is (r=%d) \n", r);
  print_matrix(v, 1, m);
  print_matrix(Q, m, m);

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return r;
}

callback.c

Listing 17.4 callback.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      callback.c

   Purpose:   To demonstrate how to use the progress
              callback.

              Compile and link the file with  MOSE, then
              use as follows:

              callback psim 25fv47.mps
              callback dsim 25fv47.mps
              callback intpnt 25fv47.mps

              The first argument tells which optimizer to use
              i.e. psim is primal simplex, dsim is dual simplex
              and intpnt is interior-point.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "mosek.h"

#define CONT_EXAMPLE

/* Type of user data used in callback functions */
typedef struct {
  double maxtime;   // Time limit for the optimizer
  int numvars;
  double *xx;       // For storing the integer solution when mixed integer optimizer is used
} cbdata;
typedef cbdata* cbdata_t;

/* Note: This function is declared using MSKAPI,
         so the correct calling convention is
         employed. */
static int MSKAPI usercallback(MSKtask_t            task,
                               MSKuserhandle_t      handle,
                               MSKcallbackcodee     caller,
                               const MSKrealt  * douinf,
                               const MSKint32t * intinf,
                               const MSKint64t * lintinf)
{
  cbdata_t data = (cbdata_t) handle;
  double maxtime = data->maxtime;
  MSKrescodee r;

  switch (caller)
  {
    case MSK_CALLBACK_BEGIN_INTPNT:
      printf("Starting interior-point optimizer\n");
      break;
    case MSK_CALLBACK_INTPNT:
      printf("Iterations: %-3d  Time: %6.2f(%.2f)  ",
             intinf[MSK_IINF_INTPNT_ITER],
             douinf[MSK_DINF_OPTIMIZER_TIME],
             douinf[MSK_DINF_INTPNT_TIME]);

      printf("Primal obj.: %-18.6e  Dual obj.: %-18.6e\n",
             douinf[MSK_DINF_INTPNT_PRIMAL_OBJ],
             douinf[MSK_DINF_INTPNT_DUAL_OBJ]);
      break;
    case MSK_CALLBACK_END_INTPNT:
      printf("Interior-point optimizer finished.\n");
      break;
    case MSK_CALLBACK_BEGIN_PRIMAL_SIMPLEX:
      printf("Primal simplex optimizer started.\n");
      break;
    case MSK_CALLBACK_UPDATE_PRIMAL_SIMPLEX:
      printf("Iterations: %-3d  ",
             intinf[MSK_IINF_SIM_PRIMAL_ITER]);
      printf("  Elapsed time: %6.2f(%.2f)\n",
             douinf[MSK_DINF_OPTIMIZER_TIME],
             douinf[MSK_DINF_SIM_TIME]);
      printf("Obj.: %-18.6e\n",
             douinf[MSK_DINF_SIM_OBJ]);
      break;
    case MSK_CALLBACK_END_PRIMAL_SIMPLEX:
      printf("Primal simplex optimizer finished.\n");
      break;
    case MSK_CALLBACK_BEGIN_DUAL_SIMPLEX:
      printf("Dual simplex optimizer started.\n");
      break;
    case MSK_CALLBACK_UPDATE_DUAL_SIMPLEX:
      printf("Iterations: %-3d  ", intinf[MSK_IINF_SIM_DUAL_ITER]);
      printf("  Elapsed time: %6.2f(%.2f)\n",
             douinf[MSK_DINF_OPTIMIZER_TIME],
             douinf[MSK_DINF_SIM_TIME]);
      printf("Obj.: %-18.6e\n", douinf[MSK_DINF_SIM_OBJ]);
      break;
    case MSK_CALLBACK_END_DUAL_SIMPLEX:
      printf("Dual simplex optimizer finished.\n");
      break;
    case MSK_CALLBACK_NEW_INT_MIO:
      printf("New integer solution has been located.\n");

      r = MSK_getxx(task, MSK_SOL_ITG, data->xx);
      if (r == MSK_RES_OK) {
        int i;
        printf("xx = ");
        for (i = 0; i < data->numvars; i++) printf("%lf ", data->xx[i]);
        printf("\nObj.: %f\n", douinf[MSK_DINF_MIO_OBJ_INT]);
      }
    default:
      break;
  }

  if (douinf[MSK_DINF_OPTIMIZER_TIME] >= maxtime)
  {
    /* mosek is spending too much time.
       Terminate it. */
    return (1);
  }

  return (0);
} /* usercallback */

static void MSKAPI printtxt(void       *info,
                            const char *buffer)
{ printf("%s", buffer);
} /* printtxt */

int main(int argc, const char *argv[])
{
  double    maxtime,
            *xx, *y;
  MSKrescodee r, trmr;
  int       j, i, numcon, numvar;
  FILE      *f;
  MSKenv_t  env;
  MSKtask_t task;
  cbdata    data;

#ifdef MIO_EXAMPLE
  /* To test the mixed-integer solution callback. */
  const char * slvr = "";
  const char * filename = "../data/milo1.lp";
#endif

#ifdef CONT_EXAMPLE
  /* To test interior-point/simplex optimizers on continuous problems */
  const char * slvr = "intpnt";
  const char * filename = "../data/25fv47.mps";
#endif

  printf("callback example\n");

  if (argc < 3)
  {
    printf("Usage: callback ( psim | dsim | intpnt ) filename\n");
    if (argc == 2)
      slvr = argv[1];
  }
  else
  {
    slvr = argv[1];
    filename = argv[2];
  }

  /* Create mosek environment. */
  /* "" should be replaced by NULL in a NON DEBUG setting. */
  r = MSK_makeenv(&env, ""); 

  /* Check the return code. */
  if (r == MSK_RES_OK)
  {
    /* Create an (empty) optimization task. */
    r = MSK_makeemptytask(env, &task);

    if (r == MSK_RES_OK)
    {
      MSK_linkfunctotaskstream(task, MSK_STREAM_MSG, NULL, printtxt);
      MSK_linkfunctotaskstream(task, MSK_STREAM_ERR, NULL, printtxt);
    }

    /* Specifies that data should be read from the file.
     */

    if (r == MSK_RES_OK)
      r = MSK_readdata(task, filename);

    if (r == MSK_RES_OK)
    {
      if (0 == strcmp(slvr, "psim"))
        MSK_putintparam(task, MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_PRIMAL_SIMPLEX);
      else  if (0 == strcmp(slvr, "dsim"))
        MSK_putintparam(task, MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_DUAL_SIMPLEX);
      else  if (0 == strcmp(slvr, "intpnt"))
        MSK_putintparam(task, MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_INTPNT);


      /* Tell mosek about the call-back function. */
      data.maxtime = 0.05;
      MSK_getnumvar(task, &data.numvars);
      data.xx = MSK_callocenv(env, data.numvars, sizeof(double));

      MSK_putcallbackfunc(task,
                          usercallback,
                          (void *) &data);

      /* Turn all MOSEK logging off. */
      MSK_putintparam(task,
                      MSK_IPAR_LOG,
                      0);

      r = MSK_optimizetrm(task, &trmr);
      printf("Return code - %d, termination code - %d\n", r, trmr);

      MSK_solutionsummary(task, MSK_STREAM_MSG);
      MSK_freeenv(env, data.xx);
    }


    {
      MSKrescodee lr = MSK_deletetask(&task);

      if ( r==MSK_RES_OK ) 
        r = lr;  
    }  
  }

  {
    MSKrescodee lr=MSK_deleteenv(&env);;

    if ( r==MSK_RES_OK ) 
      r = lr;  
  }  

  printf("Return code - %d\n", r);
  return (0);
} /* main */

ceo1.c

Listing 17.5 ceo1.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      ceo1.c

   Purpose:   To demonstrate how to solve a small conic exponential
              optimization problem using the MOSEK API.
 */
#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  MSKrescodee  r;

  const MSKint32t numvar = 3,
                  numcon = 1;
  const MSKint64t numafe = 3,
                  f_nnz  = 3;

  MSKboundkeye bkc = MSK_BK_FX;
  double       blc = 1.0;
  double       buc = 1.0;

  MSKboundkeye bkx[] = {MSK_BK_FR,
                        MSK_BK_FR,
                        MSK_BK_FR
                       };
  double       blx[] = {-MSK_INFINITY,
                        -MSK_INFINITY,
                        -MSK_INFINITY
                       };
  double       bux[] = { +MSK_INFINITY,
                         +MSK_INFINITY,
                         +MSK_INFINITY
                       };

  double       c[]   = {1.0,
                        1.0,
                        0.0
                       };

  double    a[]     = {1.0, 1.0, 1.0};
  MSKint32t asub[]  = {0, 1, 2};

  MSKint64t   afeidx[] = {0, 1, 2};
  MSKint32t   varidx[] = {0, 1, 2};
  MSKrealt    f_val[]  = {1, 1, 1};

  MSKint64t   domidx = 0;

  MSKint32t   i, j, csub[3];

  MSKenv_t    env  = NULL;
  MSKtask_t   task = NULL;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, numcon, numvar, &task);

    if (r == MSK_RES_OK)
    {
      MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      /* Append 'numcon' empty constraints.
      The constraints will initially have no bounds. */
      if (r == MSK_RES_OK)
        r = MSK_appendcons(task, numcon);

      /* Append 'numvar' variables.
      The variables will initially be fixed at zero (x=0). */
      if (r == MSK_RES_OK)
        r = MSK_appendvars(task, numvar);

      /* Append 'numafe' affine expressions.
      The affine expressions will initially be empty. */
      if (r == MSK_RES_OK)
        r = MSK_appendafes(task, numafe);


      /* Set up the linear part */
      if (r == MSK_RES_OK)
        r = MSK_putcslice(task, 0, numvar, c);
      if (r == MSK_RES_OK)
        r = MSK_putarow(task, 0, numvar, asub, a);
      if (r == MSK_RES_OK)
        r = MSK_putconbound(task, 0, bkc, blc, buc);
      if (r == MSK_RES_OK)
        r = MSK_putvarboundslice(task, 0, numvar, bkx, blx, bux);

      if (r == MSK_RES_OK)
      {
        r = MSK_putafefentrylist(task, f_nnz, afeidx, varidx, f_val);
        
        if (r == MSK_RES_OK)
          r = MSK_appendprimalexpconedomain(task, &domidx);
        
        if (r == MSK_RES_OK)
          r = MSK_appendaccseq(task, domidx, numafe, 0, NULL);
      }

      if (r == MSK_RES_OK)
      {
        MSKrescodee trmcode;

        /* Run optimizer */
        r = MSK_optimizetrm(task, &trmcode);

        /* Print a summary containing information
           about the solution for debugging purposes*/
        MSK_solutionsummary(task, MSK_STREAM_MSG);

        if (r == MSK_RES_OK)
        {
          MSKsolstae solsta;

          MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

          switch (solsta)
          {
            case MSK_SOL_STA_OPTIMAL:
              {
                double *xx = NULL;

                xx = calloc(numvar, sizeof(double));
                if (xx)
                {
                  MSK_getxx(task,
                            MSK_SOL_ITR,    /* Request the interior solution. */
                            xx);

                  printf("Optimal primal solution\n");
                  for (j = 0; j < numvar; ++j)
                    printf("x[%d]: %e\n", j, xx[j]);
                }
                else
                {
                  r = MSK_RES_ERR_SPACE;
                }
                free(xx);
              }
              break;
            case MSK_SOL_STA_DUAL_INFEAS_CER:
            case MSK_SOL_STA_PRIM_INFEAS_CER:
              printf("Primal or dual infeasibility certificate found.\n");
              break;
            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined. Termination code: %d.\n", trmcode);
              break;
            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          printf("Error while optimizing.\n");
        }
      }

      if (r != MSK_RES_OK)
      {
        /* In case of an error print error code and description. */
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];

        printf("An error occurred while optimizing.\n");
        MSK_getcodedesc(r,
                        symname,
                        desc);
        printf("Error %s - '%s'\n", symname, desc);
      }
    }
    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return (r);
} /* main */

concurrent1.cc

Listing 17.6 concurrent1.cc Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      concurrent1.cc

   Purpose: Demonstrates a simple implementation of a concurrent optimizer.
            
            The concurrent optimizer starts a few parallel optimizations
            of the same problem using different algorithms, and reports
            a solution when the first optimizer is ready.

            This example also demonstrates how to define a simple callback handler
            that stops the optimizer when requested.

            For cross-platform compatibility this is C++ code requiring C++11.
            This is due to the use of std::thread. The rest is pure C. It should be easily  
            adapted to C and any threading model by only modifying the "optimize"
            function.
*/
#include<stdio.h>
#include<stdlib.h>
#include<mosek.h>
#include<thread>

/* This function prints log output from MOSEK to the terminal. */
static void MSKAPI printstr(void       *handle,
                            const char str[])
{
  printf("%s", str);
}

/**
   Defines a Mosek callback function whose only function
   is to indicate if the optimizer should be stopped.
 */
int stop = 0;
int firstStop = -1;
MSKint32t MSKAPI callback (MSKtask_t task, MSKuserhandle_t usrptr, MSKcallbackcodee caller,
                           const MSKrealt * douinf, const MSKint32t * intinf, const MSKint64t * lintinf)
{
  return stop;
}


/** Takes a list of tasks and optimizes then in parallel. The
    response code and termination code from each optimization is
    stored in ``res`` and ``trm``.

    When one task completes with rescode == ok, others are terminated.

    Returns the index of the first optimizer which returned without error.
    Whether or not this is the task firstOK contains the
    most valuable answer, is for the caller to decide.
    Otherwise returns -1.
 */
void runTask(int         num, 
             MSKtask_t   task,
             MSKrescodee *res,
             MSKrescodee *trm) 
{
  *res = MSK_optimizetrm(task, trm);
  if (*res == MSK_RES_OK) {
    if (!stop) firstStop = num;
    stop = 1;
  }
}

int optimize(int         n,
             MSKtask_t   *tasks,
             MSKrescodee *res,
             MSKrescodee *trm)
{
  int i;
  std::thread * jobs = new std::thread[n];

  // Set a callback function and start optimization
  for(i = 0; i < n; ++i) {
    MSK_putcallbackfunc(tasks[i], callback, NULL);
    res[i] = trm[i] = MSK_RES_ERR_UNKNOWN;
    jobs[i] = std::thread(runTask, i, tasks[i], &(res[i]), &(trm[i]));
  }

  // Join all threads
  for(i = 0; i < n; ++i) jobs[i].join();
  delete[] jobs;

  // For debugging, print res and trm codes for all optimizers
  for(i = 0; i < n; ++i)
    printf("Optimizer  %d   res %d   trm %d\n", i, res[i], trm[i]);

  return firstStop;
}


/** 
    Given a continuous task, set up jobs to optimize it 
    with a list of different solvers.

    Returns an index, corresponding to the optimization
    task that is returned as winTask. This is the task
    with the best possible status of those that finished.
    If none task is considered successfull returns -1.
 */
int optimizeconcurrent(MSKtask_t            task, 
                       int                  n,
                       MSKoptimizertypee    *optimizers,
                       MSKtask_t            *winTask,
                       MSKrescodee          *winTrm,
                       MSKrescodee          *winRes)
{
  MSKtask_t   *tasks = new MSKtask_t[n];
  MSKrescodee *res   = new MSKrescodee[n];
  MSKrescodee *trm   = new MSKrescodee[n];

  // Clone tasks and choose various optimizers
  for (int i = 0; i < n; ++i)
  {
    MSK_clonetask(task, &(tasks[i]));
    MSK_putintparam(tasks[i], MSK_IPAR_OPTIMIZER, optimizers[i]);
  }

  // Solve tasks in parallel
  int firstOK = optimize(n, tasks, res, trm);

  if (firstOK >= 0) 
  {
    *winTask  = tasks[firstOK];
    *winTrm   = trm[firstOK]; 
    *winRes   = res[firstOK];
  }
  else
  {
    *winTask  = NULL; 
    *winTrm   = MSK_RES_ERR_UNKNOWN; 
    *winRes   = MSK_RES_ERR_UNKNOWN;        
  }

  // Cleanup 
  for (int i = 0; i < n; ++i)
    if (i != firstOK) MSK_deletetask(&(tasks[i]));

  delete[] tasks; delete[] res; delete[] trm;
  return firstOK;
}
 

/** 
    Given a mixed-integer task, set up jobs to optimize it 
    with different values of seed. That will lead to
    different execution paths of the optimizer.

    Returns an index, corresponding to the optimization
    task that is returned as winTask. This is the task
    with the best value of the objective function.
    If none task is considered successful returns -1.

    Typically, the input task would contain a time limit. The two
    major scenarios are:
    1. Some clone ends before time limit - then it has optimum.
    2. All clones reach time limit - pick the one with best objective.
 */
int optimizeconcurrentMIO(MSKtask_t            task, 
                          int                  n,
                          int                  *seeds,
                          MSKtask_t            *winTask,
                          MSKrescodee          *winTrm,
                          MSKrescodee          *winRes)
{
  MSKtask_t   *tasks = new MSKtask_t[n];
  MSKrescodee *res   = new MSKrescodee[n];
  MSKrescodee *trm   = new MSKrescodee[n];
  *winTask  = NULL; 
  *winTrm   = MSK_RES_ERR_UNKNOWN; 
  *winRes   = MSK_RES_ERR_UNKNOWN;        
  int bestPos = -1;

  // Clone tasks and choose various optimizers
  for (int i = 0; i < n; ++i)
  {
    MSK_clonetask(task, &(tasks[i]));
    MSK_putintparam(tasks[i], MSK_IPAR_MIO_SEED, seeds[i]);
  }

  // Solve tasks in parallel
  int firstOK = optimize(n, tasks, res, trm);

  if (firstOK >= 0) 
  {
    // Pick the task that ended with res = ok
    // and contains an integer solution with best objective value
    MSKobjsensee sense;
    double bestObj;

    MSK_getobjsense(task, &sense);
    bestObj = (sense == MSK_OBJECTIVE_SENSE_MINIMIZE) ? 1.0e+10 : -1.0e+10;

    for (int i = 0; i < n; ++i) {
      double priObj;
      MSK_getprimalobj(tasks[i], MSK_SOL_ITG, &priObj);
      printf("%d      %f\n", i, priObj);
    }

    for (int i = 0; i < n; ++i) {
      double priObj;
      MSKsolstae solsta;
      MSK_getprimalobj(tasks[i], MSK_SOL_ITG, &priObj);
      MSK_getsolsta(tasks[i], MSK_SOL_ITG, &solsta);

      if ((res[i] == MSK_RES_OK) &&
          (solsta == MSK_SOL_STA_PRIM_FEAS ||
           solsta == MSK_SOL_STA_INTEGER_OPTIMAL) &&
          ((sense == MSK_OBJECTIVE_SENSE_MINIMIZE) ? 
              (priObj < bestObj) : (priObj > bestObj) ) )
      {
        bestObj = priObj;
        bestPos = i;
      }
    }

    if (bestPos != -1)
    {
      *winTask  = tasks[bestPos];
      *winTrm   = trm[bestPos]; 
      *winRes   = res[bestPos];
    }
  }

  // Cleanup 
  for (int i = 0; i < n; ++i)
    if (i != bestPos) MSK_deletetask(&(tasks[i]));

  delete[] tasks; delete[] res; delete[] trm;
  return bestPos;
}


/** 
   This is an example of how one can use the methods
       optimizeconcurrent
       optimizeconcurrentMIO

   argv[0] : name of file with input problem
   argv[1]: (optional) time limit
 */
int main(int argc, char **argv)
{

  MSKenv_t env;
  MSKtask_t task;
  int numint;

  MSKtask_t t;
  MSKrescodee res, trm;
  int idx;

  MSK_makeenv(&env, NULL);
  MSK_makeemptytask(env, &task);
  if (argc >= 2)
    MSK_readdata(task, argv[1]);
  else
    MSK_readdata(task, "../data/25fv47.mps");

  if (argc >= 3) 
    MSK_putdouparam(task, MSK_DPAR_OPTIMIZER_MAX_TIME, atof(argv[2]));

  MSK_getnumintvar(task, &numint);

  if (numint == 0) 
  {
    /* If the problem is continuous
       optimize it with three continuous optimizers.
      (Simplex will fail for non-linear problems)
    */
    MSKoptimizertypee optimizers[3] = { 
      MSK_OPTIMIZER_CONIC,
      MSK_OPTIMIZER_DUAL_SIMPLEX,
      MSK_OPTIMIZER_PRIMAL_SIMPLEX
    };

    idx = optimizeconcurrent(task, 3, optimizers, &t, &trm, &res);
  }          
  else
  {
    /* Mixed-integer problem.
       Try various seeds.
    */
    int seeds[3] = { 42, 13, 71749373 };

    idx = optimizeconcurrentMIO(task, 3, seeds, &t, &trm, &res);
  }          

  // Check results and print the best answer
  if (idx >= 0) 
  {
    printf("Result from optimizer with index %d: res %d  trm %d\n", idx, res, trm);
    MSK_linkfunctotaskstream(t, MSK_STREAM_LOG, NULL, printstr);
    MSK_optimizersummary(t, MSK_STREAM_LOG);
    MSK_solutionsummary(t, MSK_STREAM_LOG);
    MSK_deletetask(&t);
  }
  else 
  {
    printf("All optimizers failed.");
  }

  MSK_deletetask(&task);
  MSK_deleteenv(&env);
  return 0;
}

cqo1.c

Listing 17.7 cqo1.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      cqo1.c

   Purpose:   To demonstrate how to solve a small conic quadratic
              optimization problem using the MOSEK API.
 */
#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  MSKrescodee  r;

  const MSKint32t numvar = 6,
                  numcon = 1;
  const MSKint64t numafe = 6,
                  numacc = 2,
                  f_nnz  = 6;

  MSKboundkeye bkc[] = { MSK_BK_FX };
  double       blc[] = { 1.0 };
  double       buc[] = { 1.0 };

  MSKboundkeye bkx[] = {MSK_BK_LO,
                        MSK_BK_LO,
                        MSK_BK_LO,
                        MSK_BK_FR,
                        MSK_BK_FR,
                        MSK_BK_FR
                       };
  double       blx[] = {0.0,
                        0.0,
                        0.0,
                        -MSK_INFINITY,
                        -MSK_INFINITY,
                        -MSK_INFINITY
                       };
  double       bux[] = { +MSK_INFINITY,
                         +MSK_INFINITY,
                         +MSK_INFINITY,
                         +MSK_INFINITY,
                         +MSK_INFINITY,
                         +MSK_INFINITY
                       };

  double       c[]   = {0.0,
                        0.0,
                        0.0,
                        1.0,
                        1.0,
                        1.0
                       };

  MSKint32t   aptrb[] = {0, 1, 2, 3, 3, 3},
              aptre[] = {1, 2, 3, 3, 3, 3},
              asub[]  = {0, 0, 0, 0};
  double      aval[]  = {1.0, 1.0, 2.0};

  MSKint64t   afeidx[] = {0, 1, 2, 3, 4, 5};
  MSKint32t   varidx[] = {3, 0, 1, 4, 5, 2};
  MSKrealt    f_val[]  = {1, 1, 1, 1, 1, 1};

  MSKint64t   domidx[] = {0, 0};

  MSKint32t   i, j, csub[3];

  MSKenv_t    env  = NULL;
  MSKtask_t   task = NULL;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, numcon, numvar, &task);

    if (r == MSK_RES_OK)
    {
      MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      /* Append 'numcon' empty constraints.
      The constraints will initially have no bounds. */
      if (r == MSK_RES_OK)
        r = MSK_appendcons(task, numcon);

      /* Append 'numvar' variables.
      The variables will initially be fixed at zero (x=0). */
      if (r == MSK_RES_OK)
        r = MSK_appendvars(task, numvar);

      /* Append 'numafe' affine expressions.
      The affine expressions will initially be empty. */
      if (r == MSK_RES_OK)
        r = MSK_appendafes(task, numafe);


      for (j = 0; j < numvar && r == MSK_RES_OK; ++j)
      {
        /* Set the linear term c_j in the objective.*/
        if (r == MSK_RES_OK)
          r = MSK_putcj(task, j, c[j]);

        /* Set the bounds on variable j.
        blx[j] <= x_j <= bux[j] */
        if (r == MSK_RES_OK)
          r = MSK_putvarbound(task,
                              j,           /* Index of variable.*/
                              bkx[j],      /* Bound key.*/
                              blx[j],      /* Numerical value of lower bound.*/
                              bux[j]);     /* Numerical value of upper bound.*/

        /* Input column j of A */
        if (r == MSK_RES_OK)
          r = MSK_putacol(task,
                          j,                 /* Variable (column) index.*/
                          aptre[j] - aptrb[j], /* Number of non-zeros in column j.*/
                          asub + aptrb[j],   /* Pointer to row indexes of column j.*/
                          aval + aptrb[j]);  /* Pointer to Values of column j.*/

      }

      /* Set the bounds on constraints.
       for i=1, ...,numcon : blc[i] <= constraint i <= buc[i] */
      for (i = 0; i < numcon && r == MSK_RES_OK; ++i)
        r = MSK_putconbound(task,
                            i,           /* Index of constraint.*/
                            bkc[i],      /* Bound key.*/
                            blc[i],      /* Numerical value of lower bound.*/
                            buc[i]);     /* Numerical value of upper bound.*/

      if (r == MSK_RES_OK)
      {
        /* Set the non-zero entries of the F matrix */
        r = MSK_putafefentrylist(task, f_nnz, afeidx, varidx, f_val); 

        /* Append quadratic cone domain */
        if (r == MSK_RES_OK)
          r = MSK_appendquadraticconedomain(task, 3, domidx);
        /* Append rotated quadratic cone domain */        
        if (r == MSK_RES_OK)
          r = MSK_appendrquadraticconedomain(task, 3, domidx+1);

        /* Append two ACCs made up of the AFEs and the domains defined above. */
        if (r == MSK_RES_OK)
          r = MSK_appendaccsseq(task, numacc, domidx, numafe, afeidx[0], NULL);
      }

      if (r == MSK_RES_OK)
      {
        MSKrescodee trmcode;

        /* Run optimizer */
        r = MSK_optimizetrm(task, &trmcode);


        /* Print a summary containing information
           about the solution for debugging purposes*/
        MSK_solutionsummary(task, MSK_STREAM_MSG);

        if (r == MSK_RES_OK)
        {
          MSKsolstae solsta;

          MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

          switch (solsta)
          {
            case MSK_SOL_STA_OPTIMAL:
              {
                double *xx = NULL;

                xx = calloc(numvar, sizeof(double));
                if (xx)
                {
                  MSK_getxx(task,
                            MSK_SOL_ITR,    /* Request the interior solution. */
                            xx);

                  printf("Optimal primal solution\n");
                  for (j = 0; j < numvar; ++j)
                    printf("x[%d]: %e\n", j, xx[j]);
                }
                else
                {
                  r = MSK_RES_ERR_SPACE;
                }
                free(xx);
              }
              break;
            case MSK_SOL_STA_DUAL_INFEAS_CER:
            case MSK_SOL_STA_PRIM_INFEAS_CER:
              printf("Primal or dual infeasibility certificate found.\n");
              break;
            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined. Termination code: %d.\n", trmcode);
              break;
            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          printf("Error while optimizing.\n");
        }
      }

      if (r != MSK_RES_OK)
      {
        /* In case of an error print error code and description. */
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];

        printf("An error occurred while optimizing.\n");
        MSK_getcodedesc(r,
                        symname,
                        desc);
        printf("Error %s - '%s'\n", symname, desc);
      }
    }
    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return (r);
} /* main */

djc1.c

Listing 17.8 djc1.c Click here to download.
////
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      djc1.c
//
//   Purpose: Demonstrates how to solve the problem with two disjunctions:
//
//      minimize    2x0 + x1 + 3x2 + x3
//      subject to   x0 + x1 + x2 + x3 >= -10
//                  (x0-2x1<=-1 and x2=x3=0) or (x2-3x3<=-2 and x1=x2=0)
//                  x0=2.5 or x1=2.5 or x2=2.5 or x3=2.5
////
#include <stdio.h>
#include "mosek.h"

/* This function prints log output from MOSEK to the terminal. */
static void MSKAPI printstr(void       *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  MSKenv_t           env  = NULL;
  MSKtask_t          task = NULL;
  MSKrescodee        r = MSK_RES_OK;
  MSKint32t          i, j, numvar;
  MSKint64t          k, l, numafe, numdjc;
  MSKint64t          zero1, zero2, rminus1;  // Domain indices

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, 0, 0, &task);

    if (r == MSK_RES_OK)
    {
      // Append free variables
      numvar = 4;
      r = MSK_appendvars(task, numvar);
      if (r == MSK_RES_OK) 
        MSK_putvarboundsliceconst(task, 0, numvar, MSK_BK_FR, -MSK_INFINITY, MSK_INFINITY);
    }

    if (r == MSK_RES_OK)
    {
      // The linear part: the linear constraint
      const MSKint32t idx[] = {0, 1, 2, 3};
      const MSKrealt  val[] = {1, 1, 1, 1};
      
      r = MSK_appendcons(task, 1);
      if (r == MSK_RES_OK) MSK_putarow(task, 0, 4, idx, val);
      if (r == MSK_RES_OK) MSK_putconbound(task, 0, MSK_BK_LO, -10.0, -10.0);
    }

    if (r == MSK_RES_OK)
    {
      // The linear part: objective
      const MSKint32t idx[] = {0, 1, 2, 3};
      const MSKrealt  val[] = {2, 1, 3, 1};

      r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MINIMIZE);
      if (r == MSK_RES_OK) MSK_putclist(task, 4, idx, val);
    }

    // Fill in the affine expression storage F, g
    if (r == MSK_RES_OK)
    {
      numafe = 10;
      r = MSK_appendafes(task, numafe);
    }

    if (r == MSK_RES_OK)
    {
      const MSKint64t fafeidx[] = {0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9};
      const MSKint32t fvaridx[] = {0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3};
      const MSKrealt  fval[]    = {1.0, -2.0, 1.0, -3.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
      const MSKrealt  g[]       = {1.0, 2.0, 0.0, 0.0, 0.0, 0.0, -2.5, -2.5, -2.5, -2.5};

      r = MSK_putafefentrylist(task, 12, fafeidx, fvaridx, fval);
      if (r == MSK_RES_OK)
        r = MSK_putafegslice(task, 0, numafe, g);
    }

    if (r == MSK_RES_OK)
    {
      // Create domains
      MSK_appendrzerodomain(task, 1, &zero1);
      MSK_appendrzerodomain(task, 2, &zero2);
      MSK_appendrminusdomain(task, 1, &rminus1);
    }

    if (r == MSK_RES_OK)
    {
      // Append disjunctive constraints
      numdjc = 2;
      r = MSK_appenddjcs(task, numdjc);
    }

    if (r == MSK_RES_OK)
    {
      // First disjunctive constraint
      const MSKint64t domidxlist[] = {rminus1, zero2, rminus1, zero2};
      const MSKint64t afeidxlist[] = {0, 4, 5, 1, 2, 3};
      const MSKint64t termsizelist[] = {2, 2};

      r = MSK_putdjc(task, 
                     0,                      // DJC index
                     4, domidxlist,
                     6, afeidxlist, 
                     NULL,                   // Unused
                     2, termsizelist);
    }

    if (r == MSK_RES_OK)
    {
      // Second disjunctive constraint
      const MSKint64t domidxlist[] = {zero1, zero1, zero1, zero1};
      const MSKint64t afeidxlist[] = {6, 7, 8, 9};
      const MSKint64t termsizelist[] = {1, 1, 1, 1};

      r = MSK_putdjc(task, 
                     1,                      // DJC index
                     4, domidxlist,
                     4, afeidxlist, 
                     NULL,                   // Unused
                     4, termsizelist);
    }

    // Useful for debugging
    if (r == MSK_RES_OK)
    {
      // Write a human-readable file
      r = MSK_writedata(task, "djc.ptf");
      // Directs the log task stream to the 'printstr' function. 
      if (r == MSK_RES_OK)
        r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);
    }

    // Solve the problem
    if (r == MSK_RES_OK)
    {
      MSKrescodee trmcode;

      r = MSK_optimizetrm(task, &trmcode);
 
      /* Print a summary containing information
         about the solution for debugging purposes. */
      MSK_solutionsummary(task, MSK_STREAM_LOG);

      if (r == MSK_RES_OK)
      {
        MSKsolstae solsta;

        if (r == MSK_RES_OK)
          r = MSK_getsolsta(task,
                            MSK_SOL_ITG,
                            &solsta);
        switch (solsta)
        {
          case MSK_SOL_STA_INTEGER_OPTIMAL:
            {
              double *xx = (double*) calloc(numvar, sizeof(double));
              if (xx)
              {
                MSK_getxx(task,
                          MSK_SOL_ITG,
                          xx);
               
                printf("Optimal primal solution\n");
                for (j = 0; j < numvar; ++j)
                  printf("x[%d]: %e\n", j, xx[j]);

                free(xx);
              }
              else
                r = MSK_RES_ERR_SPACE;

              break;
            }
          default:
            printf("Another solution status.\n");
            break;
        }
      }
    }

    if (r != MSK_RES_OK)
    {
      /* In case of an error print error code and description. */
      char symname[MSK_MAX_STR_LEN];
      char desc[MSK_MAX_STR_LEN];

      printf("An error occurred while optimizing.\n");
      MSK_getcodedesc(r,
                      symname,
                      desc);
      printf("Error %s - '%s'\n", symname, desc);
    }

    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return r;
}

errorreporting.c

Listing 17.9 errorreporting.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File: errorreporting.c

   Purpose:   To demonstrate how the error reporting can be customized.
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "mosek.h"

static MSKrescodee MSKAPI handleresponse(MSKuserhandle_t handle,
                                         MSKrescodee     r,
                                         const char      msg[])
/* A custom response handler. */
{
  /* Find out if we had an error or warning */
  MSKrescodetypee rct;
  MSK_getresponseclass(r, &rct);

  if (r == MSK_RES_OK)
  {
    /* Do nothing */
  }
  else if (rct == MSK_RESPONSE_WRN)
  {
    printf("MOSEK reports warning number %d: %s\n", r, msg);
    r = MSK_RES_OK;
  }
  else if (rct == MSK_RESPONSE_ERR)
  {
    printf("MOSEK reports error number %d: %s\n", r, msg);
  }
  else
  {
    printf("Unexpected response %d: %s\n", r, msg);
  }

  return (r);
} /* handlerespone */


int main(int argc, char *argv[])
{
  MSKenv_t    env;
  MSKrescodee r;
  MSKtask_t   task;

  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    r = MSK_makeemptytask(env, &task);
    if (r == MSK_RES_OK)
    {
      /*
       * Input a custom warning and error handler function.
       */

      MSK_putresponsefunc(task, handleresponse, NULL);

      /* User defined code goes here */
      /* This will provoke an error */

      if (r == MSK_RES_OK)
        r = MSK_putaij(task, 10, 10, 1.0);

    }
    MSK_deletetask(&task);
  }
  MSK_deleteenv(&env);

  printf("Return code - %d\n", r);

  if (r == MSK_RES_ERR_INDEX_IS_TOO_LARGE)
    return MSK_RES_OK;
  else
    return 1;
} /* main */

feasrepairex1.c

Listing 17.10 feasrepairex1.c Click here to download.

/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:      feasrepairex1.c

  Purpose:   To demonstrate how to use the MSK_primalrepair function to
             repair an infeasible problem.
*/

#include <math.h>
#include <stdio.h>

#include "mosek.h"

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  fputs(str, stdout);
} /* printstr */

static const char * feasrepair_lp =
"minimize                                    \n"
" obj: - 10 x1 - 9 x2                        \n"
"st                                          \n"
" c1: + 7e-01 x1 + x2 <= 630                 \n"
" c2: + 5e-01 x1 + 8.333333333e-01 x2 <= 600 \n"
" c3: + x1 + 6.6666667e-01 x2 <= 708         \n"
" c4: + 1e-01 x1 + 2.5e-01 x2 <= 135         \n"
"bounds                                      \n"
"x2 >= 650                                   \n"
"end                                         \n";

int main(int argc, const char *argv[])
{
  MSKenv_t    env;
  MSKrescodee r, trmcode;
  MSKtask_t   task;

  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
    r = MSK_makeemptytask(env, &task);

  if (r == MSK_RES_OK)
    MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

  if (r == MSK_RES_OK)
    r = MSK_readlpstring(task,feasrepair_lp); /* Read problem from string */

  if (r == MSK_RES_OK)
    r = MSK_putintparam(task, MSK_IPAR_LOG_FEAS_REPAIR, 3);

  if (r == MSK_RES_OK)
  {
    /* Weights are NULL implying all weights are 1. */
    r = MSK_primalrepair(task, NULL, NULL, NULL, NULL);
  }

  if (r == MSK_RES_OK)
  {
    double sum_viol;

    r = MSK_getdouinf(task, MSK_DINF_PRIMAL_REPAIR_PENALTY_OBJ, &sum_viol);

    if (r == MSK_RES_OK)
    {
      printf("Minimized sum of violations = %e\n", sum_viol);

      r = MSK_optimizetrm(task, &trmcode); /* Optimize the repaired task. */

      MSK_solutionsummary(task, MSK_STREAM_MSG);
    }
  }

  printf("Return code: %d\n", r);

  return (r);
}

gp1.c

Listing 17.11 gp1.c Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      gp1.c
//
//   Purpose:   Demonstrates how to solve a simple Geometric Program (GP)
//              cast into conic form with exponential cones and log-sum-exp.
//
//              Example from
//                https://gpkit.readthedocs.io/en/latest/examples.html//maximizing-the-volume-of-a-box
//
#include <stdio.h>
#include <math.h>
#include "mosek.h"

/* This function prints log output from MOSEK to the terminal. */
static void MSKAPI printstr(void       *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

// maximize     h*w*d
// subjecto to  2*(h*w + h*d) <= Awall
//              w*d <= Afloor
//              alpha <= h/w <= beta
//              gamma <= d/w <= delta
//
// Variable substitutions:  h = exp(x), w = exp(y), d = exp(z).
//
// maximize     x+y+z
// subject      log( exp(x+y+log(2/Awall)) + exp(x+z+log(2/Awall)) ) <= 0
//                              y+z <= log(Afloor)
//              log( alpha ) <= x-y <= log( beta )
//              log( gamma ) <= z-y <= log( delta )
int max_volume_box(double Aw, double Af, 
                   double alpha, double beta, double gamma, double delta,
                   double hwd[])
{
  // Basic dimensions of our problem
  const int numvar    = 3;  // Variables in original problem
  const int numcon    = 3;  // Linear constraints in original problem

  // Linear part of the problem involving x, y, z
  const double       cval[]  = {1, 1, 1};
  const int          asubi[] = {0, 0, 1, 1, 2, 2};
  const int          asubj[] = {1, 2, 0, 1, 2, 1};
  const int          alen    = 6;
  const double       aval[]  = {1.0, 1.0, 1.0, -1.0, 1.0, -1.0};
  const MSKboundkeye bkc[]   = {MSK_BK_UP, MSK_BK_RA, MSK_BK_RA};
  const double       blc[]   = {-MSK_INFINITY, log(alpha), log(gamma)};
  const double       buc[]   = {log(Af), log(beta), log(delta)};

  // Affine conic constraint data of the problem
  MSKint64t       expdomidx, rzerodomidx;
  const MSKint64t numafe = 6, f_nnz = 8;
  const MSKint64t afeidx[] = {0, 1, 2, 2, 3, 3, 5, 5};
  const MSKint32t varidx[] = {3, 4, 0, 1, 0, 2, 3, 4};
  const double     f_val[] = {1, 1, 1, 1, 1, 1, 1, 1};
  const double         g[] = {0, 0, log(2/Aw), log(2/Aw), 1, -1};

  MSKtask_t          task = NULL;
  MSKrescodee        r = MSK_RES_OK, trmcode;
  MSKsolstae         solsta;
  MSKint32t          i;
  double             *xyz = (double*) calloc(numvar, sizeof(double));

  if (r == MSK_RES_OK)
    r = MSK_maketask(NULL, 0, 0, &task);

  if (r == MSK_RES_OK)
    r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

  if (r == MSK_RES_OK)
    r = MSK_appendvars(task, numvar);

  if (r == MSK_RES_OK)
    r = MSK_appendcons(task, numcon);

  if (r == MSK_RES_OK)
    r = MSK_appendafes(task, numafe);

  // Objective is the sum of three first variables
  if (r == MSK_RES_OK)
    r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE);
  if (r == MSK_RES_OK)
    r = MSK_putcslice(task, 0, numvar, cval);
  if (r == MSK_RES_OK)
    r = MSK_putvarboundsliceconst(task, 0, numvar, MSK_BK_FR, -MSK_INFINITY, +MSK_INFINITY);

  // Add the three linear constraints
  if (r == MSK_RES_OK)
    r = MSK_putaijlist(task, alen, asubi, asubj, aval);
  if (r == MSK_RES_OK)  
    r = MSK_putconboundslice(task, 0, numvar, bkc, blc, buc);

  if (r == MSK_RES_OK)
  {
    MSKint64t acc1_afeidx[] = {0, 4, 2};
    MSKint64t acc2_afeidx[] = {1, 4, 3};
    MSKint64t acc3_afeidx[] = {5};

    // Affine expressions appearing in affine conic constraints
    // in this order:
    // u1, u2, x+y+log(2/Awall), x+z+log(2/Awall), 1.0, u1+u2-1.0
    if (r == MSK_RES_OK) 
      r = MSK_appendvars(task, 2);
    if (r == MSK_RES_OK)
      r = MSK_putvarboundsliceconst(task, numvar, numvar+2, MSK_BK_FR, -MSK_INFINITY, +MSK_INFINITY);

    if (r == MSK_RES_OK) 
      r = MSK_putafefentrylist(task, f_nnz, afeidx, varidx, f_val);
    if (r == MSK_RES_OK) 
      r = MSK_putafegslice(task, 0, numafe, g);
    
    /* Append the primal exponential cone domain */
    if (r == MSK_RES_OK)
      r = MSK_appendprimalexpconedomain(task, &expdomidx);
    
    /* (u1, 1, x+y+log(2/Awall)) \in EXP */
    if (r == MSK_RES_OK)
      r = MSK_appendacc(task, expdomidx, 3, acc1_afeidx, NULL);

    /* (u2, 1, x+z+log(2/Awall)) \in EXP */
    if (r == MSK_RES_OK)
      r = MSK_appendacc(task, expdomidx, 3, acc2_afeidx, NULL);

    /* The constraint u1+u2-1 \in \ZERO is added also as an ACC */
    if (r == MSK_RES_OK)
      r = MSK_appendrzerodomain(task, 1, &rzerodomidx);
    if (r == MSK_RES_OK)
      r = MSK_appendacc(task, rzerodomidx, 1, acc3_afeidx, NULL);
  }

  // Solve and map to original h, w, d
  if (r == MSK_RES_OK)
    r = MSK_optimizetrm(task, &trmcode);

  if (r == MSK_RES_OK)
    MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

  if (solsta == MSK_SOL_STA_OPTIMAL)
  {
    if (r == MSK_RES_OK)
      r = MSK_getxxslice(task, MSK_SOL_ITR, 0, numvar, xyz);
    for(i = 0; i < numvar; i++) hwd[i] = exp(xyz[i]);
  }
  else
  {
    printf("Solution not optimal, termination code %d.\n", trmcode);
    r = trmcode;
  }

  free(xyz);
  return r;
}
    
int main()
{
  const double Aw    = 200.0;
  const double Af    = 50.0;
  const double alpha = 2.0;
  const double beta  = 10.0;
  const double gamma = 2.0;
  const double delta = 10.0;
  MSKrescodee  r;
  double       hwd[3];

  r = max_volume_box(Aw, Af, alpha, beta, gamma, delta, hwd);

  printf("Response code: %d\n", r);
  if (r == MSK_RES_OK)
    printf("Solution h=%.4f w=%.4f d=%.4f\n", hwd[0], hwd[1], hwd[2]);

  return r;
}

helloworld.c

Listing 17.12 helloworld.c Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      helloworld.c
//
//  The most basic example of how to get started with MOSEK.

#include "mosek.h"
#include <stdio.h>

/* Error checking not included */
int main() {
  MSKrescodee      r, trmcode;
  MSKenv_t         env  = NULL;
  MSKtask_t        task = NULL;
  double           xx = 0.0;

  MSK_maketask(NULL, 0, 1, &task);      // Create task

  MSK_appendvars(task, 1);                             // 1 variable x
  MSK_putcj(task, 0, 1.0);                             // c_0 = 1.0
  MSK_putvarbound(task, 0, MSK_BK_RA, 2.0, 3.0);       // 2.0 <= x <= 3.0
  MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MINIMIZE); // Minimize

  MSK_optimizetrm(task, &trmcode);         // Optimize

  MSK_getxx(task, MSK_SOL_ITR, &xx);       // Get solution
  printf("Solution x = %f\n", xx);         // Print solution

  MSK_deletetask(&task); // Clean up task
  return 0;
}

lo1.c

Listing 17.13 lo1.c Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:      lo1.c

  Purpose:   To demonstrate how to solve a small linear
             optimization problem using the MOSEK C API,
             and handle the solver result and the problem
             solution.
*/
#include <stdio.h>
#include "mosek.h"

/* This function prints log output from MOSEK to the terminal. */
static void MSKAPI printstr(void       *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  const MSKint32t    numvar = 4,
                     numcon = 3;

  const double       c[]     = {3.0, 1.0, 5.0, 1.0};
  /* Below is the sparse representation of the A
     matrix stored by column. */
  const MSKint32t    aptrb[] = {0, 2, 5, 7},
                     aptre[] = {2, 5, 7, 9},
                     asub[]  = { 0, 1,
                                 0, 1, 2,
                                 0, 1,
                                 1, 2
                               };
  const double       aval[]  = { 3.0, 2.0,
                                 1.0, 1.0, 2.0,
                                 2.0, 3.0,
                                 1.0, 3.0
                               };

  /* Bounds on constraints. */
  const MSKboundkeye bkc[]  = {MSK_BK_FX, MSK_BK_LO,     MSK_BK_UP    };
  const double       blc[]  = {30.0,      15.0,          -MSK_INFINITY};
  const double       buc[]  = {30.0,      +MSK_INFINITY, 25.0         };
  /* Bounds on variables. */
  const MSKboundkeye bkx[]  = {MSK_BK_LO,     MSK_BK_RA, MSK_BK_LO,     MSK_BK_LO     };
  const double       blx[]  = {0.0,           0.0,       0.0,           0.0           };
  const double       bux[]  = { +MSK_INFINITY, 10.0,      +MSK_INFINITY, +MSK_INFINITY };
  MSKtask_t          task = NULL;
  MSKrescodee        r = MSK_RES_OK;
  MSKint32t          i, j;

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(NULL, numcon, numvar, &task);

    /* Directs the log task stream to the 'printstr' function. */
    if (r == MSK_RES_OK)
      r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

    /* Append 'numcon' empty constraints.
     The constraints will initially have no bounds. */
    if (r == MSK_RES_OK)
      r = MSK_appendcons(task, numcon);

    /* Append 'numvar' variables.
     The variables will initially be fixed at zero (x=0). */
    if (r == MSK_RES_OK)
      r = MSK_appendvars(task, numvar);

    for (j = 0; j < numvar && r == MSK_RES_OK; ++j)
    {
      /* Set the linear term c_j in the objective.*/
      if (r == MSK_RES_OK)
        r = MSK_putcj(task, j, c[j]);


      /* Set the bounds on variable j.
       blx[j] <= x_j <= bux[j] */
      if (r == MSK_RES_OK)
        r = MSK_putvarbound(task,
                            j,           /* Index of variable.*/
                            bkx[j],      /* Bound key.*/
                            blx[j],      /* Numerical value of lower bound.*/
                            bux[j]);     /* Numerical value of upper bound.*/

      /* Input column j of A */
      if (r == MSK_RES_OK)
        r = MSK_putacol(task,
                        j,                 /* Variable (column) index.*/
                        aptre[j] - aptrb[j], /* Number of non-zeros in column j.*/
                        asub + aptrb[j],   /* Pointer to row indexes of column j.*/
                        aval + aptrb[j]);  /* Pointer to Values of column j.*/
    }

    /* Set the bounds on constraints.
       for i=1, ...,numcon : blc[i] <= constraint i <= buc[i] */
    for (i = 0; i < numcon && r == MSK_RES_OK; ++i)
      r = MSK_putconbound(task,
                          i,           /* Index of constraint.*/
                          bkc[i],      /* Bound key.*/
                          blc[i],      /* Numerical value of lower bound.*/
                          buc[i]);     /* Numerical value of upper bound.*/

    /* Maximize objective function. */
    if (r == MSK_RES_OK)
      r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE);

    if (r == MSK_RES_OK)
    {
      MSKrescodee trmcode;

      /* Run optimizer */
      r = MSK_optimizetrm(task, &trmcode);

      /* Print a summary containing information
         about the solution for debugging purposes. */
      MSK_solutionsummary(task, MSK_STREAM_LOG);

      if (r == MSK_RES_OK)
      {
        MSKsolstae solsta;

        if (r == MSK_RES_OK)
          r = MSK_getsolsta(task,
                            MSK_SOL_BAS,
                            &solsta);
        switch (solsta)
        {
          case MSK_SOL_STA_OPTIMAL:
            {
              double *xx = (double*) calloc(numvar, sizeof(double));
              if (xx)
              {
                MSK_getxx(task,
                          MSK_SOL_BAS,    /* Request the basic solution. */
                          xx);

                printf("Optimal primal solution\n");
                for (j = 0; j < numvar; ++j)
                  printf("x[%d]: %e\n", j, xx[j]);

                free(xx);
              }
              else
                r = MSK_RES_ERR_SPACE;

              break;
            }
          case MSK_SOL_STA_DUAL_INFEAS_CER:
          case MSK_SOL_STA_PRIM_INFEAS_CER:
            printf("Primal or dual infeasibility certificate found.\n");
            break;
          case MSK_SOL_STA_UNKNOWN:
            {
              char symname[MSK_MAX_STR_LEN];
              char desc[MSK_MAX_STR_LEN];

              /* If the solutions status is unknown, print the termination code
                 indicating why the optimizer terminated prematurely. */

              MSK_getcodedesc(trmcode,
                              symname,
                              desc);

              printf("The solution status is unknown.\n");
              printf("The optimizer terminitated with code: %s\n", symname);
              break;
            }
          default:
            printf("Other solution status.\n");
            break;
        }
      }
    }

    if (r != MSK_RES_OK)
    {
      /* In case of an error print error code and description. */
      char symname[MSK_MAX_STR_LEN];
      char desc[MSK_MAX_STR_LEN];

      printf("An error occurred while optimizing.\n");
      MSK_getcodedesc(r,
                      symname,
                      desc);
      printf("Error %s - '%s'\n", symname, desc);
    }

    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  return r;
}

lo2.c

Listing 17.14 lo2.c Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:      lo2.c

  Purpose:   To demonstrate how to solve a small linear
             optimization problem using the MOSEK C API,
             and handle the solver result and the problem
             solution.
*/
#include <stdio.h>
#include "mosek.h"

/* This function prints log output from MOSEK to the terminal. */
static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  const int numvar = 4,
            numcon = 3;

  double       c[]     = {3.0, 1.0, 5.0, 1.0};
  /* Below is the sparse representation of the A
     matrix stored by row. */
  MSKint32t    aptrb[] = {0, 3, 7};
  MSKint32t    aptre[] = {3, 7, 9};
  MSKint32t    asub[]  = { 0, 1, 2,
                           0, 1, 2, 3,
                           1, 3
                         };
  double       aval[]  = { 3.0, 1.0, 2.0,
                           2.0, 1.0, 3.0, 1.0,
                           2.0, 3.0
                         };

  /* Bounds on constraints. */
  MSKboundkeye bkc[]  = {MSK_BK_FX, MSK_BK_LO,     MSK_BK_UP    };
  double       blc[]  = {30.0,      15.0,          -MSK_INFINITY};
  double       buc[]  = {30.0,      +MSK_INFINITY, 25.0         };
  /* Bounds on variables. */
  MSKboundkeye bkx[]  = {MSK_BK_LO,     MSK_BK_RA, MSK_BK_LO,     MSK_BK_LO     };
  double       blx[]  = {0.0,           0.0,       0.0,           0.0           };
  double       bux[]  = { +MSK_INFINITY, 10.0,      +MSK_INFINITY, +MSK_INFINITY };
  MSKenv_t     env  = NULL;
  MSKtask_t    task = NULL;
  MSKrescodee  r;
  MSKint32t    i, j;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, numcon, numvar, &task);

    /* Directs the log task stream to the 'printstr' function. */
    if (r == MSK_RES_OK)
      r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

    /* Append 'numcon' empty constraints.
     The constraints will initially have no bounds. */
    if (r == MSK_RES_OK)
      r = MSK_appendcons(task, numcon);

    /* Append 'numvar' variables.
     The variables will initially be fixed at zero (x=0). */
    if (r == MSK_RES_OK)
      r = MSK_appendvars(task, numvar);

    for (j = 0; j < numvar && r == MSK_RES_OK; ++j)
    {
      /* Set the linear term c_j in the objective.*/
      if (r == MSK_RES_OK)
        r = MSK_putcj(task, j, c[j]);

      /* Set the bounds on variable j.
       blx[j] <= x_j <= bux[j] */
      if (r == MSK_RES_OK)
        r = MSK_putvarbound(task,
                            j,           /* Index of variable.*/
                            bkx[j],      /* Bound key.*/
                            blx[j],      /* Numerical value of lower bound.*/
                            bux[j]);     /* Numerical value of upper bound.*/
    }

    /* Set the bounds on constraints.
       for i=1, ...,numcon : blc[i] <= constraint i <= buc[i] */
    for (i = 0; i < numcon && r == MSK_RES_OK; ++i)
    {
      r = MSK_putconbound(task,
                          i,           /* Index of constraint.*/
                          bkc[i],      /* Bound key.*/
                          blc[i],      /* Numerical value of lower bound.*/
                          buc[i]);     /* Numerical value of upper bound.*/

      /* Input row i of A */
      if (r == MSK_RES_OK)
        r = MSK_putarow(task,
                        i,                 /* Row index.*/
                        aptre[i] - aptrb[i], /* Number of non-zeros in row i.*/
                        asub + aptrb[i],   /* Pointer to column indexes of row i.*/
                        aval + aptrb[i]);  /* Pointer to values of row i.*/
    }

    /* Maximize objective function. */
    if (r == MSK_RES_OK)
      r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE);

    if (r == MSK_RES_OK)
    {
      MSKrescodee trmcode;

      /* Run optimizer */
      r = MSK_optimizetrm(task, &trmcode);

      /* Print a summary containing information
         about the solution for debugging purposes. */
      MSK_solutionsummary(task, MSK_STREAM_LOG);

      if (r == MSK_RES_OK)
      {
        MSKsolstae solsta;

        if (r == MSK_RES_OK)
          r = MSK_getsolsta(task, MSK_SOL_BAS, &solsta);
        switch (solsta)
        {
          case MSK_SOL_STA_OPTIMAL:
            {
              double *xx = (double*) calloc(numvar, sizeof(double));
              if (xx)
              {
                MSK_getxx(task,
                          MSK_SOL_BAS,    /* Request the basic solution. */
                          xx);

                printf("Optimal primal solution\n");
                for (j = 0; j < numvar; ++j)
                  printf("x[%d]: %e\n", j, xx[j]);
              }
              else
              {
                r = MSK_RES_ERR_SPACE;
              }

              free(xx);
              break;
            }
          case MSK_SOL_STA_DUAL_INFEAS_CER:
          case MSK_SOL_STA_PRIM_INFEAS_CER:
            printf("Primal or dual infeasibility certificate found.\n");
            break;
          case MSK_SOL_STA_UNKNOWN:
            {
              char symname[MSK_MAX_STR_LEN];
              char desc[MSK_MAX_STR_LEN];

              /* If the solutions status is unknown, print the termination code
                 indicating why the optimizer terminated prematurely. */

              MSK_getcodedesc(trmcode,
                              symname,
                              desc);

              printf("The solutuion status is unknown.\n");
              printf("The optimizer terminitated with code: %s\n", symname);
              break;
            }
          default:
            printf("Other solution status.\n");
            break;
        }
      }
    }

    if (r != MSK_RES_OK)
    {
      /* In case of an error print error code and description. */
      char symname[MSK_MAX_STR_LEN];
      char desc[MSK_MAX_STR_LEN];

      printf("An error occurred while optimizing.\n");
      MSK_getcodedesc(r,
                      symname,
                      desc);
      printf("Error %s - '%s'\n", symname, desc);
    }

    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return r;
}

logistic.c

Listing 17.15 logistic.c Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      logistic.c
//
// Purpose: Implements logistic regression with regulatization.
//
//          Demonstrates using the exponential cone and log-sum-exp in Optimizer API.

#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */

#define MSKCALL(x) if (res==MSK_RES_OK) res = (x);

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

const double inf = 0.0;

// Adds ACCs for t_i >= log ( 1 + exp((1-2*y[i]) * theta' * X[i]) )
// Adds auxiliary variables, AFE rows and constraints
MSKrescodee softplus(MSKtask_t task, int d, int n, MSKint32t theta, MSKint32t t, double* X, int* y)
{
  MSKint32t nvar, ncon;
  MSKint64t nafe, thetaafe, tafe, z1afe, z2afe, oneafe, expdomain;
  MSKint32t z1, z2, zcon, v1con, v2con;
  MSKint32t  *subi = (MSKint32t*) calloc(2*n, sizeof(MSKint32t));
  MSKint32t  *subj = (MSKint32t*) calloc(3*n, sizeof(MSKint32t));
  MSKrealt   *aval = (MSKrealt*) calloc(2*n, sizeof(MSKrealt));
  MSKint64t  *afeidx = (MSKint64t*) calloc(d*n+4*n, sizeof(MSKint64t));
  MSKint32t  *varidx = (MSKint32t*) calloc(d*n+4*n, sizeof(MSKint32t));
  MSKrealt   *fval   = (MSKrealt*) calloc(d*n+4*n, sizeof(MSKrealt));
  MSKint64t  idx[3];
  int        k, i, j;
  MSKrescodee res = MSK_RES_OK;

  MSKCALL(MSK_getnumvar(task, &nvar));
  MSKCALL(MSK_getnumcon(task, &ncon));
  MSKCALL(MSK_getnumafe(task, &nafe));
  MSKCALL(MSK_appendvars(task, 2*n));   // z1, z2
  MSKCALL(MSK_appendcons(task, n));     // z1 + z2 = 1
  MSKCALL(MSK_appendafes(task, 4*n));   //theta * X[i] - t[i], -t[i], z1[i], z2[i]

  z1 = nvar, z2 = nvar+n;
  zcon = ncon;
  thetaafe = nafe, tafe = nafe+n, z1afe = nafe+2*n, z2afe = nafe+3*n;
  
  // Linear constraints
  k = 0;
  for(i = 0; i < n; i++)
  {
    // z1 + z2 = 1
    subi[k] = zcon+i;  subj[k] = z1+i;  aval[k] = 1;  k++;
    subi[k] = zcon+i;  subj[k] = z2+i;  aval[k] = 1;  k++;
  }
  MSKCALL(MSK_putaijlist(task, 2*n, subi, subj, aval));
  MSKCALL(MSK_putconboundsliceconst(task, zcon, zcon+n, MSK_BK_FX, 1, 1));
  MSKCALL(MSK_putvarboundsliceconst(task, nvar, nvar+2*n, MSK_BK_FR, -inf, inf));

  // Affine conic expressions
  k = 0;

  // Thetas
  for(i = 0; i < n; i++) {
    for(j = 0; j < d; j++) {
      afeidx[k] = thetaafe + i; varidx[k] = theta + j; 
      fval[k] = ((y[i]) ? -1 : 1) * X[i*d+j];
      k++;
    }
  }

  // -t[i]
  for(i = 0; i < n; i++) {
    afeidx[k] = thetaafe + i; varidx[k] = t + i; fval[k] = -1; k++;
    afeidx[k] = tafe + i;     varidx[k] = t + i; fval[k] = -1; k++;
  }

  // z1, z2
  for(i = 0; i < n; i++) {
    afeidx[k] = z1afe + i; varidx[k] = z1 + i; fval[k] = 1; k++;
    afeidx[k] = z2afe + i; varidx[k] = z2 + i; fval[k] = 1; k++;
  }

  // Add the expressions
  MSKCALL(MSK_putafefentrylist(task, d*n+4*n, afeidx, varidx, fval));

  // Add a single row with the constant expression "1.0"
  MSKCALL(MSK_getnumafe(task, &oneafe));
  MSKCALL(MSK_appendafes(task,1));
  MSKCALL(MSK_putafeg(task, oneafe, 1.0));

  // Add an exponential cone domain
  MSKCALL(MSK_appendprimalexpconedomain(task, &expdomain));

  // Conic constraints
  for(i = 0; i < n; i++)
  {
    idx[0] = z1afe+i, idx[1] = oneafe, idx[2] = thetaafe+i;
    MSKCALL(MSK_appendacc(task, expdomain, 3, idx, NULL));
    idx[0] = z2afe+i, idx[1] = oneafe, idx[2] = tafe+i;
    MSKCALL(MSK_appendacc(task, expdomain, 3, idx, NULL));
  }
    
  free(subi); free(subj); free(aval); 
  free(afeidx); free(varidx); free(fval); 
  return res;
}

// Model logistic regression (regularized with full 2-norm of theta)
// X - n x d matrix of data points
// y - length n vector classifying training points
// lamb - regularization parameter
MSKrescodee logisticRegression(MSKenv_t       env,
                               int            n,    // num samples
                               int            d,    // dimension
                               double        *X, 
                               int           *y,
                               double         lamb,
                               double        *thetaVal)   // result
{
  MSKrescodee res = MSK_RES_OK;
  MSKrescodee trm = MSK_RES_OK;
  MSKtask_t task = NULL;
  MSKint32t nvar = 1+d+n;
  MSKint32t r = 0, theta = 1, t = 1+d;
  MSKint64t numafe, quadDom;
  int i = 0;

  MSKCALL(MSK_maketask(env, 0, 0, &task));
  MSKCALL(MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr));

  // Variables [r; theta; t]
  MSKCALL(MSK_appendvars(task, nvar));
  MSKCALL(MSK_putvarboundsliceconst(task, 0, nvar, MSK_BK_FR, -inf, inf));

  // Objective lambda*r + sum(t)
  MSKCALL(MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MINIMIZE));
  MSKCALL(MSK_putcj(task, r, lamb));
  for(i = 0; i < n && res == MSK_RES_OK; i++) 
    MSKCALL(MSK_putcj(task, t+i, 1.0));

  // Softplus function constraints
  MSKCALL(softplus(task, d, n, theta, t, X, y));

  // Regularization
  // Append a sequence of linear expressions (r, theta) to F
  MSKCALL(MSK_getnumafe(task, &numafe));
  MSKCALL(MSK_appendafes(task,1+d));
  MSKCALL(MSK_putafefentry(task, numafe, r, 1.0));
  for(i = 0; i < d; i++)
    MSKCALL(MSK_putafefentry(task, numafe + i + 1, theta + i, 1.0));

  // Add the constraint
  MSKCALL(MSK_appendquadraticconedomain(task, 1+d, &quadDom));
  MSKCALL(MSK_appendaccseq(task, quadDom, 1+d, numafe, NULL));

  // Solution
  MSKCALL(MSK_optimizetrm(task, &trm));
  MSKCALL(MSK_solutionsummary(task, MSK_STREAM_MSG));

  MSKCALL(MSK_getxxslice(task, MSK_SOL_ITR, theta, theta+d, thetaVal));
  
  return res;
}

int main()
{
  MSKenv_t env;
  MSKrescodee res = MSK_RES_OK;

  MSKCALL(MSK_makeenv(&env, NULL));
  
  // Test: detect and approximate a circle using degree 2 polynomials
  {
    int     n = 30;
    double  X[6*30*30];
    int     Y[30*30];
    int     i,j;
    double  theta[6];
  
    for(i=0; i<n; i++) 
    for(j=0; j<n; j++)
    {
      int k = i*n+j;
      double x = -1 + 2.0*i/(n-1);
      double y = -1 + 2.0*j/(n-1);
      X[6*k+0] = 1.0; X[6*k+1] = x; X[6*k+2] = y; X[6*k+3] = x*y;
      X[6*k+4] = x*x; X[6*k+5] = y*y;
      Y[k] = (x*x+y*y>=0.69) ? 1 : 0;
    }

    MSKCALL(logisticRegression(env, n*n, 6, X, Y, 0.1, theta));

    if (res == MSK_RES_OK)
      for(i=0;i<6;i++) printf("%.4f\n", theta[i]);
  }

  return res;
}

mico1.c

Listing 17.16 mico1.c Click here to download.
/*
   Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File :      mico1.c

   Purpose :   Demonstrates how to solve a small mixed
               integer conic optimization problem.

               minimize    x^2 + y^2
               subject to  x >= e^y + 3.8
                           x, y - integer
*/

#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, char *argv[])
{
  MSKint32t         numvar = 3;  /* x, y, t */

  MSKvariabletypee  vart[] = { MSK_VAR_TYPE_INT, MSK_VAR_TYPE_INT };
  MSKint32t       intsub[] = { 0, 1 }; 

  MSKint32t    i, j;

  MSKenv_t     env = NULL;
  MSKtask_t    task = NULL;
  MSKrescodee  r, trm;

  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    r = MSK_maketask(env, 0, 0, &task);

    if (r == MSK_RES_OK)
      r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

    if (r == MSK_RES_OK)
      r = MSK_appendvars(task, numvar);
    if (r == MSK_RES_OK)
      r = MSK_putvarboundsliceconst(task, 0, numvar, MSK_BK_FR, -0.0, 0.0);    

    /* Integrality constraints */
    if (r == MSK_RES_OK)
      r = MSK_putvartypelist(task, 2, intsub, vart);        

    /* Objective */
    if (r == MSK_RES_OK)
      r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MINIMIZE);
    if (r == MSK_RES_OK)
      r = MSK_putcj(task, 2, 1.0);    /* Minimize t */

    /* Conic part of the problem */
    if (r == MSK_RES_OK) 
    {
      /* Set up the affine expressions */
      /* x, x-3.8, y, t, 1.0 */
      MSKint64t afeidx[] = {0, 1, 2, 3};
      MSKint32t varidx[] = {0, 0, 1, 2};
      MSKrealt     val[] = {1.0, 1.0, 1.0, 1.0};
      MSKrealt       g[] = {0.0, -3.8, 0.0, 0.0, 1.0};
      MSKint64t domExp, domQuad;
      MSKint64t afeidxExp[]  = {1, 4, 2};
      MSKint64t afeidxQuad[] = {3, 0, 2};

      if (r == MSK_RES_OK)
        r = MSK_appendafes(task, 5);
      if (r == MSK_RES_OK)
        r = MSK_putafefentrylist(task,
                                 4,
                                 afeidx,
                                 varidx,
                                 val);
      if (r == MSK_RES_OK)
        r = MSK_putafegslice(task, 0, 5, g);

      // Add constraint (x-3.8, 1, y) \in \EXP
      if (r == MSK_RES_OK)
        r = MSK_appendprimalexpconedomain(task, &domExp);
      if (r == MSK_RES_OK)
        r = MSK_appendacc(task, domExp, 3, afeidxExp, NULL);

      // Add constraint (t, x, y) \in \QUAD
      if (r == MSK_RES_OK)
        r = MSK_appendquadraticconedomain(task, 3, &domQuad);
      if (r == MSK_RES_OK)
        r = MSK_appendacc(task, domQuad, 3, afeidxQuad, NULL);
    }

    /* Optimize the problem */
    if (r == MSK_RES_OK)
      r = MSK_optimizetrm(task, &trm);

    if (r == MSK_RES_OK)
      r = MSK_solutionsummary(task, MSK_STREAM_MSG);

    if (r == MSK_RES_OK)
    {
      MSKrealt xx[] = {0, 0};

      r = MSK_getxxslice(task, MSK_SOL_ITG, 0, 2, xx);
      
      if (r == MSK_RES_OK)
        printf("x = %.2f, y = %.2f\n", xx[0], xx[1]);
    }

    if (task) MSK_deletetask(&task);
  }

  if (env) MSK_deleteenv(&env);
  return r;
} 

milo1.c

Listing 17.17 milo1.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      milo1.c

   Purpose:   To demonstrate how to solve a small mixed
              integer linear optimization problem using
              the MOSEK API.
*/
#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, char *argv[])
{
  const MSKint32t numvar = 2,
                  numcon = 2;

  double       c[]   = {  1.0, 0.64 };
  MSKboundkeye bkc[] = { MSK_BK_UP,    MSK_BK_LO };
  double       blc[] = { -MSK_INFINITY, -4.0 };
  double       buc[] = { 250.0,        MSK_INFINITY };

  MSKboundkeye bkx[] = { MSK_BK_LO,    MSK_BK_LO };
  double       blx[] = { 0.0,          0.0 };
  double       bux[] = { MSK_INFINITY, MSK_INFINITY };


  MSKint32t    aptrb[] = { 0, 2 },
               aptre[] = { 2, 4 },
               asub[] = { 0,    1,   0,    1 };
  double       aval[] = { 50.0, 3.0, 31.0, -2.0 };
  MSKint32t    i, j;

  MSKenv_t     env = NULL;
  MSKtask_t    task = NULL;
  MSKrescodee  r;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  /* Check if return code is ok. */
  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, 0, 0, &task);

    if (r == MSK_RES_OK)
      r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

    /* Append 'numcon' empty constraints.
     The constraints will initially have no bounds. */
    if (r == MSK_RES_OK)
      r = MSK_appendcons(task, numcon);

    /* Append 'numvar' variables.
     The variables will initially be fixed at zero (x=0). */
    if (r == MSK_RES_OK)
      r = MSK_appendvars(task, numvar);

    /* Optionally add a constant term to the objective. */
    if (r == MSK_RES_OK)
      r = MSK_putcfix(task, 0.0);
    for (j = 0; j < numvar && r == MSK_RES_OK; ++j)
    {
      /* Set the linear term c_j in the objective.*/
      if (r == MSK_RES_OK)
        r = MSK_putcj(task, j, c[j]);

      /* Set the bounds on variable j.
       blx[j] <= x_j <= bux[j] */
      if (r == MSK_RES_OK)
        r = MSK_putvarbound(task,
                            j,           /* Index of variable.*/
                            bkx[j],      /* Bound key.*/
                            blx[j],      /* Numerical value of lower bound.*/
                            bux[j]);     /* Numerical value of upper bound.*/

      /* Input column j of A */
      if (r == MSK_RES_OK)
        r = MSK_putacol(task,
                        j,                 /* Variable (column) index.*/
                        aptre[j] - aptrb[j], /* Number of non-zeros in column j.*/
                        asub + aptrb[j],   /* Pointer to row indexes of column j.*/
                        aval + aptrb[j]);  /* Pointer to Values of column j.*/

    }

    /* Set the bounds on constraints.
       for i=1, ...,numcon : blc[i] <= constraint i <= buc[i] */
    for (i = 0; i < numcon && r == MSK_RES_OK; ++i)
      r = MSK_putconbound(task,
                          i,           /* Index of constraint.*/
                          bkc[i],      /* Bound key.*/
                          blc[i],      /* Numerical value of lower bound.*/
                          buc[i]);     /* Numerical value of upper bound.*/

    /* Specify integer variables. */
    for (j = 0; j < numvar && r == MSK_RES_OK; ++j)
      r = MSK_putvartype(task, j, MSK_VAR_TYPE_INT);

    if (r == MSK_RES_OK)
      r =  MSK_putobjsense(task,
                           MSK_OBJECTIVE_SENSE_MAXIMIZE);

    if (r == MSK_RES_OK)
      /* Set max solution time */
      r = MSK_putdouparam(task,
                          MSK_DPAR_MIO_MAX_TIME,
                          60.0);

    if (r == MSK_RES_OK)
    {
      MSKrescodee trmcode;

      /* Run optimizer */
      r = MSK_optimizetrm(task, &trmcode);

      /* Print a summary containing information
         about the solution for debugging purposes*/
      MSK_solutionsummary(task, MSK_STREAM_MSG);

      if (r == MSK_RES_OK)
      {
        MSKint32t  j;
        MSKsolstae solsta;
        double     *xx = NULL;

        MSK_getsolsta(task, MSK_SOL_ITG, &solsta);

        xx = calloc(numvar, sizeof(double));
        if (xx)
        {
          switch (solsta)
          {
            case MSK_SOL_STA_INTEGER_OPTIMAL:
              MSK_getxx(task,
                        MSK_SOL_ITG,    /* Request the integer solution. */
                        xx);

              printf("Optimal solution.\n");
              for (j = 0; j < numvar; ++j)
                printf("x[%d]: %e\n", j, xx[j]);
              break;
            case MSK_SOL_STA_PRIM_FEAS:
              /* A feasible but not necessarily optimal solution was located. */
              MSK_getxx(task, MSK_SOL_ITG, xx);

              printf("Feasible solution.\n");
              for (j = 0; j < numvar; ++j)
                printf("x[%d]: %e\n", j, xx[j]);
              break;
            case MSK_SOL_STA_UNKNOWN:
              {
                MSKprostae prosta;
                MSK_getprosta(task, MSK_SOL_ITG, &prosta);
                switch (prosta)
                {
                  case MSK_PRO_STA_PRIM_INFEAS_OR_UNBOUNDED:
                    printf("Problem status Infeasible or unbounded\n");
                    break;
                  case MSK_PRO_STA_PRIM_INFEAS:
                    printf("Problem status Infeasible.\n");
                    break;
                  case MSK_PRO_STA_UNKNOWN:
                    printf("Problem status unknown. Termination code %d.\n", trmcode);
                    break;
                  default:
                    printf("Other problem status.");
                    break;
                }
              }
              break;
            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          r = MSK_RES_ERR_SPACE;
        }
        free(xx);
      }
    }

    if (r != MSK_RES_OK)
    {
      /* In case of an error print error code and description. */
      char symname[MSK_MAX_STR_LEN];
      char desc[MSK_MAX_STR_LEN];

      printf("An error occurred while optimizing.\n");
      MSK_getcodedesc(r,
                      symname,
                      desc);
      printf("Error %s - '%s'\n", symname, desc);
    }

    MSK_deletetask(&task);
  }
  MSK_deleteenv(&env);

  printf("Return code: %d.\n", r);
  return (r);
} /* main */

mioinitsol.c

Listing 17.18 mioinitsol.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      mioinitsol.c

   Purpose:   To demonstrate how to solve a MIP with a start guess.

 */

#include "mosek.h"
#include <stdio.h>

static void MSKAPI printstr(void        *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */


int main(int argc, char *argv[])
{
  char         buffer[512];

  const MSKint32t numvar    = 4,
                  numcon    = 1,
                  numintvar = 3;

  MSKrescodee  r;

  MSKenv_t     env;
  MSKtask_t    task;

  double       c[] = { 7.0, 10.0, 1.0, 5.0 };

  MSKboundkeye bkc[] = {MSK_BK_UP};
  double       blc[] = { -MSK_INFINITY};
  double       buc[] = {2.5};

  MSKboundkeye bkx[] = {MSK_BK_LO, MSK_BK_LO, MSK_BK_LO, MSK_BK_LO};
  double       blx[] = {0.0,       0.0,       0.0,      0.0      };
  double       bux[] = {MSK_INFINITY, MSK_INFINITY, MSK_INFINITY, MSK_INFINITY};

  MSKint32t    ptrb[] = {0, 1, 2, 3},
               ptre[] = {1, 2, 3, 4},
               asub[] = {0,   0,   0,   0  };

  double       aval[] = {1.0, 1.0, 1.0, 1.0};
  MSKint32t    intsub[] = {0, 1, 2};
  MSKint32t    j;

  double       xx[4];

  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
    r = MSK_maketask(env, 0, 0, &task);

  if (r == MSK_RES_OK)
    r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

  if (r == MSK_RES_OK)
    r = MSK_inputdata(task,
                      numcon, numvar,
                      numcon, numvar,
                      c,
                      0.0,
                      ptrb,
                      ptre,
                      asub,
                      aval,
                      bkc,
                      blc,
                      buc,
                      bkx,
                      blx,
                      bux);

  if (r == MSK_RES_OK)
    r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE);

  for (j = 0; j < numintvar && r == MSK_RES_OK; ++j)
    r = MSK_putvartype(task, intsub[j], MSK_VAR_TYPE_INT);

  if (r == MSK_RES_OK)
  {
    /* Assign values to integer variables 
       (we only set a slice of xx) */
    double       xxInit[] = {1.0, 1.0, 0.0};
    r = MSK_putxxslice(task, MSK_SOL_ITG, 0, 3, xxInit);
  }
  if (r == MSK_RES_OK)
  {
    /* Request constructing the solution from integer variable values */
    r = MSK_putintparam(task, MSK_IPAR_MIO_CONSTRUCT_SOL, MSK_ON);
  }
  /* solve */

  if (r == MSK_RES_OK)
  {
    MSKrescodee trmcode;
    r = MSK_optimizetrm(task, &trmcode);
    MSK_solutionsummary(task, MSK_STREAM_LOG);
  }

  /* Read back the solution */
  r = MSK_getxx(task, MSK_SOL_ITG, xx);

  if (r == MSK_RES_OK) 
  {
    int constr;
    double constr_obj;

    printf("Solution:\n");
    for(j = 0; j < numvar; j++)
      printf("%f ", xx[j]);
    printf("\n");
    
    /*  Was the initial guess used?     */
    MSK_getintinf(task, MSK_IINF_MIO_CONSTRUCT_SOLUTION, &constr);
    MSK_getdouinf(task, MSK_DINF_MIO_CONSTRUCT_SOLUTION_OBJ, &constr_obj);
    printf("Construct solution utilization: %d\nConstruct solution objective: %.3f\n", constr, constr_obj);
  }
  else
  {
    /* In case of an error print error code and description. */
    char symname[MSK_MAX_STR_LEN];
    char desc[MSK_MAX_STR_LEN];

    printf("An error occurred while optimizing.\n");
    MSK_getcodedesc(r,
                    symname,
                    desc);
    printf("Error %s - '%s'\n", symname, desc);
  }

  MSK_deletetask(&task);
  MSK_deleteenv(&env);

  return (r);
}

opt_server_async.c

Listing 17.19 opt_server_async.c Click here to download.
/*
  Copyright:  Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:       opt_server_async.c

  Purpose :   Demonstrates how to use MOSEK OptServer
              to solve optimization problem asynchronously
*/
#include "mosek.h"
#ifdef _WIN32
#include "windows.h"
#else
#include "unistd.h"
#endif

static void MSKAPI printstr(void *handle, const char str[])
{
  printf("%s", str);
}

int main(int argc, char * argv[])
{

  char token[33];

  int         numpolls = 10;
  int         i = 0;

  MSKbooleant respavailable;

  MSKenv_t    env   = NULL;
  MSKtask_t   task  = NULL;

  MSKrescodee res   = MSK_RES_OK;
  MSKrescodee trm;
  MSKrescodee resp;

  const char * filename = "../data/25fv47.mps";
  const char * addr     = "solve.mosek.com:30080";
  const char * cert = NULL;

  if (argc < 4)
  {
    fprintf(stderr, "Syntax: opt_server_async filename host:port numpolls [cert]\n");
    return 0;
  }

  filename = argv[1];
  addr     = argv[2];
  numpolls = atoi(argv[3]);
  cert     = argc < 5 ? NULL : argv[4];

  res = MSK_makeenv(&env, NULL);

  if (res == MSK_RES_OK)
    res = MSK_maketask(env, 0, 0, &task);

  if (res == MSK_RES_OK)
    res = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

  if (res == MSK_RES_OK)
    res = MSK_readdata(task, filename);

  if (MSK_RES_OK == res && cert)
    res = MSK_putstrparam(task, MSK_SPAR_REMOTE_TLS_CERT_PATH,cert);

  res = MSK_asyncoptimize(task,
                          addr,
                          "",
                          token);
  MSK_deletetask(&task);
  printf("token = %s\n", token);

  if (res == MSK_RES_OK)
    res = MSK_maketask(env, 0, 0, &task);

  if (res == MSK_RES_OK)
    res = MSK_readdata(task, filename);

  if (MSK_RES_OK == res && cert)
    res = MSK_putstrparam(task, MSK_SPAR_REMOTE_TLS_CERT_PATH,cert);

  if (res == MSK_RES_OK)
    res = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

  for (i = 0; i < numpolls &&  res == MSK_RES_OK ; i++)
  {
#if __linux__
    sleep(1);
#elif defined(_WIN32)
    Sleep(1000);
#endif

    puts("+++++++++++++++++++++++++++++++++++++++++++++++++++++\n");
    printf("poll %d\n ", i);


    res = MSK_asyncpoll(task,
                        addr,
                        "",
                        token,
                        &respavailable,
                        &resp,
                        &trm);

    puts("polling done\n");
    puts("+++++++++++++++++++++++++++++++++++++++++++++++++++++\n");

    if (respavailable)
    {
      puts("solution available!");

      res = MSK_asyncgetresult(task,
                               addr,
                               "",
                               token,
                               &respavailable,
                               &resp,
                               &trm);

      MSK_solutionsummary(task, MSK_STREAM_LOG);
      break;
    }

  }


  if (i == numpolls)
  {
    printf("max num polls reached, stopping %s", addr);
    MSK_asyncstop(task, addr, "", token);
  }

  MSK_deletetask(&task);
  MSK_deleteenv(&env);

  printf("%s:%d: Result = %d\n", __FILE__, __LINE__, res); fflush(stdout);

  return res;
}

opt_server_sync.c

Listing 17.20 opt_server_sync.c Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:    opt_server_sync.c

  Purpose :   Demonstrates how to use MOSEK OptServer
              to solve optimization problem synchronously
*/
#include "mosek.h"

static void MSKAPI printstr(void *handle, const char str[])
{
  printf("%s", str);
}

int main(int argc, const char * argv[])
{
  MSKenv_t    env  = NULL;
  MSKtask_t   task = NULL;
  MSKrescodee res  = MSK_RES_OK;
  MSKrescodee trm  = MSK_RES_OK;

  if (argc <= 3)
  {
    fprintf(stderr, "Syntax: opt_server_sync inputfile addr [certpath]\n");    
    return 0;
  }
  else
  {
    const char * taskfile = argv[1];
    const char * address  = argv[2];
    const char * certfile = argc > 3 ? argv[3] : NULL;

    // Create the mosek environment.
    // The `NULL' arguments here, are used to specify customized
    // memory allocators and a memory debug file. These can
    // safely be ignored for now.
    res = MSK_makeenv(&env, NULL);

    // Create a task object linked with the environment env.
    // We create it with 0 variables and 0 constraints initially,
    // since we do not know the size of the problem.
    if (res == MSK_RES_OK)
      res = MSK_maketask(env, 0, 0, &task);

    // Direct the task log stream to a user specified function
    if (res == MSK_RES_OK)
      res = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

    // We assume that a problem file was given as the first command
    // line argument (received in `argv')
    if (res == MSK_RES_OK)
      res = MSK_readdata(task, taskfile);

    // Set OptServer URL
    if (res == MSK_RES_OK)
      res = MSK_putoptserverhost(task, address);
    
    // Path to certificate, if any
    if (MSK_RES_OK == res && certfile)
      res = MSK_putstrparam(task, MSK_SPAR_REMOTE_TLS_CERT_PATH,certfile);

    // Solve the problem remotely
    if (res == MSK_RES_OK)
      res = MSK_optimizetrm(task, &trm);

    // Print a summary of the solution.
    if (res == MSK_RES_OK)
      res = MSK_solutionsummary(task, MSK_STREAM_LOG);

    // Delete task and environment
    MSK_deletetask(&task);
    MSK_deleteenv(&env);
  }
  return res;
}

parallel.c

Listing 17.21 parallel.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      parallel.c

   Purpose: Demonstrates parallel optimization usint optimizebatch()
*/
#include<stdio.h>
#include<stdlib.h>
#include<mosek.h>


/** Example of how to use MSK_optimizebatch(). 
    Optimizes tasks whose names were read from command line.
*/
int main(int argc, char **argv)
{
  MSKenv_t env;
  int n = argc - 1;
  MSKtask_t   *tasks = NULL;
  MSKrescodee *res   = NULL;
  MSKrescodee *trm   = NULL;
  MSKrescodee r = MSK_RES_OK;
  int i;
  /* Size of thread pool available for all tasks */
  int threadpoolsize = 6; 

  tasks = calloc(n, sizeof(MSKtask_t));
  res   = calloc(n, sizeof(MSKrescodee));
  trm   = calloc(n, sizeof(MSKrescodee));

  MSK_makeenv(&env, NULL);

  /* Create an example list of tasks to optimize */
  for (i = 0; i < n; i++) {
    MSK_makeemptytask(env, &(tasks[i]));
    MSK_readdata(tasks[i], argv[i+1]);
    /* We can set the number of threads for each task */
    MSK_putintparam(tasks[i], MSK_IPAR_NUM_THREADS, 2);
  }

  /* Optimize all the given tasks in parallel */
  r = MSK_optimizebatch(env,
                        0,              // No race
                        -1.0,           // No time limit
                        threadpoolsize,
                        n, 
                        tasks,          // Array of tasks to optimize
                        trm,            
                        res);

  for(i = 0; i < n; i++) {
    double obj, tm;
    MSK_getdouinf(tasks[i], MSK_DINF_INTPNT_PRIMAL_OBJ, &obj);
    MSK_getdouinf(tasks[i], MSK_DINF_OPTIMIZER_TIME, &tm);

    printf("Task  %d  res %d   trm %d   obj_val  %.5f  time %.5f\n", 
            i, 
            res[i], 
            trm[i],  
            obj,
            tm);
  }

  for(i = 0; i < n; i++)
    MSK_deletetask(&(tasks[i]));
  free(tasks);
  free(trm);
  free(res);
  MSK_deleteenv(&env);
  return 0;
}

parameters.c

Listing 17.22 parameters.c Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:      parameters.c

 Purpose: Demonstrates a very simple example about how to get/set
 parameters with MOSEK C API
*/

#include "mosek.h"
#include "stdio.h"

int main()
{

  MSKenv_t env = NULL;
  MSKtask_t task = NULL;
  MSKrescodee res, trmcode;

  double param = 0.0;
  double      tm;
  MSKint32t   iter;

  /* Create an environment */
  res = MSK_makeenv(&env, NULL);

  if (res == MSK_RES_OK)
  {
    /* Create a task */
    res = MSK_maketask(env, 0, 0, &task);

    if (res == MSK_RES_OK)
    {
      printf("Test MOSEK parameter get/set functions\n");

      // Set log level (integer parameter)
      res = MSK_putintparam(task, MSK_IPAR_LOG, 1);
      // Select interior-point optimizer... (integer parameter)
      res = MSK_putintparam(task, MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_INTPNT);
      // ... without basis identification (integer parameter)
      res = MSK_putintparam(task, MSK_IPAR_INTPNT_BASIS, MSK_BI_NEVER);
      // Set relative gap tolerance (double parameter)
      res = MSK_putdouparam(task, MSK_DPAR_INTPNT_CO_TOL_REL_GAP, 1.0e-7);

      // The same using explicit string names
      res = MSK_putparam(task, "MSK_DPAR_INTPNT_CO_TOL_REL_GAP", "1.0e-7");
      res = MSK_putnadouparam(task, "MSK_DPAR_INTPNT_CO_TOL_REL_GAP",  1.0e-7);

      // Incorrect value
      res = MSK_putdouparam(task, MSK_DPAR_INTPNT_CO_TOL_REL_GAP,  -1.0);
      if (res != MSK_RES_OK)
        printf("Wrong parameter value\n");

      res = MSK_getdouparam(task, MSK_DPAR_INTPNT_CO_TOL_REL_GAP, &param);
      printf("Current value for parameter MSK_DPAR_INTPNT_CO_TOL_REL_GAP = %e\n", param);

      /* Define and solve an optimization problem here */
      /* MSK_optimizetrm(task, &trmcode) */
      /* After optimization: */

      printf("Get MOSEK information items\n");

      res = MSK_getdouinf(task, MSK_DINF_OPTIMIZER_TIME,     &tm);
      res = MSK_getintinf(task, MSK_IINF_INTPNT_ITER,        &iter);

      printf("Time: %f\nIterations: %d\n", tm, iter);

      MSK_deletetask(&task);
    }
    MSK_deleteenv(&env);
  }
  return 0;
}

pinfeas.c

Listing 17.23 pinfeas.c Click here to download.
//  File : pinfeas.c
//
//  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  Purpose: Demonstrates how to fetch a primal infeasibility certificate
//           for a linear problem
//
#include <stdio.h>
#include "mosek.h"

/* This function prints log output from MOSEK to the terminal. */
static void MSKAPI printstr(void       *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

// Set up a simple linear problem from the manual for test purposes
static MSKrescodee testProblem(MSKtask_t *task) {
  MSKrescodee r = MSK_RES_OK;
  const char data[] = "Task ''\n\
Objective ''\n\
    Minimize + @x0 + 2 @x1 + 5 @x2 + 2 @x3 + @x4 + 2 @x5 + @x6\n\
Constraints\n\
    @c0 [-inf;200] + @x0 + @x1\n\
    @c1 [-inf;1000] + @x2 + @x3\n\
    @c2 [-inf;1000] + @x4 + @x5 + @x6\n\
    @c3 [1100] + @x0 + @x4\n\
    @c4 [200] + @x1\n\
    @c5 [500] + @x2 + @x5\n\
    @c6 [500] + @x3 + @x6\n\
Variables\n\
    @x0 [0;+inf]\n\
    @x1 [0;+inf]\n\
    @x2 [0;+inf]\n\
    @x3 [0;+inf]\n\
    @x4 [0;+inf]\n\
    @x5 [0;+inf]\n\
    @x6 [0;+inf]\n";

  r = MSK_makeemptytask(NULL, task);
  if (r == MSK_RES_OK)
    r = MSK_readptfstring(*task, data);
  return r;
}

// Analyzes and prints infeasibility contributing elements
// n - length of arrays
// sl - dual values for lower bounds
// su - dual values for upper bounds
// eps - tolerance for when a nunzero dual value is significant
static void analyzeCertificate(MSKint32t n, MSKrealt *sl, MSKrealt *su, double eps) {
  MSKint32t i;
  for(i = 0; i < n; i++) {
    if (abs(sl[i]) > eps)
      printf("#%d, lower,  dual = %e\n", i, sl[i]);
    if (abs(su[i]) > eps)
      printf("#%d, upper,  dual = %e\n", i, su[i]);
  }
}

int main(int argc, const char *argv[])
{
  MSKtask_t          task;
  MSKrescodee        r = MSK_RES_OK;
  double             inf = 0.0;

  // In this example we set up a simple problem
  // One could use any task or a task read from a file
  r = testProblem(&task);

  // Useful for debugging
  if (r == MSK_RES_OK)  
    r = MSK_writedata(task, "pinfeas.ptf");                          // Write file in human-readable format
  // Directs the log task stream to the 'printstr' function. 
  if (r == MSK_RES_OK)
    r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

  // Perform the optimization.
  if (r == MSK_RES_OK)
    r = MSK_optimize(task);
  if (r == MSK_RES_OK)
    r = MSK_solutionsummary(task, MSK_STREAM_LOG);

  // Check problem status, we use the interior point solution
  {
    MSKprostae prosta;
    if (r == MSK_RES_OK)
      r = MSK_getprosta(task, MSK_SOL_ITR, &prosta);

    if (r == MSK_RES_OK && prosta == MSK_PRO_STA_PRIM_INFEAS) {
      // Set the tolerance at which we consider a dual value as essential
      double eps = 1e-7;

      printf("Variable bounds important for infeasibility: \n");
      if (r == MSK_RES_OK) {
        MSKint32t m;
        r = MSK_getnumvar(task, &m);
        if (r == MSK_RES_OK) {
          MSKrealt *slx = NULL, *sux = NULL;
          slx = (MSKrealt*) calloc(m, sizeof(MSKrealt));
          sux = (MSKrealt*) calloc(m, sizeof(MSKrealt));
          
          if (r == MSK_RES_OK) r = MSK_getslx(task, MSK_SOL_ITR, slx);
          if (r == MSK_RES_OK) r = MSK_getsux(task, MSK_SOL_ITR, sux);

          if (r == MSK_RES_OK)
            analyzeCertificate(m, slx, sux, eps);

          free(sux);
          free(slx);
        }
      }

      printf("Constraint bounds important for infeasibility: \n");
      if (r == MSK_RES_OK) {
        MSKint32t n;
        r = MSK_getnumcon(task, &n);
        if (r == MSK_RES_OK) {
          MSKrealt *slc = NULL, *suc = NULL;
          slc = (MSKrealt*) calloc(n, sizeof(MSKrealt));
          suc = (MSKrealt*) calloc(n, sizeof(MSKrealt));
          
          if (r == MSK_RES_OK) r = MSK_getslc(task, MSK_SOL_ITR, slc);
          if (r == MSK_RES_OK) r = MSK_getsuc(task, MSK_SOL_ITR, suc);

          if (r == MSK_RES_OK)
            analyzeCertificate(n, slc, suc, eps);

          free(suc);
          free(slc);
        }
      }
    }
    else {
      printf("The problem is not primal infeasible, no certificate to show.\n");
    }
  }

  MSK_deletetask(&task);
  return r;
}

portfolio_1_basic.c

Listing 17.24 portfolio_1_basic.c Click here to download.
/*
  File : portfolio_1_basic.c

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model.
*/
#include <math.h>
#include <stdio.h>
#include "mosek.h"

#define MOSEKCALL(_r,_call)  if ( (_r)==MSK_RES_OK ) (_r) = (_call)

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char **argv)
{
  char            buf[128];

  double          expret  = 0.0,
                  stddev  = 0.0,
                  xj;

  const MSKint32t n       = 8;
  const MSKrealt  gamma   = 36.0;
  const MSKrealt  mu[]    = {0.07197349, 0.15518171, 0.17535435, 0.0898094 , 0.42895777, 0.39291844, 0.32170722, 0.18378628};
  // GT must have size n rows
  const MSKrealt  GT[][8] = {
    {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
    {0.0,     0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
    {0.0,     0.0,     0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
    {0.0,     0.0,     0.0,     0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
    {0.0,     0.0,     0.0,     0.0,     0.36096, 0.12574, 0.10157, 0.0571 },
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.21552, 0.05663, 0.06187},
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.22514, 0.03327},
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.2202 }
  };
  const MSKint32t k       = sizeof(GT) / (n * sizeof(MSKrealt));
  const MSKrealt  x0[]    = {8.0, 5.0, 3.0, 5.0, 2.0, 9.0, 3.0, 6.0};
  const MSKrealt  w       = 59;
  MSKrealt        totalBudget;

  MSKenv_t        env;
  MSKint32t       i, j, *sub;
  MSKrescodee     res = MSK_RES_OK, trmcode;
  MSKtask_t       task;

  //Offset of variables into the API variable.
  MSKint32t numvar = n;
  MSKint32t voff_x = 0;

  // Constraints offsets
  MSKint32t numcon = 1;
  MSKint32t coff_bud = 0;

  /* Initial setup. */
  env  = NULL;
  task = NULL;
  MOSEKCALL(res, MSK_makeenv(&env, NULL));
  MOSEKCALL(res, MSK_maketask(env, 0, 0, &task));
  MOSEKCALL(res, MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr));

  // Holding variable x of length n
  // No other auxiliary variables are needed in this formulation
  MOSEKCALL(res, MSK_appendvars(task, numvar));
  // Setting up variable x 
  for (j = 0; j < n; ++j)
  {
    /* Optionally we can give the variables names */
    sprintf(buf, "x[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_x + j, buf));
    /* No short-selling - x^l = 0, x^u = inf */
    MOSEKCALL(res, MSK_putvarbound(task, voff_x + j, MSK_BK_LO, 0.0, MSK_INFINITY));
  }

  // One linear constraint: total budget
  MOSEKCALL(res, MSK_appendcons(task, 1));
  sprintf(buf, "%s", "budget");
  MOSEKCALL(res, MSK_putconname(task, coff_bud, buf));
  for (j = 0; j < n; ++j)
  {
    /* Coefficients in the first row of A */
    MOSEKCALL(res, MSK_putaij(task, coff_bud, voff_x + j, 1.0));
  }
  totalBudget = w;
  for (i = 0; i < n; ++i)
  {
    totalBudget += x0[i];
  }
  MOSEKCALL(res, MSK_putconbound(task, coff_bud, MSK_BK_FX, totalBudget, totalBudget));

  // Input (gamma, GTx) in the AFE (affine expression) storage
  // We need k+1 rows
  MOSEKCALL(res, MSK_appendafes(task, k + 1));
  // The first affine expression = gamma
  MOSEKCALL(res, MSK_putafeg(task, 0, gamma));
  // The remaining k expressions comprise GT*x, we add them row by row
  // In more realisic scenarios it would be better to extract nonzeros and input in sparse form
  MSKint32t* vslice_x = (MSKint32t*) malloc(n * sizeof(MSKint32t));
  for (i = 0; i < n; ++i)
  {
    vslice_x[i] = voff_x + i;
  } 
  for (i = 0; i < k; ++i)
  {
    MOSEKCALL(res, MSK_putafefrow(task, i + 1, n, vslice_x, GT[i]));
  }
  free(vslice_x);

  // Input the affine conic constraint (gamma, GT*x) \in QCone
  // Add the quadratic domain of dimension k+1
  MSKint64t qdom;
  MOSEKCALL(res, MSK_appendquadraticconedomain(task, k + 1, &qdom));
  // Add the constraint
  MOSEKCALL(res, MSK_appendaccseq(task, qdom, k + 1, 0, NULL));
  sprintf(buf, "%s", "risk");            
  MOSEKCALL(res, MSK_putaccname(task, 0, buf));

  // Objective: maximize expected return mu^T x
  for (j = 0; j < n; ++j)
  {
    MOSEKCALL(res, MSK_putcj(task, voff_x + j, mu[j]));
  }
  MOSEKCALL(res, MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE));

#if 0
  /* No log output */
  MOSEKCALL(res, MSK_putintparam(task, MSK_IPAR_LOG, 0));
#endif

#if 0
  /* Dump the problem to a human readable PTF file. */
  MOSEKCALL(res, MSK_writedata(task, "dump.ptf"));
#endif

  MOSEKCALL(res, MSK_optimizetrm(task, &trmcode));

#if 1
  /* Display the solution summary for quick inspection of results. */
  MSK_solutionsummary(task, MSK_STREAM_MSG);
#endif

  if ( res == MSK_RES_OK )
  {
    expret = 0.0;
    stddev = 0.0;

    /* Read the x variables one by one and compute expected return. */
    /* Can also be obtained as value of the objective. */
    for (j = 0; j < n; ++j)
    {
      MOSEKCALL(res, MSK_getxxslice(task, MSK_SOL_ITR, voff_x + j, voff_x + j + 1, &xj));
      expret += mu[j] * xj;
    }

    /* Read the value of s. This should be gamma. */
    printf("\nExpected return %e for gamma %e\n", expret, gamma);
  }

  if ( task != NULL )
    MSK_deletetask(&task);

  if ( env != NULL )
    MSK_deleteenv(&env);

  return ( 0 );
}

portfolio_2_frontier.c

Listing 17.25 portfolio_2_frontier.c Click here to download.
/*
  File : portfolio_2_frontier.c

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model.
                 Determines points on the efficient frontier.
 */
#include <math.h>
#include <stdio.h>

#include "mosek.h"

#define MOSEKCALL(_r,_call)  if ( (_r)==MSK_RES_OK ) (_r) = (_call)
#define LOGLEVEL             0

static void MSKAPI printstr(void       *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */


int main(int argc, const char **argv)
{
  char buf[128];
  const MSKint32t n       = 8;
  const MSKrealt  mu[]    = {0.07197349, 0.15518171, 0.17535435, 0.0898094 , 0.42895777, 0.39291844, 0.32170722, 0.18378628};
  // GT must have size n rows
  const MSKrealt  GT[][8] = {
    {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
    {0.0,     0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
    {0.0,     0.0,     0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
    {0.0,     0.0,     0.0,     0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
    {0.0,     0.0,     0.0,     0.0,     0.36096, 0.12574, 0.10157, 0.0571 },
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.21552, 0.05663, 0.06187},
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.22514, 0.03327},
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.2202 }
  };
  const MSKint32t k       = sizeof(GT) / (n * sizeof(MSKrealt));
  const MSKrealt  x0[]    = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  const MSKrealt  w       = 1.0;
  const MSKrealt alphas[] = {0.0, 0.01, 0.1, 0.25, 0.30, 0.35, 0.4, 0.45, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 10.0};
  const MSKint32t numalpha = 15;
  MSKrealt totalBudget;

  MSKenv_t        env;
  MSKint32t       i, j;
  MSKrescodee     res = MSK_RES_OK, lres;
  MSKtask_t       task;
  MSKrealt        xj;
  MSKsolstae      solsta;

  //Offset of variables into the API variable.
  MSKint32t numvar = n + 1;
  MSKint32t voff_x = 0;
  MSKint32t voff_s = n;

  // Constraints offsets
  MSKint32t numcon = 1;
  MSKint32t coff_bud = 0;

  /* Initial setup. */
  env  = NULL;
  task = NULL;

  /* Replace "" with NULL in production. */
  MOSEKCALL(res, MSK_makeenv(&env, "")); 

  MOSEKCALL(res, MSK_maketask(env, 0, 0, &task));
  MOSEKCALL(res, MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr));

  // Holding variable x of length n
  MOSEKCALL(res, MSK_appendvars(task, numvar));
  // Setting up variable x 
  for (j = 0; j < n; ++j)
  {
    /* Optionally we can give the variables names */
    sprintf(buf, "x[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_x + j, buf));
    /* No short-selling - x^l = 0, x^u = inf */
    MOSEKCALL(res, MSK_putvarbound(task, voff_x + j, MSK_BK_LO, 0.0, MSK_INFINITY));
  }
  sprintf(buf, "s");
  MOSEKCALL(res, MSK_putvarname(task, voff_s, buf));
  MOSEKCALL(res, MSK_putvarbound(task, voff_s, MSK_BK_FR, -MSK_INFINITY, MSK_INFINITY));

  // One linear constraint: total budget
  MOSEKCALL(res, MSK_appendcons(task, 1));
  sprintf(buf, "%s", "budget");
  MOSEKCALL(res, MSK_putconname(task, coff_bud, buf));
  for (j = 0; j < n; ++j)
  {
    /* Coefficients in the first row of A */
    MOSEKCALL(res, MSK_putaij(task, coff_bud, voff_x + j, 1.0));
  }
  totalBudget = w;
  for (i = 0; i < n; ++i)
  {
    totalBudget += x0[i];
  }
  MOSEKCALL(res, MSK_putconbound(task, coff_bud, MSK_BK_FX, totalBudget, totalBudget));

  // Input (gamma, GTx) in the AFE (affine expression) storage
  // We build the following F and g for variables [x, s]:
  //     [0, 1]      [0  ]
  // F = [0, 0], g = [0.5]
  //     [GT,0]      [0  ]
  // We need k+2 rows
  MOSEKCALL(res, MSK_appendafes(task, k + 2));
  // The first affine expression is variable s (last variable, index n)
  MOSEKCALL(res, MSK_putafefentry(task, 0, n, 1.0));
  // The second affine expression is constant 0.5
  MOSEKCALL(res, MSK_putafeg(task, 1, 0.5));
  // The remaining k expressions comprise GT*x, we add them row by row
  // In more realisic scenarios it would be better to extract nonzeros and input in sparse form
  MSKint32t* vslice_x = (MSKint32t*) malloc(n * sizeof(MSKint32t));
  for (i = 0; i < n; ++i)
  {
    vslice_x[i] = voff_x + i;
  } 
  for (i = 0; i < k; ++i)
  {
    MOSEKCALL(res, MSK_putafefrow(task, i + 2, n, vslice_x, GT[i]));
  }
  free(vslice_x);

  // Input the affine conic constraint (gamma, GT*x) \in QCone
  // Add the quadratic domain of dimension k+1
  MSKint64t rqdom;
  MOSEKCALL(res, MSK_appendrquadraticconedomain(task, k + 2, &rqdom));
  // Add the constraint
  MOSEKCALL(res, MSK_appendaccseq(task, rqdom, k + 2, 0, NULL));
  sprintf(buf, "%s", "risk");            
  MOSEKCALL(res, MSK_putaccname(task, 0, buf));

  // Objective: maximize expected return mu^T x
  for (j = 0; j < n; ++j)
  {
    MOSEKCALL(res, MSK_putcj(task, voff_x + j, mu[j]));
  }
  MOSEKCALL(res, MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE));

  /* Set the log level */
  MOSEKCALL(res, MSK_putintparam(task, MSK_IPAR_LOG, LOGLEVEL));

  printf("%-12s  %-12s  %-12s\n", "alpha", "exp ret", "std. dev");

  for (i = 0; i < numalpha && res==MSK_RES_OK; ++i)
  {
    const double alpha = alphas[i];
    MSKrescodee  trmcode;

    /* Sets the objective function coefficient for s. */
    MOSEKCALL(res, MSK_putcj(task, voff_s + 0, -alpha));

    MOSEKCALL(res, MSK_optimizetrm(task, &trmcode));

    MOSEKCALL(res, MSK_getsolsta(task, MSK_SOL_ITR, &solsta));

    if (solsta == MSK_SOL_STA_OPTIMAL)
    {
      double expret = 0.0,
             stddev;

      for (j = 0; j < n; ++j)
      {
        MOSEKCALL(res, MSK_getxxslice(task, MSK_SOL_ITR, voff_x + j, voff_x + j + 1, &xj));
        expret += mu[j] * xj;
      }

      MOSEKCALL(res, MSK_getxxslice(task, MSK_SOL_ITR, voff_s + 0, voff_s + 1, &stddev));

      printf("%-12.3e  %-12.3e  %-12.3e\n", alpha, expret, sqrt(stddev));
    }
    else
    {
      printf("An error occurred when solving for alpha=%e\n", alpha);
    }
  }
  lres = MSK_deletetask(&task);
  res  = res==MSK_RES_OK ? lres : res;
  
  lres = MSK_deleteenv(&env);
  res  = res==MSK_RES_OK ? lres : res;

  return ( res );
}

portfolio_3_impact.c

Listing 17.26 portfolio_3_impact.c Click here to download.
/*
  File : portfolio_3_impact.c

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model
                 with transaction costs of order x^(3/2).
 */
#include <math.h>
#include <stdio.h>

#include "mosek.h"

#define MOSEKCALL(_r,_call)  if ( (_r)==MSK_RES_OK ) (_r) = (_call)

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char **argv)
{
  char            buf[128];
  const MSKint32t n       = 8;
  const MSKrealt  mu[]    = {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379};
  // GT must have size n rows
  const MSKrealt  GT[][8] = {
    {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
    {0.0,     0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
    {0.0,     0.0,     0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
    {0.0,     0.0,     0.0,     0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
    {0.0,     0.0,     0.0,     0.0,     0.36096, 0.12574, 0.10157, 0.0571 },
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.21552, 0.05663, 0.06187},
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.22514, 0.03327},
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.2202 }
  };
  const MSKint32t k       = sizeof(GT) / (n * sizeof(MSKrealt));
  const MSKrealt  x0[]    = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  const MSKrealt  w       = 1.0;
  const MSKrealt  gamma   = 0.36;
  MSKrealt totalBudget;
  MSKint32t       i, j;

  MSKrealt* m = (MSKrealt*) malloc(n * sizeof(MSKrealt));
  for (i = 0; i < n; ++i)
  {
    m[i] = 0.01;
  }

  // Offset of variables into the API variable.
  MSKint32t numvar = 3 * n;
  MSKint32t voff_x = 0;
  MSKint32t voff_c = n;
  MSKint32t voff_z = 2 * n;

  // Offset of constraints.
  MSKint32t numcon = 2 * n + 1;
  MSKint32t coff_bud = 0;
  MSKint32t coff_abs1 = 1; 
  MSKint32t coff_abs2 = 1 + n;

  double          expret,
                  stddev,
                  xj;
  MSKenv_t        env;
  MSKrescodee     res = MSK_RES_OK, trmcode;
  MSKtask_t       task;

  /* Initial setup. */
  env  = NULL;
  task = NULL;
  MOSEKCALL(res, MSK_makeenv(&env, NULL));
  MOSEKCALL(res, MSK_maketask(env, 0, 0, &task));
  MOSEKCALL(res, MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr));

  // Variables (vector of x, c, z)
  MOSEKCALL(res, MSK_appendvars(task, numvar));
  for (j = 0; j < n; ++j)
  {
    /* Optionally we can give the variables names */
    sprintf(buf, "x[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_x + j, buf));
    sprintf(buf, "c[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_c + j, buf));
    sprintf(buf, "z[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_z + j, buf));
    /* Apply variable bounds (x >= 0, c and z free) */
    MOSEKCALL(res, MSK_putvarbound(task, voff_x + j, MSK_BK_LO, 0.0, MSK_INFINITY));
    MOSEKCALL(res, MSK_putvarbound(task, voff_c + j, MSK_BK_FR, -MSK_INFINITY, MSK_INFINITY));
    MOSEKCALL(res, MSK_putvarbound(task, voff_z + j, MSK_BK_FR, -MSK_INFINITY, MSK_INFINITY));
  }

  // Linear constraints
  // - Total budget
  MOSEKCALL(res, MSK_appendcons(task, 1));
  sprintf(buf, "%s", "budget");
  MOSEKCALL(res, MSK_putconname(task, coff_bud, buf));
  for (j = 0; j < n; ++j)
  {
    /* Coefficients in the first row of A */
    MOSEKCALL(res, MSK_putaij(task, coff_bud, voff_x + j, 1.0));
    MOSEKCALL(res, MSK_putaij(task, coff_bud, voff_c + j, m[j]));
  }
  totalBudget = w;
  for (i = 0; i < n; ++i)
  {
    totalBudget += x0[i];
  }
  MOSEKCALL(res, MSK_putconbound(task, coff_bud, MSK_BK_FX, totalBudget, totalBudget));
  
  // - Absolute value
  MOSEKCALL(res, MSK_appendcons(task, 2 * n));
  for (i = 0; i < n; ++i)
  {
    sprintf(buf, "zabs1[%d]", 1 + i);
    MOSEKCALL(res, MSK_putconname(task, coff_abs1 + i, buf));
    MOSEKCALL(res, MSK_putaij(task, coff_abs1 + i, voff_x + i, -1.0));
    MOSEKCALL(res, MSK_putaij(task, coff_abs1 + i, voff_z + i, 1.0));
    MOSEKCALL(res, MSK_putconbound(task, coff_abs1 + i, MSK_BK_LO, -x0[i], MSK_INFINITY));
    sprintf(buf, "zabs2[%d]", 1 + i);
    MOSEKCALL(res, MSK_putconname(task, coff_abs2 + i, buf));
    MOSEKCALL(res, MSK_putaij(task, coff_abs2 + i, voff_x + i, 1.0));
    MOSEKCALL(res, MSK_putaij(task, coff_abs2 + i, voff_z + i, 1.0));
    MOSEKCALL(res, MSK_putconbound(task, coff_abs2 + i, MSK_BK_LO, x0[i], MSK_INFINITY));          
  } 

  // ACCs
  MSKint32t aoff_q = 0;
  MSKint32t aoff_pow = k + 1;
  // - (gamma, GTx) in Q(k+1)
  // The part of F and g for variable x:
  //     [0,  0, 0]      [gamma]
  // F = [GT, 0, 0], g = [0    ]    
  MOSEKCALL(res, MSK_appendafes(task, k + 1));
  MOSEKCALL(res, MSK_putafeg(task, aoff_q, gamma));
  MSKint32t* vslice_x = (MSKint32t*) malloc(n * sizeof(MSKint32t));
  for (i = 0; i < n; ++i)
  {
    vslice_x[i] = voff_x + i;
  } 
  for (i = 0; i < k; ++i)
  {
    MOSEKCALL(res, MSK_putafefrow(task, aoff_q + i + 1, n, vslice_x, GT[i]));
  }
  free(vslice_x);
  MSKint64t qdom;
  MOSEKCALL(res, MSK_appendquadraticconedomain(task, k + 1, &qdom));
  MOSEKCALL(res, MSK_appendaccseq(task, qdom, k + 1, aoff_q, NULL));
  sprintf(buf, "%s", "risk");            
  MOSEKCALL(res, MSK_putaccname(task, aoff_q, buf));

  // - (c_j, 1, z_j) in P3(2/3, 1/3)
  // The part of F and g for variables [c, z]:
  //     [0, I, 0]      [0]
  // F = [0, 0, I], g = [0]
  //     [0, 0, 0]      [1]
  MOSEKCALL(res, MSK_appendafes(task, 2 * n + 1));
  for (i = 0; i < n; ++i)
  {
    MOSEKCALL(res, MSK_putafefentry(task, aoff_pow + i, voff_c + i, 1.0));
    MOSEKCALL(res, MSK_putafefentry(task, aoff_pow + n + i, voff_z + i, 1.0));
  } 
  MOSEKCALL(res, MSK_putafeg(task, aoff_pow + 2 * n, 1.0));
  // We use one row from F and g for both c_j and z_j, and the last row of F and g for the constant 1.
  // NOTE: Here we reuse the last AFE and the power cone n times, but we store them only once.
  MSKrealt exponents[] = {2, 1};
  MSKint64t powdom;
  MOSEKCALL(res, MSK_appendprimalpowerconedomain(task, 3, 2, exponents, &powdom));
  MSKint64t* flat_afe_list = (MSKint64t*) malloc(3 * n * sizeof(MSKint64t));
  MSKint64t* dom_list = (MSKint64t*) malloc(n * sizeof(MSKint64t));
  for (i = 0; i < n; ++i)
  {
    flat_afe_list[3 * i + 0] = aoff_pow + i;
    flat_afe_list[3 * i + 1] = aoff_pow + 2 * n;
    flat_afe_list[3 * i + 2] = aoff_pow + n + i;
    dom_list[i] = powdom;
  }
  MOSEKCALL(res, MSK_appendaccs(task, n, dom_list, 3 * n, flat_afe_list, NULL));
  free(flat_afe_list);
  free(dom_list);
  for (i = 0; i < n; ++i)
  {
    sprintf(buf, "market_impact[%d]", i);            
    MOSEKCALL(res, MSK_putaccname(task, i + 1, buf));
  }             

  // Objective: maximize expected return mu^T x
  for (j = 0; j < n; ++j)
  {
    MOSEKCALL(res, MSK_putcj(task, voff_x + j, mu[j]));
  }
  MOSEKCALL(res, MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE));

  #if 0
    /* no log output. */
    MOSEKCALL(res, MSK_putintparam(task, MSK_IPAR_LOG, 0));
  #endif

  #if 0
    /* Dump the problem to a human readable OPF file. */
    MOSEKCALL(res, MSK_writedata(task, "dump.ptf"));
  #endif

    MOSEKCALL(res, MSK_optimizetrm(task, &trmcode));

    /* Display the solution summary for quick inspection of results. */
  #if 1
    MSK_solutionsummary(task, MSK_STREAM_MSG);
  #endif

  if (res == MSK_RES_OK)
  {
    expret = 0.0;
    stddev = 0.0;

    for (j = 0; j < n; ++j)
    {
      MOSEKCALL(res, MSK_getxxslice(task, MSK_SOL_ITR, voff_x + j, voff_x + j + 1, &xj));
      expret += mu[j] * xj;
    }

    printf("\nExpected return %e for gamma %e\n", expret, gamma);
  }

  MSK_deletetask(&task);
  MSK_deleteenv(&env);

  free(m);

  return (0);
}

portfolio_4_transcost.c

Listing 17.27 portfolio_4_transcost.c Click here to download.
/*
  File : portfolio_4_transcost.c

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model
                 with fixed setup costs and transaction costs
                 as a mixed-integer problem.
 */
#include <math.h>
#include <stdio.h>

#include "mosek.h"

#define MOSEKCALL(_r,_call)  if ( (_r)==MSK_RES_OK ) (_r) = (_call)

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char **argv)
{
  char            buf[128];
  const MSKint32t n       = 8;
  const MSKrealt  mu[]    = {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379};
  // GT must have size n rows
  const MSKrealt  GT[][8] = {
    {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
    {0.0,     0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
    {0.0,     0.0,     0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
    {0.0,     0.0,     0.0,     0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
    {0.0,     0.0,     0.0,     0.0,     0.36096, 0.12574, 0.10157, 0.0571 },
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.21552, 0.05663, 0.06187},
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.22514, 0.03327},
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.2202 }
  };
  const MSKint32t k       = sizeof(GT) / (n * sizeof(MSKrealt));
  const MSKrealt  x0[]    = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  const MSKrealt  w       = 1.0;
  const MSKrealt  gamma   = 0.36;

  MSKint32t       i, j;  
  MSKrealt* f = (MSKrealt*) malloc(n * sizeof(MSKrealt));
  MSKrealt* g = (MSKrealt*) malloc(n * sizeof(MSKrealt));
  for (i = 0; i < n; ++i)
  {
    f[i] = 0.01;
    g[i] = 0.001;
  }

  // Offset of variables.
  MSKint32t numvar = 3 * n;
  MSKint32t voff_x = 0;
  MSKint32t voff_z = n;
  MSKint32t voff_y = 2 * n;

  // Offset of constraints.
  MSKint32t numcon = 3 * n + 1;
  MSKint32t coff_bud = 0;
  MSKint32t coff_abs1 = 1;
  MSKint32t coff_abs2 = 1 + n;
  MSKint32t coff_swi = 1 + 2 * n; 

  double          expret,
                  xj;
  MSKenv_t        env;
  MSKrescodee     res = MSK_RES_OK, trmcode;
  MSKtask_t       task;

  /* Initial setup. */
  env  = NULL;
  task = NULL;
  MOSEKCALL(res, MSK_makeenv(&env, NULL));
  MOSEKCALL(res, MSK_maketask(env, 0, 0, &task));
  MOSEKCALL(res, MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr));

  // Variables (vector of x, z, y)
  MOSEKCALL(res, MSK_appendvars(task, numvar));
  for (j = 0; j < n; ++j)
  {
    /* Optionally we can give the variables names */
    sprintf(buf, "x[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_x + j, buf));
    sprintf(buf, "z[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_z + j, buf));
    sprintf(buf, "y[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_y + j, buf));
    /* Apply variable bounds (x >= 0, z free, y binary) */
    MOSEKCALL(res, MSK_putvarbound(task, voff_x + j, MSK_BK_LO, 0.0, MSK_INFINITY));
    MOSEKCALL(res, MSK_putvarbound(task, voff_z + j, MSK_BK_FR, -MSK_INFINITY, MSK_INFINITY));
    MOSEKCALL(res, MSK_putvarbound(task, voff_y + j, MSK_BK_RA, 0.0, 1.0));
    MOSEKCALL(res, MSK_putvartype(task, voff_y + j, MSK_VAR_TYPE_INT));
  }

  // Linear constraints
  // - Total budget
  MOSEKCALL(res, MSK_appendcons(task, 1));
  sprintf(buf, "%s", "budget");
  MOSEKCALL(res, MSK_putconname(task, coff_bud, buf));
  for (j = 0; j < n; ++j)
  {
    /* Coefficients in the first row of A */
    MOSEKCALL(res, MSK_putaij(task, coff_bud, voff_x + j, 1.0));
    MOSEKCALL(res, MSK_putaij(task, coff_bud, voff_z + j, g[j]));
    MOSEKCALL(res, MSK_putaij(task, coff_bud, voff_y + j, f[j]));
  }
  MSKrealt U = w;
  for (i = 0; i < n; ++i)
  {
    U += x0[i];
  }
  MOSEKCALL(res, MSK_putconbound(task, coff_bud, MSK_BK_FX, U, U));
  
  // - Absolute value
  MOSEKCALL(res, MSK_appendcons(task, 2 * n));
  for (i = 0; i < n; ++i)
  {
    sprintf(buf, "zabs1[%d]", 1 + i);
    MOSEKCALL(res, MSK_putconname(task, coff_abs1 + i, buf));
    MOSEKCALL(res, MSK_putaij(task, coff_abs1 + i, voff_x + i, -1.0));
    MOSEKCALL(res, MSK_putaij(task, coff_abs1 + i, voff_z + i, 1.0));
    MOSEKCALL(res, MSK_putconbound(task, coff_abs1 + i, MSK_BK_LO, -x0[i], MSK_INFINITY));
    sprintf(buf, "zabs2[%d]", 1 + i);
    MOSEKCALL(res, MSK_putconname(task, coff_abs2 + i, buf));
    MOSEKCALL(res, MSK_putaij(task, coff_abs2 + i, voff_x + i, 1.0));
    MOSEKCALL(res, MSK_putaij(task, coff_abs2 + i, voff_z + i, 1.0));
    MOSEKCALL(res, MSK_putconbound(task, coff_abs2 + i, MSK_BK_LO, x0[i], MSK_INFINITY));
  }

  // - Switch 
  MOSEKCALL(res, MSK_appendcons(task, n));
  for (i = 0; i < n; ++i)
  {
    sprintf(buf, "switch[%d]", 1 + i);
    MOSEKCALL(res, MSK_putconname(task, coff_swi + i, buf));
    MOSEKCALL(res, MSK_putaij(task, coff_swi + i, voff_z + i, 1.0));
    MOSEKCALL(res, MSK_putaij(task, coff_swi + i, voff_y + i, -U));
    MOSEKCALL(res, MSK_putconbound(task, coff_swi + i, MSK_BK_UP, -MSK_INFINITY, 0.0));
  }      

  // ACCs
  MSKint32t aoff_q = 0;
  // - (gamma, GTx) in Q(k+1)
  // The part of F and g for variable x:
  //     [0,  0, 0]      [gamma]
  // F = [GT, 0, 0], g = [0    ]    
  MOSEKCALL(res, MSK_appendafes(task, k + 1));
  MOSEKCALL(res, MSK_putafeg(task, aoff_q, gamma));
  MSKint32t* vslice_x = (MSKint32t*) malloc(n * sizeof(MSKint32t));
  for (i = 0; i < n; ++i)
  {
    vslice_x[i] = voff_x + i;
  } 
  for (i = 0; i < k; ++i)
  {
    MOSEKCALL(res, MSK_putafefrow(task, aoff_q + i + 1, n, vslice_x, GT[i]));
  }
  free(vslice_x);
  MSKint64t qdom;
  MOSEKCALL(res, MSK_appendquadraticconedomain(task, k + 1, &qdom));
  MOSEKCALL(res, MSK_appendaccseq(task, qdom, k + 1, aoff_q, NULL));
  sprintf(buf, "%s", "risk");            
  MOSEKCALL(res, MSK_putaccname(task, aoff_q, buf));

  // Objective: maximize expected return mu^T x
  for (j = 0; j < n; ++j)
  {
    MOSEKCALL(res, MSK_putcj(task, voff_x + j, mu[j]));
  }
  MOSEKCALL(res, MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE));

#if 0
  /* no log output. */
  MOSEKCALL(res, MSK_putintparam(task, MSK_IPAR_LOG, 0));
#endif

#if 0
  /* Dump the problem to a human readable OPF file. */
  MOSEKCALL(res, MSK_writedata(task, "dump.ptf"));
#endif

  MOSEKCALL(res, MSK_optimizetrm(task, &trmcode));

#if 1
  /* Display the solution summary for quick inspection of results. */
  MSK_solutionsummary(task, MSK_STREAM_MSG);
#endif
  
  if (res == MSK_RES_OK)
  {
    expret = 0.0;

    for (j = 0; j < n; ++j)
    {
      MOSEKCALL(res, MSK_getxxslice(task, MSK_SOL_ITG, voff_x + j, voff_x + j + 1, &xj));
      expret += mu[j] * xj;
    }

    printf("\nExpected return %e for gamma %e\n", expret, gamma);
  }

  MSK_deletetask(&task);
  MSK_deleteenv(&env);

  free(f);
  free(g);

  return (0);
}

portfolio_5_card.c

Listing 17.28 portfolio_5_card.c Click here to download.
/*
  File : portfolio_5_card.c

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model
                 with cardinality constraints on number of assets traded.
 */
#include <math.h>
#include <stdio.h>

#include "mosek.h"

#define MOSEKCALL(_r,_call)  if ( (_r)==MSK_RES_OK ) (_r) = (_call)

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */
  
MSKrescodee markowitz_with_card(const int        n,
                                const int        k,
                                const double     x0[],
                                const double     w,
                                const double     gamma,
                                const double     mu[],
                                const double     GT[][8],
                                const int        K,
                                      double     *xx)
{

  // Offset of variables.
  MSKint32t numvar = 3 * n;
  MSKint32t voff_x = 0; 
  MSKint32t voff_z = n;
  MSKint32t voff_y = 2 * n;

  // Offset of constraints.
  MSKint32t numcon = 3 * n + 2;
  MSKint32t coff_bud = 0;
  MSKint32t coff_abs1 = 1;
  MSKint32t coff_abs2 = 1 + n;
  MSKint32t coff_swi = 1 + 2 * n;
  MSKint32t coff_card = 1 + 3 * n;  

  char            buf[128];
  MSKenv_t        env;
  MSKint32t       i, j;
  MSKrescodee     res = MSK_RES_OK, trmcode;
  MSKtask_t       task;
  
  /* Initial setup. */
  env  = NULL;
  task = NULL;
  MOSEKCALL(res, MSK_makeenv(&env, NULL));
  MOSEKCALL(res, MSK_maketask(env, 0, 0, &task));
  MOSEKCALL(res, MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr));

  // Variables (vector of x, z, y)
  MOSEKCALL(res, MSK_appendvars(task, numvar));
  for (j = 0; j < n; ++j)
  {
    /* Optionally we can give the variables names */
    sprintf(buf, "x[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_x + j, buf));
    sprintf(buf, "z[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_z + j, buf));
    sprintf(buf, "y[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_y + j, buf));
    /* Apply variable bounds (x >= 0, z free, y binary) */
    MOSEKCALL(res, MSK_putvarbound(task, voff_x + j, MSK_BK_LO, 0.0, MSK_INFINITY));
    MOSEKCALL(res, MSK_putvarbound(task, voff_z + j, MSK_BK_FR, -MSK_INFINITY, MSK_INFINITY));
    MOSEKCALL(res, MSK_putvarbound(task, voff_y + j, MSK_BK_RA, 0.0, 1.0));
    MOSEKCALL(res, MSK_putvartype(task, voff_y + j, MSK_VAR_TYPE_INT));
  }

  // Linear constraints
  // - Total budget
  MOSEKCALL(res, MSK_appendcons(task, 1));
  sprintf(buf, "%s", "budget");
  MOSEKCALL(res, MSK_putconname(task, coff_bud, buf));
  for (j = 0; j < n; ++j)
  {
    /* Coefficients in the first row of A */
    MOSEKCALL(res, MSK_putaij(task, coff_bud, voff_x + j, 1.0));
  }
  MSKrealt U = w;
  for (i = 0; i < n; ++i)
  {
    U += x0[i];
  }
  MOSEKCALL(res, MSK_putconbound(task, coff_bud, MSK_BK_FX, U, U));
  
  // - Absolute value
  MOSEKCALL(res, MSK_appendcons(task, 2 * n));
  for (i = 0; i < n; ++i)
  {
    sprintf(buf, "zabs1[%d]", 1 + i);
    MOSEKCALL(res, MSK_putconname(task, coff_abs1 + i, buf));
    MOSEKCALL(res, MSK_putaij(task, coff_abs1 + i, voff_x + i, -1.0));
    MOSEKCALL(res, MSK_putaij(task, coff_abs1 + i, voff_z + i, 1.0));
    MOSEKCALL(res, MSK_putconbound(task, coff_abs1 + i, MSK_BK_LO, -x0[i], MSK_INFINITY));
    sprintf(buf, "zabs2[%d]", 1 + i);
    MOSEKCALL(res, MSK_putconname(task, coff_abs2 + i, buf));
    MOSEKCALL(res, MSK_putaij(task, coff_abs2 + i, voff_x + i, 1.0));
    MOSEKCALL(res, MSK_putaij(task, coff_abs2 + i, voff_z + i, 1.0));
    MOSEKCALL(res, MSK_putconbound(task, coff_abs2 + i, MSK_BK_LO, x0[i], MSK_INFINITY));
  }

  // - Switch 
  MOSEKCALL(res, MSK_appendcons(task, n));
  for (i = 0; i < n; ++i)
  {
    sprintf(buf, "switch[%d]", 1 + i);
    MOSEKCALL(res, MSK_putconname(task, coff_swi + i, buf));
    MOSEKCALL(res, MSK_putaij(task, coff_swi + i, voff_z + i, 1.0));
    MOSEKCALL(res, MSK_putaij(task, coff_swi + i, voff_y + i, -U));
    MOSEKCALL(res, MSK_putconbound(task, coff_swi + i, MSK_BK_UP, -MSK_INFINITY, 0.0));
  }    

  // - Cardinality
  MOSEKCALL(res, MSK_appendcons(task, 1));
  sprintf(buf, "cardinality");
  MOSEKCALL(res, MSK_putconname(task, coff_card, buf));  
  for (i = 0; i < n; ++i)
  {
    MOSEKCALL(res, MSK_putaij(task, coff_card, voff_y + i, 1.0));
    MOSEKCALL(res, MSK_putconbound(task, coff_card, MSK_BK_UP, -MSK_INFINITY, K));
  }

  // ACCs
  MSKint32t aoff_q = 0;
  // - (gamma, GTx) in Q(k+1)
  // The part of F and g for variable x:
  //     [0,  0, 0]      [gamma]
  // F = [GT, 0, 0], g = [0    ]    
  MOSEKCALL(res, MSK_appendafes(task, k + 1));
  MOSEKCALL(res, MSK_putafeg(task, aoff_q, gamma));
  MSKint32t* vslice_x = (MSKint32t*) malloc(n * sizeof(MSKint32t));
  for (i = 0; i < n; ++i)
  {
    vslice_x[i] = voff_x + i;
  } 
  for (i = 0; i < k; ++i)
  {
    MOSEKCALL(res, MSK_putafefrow(task, aoff_q + i + 1, n, vslice_x, GT[i]));
  }
  free(vslice_x);
  MSKint64t qdom;
  MOSEKCALL(res, MSK_appendquadraticconedomain(task, k + 1, &qdom));
  MOSEKCALL(res, MSK_appendaccseq(task, qdom, k + 1, aoff_q, NULL));
  sprintf(buf, "%s", "risk");            
  MOSEKCALL(res, MSK_putaccname(task, aoff_q, buf));

  // Objective: maximize expected return mu^T x
  for (j = 0; j < n; ++j)
  {
    MOSEKCALL(res, MSK_putcj(task, voff_x + j, mu[j]));
  }
  MOSEKCALL(res, MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE));

#if 0
  /* no log output. */
  MOSEKCALL(res, MSK_putintparam(task, MSK_IPAR_LOG, 0));
#endif

#if 0
  /* Dump the problem to a human readable OPF file. */
  MOSEKCALL(res, MSK_writedata(task, "dump.ptf"));
#endif

  MOSEKCALL(res, MSK_optimizetrm(task, &trmcode));

#if 1
  /* Display the solution summary for quick inspection of results. */
  MSK_solutionsummary(task, MSK_STREAM_MSG);
#endif
  
  if (res == MSK_RES_OK)
    MOSEKCALL(res, MSK_getxxslice(task, MSK_SOL_ITG, voff_x, voff_x + n, xx));

  MSK_deletetask(&task);
  MSK_deleteenv(&env);

  return res;
}

int main(int argc, const char **argv)
{
  int i;
  double expret;
  MSKrescodee res;

  const MSKint32t n       = 8;
  const MSKrealt  mu[]    = {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379};
  // GT must have size n rows
  const MSKrealt  GT[][8] = {
    {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
    {0.0,     0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
    {0.0,     0.0,     0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
    {0.0,     0.0,     0.0,     0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
    {0.0,     0.0,     0.0,     0.0,     0.36096, 0.12574, 0.10157, 0.0571 },
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.21552, 0.05663, 0.06187},
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.22514, 0.03327},
    {0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.0,     0.2202 }
  };
  const MSKint32t k       = sizeof(GT) / (n * sizeof(MSKrealt));
  const MSKrealt  x0[]    = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  const MSKrealt  w       = 1.0;
  const MSKrealt  gamma   = 0.25;
  double xx[]    = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  int   K, j;

  for (K = 1; K <= n; K++)
  {
    res = markowitz_with_card(n, k, x0, w, gamma, mu, GT, K, xx);
    if (res == MSK_RES_OK)
    {
      expret = 0;
      printf("Bound %d:  x = ", K);
      for(i=0; i<n; i++)
      { 
        printf("%.5f ", xx[i]);
        expret += xx[i]*mu[i];
      }
      printf("  Return:  %.5f\n", expret);
    }
  }

  return (0);
}

portfolio_6_factor.c

Listing 17.29 portfolio_6_factor.c Click here to download.
/*
  File : portfolio_6_factor.c

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a portfolio optimization model using factor model.
*/
#include <math.h>
#include <stdio.h>
#include "mosek.h"

#define MOSEKCALL(_r,_call)  if ( (_r)==MSK_RES_OK ) (_r) = (_call)

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

typedef struct 
{
  MSKrealt** m;
  int nr;
  int nc;
} matrix;

static void matrix_print(matrix* m)
{
  int i,j;
  for (i = 0; i < m->nr; ++i)
  {
    printf("[");
    for (j = 0; j < m->nc; ++j)
    {
      printf("%f, ", m->m[i][j]);    
    }
    printf("\b\b]\n");
  }
}

static void array_print(MSKrealt* a, int len)
{
  int j;
  printf("[");
  for (j = 0; j < len; ++j)
  {
    printf("%f, ", a[j]);    
  }
  printf("\b\b]\n");
}

static matrix* matrix_alloc(int dim1, int dim2)
{
  int i;
  matrix* m = (matrix*) malloc(sizeof(matrix));
  m->nr = dim1;
  m->nc = dim2;
  m->m = (MSKrealt**) malloc(dim1 * sizeof(MSKrealt*));
  for (i = 0; i < dim1; ++i)
  {
    m->m[i] = (MSKrealt*) malloc(dim2 * sizeof(MSKrealt));
  }
  return m;
}

static void matrix_free(matrix* m)
{
  int i;
  for (i = 0; i < m->nr; ++i)
  {
    free(m->m[i]);
  }
  free(m->m);
  free(m);
}

static MSKrealt* vector_alloc(int dim)
{
  MSKrealt* v = (MSKrealt*) malloc(dim * sizeof(MSKrealt));
  return v;
}

static void vector_free(MSKrealt* v)
{
  free(v);
}

static MSKrealt sum(MSKrealt* x, int n) 
{
  int i;
  MSKrealt r = 0.0;
  for (i = 0; i < n; ++i) r += x[i];
  return r;
}

static MSKrealt dot(MSKrealt* x, MSKrealt* y, int n) 
{
  int i;
  MSKrealt r = 0.0;
  for (i = 0; i < n; ++i) r += x[i] * y[i];
  return r;
}

// Vectorize matrix (column-major order)
static MSKrealt* mat_to_vec_c(matrix* m) 
{  
  int ni = m->nr;
  int nj = m->nc;
  int i,j;

  MSKrealt* c = vector_alloc(ni * nj);
  for (j = 0; j < nj; ++j) 
  {
    for (i = 0; i < ni; ++i) 
    {
      c[j * ni + i] = m->m[i][j];
    }
  }
  return c;
}

// Reshape vector to matrix (column-major order)
static matrix* vec_to_mat_c(MSKrealt* c, int ni, int nj) 
{  
  int i,j;
  matrix* m = matrix_alloc(ni, nj);
  for (j = 0; j < nj; ++j) 
  {
    for (i = 0; i < ni; ++i) 
    {
      m->m[i][j] = c[j * ni + i];
    }
  }
  return m;
}

// Reshape vector to matrix (row-major order)
static matrix* vec_to_mat_r(MSKrealt* r, int ni, int nj) 
{  
  int i,j;
  matrix* m = matrix_alloc(ni, nj);
  for (i = 0; i < ni; ++i) 
  {
    for (j = 0; j < nj; ++j) 
    {
      m->m[i][j] = r[i * nj + j];
    }
  }
  return m;
}

static matrix* cholesky(MSKenv_t* env, matrix* m) 
{
  MSKrescodee res = MSK_RES_OK;
  int i,j;
  int n = m->nr;  
  MSKrealt* vecs = mat_to_vec_c(m);
  MOSEKCALL(res, MSK_potrf(*env, MSK_UPLO_LO, n, vecs));
  matrix* s = vec_to_mat_c(vecs, n, n);
  vector_free(vecs);

  // Zero out upper triangular part (MSK_potrf does not use it, original matrix values remain there)
  for (i = 0; i < n; ++i) 
  {
    for (j = i+1; j < n; ++j) 
    {
      s->m[i][j] = 0.0;
    }
  }
  return s;
}

// Matrix multiplication
static matrix* matrix_mul(MSKenv_t* env, matrix* a, matrix* b) 
{
  MSKrescodee res = MSK_RES_OK;
  int na = a->nr;
  int nb = b->nc;
  int k = b->nr;
  int i;

  MSKrealt* vecm = vector_alloc(na * nb);
  for (i = 0; i < na * nb; ++i) vecm[i] = 0.0;
  MSKrealt* veca = mat_to_vec_c(a);
  MSKrealt* vecb = mat_to_vec_c(b);
  MOSEKCALL(res, MSK_gemm(*env, MSK_TRANSPOSE_NO, MSK_TRANSPOSE_NO, na, nb, k, 1.0, veca, vecb, 1.0, vecm));
  vector_free(veca);
  vector_free(vecb);
  matrix* ab = vec_to_mat_c(vecm, na, nb);
  vector_free(vecm);
  return ab;
}

int main(int argc, const char **argv)
{
  char            buf[128];

  MSKrealt        expret  = 0.0,
                  xj;

  MSKrealt        rtemp;
  MSKenv_t        env;
  MSKint32t       i, j, *sub;
  MSKrescodee     res = MSK_RES_OK, trmcode;
  MSKtask_t       task;

  MSKint32t  n      = 8;
  MSKrealt   w      = 1.0;
  MSKrealt   mu[]   = {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379};
  MSKrealt   x0[]   = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

  /* Initial setup. */
  env  = NULL;
  task = NULL;
  MOSEKCALL(res, MSK_makeenv(&env, NULL));
  MOSEKCALL(res, MSK_maketask(env, 0, 0, &task));
  MOSEKCALL(res, MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr));

  // NOTE: Here we specify matrices as vectors (row major order) to avoid having 
  // to initialize them as double(*)[] type, which is incompatible with double**. 
  
  // Factor exposure matrix
  MSKrealt vecB[] = 
  {
    0.4256, 0.1869,
    0.2413, 0.3877,
    0.2235, 0.3697,
    0.1503, 0.4612,
    1.5325, -0.2633,
    1.2741, -0.2613,
    0.6939, 0.2372,
    0.5425, 0.2116
  };
  matrix* B = vec_to_mat_r(vecB, n, 2); 
  // Factor covariance matrix
  MSKrealt vecS_F[] = 
  {
    0.0620, 0.0577,
    0.0577, 0.0908
  };
  matrix* S_F = vec_to_mat_r(vecS_F, 2, 2);

  // Specific risk components
  MSKrealt theta[] = {0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459};

  matrix* P = cholesky(&env, S_F);
  matrix* G_factor = matrix_mul(&env, B, P);
  matrix_free(B);
  matrix_free(S_F);
  matrix_free(P);

  MSKint32t k = G_factor->nc;
  MSKrealt gammas[] = {0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48};
  int num_gammas = 7;
  MSKrealt   totalBudget;

  //Offset of variables into the API variable.
  MSKint32t numvar = n;
  MSKint32t voff_x = 0;

  // Constraint offset
  MSKint32t coff_bud = 0;

  // Holding variable x of length n
  // No other auxiliary variables are needed in this formulation
  MOSEKCALL(res, MSK_appendvars(task, numvar));
  // Setting up variable x 
  for (j = 0; j < n; ++j)
  {
    /* Optionally we can give the variables names */
    sprintf(buf, "x[%d]", 1 + j);
    MOSEKCALL(res, MSK_putvarname(task, voff_x + j, buf));
    /* No short-selling - x^l = 0, x^u = inf */
    MOSEKCALL(res, MSK_putvarbound(task, voff_x + j, MSK_BK_LO, 0.0, MSK_INFINITY));
  }

  // One linear constraint: total budget
  MOSEKCALL(res, MSK_appendcons(task, 1));
  sprintf(buf, "%s", "budget");
  MOSEKCALL(res, MSK_putconname(task, coff_bud, buf));
  for (j = 0; j < n; ++j)
  {
    /* Coefficients in the first row of A */
    MOSEKCALL(res, MSK_putaij(task, coff_bud, voff_x + j, 1.0));
  }
  totalBudget = w;
  for (i = 0; i < n; ++i)
  {
    totalBudget += x0[i];
  }
  MOSEKCALL(res, MSK_putconbound(task, coff_bud, MSK_BK_FX, totalBudget, totalBudget));

  // Input (gamma, G_factor_T x, diag(sqrt(theta))*x) in the AFE (affine expression) storage
  // We need k+n+1 rows and we fill them in in three parts
  MOSEKCALL(res, MSK_appendafes(task, k + n + 1));
  // 1. The first affine expression = gamma, will be specified later
  // 2. The next k expressions comprise G_factor_T*x, we add them column by column since
  //    G_factor is stored row-wise and we transpose on the fly 
  MSKint64t* afeidxs = (MSKint64t*) malloc(k * sizeof(MSKint64t));
  for (i = 0; i < k; ++i)
  {
    afeidxs[i] = i + 1;
  } 
  for (i = 0; i < n; ++i)
  {
    MOSEKCALL(res, MSK_putafefcol(task, i, k, afeidxs, G_factor->m[i])); //i-th row of G_factor goes in i-th column of F
  }
  free(afeidxs);
  // 3. The remaining n rows contain sqrt(theta) on the diagonal
  for (i = 0; i < n; ++i)
  {
    MOSEKCALL(res, MSK_putafefentry(task, k + 1 + i, voff_x + i, sqrt(theta[i])));
  }

  // Input the affine conic constraint (gamma, G_factor_T x, diag(sqrt(theta))*x) \in QCone
  // Add the quadratic domain of dimension k+n+1
  MSKint64t qdom;
  MOSEKCALL(res, MSK_appendquadraticconedomain(task, k + n + 1, &qdom));
  // Add the constraint
  MOSEKCALL(res, MSK_appendaccseq(task, qdom, k + n + 1, 0, NULL));
  sprintf(buf, "%s", "risk");            
  MOSEKCALL(res, MSK_putaccname(task, 0, buf));

  // Objective: maximize expected return mu^T x
  for (j = 0; j < n; ++j)
  {
    MOSEKCALL(res, MSK_putcj(task, voff_x + j, mu[j]));
  }
  MOSEKCALL(res, MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE));

#if 0
  /* No log output */
  MOSEKCALL(res, MSK_putintparam(task, MSK_IPAR_LOG, 0));
#endif

  for (i = 0; i < num_gammas; i++)
  { 
    MSKrealt gamma = gammas[i];

    // Specify gamma in ACC
    MOSEKCALL(res, MSK_putafeg(task, 0, gamma));

#if 0
    /* Dump the problem to a human readable PTF file. */
    MOSEKCALL(res, MSK_writedata(task, "dump.ptf"));
#endif

    MOSEKCALL(res, MSK_optimizetrm(task, &trmcode));

#if 1
    /* Display the solution summary for quick inspection of results. */
    MSK_solutionsummary(task, MSK_STREAM_MSG);
#endif

    if ( res == MSK_RES_OK )
    {
      expret = 0.0;

      /* Read the x variables one by one and compute expected return. */
      /* Can also be obtained as value of the objective. */
      for (j = 0; j < n; ++j)
      {
        MOSEKCALL(res, MSK_getxxslice(task, MSK_SOL_ITR, voff_x + j, voff_x + j + 1, &xj));
        expret += mu[j] * xj;
      }

      printf("\nExpected return %e for gamma %e\n", expret, gamma);
    }
  }

  if ( task != NULL )
    MSK_deletetask(&task);

  if ( env != NULL )
    MSK_deleteenv(&env);

  matrix_free(G_factor);

  return ( 0 );
}

pow1.c

Listing 17.30 pow1.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      pow1.c

  Purpose: Demonstrates how to solve the problem

    maximize x^0.2*y^0.8 + z^0.4 - x
          st x + y + 0.5z = 2
             x,y,z >= 0
 */
#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  MSKrescodee  r;

  const MSKint32t numvar = 5,
                  numcon = 1;
  const MSKint64t numafe = 6,
                  numacc = 2,
                  f_nnz  = 5;

  MSKboundkeye bkx[5];
  double       blx[5], bux[5];

  double       val[] = { 1.0, 1.0, -1.0 };
  MSKint32t    sub[] = { 3, 4, 0 };   

  double       aval[] = { 1.0, 1.0, 0.5 };
  MSKint32t    asub[] = { 0, 1, 2 };  

  MSKint64t  afeidx[] = {0, 1, 2, 3, 5};
  MSKint32t  varidx[] = {0, 1, 3, 2, 4};
  MSKrealt    f_val[] = {1, 1, 1, 1, 1},
                    g = 1;

  const MSKrealt alpha_1[]= {0.2, 0.8}, 
                 alpha_2[]= {0.4, 0.6};

  MSKint64t   domidx[] = {0, 0}; 

  MSKint32t   i, j;

  MSKenv_t    env  = NULL;
  MSKtask_t   task = NULL;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, numcon, numvar, &task);

    if (r == MSK_RES_OK)
    {
      MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      /* Append 'numcon' empty constraints.
      The constraints will initially have no bounds. */
      if (r == MSK_RES_OK)
        r = MSK_appendcons(task, numcon);

      /* Append 'numvar' variables.
      The variables will initially be fixed at zero (x=0). */
      if (r == MSK_RES_OK)
        r = MSK_appendvars(task, numvar);

      /* Append 'numafe' affine expressions.
      The affine expressions will initially be empty. */
      if (r == MSK_RES_OK)
        r = MSK_appendafes(task, numafe);


      /* Set up the linear part */
      MSK_putclist(task, 3, sub, val);
      MSK_putarow(task, 0, 3, asub, aval);
      MSK_putconbound(task, 0, MSK_BK_FX, 2.0, 2.0);
      for(i=0;i<5;i++)
        bkx[i] = MSK_BK_FR, blx[i] = -MSK_INFINITY, bux[i] = MSK_INFINITY;
      MSK_putvarboundslice(task, 0, numvar, bkx, blx, bux);

      if (r == MSK_RES_OK)
      {
        /* Set the non-zero entries of the F matrix */
        r = MSK_putafefentrylist(task, f_nnz, afeidx, varidx, f_val);
        /* Set the non-zero element of the g vector */
        if (r == MSK_RES_OK)
          r = MSK_putafeg(task, 4, g);

        /* Append the primal power cone domains with their respective parameter values. */
        if (r == MSK_RES_OK)
          r = MSK_appendprimalpowerconedomain(task, 3, 2, alpha_1, domidx);
        if (r == MSK_RES_OK)
          r = MSK_appendprimalpowerconedomain(task, 3, 2, alpha_2, domidx+1);

        /* Append two ACCs made up of the AFEs and the domains defined above. */
        if (r == MSK_RES_OK)
          r = MSK_appendaccsseq(task, numacc, domidx, numafe, afeidx[0], NULL);
      }

      MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE);

      if (r == MSK_RES_OK)
      {
        MSKrescodee trmcode;

        /* Run optimizer */
        r = MSK_optimizetrm(task, &trmcode);

        /* Print a summary containing information
           about the solution for debugging purposes*/
        MSK_solutionsummary(task, MSK_STREAM_MSG);

        if (r == MSK_RES_OK)
        {
          MSKsolstae solsta;

          MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

          switch (solsta)
          {
            case MSK_SOL_STA_OPTIMAL:
              {
                double *xx = NULL;

                xx = calloc(numvar, sizeof(double));
                if (xx)
                {
                  MSK_getxx(task,
                            MSK_SOL_ITR,    /* Request the interior solution. */
                            xx);

                  printf("Optimal primal solution\n");
                  for (j = 0; j < 3; ++j)
                    printf("x[%d]: %e\n", j, xx[j]);
                }
                else
                {
                  r = MSK_RES_ERR_SPACE;
                }
                free(xx);
              }
              break;
            case MSK_SOL_STA_DUAL_INFEAS_CER:
            case MSK_SOL_STA_PRIM_INFEAS_CER:
              printf("Primal or dual infeasibility certificate found.\n");
              break;
            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined. Termination code: %d.\n", trmcode);
              break;
            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          printf("Error while optimizing.\n");
        }
      }

      if (r != MSK_RES_OK)
      {
        /* In case of an error print error code and description. */
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];

        printf("An error occurred while optimizing.\n");
        MSK_getcodedesc(r,
                        symname,
                        desc);
        printf("Error %s - '%s'\n", symname, desc);
      }
    }
    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return (r);
} /* main */

qcqo1.c

Listing 17.31 qcqo1.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      qcqo1.c

   Purpose:   To demonstrate how to solve a quadratic
              optimization problem using the MOSEK API.

              minimize  x_1^2 + 0.1 x_2^2 +  x_3^2 - x_1 x_3 - x_2
              s.t 1 <=  x_1 + x_2 + x_3 - x_1^2 - x_2^2 - 0.1 x_3^2 + 0.2 x_1 x_3
              x >= 0

 */
#include <stdio.h>

#include "mosek.h" /* Include the MOSEK definition file. */

#define NUMCON 1   /* Number of constraints.             */
#define NUMVAR 3   /* Number of variables.               */
#define NUMANZ 3   /* Number of non-zeros in A.          */
#define NUMQNZ 4   /* Number of non-zeros in Q.          */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{

  MSKrescodee  r;

  double       c[]    = {0.0, -1.0, 0.0};

  MSKboundkeye bkc[]  = {MSK_BK_LO};
  double       blc[]  = {1.0};
  double       buc[]  = { +MSK_INFINITY};

  MSKboundkeye bkx[]  = {MSK_BK_LO,
                         MSK_BK_LO,
                         MSK_BK_LO
                        };
  double       blx[]  = {0.0,
                         0.0,
                         0.0
                        };
  double       bux[]  = { +MSK_INFINITY,
                          +MSK_INFINITY,
                          +MSK_INFINITY
                        };

  MSKint32t    aptrb[] = {0, 1, 2 },
               aptre[] = {1, 2, 3},
               asub[] = { 0,   0,   0};
  double       aval[] = { 1.0, 1.0, 1.0};
  MSKint32t    qsubi[NUMQNZ],
               qsubj[NUMQNZ];
  double       qval[NUMQNZ];

  MSKint32t    j, i;
  double       xx[NUMVAR];
  MSKenv_t     env;
  MSKtask_t    task;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, NUMCON, NUMVAR, &task);

    if (r == MSK_RES_OK)
    {
      r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      /* Append 'NUMCON' empty constraints.
       The constraints will initially have no bounds. */
      if (r == MSK_RES_OK)
        r = MSK_appendcons(task, NUMCON);

      /* Append 'NUMVAR' variables.
       The variables will initially be fixed at zero (x=0). */
      if (r == MSK_RES_OK)
        r = MSK_appendvars(task, NUMVAR);

      /* Optionally add a constant term to the objective. */
      if (r == MSK_RES_OK)
        r = MSK_putcfix(task, 0.0);
      for (j = 0; j < NUMVAR && r == MSK_RES_OK; ++j)
      {
        /* Set the linear term c_j in the objective.*/
        if (r == MSK_RES_OK)
          r = MSK_putcj(task, j, c[j]);

        /* Set the bounds on variable j.
         blx[j] <= x_j <= bux[j] */
        if (r == MSK_RES_OK)
          r = MSK_putvarbound(task,
                              j,           /* Index of variable.*/
                              bkx[j],      /* Bound key.*/
                              blx[j],      /* Numerical value of lower bound.*/
                              bux[j]);     /* Numerical value of upper bound.*/

        /* Input column j of A */
        if (r == MSK_RES_OK)
          r = MSK_putacol(task,
                          j,                 /* Variable (column) index.*/
                          aptre[j] - aptrb[j], /* Number of non-zeros in column j.*/
                          asub + aptrb[j],   /* Pointer to row indexes of column j.*/
                          aval + aptrb[j]);  /* Pointer to Values of column j.*/

      }

      /* Set the bounds on constraints.
         for i=1, ...,NUMCON : blc[i] <= constraint i <= buc[i] */
      for (i = 0; i < NUMCON && r == MSK_RES_OK; ++i)
        r = MSK_putconbound(task,
                            i,           /* Index of constraint.*/
                            bkc[i],      /* Bound key.*/
                            blc[i],      /* Numerical value of lower bound.*/
                            buc[i]);     /* Numerical value of upper bound.*/

      if (r == MSK_RES_OK)
      {
        /*
         * The lower triangular part of the Q^o
         * matrix in the objective is specified.
         */

        qsubi[0] = 0;   qsubj[0] = 0;  qval[0] = 2.0;
        qsubi[1] = 1;   qsubj[1] = 1;  qval[1] = 0.2;
        qsubi[2] = 2;   qsubj[2] = 0;  qval[2] = -1.0;
        qsubi[3] = 2;   qsubj[3] = 2;  qval[3] = 2.0;

        /* Input the Q^o for the objective. */

        r = MSK_putqobj(task, NUMQNZ, qsubi, qsubj, qval);
      }

      if (r == MSK_RES_OK)
      {
        /*
        * The lower triangular part of the Q^0
        * matrix in the first constraint is specified.
        This corresponds to adding the term
        - x_1^2 - x_2^2 - 0.1 x_3^2 + 0.2 x_1 x_3
        */

        qsubi[0] = 0;   qsubj[0] = 0;  qval[0] = -2.0;
        qsubi[1] = 1;   qsubj[1] = 1;  qval[1] = -2.0;
        qsubi[2] = 2;   qsubj[2] = 2;  qval[2] = -0.2;
        qsubi[3] = 2;   qsubj[3] = 0;  qval[3] = 0.2;

        /* Put Q^0 in constraint with index 0. */

        r = MSK_putqconk(task,
                         0,
                         4,
                         qsubi,
                         qsubj,
                         qval);
      }

      if (r == MSK_RES_OK)
        r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MINIMIZE);

      if (r == MSK_RES_OK)
      {
        MSKrescodee trmcode;

        /* Run optimizer */
        r = MSK_optimizetrm(task, &trmcode);

        /* Print a summary containing information
           about the solution for debugging purposes*/
        MSK_solutionsummary(task, MSK_STREAM_LOG);

        if (r == MSK_RES_OK)
        {
          MSKsolstae solsta;
          int j;

          MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

          switch (solsta)
          {
            case MSK_SOL_STA_OPTIMAL:
              MSK_getxx(task,
                        MSK_SOL_ITR,    /* Request the interior solution. */
                        xx);

              printf("Optimal primal solution\n");
              for (j = 0; j < NUMVAR; ++j)
                printf("x[%d]: %e\n", j, xx[j]);

              break;

            case MSK_SOL_STA_DUAL_INFEAS_CER:
            case MSK_SOL_STA_PRIM_INFEAS_CER:
              printf("Primal or dual infeasibility certificate found.\n");
              break;

            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined. Termination code: %d.\n", trmcode);
              break;

            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          printf("Error while optimizing.\n");
        }
      }

      if (r != MSK_RES_OK)
      {
        /* In case of an error print error code and description. */
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];

        printf("An error occurred while optimizing.\n");
        MSK_getcodedesc(r,
                        symname,
                        desc);
        printf("Error %s - '%s'\n", symname, desc);
      }
    }

    MSK_deletetask(&task);
  }
  MSK_deleteenv(&env);

  return (r);
} /* main */

qo1.c

Listing 17.32 qo1.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      qo1.c

   Purpose: To demonstrate how to solve a quadratic optimization
              problem using the MOSEK API.
 */

#include <stdio.h>

#include "mosek.h" /* Include the MOSEK definition file. */

#define NUMCON 1   /* Number of constraints.             */
#define NUMVAR 3   /* Number of variables.               */
#define NUMANZ 3   /* Number of non-zeros in A.           */
#define NUMQNZ 4   /* Number of non-zeros in Q.           */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{

  double        c[]   = {0.0, -1.0, 0.0};

  MSKboundkeye  bkc[] = {MSK_BK_LO};
  double        blc[] = {1.0};
  double        buc[] = { +MSK_INFINITY };

  MSKboundkeye  bkx[] = {MSK_BK_LO,
                         MSK_BK_LO,
                         MSK_BK_LO
                        };
  double        blx[] = {0.0,
                         0.0,
                         0.0
                        };
  double        bux[] = { +MSK_INFINITY,
                          +MSK_INFINITY,
                          +MSK_INFINITY
                        };

  MSKint32t     aptrb[] = {0,   1,   2},
                aptre[] = {1,   2,   3},
                asub[]  = {0,   0,   0};
  double        aval[]  = {1.0, 1.0, 1.0};

  MSKint32t     qsubi[NUMQNZ];
  MSKint32t     qsubj[NUMQNZ];
  double        qval[NUMQNZ];

  MSKint32t     i, j;
  double        xx[NUMVAR];

  MSKenv_t      env = NULL;
  MSKtask_t     task = NULL;
  MSKrescodee   r;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, NUMCON, NUMVAR, &task);

    if (r == MSK_RES_OK)
    {
      r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      /* Append 'NUMCON' empty constraints.
       The constraints will initially have no bounds. */
      if (r == MSK_RES_OK)
        r = MSK_appendcons(task, NUMCON);

      /* Append 'NUMVAR' variables.
       The variables will initially be fixed at zero (x=0). */
      if (r == MSK_RES_OK)
        r = MSK_appendvars(task, NUMVAR);

      /* Optionally add a constant term to the objective. */
      if (r == MSK_RES_OK)
        r = MSK_putcfix(task, 0.0);
      for (j = 0; j < NUMVAR && r == MSK_RES_OK; ++j)
      {
        /* Set the linear term c_j in the objective.*/
        if (r == MSK_RES_OK)
          r = MSK_putcj(task, j, c[j]);

        /* Set the bounds on variable j.
         blx[j] <= x_j <= bux[j] */
        if (r == MSK_RES_OK)
          r = MSK_putvarbound(task,
                              j,           /* Index of variable.*/
                              bkx[j],      /* Bound key.*/
                              blx[j],      /* Numerical value of lower bound.*/
                              bux[j]);     /* Numerical value of upper bound.*/

        /* Input column j of A */
        if (r == MSK_RES_OK)
          r = MSK_putacol(task,
                          j,                 /* Variable (column) index.*/
                          aptre[j] - aptrb[j], /* Number of non-zeros in column j.*/
                          asub + aptrb[j],   /* Pointer to row indexes of column j.*/
                          aval + aptrb[j]);  /* Pointer to Values of column j.*/

      }

      /* Set the bounds on constraints.
         for i=1, ...,NUMCON : blc[i] <= constraint i <= buc[i] */
      for (i = 0; i < NUMCON && r == MSK_RES_OK; ++i)
        r = MSK_putconbound(task,
                            i,           /* Index of constraint.*/
                            bkc[i],      /* Bound key.*/
                            blc[i],      /* Numerical value of lower bound.*/
                            buc[i]);     /* Numerical value of upper bound.*/

      if (r == MSK_RES_OK)
      {
        /*
         * The lower triangular part of the Q
         * matrix in the objective is specified.
         */

        qsubi[0] = 0;   qsubj[0] = 0;  qval[0] = 2.0;
        qsubi[1] = 1;   qsubj[1] = 1;  qval[1] = 0.2;
        qsubi[2] = 2;   qsubj[2] = 0;  qval[2] = -1.0;
        qsubi[3] = 2;   qsubj[3] = 2;  qval[3] = 2.0;

        /* Input the Q for the objective. */

        r = MSK_putqobj(task, NUMQNZ, qsubi, qsubj, qval);
      }

      if (r == MSK_RES_OK)
      {
        MSKrescodee trmcode;

        /* Run optimizer */
        r = MSK_optimizetrm(task, &trmcode);

        /* Print a summary containing information
           about the solution for debugging purposes*/
        MSK_solutionsummary(task, MSK_STREAM_MSG);

        if (r == MSK_RES_OK)
        {
          MSKsolstae solsta;
          int j;

          MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

          switch (solsta)
          {
            case MSK_SOL_STA_OPTIMAL:
              MSK_getxx(task,
                        MSK_SOL_ITR,    /* Request the interior solution. */
                        xx);

              printf("Optimal primal solution\n");
              for (j = 0; j < NUMVAR; ++j)
                printf("x[%d]: %e\n", j, xx[j]);

              break;

            case MSK_SOL_STA_DUAL_INFEAS_CER:
            case MSK_SOL_STA_PRIM_INFEAS_CER:
              printf("Primal or dual infeasibility certificate found.\n");
              break;

            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined. Termination code: %d.\n", trmcode);
              break;

            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          printf("Error while optimizing.\n");
        }
      }

      if (r != MSK_RES_OK)
      {
        /* In case of an error print error code and description. */
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];

        printf("An error occurred while optimizing.\n");
        MSK_getcodedesc(r,
                        symname,
                        desc);
        printf("Error %s - '%s'\n", symname, desc);
      }
    }
    MSK_deletetask(&task);
  }
  MSK_deleteenv(&env);

  return (r);
} /* main */

reoptimization.c

Listing 17.33 reoptimization.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      reoptimization.c

   Purpose:   To demonstrate how to solve a  linear
              optimization problem using the MOSEK API
              and modify and re-optimize the problem.
 */
#include <stdio.h>

#define DEBUG 1    /* Set to zero to turn debugging off. */

#include "mosek.h" /* Include the MOSEK definition file. */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

void printres(int n, const double x[]) {
  int i;
  for (i = 0; i < n; i++) printf("%.3f ", x[i]);
  printf("\n");
}

int main(int argc, char *argv[])
{
  MSKint32t       numvar = 3,
                  numcon = 3;
  MSKint32t       i, j;
  double          c[]    = {1.5, 2.5, 3.0};
  MSKint32t       ptrb[] = {0, 3, 6},
                  ptre[] = {3, 6, 9},
                  asub[] = { 0, 1, 2,
                             0, 1, 2,
                             0, 1, 2
                           };

  double          aval[] = { 2.0, 3.0, 2.0,
                             4.0, 2.0, 3.0,
                             3.0, 3.0, 2.0
                           };

  MSKboundkeye    bkc[]  = {MSK_BK_UP, MSK_BK_UP, MSK_BK_UP    };
  double          blc[]  = { -MSK_INFINITY, -MSK_INFINITY, -MSK_INFINITY};
  double          buc[]  = {100000, 50000, 60000};

  MSKboundkeye    bkx[]  = {MSK_BK_LO,     MSK_BK_LO,    MSK_BK_LO};
  double          blx[]  = {0.0,           0.0,          0.0,};
  double          bux[]  = { +MSK_INFINITY, +MSK_INFINITY, +MSK_INFINITY};

  double          *xx = NULL;
  MSKenv_t        env;
  MSKtask_t       task;
  MSKint32t       varidx, conidx;
  MSKrescodee     r,lr;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env,DEBUG ? "" :  NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, numcon, numvar, &task);

    /* Append the constraints. */
    if (r == MSK_RES_OK)
      r = MSK_appendcons(task, numcon);

    /* Append the variables. */
    if (r == MSK_RES_OK)
      r = MSK_appendvars(task, numvar);

    /* Put C. */
    if (r == MSK_RES_OK)
      r = MSK_putcfix(task, 0.0);

    if (r == MSK_RES_OK)
      for (j = 0; j < numvar; ++j)
        r = MSK_putcj(task, j, c[j]);

    /* Put constraint bounds. */
    if (r == MSK_RES_OK)
      for (i = 0; i < numcon; ++i)
        r = MSK_putconbound(task, i, bkc[i], blc[i], buc[i]);

    /* Put variable bounds. */
    if (r == MSK_RES_OK)
      for (j = 0; j < numvar; ++j)
        r = MSK_putvarbound(task, j, bkx[j], blx[j], bux[j]);

    /* Put A. */
    if (r == MSK_RES_OK)
      if (numcon > 0)
        for (j = 0; j < numvar; ++j)
          r = MSK_putacol(task,
                          j,
                          ptre[j] - ptrb[j],
                          asub + ptrb[j],
                          aval + ptrb[j]);

    if (r == MSK_RES_OK)
      r = MSK_putobjsense(task,
                          MSK_OBJECTIVE_SENSE_MAXIMIZE);

    if (r == MSK_RES_OK)
      r = MSK_optimizetrm(task, NULL);

    if (r == MSK_RES_OK)
    {
      xx = calloc(numvar, sizeof(double));
      if (!xx)
        r = MSK_RES_ERR_SPACE;
    }

    if (r == MSK_RES_OK)
      r = MSK_getxx(task,
                    MSK_SOL_BAS,       /* Basic solution.       */
                    xx);

    printres(numvar, xx);

    /******************** Make a change to the A matrix **********/
    if (r == MSK_RES_OK)
      r = MSK_putaij(task, 0, 0, 3.0);
    if (r == MSK_RES_OK)
      r = MSK_optimizetrm(task, NULL);
    if (r == MSK_RES_OK)
      r = MSK_getxx(task,
                    MSK_SOL_BAS,       /* Basic solution.       */
                    xx);
    printres(numvar, xx);
    if (xx)
      free(xx);

    /*********************** Add a new variable ******************/
    /* Get index of new variable, this should be 3 */
    if (r == MSK_RES_OK)
      r = MSK_getnumvar(task, &varidx);
    /* Append a new variable x_3 to the problem */
    if (r == MSK_RES_OK)
    {
      r = MSK_appendvars(task, 1);
      numvar++;
    }
    /* Set bounds on new variable */
    if (r == MSK_RES_OK)
      r = MSK_putvarbound(task,
                          varidx,
                          MSK_BK_LO,
                          0,
                          +MSK_INFINITY);

    /* Change objective */
    if (r == MSK_RES_OK)
      r = MSK_putcj(task, varidx, 1.0);

    /* Put new values in the A matrix */
    if (r == MSK_RES_OK)
    {
      MSKint32t acolsub[] = {0,   2};
      double    acolval[] =  {4.0, 1.0};

      r = MSK_putacol(task,
                      varidx, /* column index */
                      2, /* num nz in column*/
                      acolsub,
                      acolval);
    }

    /* Change optimizer to free simplex and reoptimize */
    if (r == MSK_RES_OK)
      r = MSK_putintparam(task, MSK_IPAR_OPTIMIZER, MSK_OPTIMIZER_FREE_SIMPLEX);

    if (r == MSK_RES_OK)
      r = MSK_optimizetrm(task, NULL);
    if (r == MSK_RES_OK)
    {
      xx = calloc(numvar, sizeof(double));
      if (!xx)
        r = MSK_RES_ERR_SPACE;
    }

    if (r == MSK_RES_OK)
      r = MSK_getxx(task,
                    MSK_SOL_BAS,       /* Basic solution.       */
                    xx);

    printres(numvar, xx);

    /* **************** Add a new constraint ******************* */
    /* Get index of new constraint*/
    if (r == MSK_RES_OK)
      r = MSK_getnumcon(task, &conidx);

    /* Append a new constraint */
    if (r == MSK_RES_OK)
    {
      r = MSK_appendcons(task, 1);
      numcon++;
    }

    /* Set bounds on new constraint */
    if (r == MSK_RES_OK)
      r = MSK_putconbound(task,
                          conidx,
                          MSK_BK_UP,
                          -MSK_INFINITY,
                          30000);

    /* Put new values in the A matrix */
    if (r == MSK_RES_OK)
    {
      MSKint32t arowsub[] = {0,   1,   2,   3  };
      double    arowval[] =  {1.0, 2.0, 1.0, 1.0};

      r = MSK_putarow(task,
                      conidx, /* row index */
                      4,      /* num nz in row*/
                      arowsub,
                      arowval);
    }
    if (r == MSK_RES_OK)
      r = MSK_optimizetrm(task, NULL);

    if (r == MSK_RES_OK)
      r = MSK_getxx(task,
                    MSK_SOL_BAS,       /* Basic solution.       */
                    xx);

    printres(numvar, xx);


    /* **************** Change constraint bounds ******************* */
    if (r == MSK_RES_OK)
    {
      MSKboundkeye newbkc[]  = { MSK_BK_UP, MSK_BK_UP, MSK_BK_UP, MSK_BK_UP };
      double       newblc[]  = { -MSK_INFINITY, -MSK_INFINITY, -MSK_INFINITY, -MSK_INFINITY };
      double       newbuc[]  = { 80000, 40000, 50000, 22000 };

      r = MSK_putconboundslice(task, 0, numcon, newbkc, newblc, newbuc);
    }
    
    if (r == MSK_RES_OK)
      r = MSK_optimizetrm(task, NULL);

    if (r == MSK_RES_OK)
      r = MSK_getxx(task,
                    MSK_SOL_BAS,       /* Basic solution.       */
                    xx);

    printres(numvar, xx);

    if (xx)
      free(xx);
    
    lr = MSK_deletetask(&task);
    if (r == MSK_RES_OK)
      r = lr;

  }
  lr = MSK_deleteenv(&env);
  if (r == MSK_RES_OK)
    r = lr; 

  printf("Return code: %d (0 means no error occured.)\n", r);

  return (r);
} /* main */

response.c

Listing 17.34 response.c Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:      response.c

  Purpose:   This example demonstrates proper response handling
             for problems solved with the interior-point optimizers.
*/

#include <stdio.h>
#include "mosek.h"

/* Log handler */
void MSKAPI printlog(void *ptr,
                     const char s[])
{
  printf("%s", s);
}

int main(int argc, char const *argv[])
{
  MSKenv_t    env;
  MSKtask_t   task;
  MSKrescodee r;
  char        symname[MSK_MAX_STR_LEN];
  char        desc[MSK_MAX_STR_LEN];
  int         i, numvar;
  double      *xx = NULL;
  const char  *filename;

  if (argc >= 2) filename = argv[1];
  else           filename = "../data/cqo1.mps";

  // Create the environment
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    // Create the task
    r = MSK_makeemptytask(env, &task);

    // (Optionally) attach the log handler to receive log information
    // if ( r == MSK_RES_OK ) MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printlog);

    // (Optionally) uncomment this line to most likely see solution status Unknown
    // MSK_putintparam(task, MSK_IPAR_INTPNT_MAX_ITERATIONS, 1);

    // In this example we read an optimization problem from a file
    r = MSK_readdata(task, filename);

    if (r == MSK_RES_OK)
    {
      MSKrescodee trmcode;
      MSKsolstae  solsta;

      // Do the optimization, and exit in case of error
      r = MSK_optimizetrm(task, &trmcode);

      if (r != MSK_RES_OK) {
        MSK_getcodedesc(r, symname, desc);
        printf("Error during optimization: %s %s\n", symname, desc);
        exit(r);
      }

      MSK_solutionsummary(task, MSK_STREAM_LOG);

      /* Expected result: The solution status of the interiot-point solution is optimal. */

      if (MSK_RES_OK == MSK_getsolsta(task, MSK_SOL_ITR, &solsta))
      {
        switch (solsta)
        {
          case MSK_SOL_STA_OPTIMAL:
            printf("An optimal interior-point solution is located.\n");

            /* Read and print the variable values in the solution */
            MSK_getnumvar(task, &numvar);
            xx = calloc(numvar, sizeof(double));
            MSK_getxx(task, MSK_SOL_ITR, xx);
            for (i = 0; i < numvar; i++)
              printf("xx[%d] = %.4lf\n", i, xx[i]);
            free(xx);
            break;

          case MSK_SOL_STA_DUAL_INFEAS_CER:
            printf("Dual infeasibility certificate found.\n");
            break;

          case MSK_SOL_STA_PRIM_INFEAS_CER:
            printf("Primal infeasibility certificate found.\n");
            break;

          case MSK_SOL_STA_UNKNOWN:
            /* The solutions status is unknown. The termination code
               indicating why the optimizer terminated prematurely. */
            printf("The solution status is unknown.\n");
            /* No-error cause of termination e.g. an iteration limit is reached.  */
            MSK_getcodedesc(trmcode, symname, desc);
            printf("  Termination code: %s %s\n", symname, desc);
            break;

          default:
            MSK_solstatostr(task, solsta, desc);
            printf("An unexpected solution status %s with code %d is obtained.\n", desc, solsta);
            break;
        }
      }
      else
        printf("Could not obtain the solution status for the requested solution.\n");
    }
    else {
      MSK_getcodedesc(r, symname, desc);
      printf("Optimization was not started because of error %s(%d): %s\n", symname, r, desc);
    }

    MSK_deletetask(&task);
  }

  MSK_deleteenv(&env);
  return r;
}

sdo1.c

Listing 17.35 sdo1.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      sdo1.c

   Purpose:   Solves the following small semidefinite optimization problem
              using the MOSEK API.

     minimize    Tr [2, 1, 0; 1, 2, 1; 0, 1, 2]*X + x0

     subject to  Tr [1, 0, 0; 0, 1, 0; 0, 0, 1]*X + x0           = 1
                 Tr [1, 1, 1; 1, 1, 1; 1, 1, 1]*X      + x1 + x2 = 0.5
                 (x0,x1,x2) \in Q,  X \in PSD
*/

#include <stdio.h>

#include "mosek.h"    /* Include the MOSEK definition file.  */

#define NUMCON    2   /* Number of constraints.              */
#define NUMVAR    3   /* Number of conic quadratic variables */
#define NUMANZ    3   /* Number of non-zeros in A            */
#define NUMAFE    3   /* Number of affine expressions        */
#define NUMFNZ    3   /* Number of non-zeros in F            */
#define NUMBARVAR 1   /* Number of semidefinite variables    */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  MSKrescodee  r;

  MSKint32t    DIMBARVAR[] = {3};         /* Dimension of semidefinite cone */
  MSKint64t    LENBARVAR[] = {3 * (3 + 1) / 2}; /* Number of scalar SD variables  */

  MSKboundkeye bkc[] = { MSK_BK_FX, MSK_BK_FX };
  double       blc[] = { 1.0, 0.5 };
  double       buc[] = { 1.0, 0.5 };

  MSKint32t    barc_i[] = {0, 1, 1, 2, 2},
               barc_j[] = {0, 0, 1, 1, 2};
  double       barc_v[] = {2.0, 1.0, 2.0, 1.0, 2.0};

  MSKint32t    aptrb[]  = {0, 1},
               aptre[]  = {1, 3},
               asub[]   = {0, 1, 2}; /* column subscripts of A */
  double       aval[]   = {1.0, 1.0, 1.0};

  MSKint32t    bara_i[] = {0, 1, 2, 0, 1, 2, 1, 2, 2},
               bara_j[] = {0, 1, 2, 0, 0, 0, 1, 1, 2};
  double       bara_v[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
  MSKint32t    conesub[] = {0, 1, 2};

  MSKint64t    afeidx[] = {0, 1, 2};
  int          varidx[] = {0, 1, 2};
  double       f_val[]  = {1, 1, 1};

  MSKint32t    i, j;
  MSKint64t    idx;
  double       falpha = 1.0;

  MSKrealt     *xx;
  MSKrealt     *barx;
  MSKenv_t     env = NULL;
  MSKtask_t    task = NULL;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, NUMCON, 0, &task);

    if (r == MSK_RES_OK)
    {
      MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      /* Append 'NUMCON' empty constraints.
       The constraints will initially have no bounds. */
      if (r == MSK_RES_OK)
        r = MSK_appendcons(task, NUMCON);

      /* Append 'NUMVAR' variables.
      The variables will initially be fixed at zero (x=0). */
      if (r == MSK_RES_OK)
        r = MSK_appendvars(task, NUMVAR);

      /* Append 'NUMAFE' affine expressions.*/
      if (r == MSK_RES_OK)
        r = MSK_appendafes(task, NUMAFE);

      /* Append 'NUMBARVAR' semidefinite variables. */
      if (r == MSK_RES_OK) {
        r = MSK_appendbarvars(task, NUMBARVAR, DIMBARVAR);
      }

      /* Optionally add a constant term to the objective. */
      if (r == MSK_RES_OK)
        r = MSK_putcfix(task, 0.0);

      /* Set the linear term c_j in the objective.*/
      if (r == MSK_RES_OK)
        r = MSK_putcj(task, 0, 1.0);

      for (j = 0; j < NUMVAR && r == MSK_RES_OK; ++j)
        r = MSK_putvarbound(task,
                            j,
                            MSK_BK_FR,
                            -MSK_INFINITY,
                            MSK_INFINITY);

      /* Set the linear term barc_j in the objective.*/
      if (r == MSK_RES_OK)
        r = MSK_appendsparsesymmat(task,
                                   DIMBARVAR[0],
                                   5,
                                   barc_i,
                                   barc_j,
                                   barc_v,
                                   &idx);

      if (r == MSK_RES_OK)
        r = MSK_putbarcj(task, 0, 1, &idx, &falpha);

      /* Set the bounds on constraints.
        for i=1, ...,NUMCON : blc[i] <= constraint i <= buc[i] */
      for (i = 0; i < NUMCON && r == MSK_RES_OK; ++i)
        r = MSK_putconbound(task,
                            i,           /* Index of constraint.*/
                            bkc[i],      /* Bound key.*/
                            blc[i],      /* Numerical value of lower bound.*/
                            buc[i]);     /* Numerical value of upper bound.*/

      /* Input A row by row */
      for (i = 0; i < NUMCON && r == MSK_RES_OK; ++i)
        r = MSK_putarow(task,
                        i,
                        aptre[i] - aptrb[i],
                        asub     + aptrb[i],
                        aval     + aptrb[i]);

      /* Append the affine conic constraint with quadratic cone */
      if (r == MSK_RES_OK)
      {
        r = MSK_putafefentrylist(task, NUMFNZ, afeidx, varidx, f_val);
        if (r == MSK_RES_OK)
          r = MSK_appendquadraticconedomain(task, 3, NULL);
        if (r == MSK_RES_OK)
          r = MSK_appendacc(task, 0, 3, afeidx, NULL);
      }

      /* Add the first row of barA */
      if (r == MSK_RES_OK)
        r = MSK_appendsparsesymmat(task,
                                   DIMBARVAR[0],
                                   3,
                                   bara_i,
                                   bara_j,
                                   bara_v,
                                   &idx);

      if (r == MSK_RES_OK)
        r = MSK_putbaraij(task, 0, 0, 1, &idx, &falpha);

      /* Add the second row of barA */
      if (r == MSK_RES_OK)
        r = MSK_appendsparsesymmat(task,
                                   DIMBARVAR[0],
                                   6,
                                   bara_i + 3,
                                   bara_j + 3,
                                   bara_v + 3,
                                   &idx);

      if (r == MSK_RES_OK)
        r = MSK_putbaraij(task, 1, 0, 1, &idx, &falpha);

      if (r == MSK_RES_OK)
      {
        MSKrescodee trmcode;

        /* Run optimizer */
        r = MSK_optimizetrm(task, &trmcode);

        /* Print a summary containing information
           about the solution for debugging purposes*/
        MSK_solutionsummary(task, MSK_STREAM_MSG);

        if (r == MSK_RES_OK)
        {
          MSKsolstae solsta;

          MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

          switch (solsta)
          {
            case MSK_SOL_STA_OPTIMAL:
              xx   = (MSKrealt *) MSK_calloctask(task, NUMVAR, sizeof(MSKrealt));
              barx = (MSKrealt *) MSK_calloctask(task, LENBARVAR[0], sizeof(MSKrealt));

              MSK_getxx(task,
                        MSK_SOL_ITR,
                        xx);
              MSK_getbarxj(task,
                           MSK_SOL_ITR,    /* Request the interior solution. */
                           0,
                           barx);

              printf("Optimal primal solution\n");
              for (i = 0; i < NUMVAR; ++i)
                printf("x[%d]   : % e\n", i, xx[i]);

              for (i = 0; i < LENBARVAR[0]; ++i)
                printf("barx[%d]: % e\n", i, barx[i]);

              MSK_freetask(task, xx);
              MSK_freetask(task, barx);

              break;

            case MSK_SOL_STA_DUAL_INFEAS_CER:
            case MSK_SOL_STA_PRIM_INFEAS_CER:
              printf("Primal or dual infeasibility certificate found.\n");
              break;

            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined. Termination code: %d.\n", trmcode);
              break;

            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          printf("Error while optimizing.\n");
        }
      }

      if (r != MSK_RES_OK)
      {
        /* In case of an error print error code and description. */
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];

        printf("An error occurred while optimizing.\n");
        MSK_getcodedesc(r,
                        symname,
                        desc);
        printf("Error %s - '%s'\n", symname, desc);
      }
    }
    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return (r);
} /* main */

sdo2.c

Listing 17.36 sdo2.c Click here to download.
/*
  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File :      sdo2.c

  Purpose :   Solves the semidefinite problem with two symmetric variables:

                 min   <C1,X1> + <C2,X2>
                 st.   <A1,X1> + <A2,X2> = b
                             (X2)_{1,2} <= k
                
                 where X1, X2 are symmetric positive semidefinite,

                 C1, C2, A1, A2 are assumed to be constant symmetric matrices,
                 and b, k are constants.
*/

#include <stdio.h>
#include "mosek.h"    /* Include the MOSEK definition file.  */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  MSKrescodee  r;

  /* Input data */
  MSKint32t    numbarvar = 2;
  MSKint32t    dimbarvar[] = {3, 4};          /* Dimension of semidefinite variables */

  /* Objective coefficients concatenated */
  MSKint32t    Cj[] = { 0, 0, 1, 1, 1, 1 };   /* Which symmetric variable (j) */
  MSKint32t    Ck[] = { 0, 2, 0, 1, 1, 2 };   /* Which entry (k,l)->v */
  MSKint32t    Cl[] = { 0, 2, 0, 0, 1, 2 };
  MSKrealt     Cv[] = { 1.0, 6.0, 1.0, -3.0, 2.0, 1.0 };

  /* Equality constraints coefficients concatenated */
  MSKint32t    Ai[] = { 0, 0, 0, 0, 0, 0 };   /* Which constraint (i = 0) */
  MSKint32t    Aj[] = { 0, 0, 0, 1, 1, 1 };   /* Which symmetric variable (j) */
  MSKint32t    Ak[] = { 0, 2, 2, 1, 1, 3 };   /* Which entry (k,l)->v */
  MSKint32t    Al[] = { 0, 0, 2, 0, 1, 3 };
  MSKrealt     Av[] = { 1.0, 1.0, 2.0, 1.0, -1.0, -3.0 };

  /* The second constraint - one-term inequality */
  MSKint32t    A2i = 1;                        /* Which constraint (i = 1) */
  MSKint32t    A2j = 1;                        /* Which symmetric variable (j = 1) */
  MSKint32t    A2k = 1;                        /* Which entry A(1,0) = A(0,1) = 0.5 */
  MSKint32t    A2l = 0;
  MSKrealt     A2v = 0.5;

  /* Constraint bounds and values */
  MSKint32t    numcon = 2;
  MSKboundkeye bkc[] = { MSK_BK_FX, MSK_BK_UP };
  double       blc[] = { 23.0, -MSK_INFINITY };
  double       buc[] = { 23.0, -3.0 };

  MSKint32t    i, j, dim;
  MSKrealt     *barx;
  MSKenv_t     env = NULL;
  MSKtask_t    task = NULL;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, 0, 0, &task);

    if (r == MSK_RES_OK)
    {
      MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      /* Append empty constraints.
       The constraints will initially have no bounds. */
      if (r == MSK_RES_OK)
        r = MSK_appendcons(task, numcon);

      /* Append semidefinite variables. */
      if (r == MSK_RES_OK)
        r = MSK_appendbarvars(task, numbarvar, dimbarvar);

      /* Set objective (6 nonzeros).*/
      if (r == MSK_RES_OK)
        r = MSK_putbarcblocktriplet(task, 6, Cj, Ck, Cl, Cv);

      /* Set the equality constraint (6 nonzeros).*/
      if (r == MSK_RES_OK)
        r = MSK_putbarablocktriplet(task, 6, Ai, Aj, Ak, Al, Av);

      /* Set the inequality constraint (1 nonzero).*/
      if (r == MSK_RES_OK)
        r = MSK_putbarablocktriplet(task, 1, &A2i, &A2j, &A2k, &A2l, &A2v);

      /* Set constraint bounds */
      if (r == MSK_RES_OK)
        r = MSK_putconboundslice(task, 0, 2, bkc, blc, buc);

      if (r == MSK_RES_OK)
      {
        MSKrescodee trmcode;

        /* Run optimizer */
        r = MSK_optimizetrm(task, &trmcode);
        MSK_solutionsummary(task, MSK_STREAM_MSG);

        if (r == MSK_RES_OK)
        {
          MSKsolstae solsta;

          MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

          switch (solsta)
          {
            case MSK_SOL_STA_OPTIMAL:

              /* Retrieve the soution for all symmetric variables */
              printf("Solution (lower triangular part vectorized):\n");
              for(i = 0; i < numbarvar; i++) {
                dim = dimbarvar[i] * (dimbarvar[i] + 1) / 2;
                barx = (MSKrealt *) MSK_calloctask(task, dim, sizeof(MSKrealt));

                MSK_getbarxj(task, MSK_SOL_ITR, i, barx);
  
                printf("X%d: ", i + 1);
                for (j = 0; j < dim; ++j)
                  printf("%.3f ", barx[j]);
                printf("\n");

                MSK_freetask(task, barx);                
              }

              break;

            case MSK_SOL_STA_DUAL_INFEAS_CER:
            case MSK_SOL_STA_PRIM_INFEAS_CER:
              printf("Primal or dual infeasibility certificate found.\n");
              break;

            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined. Termination code: %d.\n", trmcode);
              break;

            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          printf("Error while optimizing.\n");
        }
      }

      if (r != MSK_RES_OK)
      {
        /* In case of an error print error code and description. */
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];

        printf("An error occurred while optimizing.\n");
        MSK_getcodedesc(r,
                        symname,
                        desc);
        printf("Error %s - '%s'\n", symname, desc);
      }
    }
    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return (r);
} /* main */

sdo_lmi.c

Listing 17.37 sdo_lmi.c Click here to download.
/*
    Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
 
    File :      sdo_lmi.c
 
    Purpose :   To solve a problem with an LMI and an affine conic constrained problem with a PSD term
 
                 minimize    Tr [1, 0; 0, 1]*X + x(1) + x(2) + 1

                 subject to  Tr [0, 1; 1, 0]*X - x(1) - x(2) >= 0
                             x(1) [0, 1; 1, 3] + x(2) [3, 1; 1, 0] - [1, 0; 0, 1] >> 0
                             X >> 0
*/

#include <stdio.h>
#include <math.h>

#include "mosek.h"    /* Include the MOSEK definition file.  */

#define NUMVAR    2   /* Number of scalar variables */
#define NUMAFE    4   /* Number of affine expressions        */
#define NUMFNZ    6   /* Number of non-zeros in F            */
#define NUMBARVAR 1   /* Number of semidefinite variables    */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
    printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
    MSKrescodee  r;

    const MSKint32t    DIMBARVAR[] = {2};         /* Dimension of semidefinite cone */
          MSKint64t    LENBARVAR[] = {2 * (2 + 1) / 2}; /* Number of scalar SD variables  */

    const MSKint32t   barc_j[] = {0, 0},
                      barc_k[] = {0, 1},
                      barc_l[] = {0, 1};
    const MSKrealt    barc_v[] = {1, 1};

    const MSKint64t   barf_i[] = {0,0};
    const MSKint32t   barf_j[] = {0,0},
                      barf_k[] = {0,1},
                      barf_l[] = {0,0};
    const MSKrealt    barf_v[] = {0,1};

    const MSKint64t   afeidx[] = {0, 0, 1, 2, 2, 3};
    const MSKint32t   varidx[] = {0, 1, 1, 0, 1, 0};
    const MSKrealt     f_val[] = {-1, -1, 3, sqrt(2), sqrt(2), 3},
                           g[] = {0, -1, 0, -1};

    MSKrealt *xx, *barx;
    MSKint32t i, j;
    
    MSKenv_t     env = NULL;
    MSKtask_t    task = NULL;

    /* Create the mosek environment. */
    r = MSK_makeenv(&env, NULL);

    if (r == MSK_RES_OK)
    {
        /* Create the optimization task. */
        r = MSK_maketask(env, 0, 0, &task);

        if (r == MSK_RES_OK)
        {
            r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

            /* Append 'NUMAFE' empty affine expressions. */
            if (r == MSK_RES_OK)
                r = MSK_appendafes(task, NUMAFE);

            /* Append 'NUMVAR' scalar variables.
            The variables will initially be fixed at zero (x=0). */
            if (r == MSK_RES_OK)
                r = MSK_appendvars(task, NUMVAR);

            /* Append 'NUMBARVAR' semidefinite variables. */
            if (r == MSK_RES_OK) 
            {
                r = MSK_appendbarvars(task, NUMBARVAR, DIMBARVAR);
            }

            /* Set the constant term in the objective. */
            if (r == MSK_RES_OK)
                r = MSK_putcfix(task, 1.0);

            /* Set c_j and the bounds for each scalar variable*/
            for (j = 0; j < NUMVAR && r == MSK_RES_OK; ++j)
            {
                r = MSK_putcj(task, j, 1.0);
                if (r==MSK_RES_OK)
                    r = MSK_putvarbound(task, j, MSK_BK_FR, -MSK_INFINITY, MSK_INFINITY);
            }

            /* Set the linear term barc_j in the objective.*/
            if (r == MSK_RES_OK)
                r = MSK_putbarcblocktriplet(task,
                                            2,
                                            barc_j, 
                                            barc_k,
                                            barc_l,
                                            barc_v);

            /* Set the F matrix */
            if (r == MSK_RES_OK)
                r = MSK_putafefentrylist(task, NUMFNZ, afeidx, varidx, f_val);
            /* Set the g vector */
            if (r == MSK_RES_OK)
                r = MSK_putafegslice(task, 0, NUMAFE, g);

            /* Set the barF matrix */
            if (r == MSK_RES_OK)
                r = MSK_putafebarfblocktriplet(task, 
                                               2, 
                                               barf_i, 
                                               barf_j, 
                                               barf_k,
                                               barf_l,
                                               barf_v);

            /* Append R+ domain and the corresponding ACC */
            MSKint64t acc1_afeidx[] = {0};
            if (r == MSK_RES_OK)
                r = MSK_appendrplusdomain(task, 1, NULL);
            if (r == MSK_RES_OK)
                r = MSK_appendacc(task, 0, 1, acc1_afeidx, NULL);
            
            /* Append the SVEC_PSD domain and the corresponding ACC */
            MSKint64t acc2_afeidx[] = {1, 2, 3};
            if (r == MSK_RES_OK)
                r = MSK_appendsvecpsdconedomain(task, 3, NULL);
            if (r == MSK_RES_OK)
                r = MSK_appendacc(task, 1, 3, acc2_afeidx, NULL);

            if (r == MSK_RES_OK)
            {
                MSKrescodee trmcode;

                /* Run optimizer */
                r = MSK_optimizetrm(task, &trmcode);

                /* Print a summary containing information
                about the solution for debugging purposes*/
                MSK_solutionsummary(task, MSK_STREAM_MSG);

                if (r == MSK_RES_OK)
                {
                    MSKsolstae solsta;
                    MSK_getsolsta(task, MSK_SOL_ITR, &solsta);

                    switch (solsta)
                    {
                        case MSK_SOL_STA_OPTIMAL:
                        xx   = (MSKrealt *) MSK_calloctask(task, NUMVAR, sizeof(MSKrealt));
                        barx = (MSKrealt *) MSK_calloctask(task, LENBARVAR[0], sizeof(MSKrealt));

                        MSK_getxx(task,MSK_SOL_ITR,xx);
                        MSK_getbarxj(task, MSK_SOL_ITR, 0, barx);

                        printf("Optimal primal solution\n");
                        for (i = 0; i < NUMVAR; ++i)
                            printf("x[%d]   : % e\n", i, xx[i]);

                        for (i = 0; i < LENBARVAR[0]; ++i)
                            printf("barx[%d]: % e\n", i, barx[i]);

                        MSK_freetask(task, xx);
                        MSK_freetask(task, barx);
                        break;

                        case MSK_SOL_STA_DUAL_INFEAS_CER:
                        case MSK_SOL_STA_PRIM_INFEAS_CER:
                            printf("Primal or dual infeasibility certificate found.\n");
                            break;
                        case MSK_SOL_STA_UNKNOWN:
                            printf("The status of the solution could not be determined. Termination code: %d.\n", trmcode);
                            break;
                        default:
                            printf("Other solution status.");
                            break;
                    }
                }
                else
                    printf("Error while optimizing.\n");
            }
            if (r != MSK_RES_OK)
            {
                /* In case of an error print error code and description. */
                char symname[MSK_MAX_STR_LEN];
                char desc[MSK_MAX_STR_LEN];
                printf("An error occurred while optimizing.\n");
                MSK_getcodedesc(r, symname, desc);
                printf("Error %s - '%s'\n", symname, desc);
            }
        }
        /* Delete the task and the associated data. */
        MSK_deletetask(&task);
    }

    /* Delete the environment and the associated data. */
    MSK_deleteenv(&env);

    return (r);
} /* main */

sensitivity.c

Listing 17.38 sensitivity.c Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:      sensitivity.c

  Purpose:   To demonstrate how to perform sensitivity
             analysis from the API on a small problem:

  minimize

  obj: +1 x11 + 2 x12 + 5 x23 + 2 x24 + 1 x31 + 2 x33 + 1 x34
  st
  c1:   +  x11 +   x12                                           <= 400
  c2:                  +   x23 +   x24                           <= 1200
  c3:                                  +   x31 +   x33 +   x34   <= 1000
  c4:   +  x11                         +   x31                   = 800
  c5:          +   x12                                           = 100
  c6:                  +   x23                 +   x33           = 500
  c7:                          +   x24                 +   x34   = 500

  The example uses basis type sensitivity analysis.
*/
#include <stdio.h>

#include "mosek.h" /* Include the MOSEK definition file. */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, char *argv[])
{
  const MSKint32t numcon = 7,
                  numvar = 7;
  MSKint32t       i, j;
  MSKboundkeye    bkc[] = {MSK_BK_UP, MSK_BK_UP, MSK_BK_UP, MSK_BK_FX,
                           MSK_BK_FX, MSK_BK_FX, MSK_BK_FX
                          };
  MSKboundkeye    bkx[] = {MSK_BK_LO, MSK_BK_LO, MSK_BK_LO,
                           MSK_BK_LO, MSK_BK_LO, MSK_BK_LO, MSK_BK_LO
                          };
  MSKint32t       ptrb[] = {0, 2, 4, 6, 8, 10, 12};
  MSKint32t       ptre[] = {2, 4, 6, 8, 10, 12, 14};
  MSKint32t       sub[] = {0, 3, 0, 4, 1, 5, 1, 6, 2, 3, 2, 5, 2, 6};
  MSKrealt        blc[] = { -MSK_INFINITY, -MSK_INFINITY, -MSK_INFINITY, 800, 100, 500, 500};
  MSKrealt        buc[] = {400,          1200,         1000,         800, 100, 500, 500};
  MSKrealt        c[]   = {1.0, 2.0, 5.0, 2.0, 1.0, 2.0, 1.0};
  MSKrealt        blx[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  MSKrealt        bux[] = {MSK_INFINITY, MSK_INFINITY, MSK_INFINITY, MSK_INFINITY,
                           MSK_INFINITY, MSK_INFINITY, MSK_INFINITY
                          };
  MSKrealt        val[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
  MSKrescodee     r;

  MSKenv_t        env;
  MSKtask_t       task;

  /* Create mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    /* Make the optimization task. */
    r = MSK_makeemptytask(env, &task);

    if (r == MSK_RES_OK)
    {
      /* Directs the log task stream to the user
         specified procedure 'printstr'. */

      MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      MSK_echotask(task,
                   MSK_STREAM_MSG,
                   "Defining the problem data.\n");
    }

    /* Append the constraints. */
    if (r == MSK_RES_OK)
      r = MSK_appendcons(task, numcon);

    /* Append the variables. */
    if (r == MSK_RES_OK)
      r = MSK_appendvars(task, numvar);

    /* Put C. */
    if (r == MSK_RES_OK)
      r = MSK_putcfix(task, 0.0);

    if (r == MSK_RES_OK)
      r = MSK_putcslice(task, 0, numvar, c);

    /* Put constraint bounds. */
    if (r == MSK_RES_OK)
      r = MSK_putconboundslice(task, 0, numcon, bkc, blc, buc);

    /* Put variable bounds. */
    if (r == MSK_RES_OK)
      r = MSK_putvarboundslice(task, 0, numvar, bkx, blx, bux);

    /* Put A. */
    if (r == MSK_RES_OK)
      r = MSK_putacolslice(task, 0, numvar, ptrb, ptre, sub, val);

    if (r == MSK_RES_OK)
      r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MINIMIZE);

    if (r == MSK_RES_OK)
      r = MSK_optimizetrm(task, NULL);

    if (r == MSK_RES_OK)
    {
      /* Analyze upper bound on c1 and the equality constraint on c4 */
      MSKint32t subi[] = {0, 3};
      MSKmarke marki[] = {MSK_MARK_UP, MSK_MARK_UP};

      /* Analyze lower bound on the variables x12 and x31 */
      MSKint32t subj[] = {1, 4};
      MSKmarke markj[] = {MSK_MARK_LO, MSK_MARK_LO};

      MSKrealt leftpricei[2];
      MSKrealt rightpricei[2];
      MSKrealt leftrangei[2];
      MSKrealt rightrangei[2];
      MSKrealt leftpricej[2];
      MSKrealt rightpricej[2];
      MSKrealt leftrangej[2];
      MSKrealt rightrangej[2];

      r = MSK_primalsensitivity(task,
                                2,
                                subi,
                                marki,
                                2,
                                subj,
                                markj,
                                leftpricei,
                                rightpricei,
                                leftrangei,
                                rightrangei,
                                leftpricej,
                                rightpricej,
                                leftrangej,
                                rightrangej);

      printf("Results from sensitivity analysis on bounds:\n");

      printf("For constraints:\n");
      for (i = 0; i < 2; ++i)
        printf("leftprice = %e, rightprice = %e,leftrange = %e, rightrange =%e\n",
               leftpricei[i], rightpricei[i], leftrangei[i], rightrangei[i]);

      printf("For variables:\n");
      for (i = 0; i < 2; ++i)
        printf("leftprice = %e, rightprice = %e,leftrange = %e, rightrange =%e\n",
               leftpricej[i], rightpricej[i], leftrangej[i], rightrangej[i]);
    }

    if (r == MSK_RES_OK)
    {
      MSKint32t subj[]  = {2, 5};
      MSKrealt  leftprice[2];
      MSKrealt  rightprice[2];
      MSKrealt  leftrange[2];
      MSKrealt  rightrange[2];

      r = MSK_dualsensitivity(task,
                              2,
                              subj,
                              leftprice,
                              rightprice,
                              leftrange,
                              rightrange
                             );

      printf("Results from sensitivity analysis on objective coefficients:\n");

      for (i = 0; i < 2; ++i)
        printf("leftprice = %e, rightprice = %e,leftrange = %e, rightrange =%e\n",
               leftprice[i], rightprice[i], leftrange[i], rightrange[i]);
    }

    MSK_deletetask(&task);
  }

  MSK_deleteenv(&env);

  printf("Return code: %d (0 means no error occured.)\n", r);
  return (r);
} /* main */

simple.c

Listing 17.39 simple.c Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:    simple.c

  Purpose: To demonstrate a very simple example using MOSEK by
           reading a problem file, solving the problem and
           writing the solution to a file.
*/
#include "mosek.h"

static void MSKAPI printstr(void *handle, const char str[])
{
  printf("%s", str);
}

int main(int argc, const char * argv[])
{
  MSKenv_t    env  = NULL;
  MSKtask_t   task = NULL;
  MSKrescodee res  = MSK_RES_OK, trmcode;

  if (argc <= 1)
  {
    printf("Missing argument, syntax is:\n");
    printf("  simple inputfile [ solutionfile ]\n");
  }
  else
  {
    // Create the mosek environment.
    // The `NULL' arguments here, are used to specify customized
    // memory allocators and a memory debug file. These can
    // safely be ignored for now.
    res = MSK_makeenv(&env, NULL);

    // Create a task object linked with the environment env.
    // We create it with 0 variables and 0 constraints initially,
    // since we do not know the size of the problem.
    if (res == MSK_RES_OK)
      res = MSK_maketask(env, 0, 0, &task);

    // Direct the task log stream to a user specified function
    if (res == MSK_RES_OK)
      res = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

    // We assume that a problem file was given as the first command
    // line argument (received in `argv')
    if (res == MSK_RES_OK)
      res = MSK_readdata(task, argv[1]);

    // Solve the problem
    if (res == MSK_RES_OK)
      res = MSK_optimizetrm(task, &trmcode);

    // Print a summary of the solution.
    if (res == MSK_RES_OK)
      res = MSK_solutionsummary(task, MSK_STREAM_LOG);

    // If an output file was specified, save problem to file
    if (res == MSK_RES_OK && argc >= 3)
    {
      // If using OPF format, these parameters will specify what to include in output
      MSK_putintparam(task, MSK_IPAR_OPF_WRITE_PROBLEM,    MSK_ON);
      MSK_putintparam(task, MSK_IPAR_OPF_WRITE_SOLUTIONS,  MSK_ON);
      MSK_putintparam(task, MSK_IPAR_OPF_WRITE_HINTS,      MSK_OFF);
      MSK_putintparam(task, MSK_IPAR_OPF_WRITE_PARAMETERS, MSK_OFF);

      res = MSK_writedata(task, argv[2]);
    }

    // Delete task and environment
    MSK_deletetask(&task);
    MSK_deleteenv(&env);
  }
  return res;
}

solutionquality.c

Listing 17.40 solutionquality.c Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:    solutionquality.c

  Purpose: To demonstrate how to examine the quality of a solution.
*/

#include <math.h>
#include "mosek.h"

static double dmin(double x,
                   double y)
{
  return (x <= y) ? (x) : (y);
} /* dmin */

static double dmax(double x,
                   double y)
{
  return (x >= y) ? (x) : (y);
} /* dmax */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */


int main(int argc, const char * argv[])
{
  double max_primal_viol, /* maximal primal violation */
         max_dual_viol,  /* maximal dual violation */
         abs_obj_gap,
         rel_obj_gap;

  MSKenv_t    env      = NULL;

  MSKint32t   numvar, j;

  MSKsolstae  solsta;
  MSKsoltypee whichsol = MSK_SOL_BAS;

  MSKrealt    primalobj, pviolcon, pviolvar, pviolbarvar, pviolcones, pviolitg,
              dualobj, dviolcon, dviolvar, dviolbarvar, dviolcones, xj;

  MSKrescodee r        = MSK_RES_OK;
  MSKrescodee trmcode;

  MSKtask_t   task     = NULL;

  int         accepted = 0;

  if (argc <= 1)
  {
    printf("Missing argument. The syntax is:\n");
    printf(" solutionquality inputfile\n");
  }
  else
  {

    r = MSK_makeenv(&env, NULL);

    if (r == MSK_RES_OK)
      r = MSK_makeemptytask(env, &task);

    if (r == MSK_RES_OK)
      MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

    /* We assume that a problem file was given as the first command
       line argument (received in `argv'). */
    if (r == MSK_RES_OK)
      r = MSK_readdata(task, argv[1]);

    /* Solve the problem */
    if (r == MSK_RES_OK)
    {
      r = MSK_optimizetrm(task, &trmcode);
    }

    /* Print a summary of the solution. */
    MSK_solutionsummary(task, MSK_STREAM_MSG);

    if (r == MSK_RES_OK)
    {

      MSK_getsolsta(task, whichsol, &solsta);

      r = MSK_getsolutioninfo(task,
                              whichsol,
                              &primalobj,
                              &pviolcon,
                              &pviolvar,
                              &pviolbarvar,
                              &pviolcones,
                              &pviolitg,
                              &dualobj,
                              &dviolcon,
                              &dviolvar,
                              &dviolbarvar,
                              &dviolcones);

      switch (solsta)
      {
        case MSK_SOL_STA_OPTIMAL:
          {

            abs_obj_gap     = fabs(dualobj - primalobj);
            rel_obj_gap     = abs_obj_gap / (1.0 + dmin(fabs(primalobj), fabs(dualobj)));
            max_primal_viol = dmax(pviolcon, pviolvar);
            max_primal_viol = dmax(max_primal_viol, pviolbarvar);
            max_primal_viol = dmax(max_primal_viol, pviolcones);

            max_dual_viol   = dmax(dviolcon, dviolvar);
            max_dual_viol   = dmax(max_dual_viol, dviolbarvar);
            max_dual_viol   = dmax(max_dual_viol, dviolcones);

            /* Assume the application needs the solution to be within
               1e-6 optimality in an absolute sense. Another approach
               would be looking at the relative objective gap */

            printf("\n\n");
            printf("Customized solution information.\n");
            printf("  Absolute objective gap: %e\n", abs_obj_gap);
            printf("  Relative objective gap: %e\n", rel_obj_gap);
            printf("  Max primal violation  : %e\n", max_primal_viol);
            printf("  Max dual violation    : %e\n", max_dual_viol);

            if (rel_obj_gap > 1e-6)
            {
              printf("Warning: The relative objective gap is LARGE.\n");
              accepted = 0;
            }

            /* We will accept a primal infeasibility of 1e-8 and
               dual infeasibility of 1e-6. These number should chosen problem
               dependent.
             */

            if (max_primal_viol > 1e-8)
            {
              printf("Warning: Primal violation is too LARGE.\n");
              accepted = 0;
            }

            if (max_dual_viol > 1e-6)
            {
              printf("Warning: Dual violation is too LARGE.\n");
              accepted = 0;
            }

            if (accepted)
            {

              if (MSK_RES_OK == MSK_getnumvar(task, &numvar))
              {
                printf("Optimal primal solution\n");
                for (j = 0; j < numvar && r == MSK_RES_OK; ++j)
                {
                  r = MSK_getxxslice(task, whichsol, j, j + 1, &xj);
                  if (r == MSK_RES_OK)
                    printf("x[%d]: %e\n", j, xj);
                }
              }
            }
            else if (r == MSK_RES_OK)
            {
              /* Print detailed information about the solution */
              r = MSK_analyzesolution(task, MSK_STREAM_LOG, whichsol);
            }
            break;
          }
        case MSK_SOL_STA_DUAL_INFEAS_CER:
        case MSK_SOL_STA_PRIM_INFEAS_CER:
          printf("Primal or dual infeasibility certificate found.\n");
          break;
        case MSK_SOL_STA_UNKNOWN:
          printf("The status of the solution is unknown. Termination code: %d.\n", trmcode);
          break;
        default:
          printf("Other solution status");
          break;
      }
    }

    MSK_deletetask(&task);
    MSK_deleteenv(&env);
  }
  return (r);
}

solvebasis.c

Listing 17.41 solvebasis.c Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:      solvebasis.c

  Purpose:   To demonstrate the usage of
             MSK_solvewithbasis on the problem:

             maximize  x0 + x1
             st.
                     x0 + 2.0 x1 <= 2
                     x0  +    x1 <= 6
                     x0 >= 0, x1>= 0

               The problem has the slack variables
               xc0, xc1 on the constraints
               and the variables x0 and x1.

               maximize  x0 + x1
               st.
                  x0 + 2.0 x1 -xc1       = 2
                  x0  +    x1       -xc2 = 6
                  x0 >= 0, x1>= 0,
                  xc1 <=  0 , xc2 <= 0

  Syntax:    solvebasis

 */
#include "mosek.h"

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  MSKenv_t     env;
  MSKtask_t    task;
  MSKint32t    numcon = 2, numvar = 2;
  double       c[]    = {1.0, 1.0};
  MSKint32t    ptrb[] = {0, 2},
                        ptre[] = {2, 3};
  MSKint32t    asub[] = {0, 1, 0, 1};
  double       aval[] = {1.0, 1.0, 2.0, 1.0};
  MSKboundkeye bkc[]  = { MSK_BK_UP, MSK_BK_UP };
  double       blc[]  = { -MSK_INFINITY, -MSK_INFINITY };
  double       buc[]  = {2.0, 6.0};

  MSKboundkeye bkx[]  = { MSK_BK_LO, MSK_BK_LO };
  double       blx[]  = {0.0, 0.0};
  double       bux[]  = { +MSK_INFINITY, +MSK_INFINITY};
  MSKrescodee  r      = MSK_RES_OK;
  MSKint32t    i, nz;
  double       w[]    = {2.0, 6.0};
  MSKint32t    sub[]  = {0, 1};
  MSKint32t    *basis;

  if (r == MSK_RES_OK)
    r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
    r = MSK_makeemptytask(env, &task);

  if (r == MSK_RES_OK)
    MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

  if (r == MSK_RES_OK)
    r = MSK_inputdata(task, numcon, numvar, numcon, numvar,
                      c, 0.0,
                      ptrb, ptre, asub, aval, bkc, blc, buc, bkx, blx, bux);

  if (r == MSK_RES_OK)
    r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE);

  if (r == MSK_RES_OK)
    r = MSK_optimizetrm(task, NULL);

  if (r == MSK_RES_OK)
    basis = MSK_calloctask(task, numcon, sizeof(MSKint32t));

  if (r == MSK_RES_OK)
    r = MSK_initbasissolve(task, basis);

  /* List basis variables corresponding to columns of B */
  for (i = 0; i < numcon && r == MSK_RES_OK; ++i)
  {
    printf("basis[%d] = %d\n", i, basis[i]);
    if (basis[sub[i]] < numcon)
      printf("Basis variable no %d is xc%d.\n", i, basis[i]);
    else
      printf("Basis variable no %d is x%d.\n", i, basis[i] - numcon);
  }

  nz = 2;
  /* solve Bx = w */
  /* sub contains index of non-zeros in w.
     On return w contains the solution x and sub
     the index of the non-zeros in x.
   */
  if (r == MSK_RES_OK)
    r = MSK_solvewithbasis(task, 0, nz, sub, w, &nz);

  if (r == MSK_RES_OK)
  {
    printf("\nSolution to Bx = w:\n\n");

    /* Print solution and b. */

    for (i = 0; i < nz; ++i)
    {
      if (basis[sub[i]] < numcon)
        printf("xc%d = %e\n", basis[sub[i]] , w[sub[i]]);
      else
        printf("x%d = %e\n", basis[sub[i]] - numcon , w[sub[i]]);
    }
  }

  /* Solve B^T y = w */
  nz     = 1;     /* Only one element in sub is nonzero. */
  sub[0] = 1;     /* Only w[1] is nonzero. */
  w[0]   = 0.0;
  w[1]   = 1.0;

  if (r == MSK_RES_OK)
    r = MSK_solvewithbasis(task, 1, nz, sub, w, &nz);

  if (r == MSK_RES_OK)
  {
    printf("\nSolution to B^T y = w:\n\n");
    /* Print solution and y. */
    for (i = 0; i < nz; ++i)
      printf("y%d = %e\n", sub[i], w[sub[i]]);
  }

  return (r);
}/* main */

solvelinear.c

Listing 17.42 solvelinear.c Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File     :  solvelinear.c

   Purpose  :  To demonstrate the usage of MSK_solvewithbasis
               to solve the linear system:

                1.0  x1             = b1
               -1.0  x0  +  1.0  x1 = b2

               with two different right hand sides

               b = (1.0, -2.0)

               and

               b = (7.0, 0.0)
 */

#include "mosek.h"

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */


static MSKrescodee setup(MSKtask_t task,
                         double    *aval,
                         MSKint32t *asub,
                         MSKint32t *ptrb,
                         MSKint32t *ptre,
                         MSKint32t numvar,
                         MSKint32t *basis)

{
  MSKint32t   i, j;
  MSKrescodee r = MSK_RES_OK;
  MSKstakeye *skx = NULL , *skc = NULL;

  skx = (MSKstakeye *) calloc(numvar, sizeof(MSKstakeye));
  if (skx == NULL && numvar)
    r = MSK_RES_ERR_SPACE;

  skc = (MSKstakeye *) calloc(numvar, sizeof(MSKstakeye));
  if (skc == NULL && numvar)
    r = MSK_RES_ERR_SPACE;

  for (i = 0; i < numvar && r == MSK_RES_OK; ++i)
  {
    skx[i] = MSK_SK_BAS;
    skc[i] = MSK_SK_FIX;
  }

  /* Create a coefficient matrix and right hand
     side with the data from the linear system */
  if (r == MSK_RES_OK)
    r = MSK_appendvars(task, numvar);

  if (r == MSK_RES_OK)
    r = MSK_appendcons(task, numvar);

  for (i = 0; i < numvar && r == MSK_RES_OK; ++i)
    r = MSK_putacol(task, i, ptre[i] - ptrb[i], asub + ptrb[i], aval + ptrb[i]);

  for (i = 0; i < numvar && r == MSK_RES_OK; ++i)
    r = MSK_putconbound(task, i, MSK_BK_FX, 0, 0);

  for (i = 0; i < numvar && r == MSK_RES_OK; ++i)
    r = MSK_putvarbound(task, i, MSK_BK_FR, -MSK_INFINITY, MSK_INFINITY);

  /* Allocate space for the solution and set status to unknown */

  if (r == MSK_RES_OK)
    r = MSK_deletesolution(task, MSK_SOL_BAS);

  /* Setup status keys. That is all that is needed. */
  if (r == MSK_RES_OK)
    r = MSK_putskcslice(task, MSK_SOL_BAS, 0, numvar, skc);

  if (r == MSK_RES_OK)
    r = MSK_putskxslice(task, MSK_SOL_BAS, 0, numvar, skx);

  if (r == MSK_RES_OK)
    r = MSK_initbasissolve(task, basis);

  free(skx);
  free(skc);

  return (r);

}

#define NUMCON 2
#define NUMVAR 2

int main(int argc, const char *argv[])
{
  const MSKint32t numvar = NUMCON,
                  numcon = NUMVAR;   /* we must have numvar == numcon */
  MSKenv_t        env;
  MSKtask_t       task;
  MSKrescodee     r = MSK_RES_OK;
  MSKint32t       i, nz;
  double          aval[] = { -1.0, 1.0, 1.0};
  MSKint32t       asub[] = {1, 0, 1};
  MSKint32t       ptrb[] = {0, 1};
  MSKint32t       ptre[] = {1, 3};
  MSKint32t       bsub[NUMCON];
  double          b[NUMCON];
  MSKint32t       *basis = NULL;

  if (r == MSK_RES_OK)
    r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
    r = MSK_makeemptytask(env, &task);

  if (r == MSK_RES_OK)