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"
#include <math.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 (fabs(sl[i]) > eps)
      printf("#%d, lower,  dual = %e\n", i, sl[i]);
    if (fabs(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

// Check if the interior point solution is an optimal point
MSKsolstae solsta;
MOSEKCALL(res, MSK_getsolsta(task, MSK_SOL_ITR, &solsta));
if (solsta != MSK_SOL_STA_OPTIMAL)
{
  // See https://docs.mosek.com/latest/capi/accessing-solution.html about handling solution statuses.
  fprintf(stderr, "Unexpected solution status: %d", solsta);
  return MSK_RES_ERR_UNHANDLED_SOLUTION_STATUS;
}

  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));

    // Check if the interior point solution is an optimal point
    MOSEKCALL(res, MSK_getsolsta(task, MSK_SOL_ITR, &solsta));
    if (solsta != MSK_SOL_STA_OPTIMAL)
    {
      // See https://docs.mosek.com/latest/capi/accessing-solution.html about handling solution statuses.
      fprintf(stderr, "Unexpected solution status: %d", solsta);
      return MSK_RES_ERR_UNHANDLED_SOLUTION_STATUS;
    }

    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

  // Check if the interior point solution is an optimal point
  MSKsolstae solsta;
  MOSEKCALL(res, MSK_getsolsta(task, MSK_SOL_ITR, &solsta));
  if (solsta != MSK_SOL_STA_OPTIMAL)
  {
    // See https://docs.mosek.com/latest/capi/accessing-solution.html about handling solution statuses.
    fprintf(stderr, "Unexpected solution status: %d", solsta);
    return MSK_RES_ERR_UNHANDLED_SOLUTION_STATUS;
  }

  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
  
// Check if the interior point solution is an optimal point
MSKsolstae solsta;
MOSEKCALL(res, MSK_getsolsta(task, MSK_SOL_ITG, &solsta));
if (solsta != MSK_SOL_STA_INTEGER_OPTIMAL)
{
  // See https://docs.mosek.com/latest/capi/accessing-solution.html about handling solution statuses.
  fprintf(stderr, "Unexpected solution status: %d", solsta);
  return MSK_RES_ERR_UNHANDLED_SOLUTION_STATUS;
}


  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

// Check if the interior point solution is an optimal point
MSKsolstae solsta;
MOSEKCALL(res, MSK_getsolsta(task, MSK_SOL_ITG, &solsta));
if (solsta != MSK_SOL_STA_INTEGER_OPTIMAL)
{
  // See https://docs.mosek.com/latest/capi/accessing-solution.html about handling solution statuses.
  fprintf(stderr, "Unexpected solution status: %d", solsta);
  return MSK_RES_ERR_UNHANDLED_SOLUTION_STATUS;
}
  
  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);