solvermosek.hpp

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>

#include "declarations.h"

//#include<iostream>
//#include<fstream>
//#include<iomanip>

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

//#define NUMCON    14   /* Number of constraints.              */
//#define NUMBARVAR 4   /* Number of semidefinite variables    *//


#include "../arithm/Scl.hpp"
#include "../arithm/Scl_mpf.hpp"


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

template<typename T>
int solvermosek1(int n,int k,T *a, constraintmatrixv<T> *constraints,int *listsize,int dimcon,double *&solmosek,double *&solprimal, int &solst,double *&soldual,int &dimsoldual)
{


    MSKrescodee  r;


      MSKint32t    NUMCON;
      MSKint32t    NUMBARVAR;
      MSKint32t    DIMBARVAR[]={};         /* Dimension of semidefinite cone */
      MSKint64t    LENBARVAR[]={}; /* Number of scalar SD variables  */

      struct sparseblockv<T> *p,*q;

//      MSKboundkeye bkc[] = { MSK_BK_FX, MSK_BK_FX, MSK_BK_FX, MSK_BK_FX,MSK_BK_FX, MSK_BK_FX, MSK_BK_FX, MSK_BK_FX };
//      double       blc[] = { 0,0,0,0,0,0,1,0};
//      double       buc[] = { 0,0,0,0,0,0,1,0};

//

//        double       blc[] = { 0,0,0,1.5,1.0,1.5,1.5,1.5};
//        double       buc[] = { 0,0,0,1.5,1.0,1.5,1.5,1.5};

        MSKint32t    *barc_i,
                     *barc_j;
        double       *barc_v;



        MSKint32t    *bara_i,
                     *bara_j;
        double       *bara_v;

      MSKint32t    i,j,l,m,numb;
      MSKint64t    idx,num=0,numa;

      double       falpha=1.0;


      double       *barx1;
      double       *bars1;
      double       *yy,*primalobj;




      cout<<"dentro de mosek"<<endl;
      NUMCON=k-1;
      NUMBARVAR=1;
      DIMBARVAR[0]=n;
      cout<<"n"<<n<<endl;
      for(i=0;i<dimcon;i++) DIMBARVAR[0]+=listsize[i];
      LENBARVAR[0]=DIMBARVAR[0]*(DIMBARVAR[0]+1)/2;

      cout<<"dentro de mosek avant task"<<endl;
      cout<<"NUMCON"<<NUMCON<<endl;
      cout<<"DIMBARVAR[0]"<<DIMBARVAR[0]<<endl;
      cout<<"LENBARVAR[0]"<<LENBARVAR[0]<<endl;
      MSKenv_t     env = NULL;
      MSKtask_t    task = NULL;

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

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

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

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



          /* Append 'NUMBARVAR' semidefinite variables. */
          if ( r == MSK_RES_OK ) {
              r = MSK_appendbarvars(task, NUMBARVAR,DIMBARVAR);
          }

          /* Optionally add a constant term to the objective. */
          if ( r ==MSK_RES_OK )
            r = MSK_putcfix(task,0.0);


        cout<<"dentro de mosek avant task2"<<endl;
          /* Set the linear term barc_j in the objective.*/
          q=constraints[k].blocks;
          num=q->numentries;

          barc_i = (MSKint32t*) MSK_calloctask(task,num,sizeof(MSKintt));
          barc_j = (MSKint32t*) MSK_calloctask(task,num,sizeof(MSKintt));
          barc_v= (double*) MSK_calloctask(task,num,sizeof(MSKrealt));

          for (l=0;l<num;l++)
          {
              barc_i[l]=q->jindices[l+1]-1;
              barc_j[l]=q->iindices[l+1]-1;
            //put this cast from float to double (doesn't work cast simple to double)

              barc_v[l]= Scl<MPF>::castDouble(q->entries[l+1])*(-1);

//              cout<<"barc_i[l]"<<barc_i[l]<<endl;
//              cout<<"barc_j[l]"<<barc_j[l]<<endl;
      //        cout<<"barc_v[l]"<<barc_v[l]<<endl;



          }

       cout<<"dentro de mosek avant task2bis"<<"num"<<num<<endl;

           if ( r == MSK_RES_OK )
            r = MSK_appendsparsesymmat(task,
                                       DIMBARVAR[0],
                                       num,
                                       barc_i,
                                       barc_j,
                                       barc_v,
                                       &idx);

          if ( r == MSK_RES_OK )
            r = MSK_putbarcj(task, 0, 1, &idx, &falpha);


     cout<<"dentro de mosek avant task3"<<endl;

    /* 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)
          {
       //      cout<<"a[i+1]"<<a[i+1]<<endl;
            r = MSK_putconbound(task,
                               i,           /* Index of constraint.*/
                               MSK_BK_FX,      /* Bound key.*/
                               Scl<MPF>::castDouble(a[i+1]),      /* Numerical value of lower bound.*/
                               Scl<MPF>::castDouble(a[i+1]));     /* Numerical value of upper bound.*/
          }

           cout<<"dentro de mosek avant task4"<<endl;
          /* Add the first constraint of barA */
         numa=0;

         for (i=1;i<=NUMCON;i++)
          {
           p=constraints[i].blocks;

           numa+=p->numentries;
          }

         bara_i = (MSKint32t*) MSK_calloctask(task,numa,sizeof(MSKintt));
         bara_j = (MSKint32t*) MSK_calloctask(task,numa,sizeof(MSKintt));
         bara_v= (double*) MSK_calloctask(task,numa,sizeof(MSKrealt));
          m=0;
          for (i=1;i<=NUMCON;i++)
          {
            p=constraints[i].blocks;


          for (l=0;l<p->numentries;l++)
           {

              bara_i[l+m]=p->jindices[l+1]-1;
              bara_j[l+m]=p->iindices[l+1]-1;
              bara_v[l+m]=Scl<MPF>::castDouble(p->entries[l+1]);


//              cout<<"bara_i[l+m]"<<bara_i[l+m]<<endl;
//              cout<<"bara_j[l+m]"<<bara_j[l+m]<<endl;
//              cout<<"bara_v[l+m]"<<bara_v[l+m]<<endl;



           }
           m+=p->numentries;
          }
          numb=0;
          for (i=1;i<=NUMCON;i++)
          {
            p=constraints[i].blocks;

             if ( r==MSK_RES_OK )
               r = MSK_appendsparsesymmat(task,
                                          DIMBARVAR[0],
                                  p->numentries,
                                  bara_i+numb,
                                  bara_j+numb,
                                  bara_v+numb,
                                  &idx);

            if ( r==MSK_RES_OK )
            r = MSK_putbaraij(task, i-1, 0, 1, &idx, &falpha);

            numb+=p->numentries;



          }

        cout<<"dentro de mosek avant task5"<<endl;







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

          if ( r==MSK_RES_OK )
          {

            cout<<"gap"<<MSK_putdouparam(task,MSK_DPAR_INTPNT_CO_TOL_REL_GAP,1.0e-6)<<endl;
//           // MSK_putdouparam(task,MSK_DPAR_INTPNT_CO_TOL_NEAR_REL,1.0);
//            MSK_putdouparam(task,MSK_DPAR_INTPNT_CO_TOL_INFEAS,1.0e-3);
//            MSK_putdouparam(task,MSK_DPAR_INTPNT_CO_TOL_DFEAS,1.0e-3);
           // MSK_putdouparam(task,MSK_DPAR_INTPNT_CO_TOL_NEAR_REL,1.0e+100);
           // MSK_putdouparam(task,MSK_DPAR_INTPNT_CO_TOL_PFEAS,1.0e-100);
            //cout<<"quelle type de probleme"<<MSK_putintparam(task,MSK_IPAR_INTPNT_SOLVE_FORM,MSK_SOLVE_PRIMAL)<<endl;
            //cout<<"pfeas"<<MSK_putdouparam(task,MSK_DPAR_INTPNT_CO_TOL_PFEAS,1.0e-5)<<endl;
            //cout<<"dfeas"<<MSK_putdouparam(task,MSK_DPAR_INTPNT_CO_TOL_DFEAS,1.0e-5)<<endl;

            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);
              //int solst;
              solst=solsta;
              cout<<"solsta"<<solst<<endl;
              dimsoldual= LENBARVAR[0];

              switch(solsta)
              {
                case MSK_SOL_STA_OPTIMAL:
                case MSK_SOL_STA_NEAR_OPTIMAL:
                case MSK_SOL_STA_UNKNOWN:

                  barx1 = (double*) MSK_calloctask(task,LENBARVAR[0],sizeof(MSKrealt));
                  bars1 = (double*) MSK_calloctask(task,LENBARVAR[0],sizeof(MSKrealt));
                  yy = (double*) MSK_calloctask(task,NUMCON,sizeof(MSKrealt));
                  solmosek=(double*) MSK_calloctask(task,NUMCON+1,sizeof(MSKrealt));
                  primalobj=(double*) MSK_calloctask(task,1,sizeof(MSKrealt));
                  solprimal=(double*) MSK_calloctask(task,1,sizeof(MSKrealt));
                  soldual = (double*) MSK_calloctask(task,LENBARVAR[0],sizeof(MSKrealt));

                  // obtain primal solution, semidefinite matrix
                  MSK_getbarxj(task,
                               MSK_SOL_ITR,    /* Request the interior solution. */
                               0,
                               barx1);

                  // obtain dual solution, semidefinite matrix
                  MSK_getbarsj(task,
                               MSK_SOL_ITR,    /* Request the interior solution. */
                               0,
                               bars1);
                  MSK_gety(task,
                               MSK_SOL_ITR,    /* Request the interior solution. */
                               yy);

                  MSK_getprimalobj(task,
                                  MSK_SOL_ITR,
                                  primalobj);

                  printf("Optimal primal solution\n");

//                  for(i=0; i<LENBARVAR[0]; ++i)
//                    printf("barx1[%d]: % e\n",i,barx1[i]);

                  printf("Optimal dual solution\n");
                  for(i=0; i<LENBARVAR[0]; ++i)
                   {
                   // printf("bars1[%d]: % e\n",i,-bars1[i]);
                    soldual[i]=-bars1[i];
                  }
                  solprimal[0]=primalobj[0];

                  cout<<"Primal Solution "<<solprimal[0]<<endl;

                  for(i=0; i<NUMCON; ++i)
                   {
  //                  printf("yy[%d]: % e\n",i,yy[i]);
                    solmosek[i]=yy[i];
                  }
                  solmosek[NUMCON]=1.0;

                  MSK_freetask(task,barx1);

                  MSK_freetask(task,bars1);
                  MSK_freetask(task,yy);

                  break;

                case MSK_SOL_STA_DUAL_INFEAS_CER:
                case MSK_SOL_STA_PRIM_INFEAS_CER:
                case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
                case MSK_SOL_STA_NEAR_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.\n");
                 // 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 */