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

#include "mex.h"
#include "loqo.h"
#include "myalloc.h"

double _HUGE=1.0e+300;

void nlinit( void *vlp );
int nlupdate( void *vlp );

double objval( double *x );
void objgrad( double *c, double *x );
void hessian( double *Q, double *x, double *y );
void conval( double *h, double *x );
void congrad( double *A, double *At, double *x );

int rd_opns();
void *binsearch(char **sp);
// void set_opns( void );

void cleanup( void );

static int NCellElems, M, *Cm, *Cn, **iC, **kC, m, n, *iA, *kA, *iQ, *kQ;
static double *c, **C, *A, *b, *Q, *l, *u, *x0, *x, *y_lin, *y_cone; 
static int neq;
static LOQO *lp=NULL;
static char *cmdstr=NULL;
static char *opnstr=NULL;
static mxArray *xk;
static int bad_opns;

typedef struct keyword keyword;
struct keyword {
					char	*name;
					int		type;
					int		ivalue;
					double	dvalue;
				};
#define KW(a,b,c,d) {a,b,c,d}
static keyword keywds[] = {
	KW("bndpush", 2, 0, -1.0),			KW("convex", 0, 0, 0.0),
	KW("dense", 1, -1, 0.0),			KW("dual", 0, 0, 0.0),
	KW("epsdiag", 2, 0, 1.0e-10),		KW("epsnum", 2, 0, 0.0),
	KW("epssol", 2, 0, 1.0e-6),			KW("honor_bnds", 0, -1, 0.0),
	KW("honor_bnds_init", 0, 0, 0.0),	
	KW("inftol", 2, 0, 1e-6),			KW("inftol2", 2, 0, 1e+5),
	KW("iterlim", 1, 200, 0.0),			
	KW("max", 0, 0, 0.0),				KW("maximize", 0, 0, 0.0),
	KW("maxit", 1, 200, 0.0),			KW("min", 0, 0, 0.0),
	KW("mindeg", 0, 0, 0.0),			KW("minimize", 0, 0, 0.0),
	KW("minlocfill", 0, 0, 0.0),		
	KW("mufactor", 2, 0, -1.0),			KW("noreord", 0, 0, 0.0),
	KW("outlev", 1, 0, 0.0),
	KW("pred_corr", 1, -1, 0.0),		KW("primal", 0, 0, 0.0),
	KW("sdp", 0, 0, 0.0),				KW("sigfig", 1, 8, 0.0),
	KW("stablty", 2, 0, 1.0),			KW("steplen", 2, 0, 0.95),
	KW("timlim", 2, 0, -1.0),
	KW("verbose", 1, 0, 0.0)
};
	
int mloqo() 
{
    int i, j, k, ii, jj, kk, row, nz, qnz;
    int status;
    double *r=NULL;

    if (kA != NULL) {nz  = kA[n];} else { nz = 0;}
    /* if (kQ != NULL) {qnz = kQ[n];} else {qnz = 0;} */
    qnz = 0;    /* It's a pain to overlay Q with Hessian of SOCP constraints */

    if (kA == NULL) { CALLOC(kA, n+1, int); }

    if (lp == NULL) {lp = openlp();}

    M = 0;
    for (i=0; i<NCellElems; i++) M += Cm[i];

    lp->m   = m + M + NCellElems;
    lp->n   = n + M; 
    lp->nz  = nz + 2*M;  
        for (i=0; i<NCellElems; i++) {
                lp->nz += kC[i][Cn[i]]-kC[i][1];
        }
    lp->qnz = qnz + 3*M - 2*NCellElems;

    CALLOC( lp->c, lp->n, double );     
    for (j=M, jj=0; jj<n; j++, jj++) lp->c[j] = c[jj];

    MALLOC( lp->A, lp->nz, double );    /* numeric values updated in nlupdate*/
    MALLOC( lp->iA,lp->nz, int    );    
    MALLOC( lp->kA,lp->n+1,int    );    
    j=0; k=0; lp->kA[0] = 0;
    for (i=0; i<NCellElems; i++) {
        for (kk=0; kk<Cm[i]; kk++) {
            lp->iA[k] = i; 
            k++;                /* lp->A[k] will be given in nlupdate
*/

            lp->iA[k] = j+NCellElems;
            lp->A[k]  = 1;
            k++;

            j++;
            lp->kA[j] = k;
        }
    }
    for (jj=0; jj<n; jj++) {
        ii = NCellElems;
        for (i=0; i<NCellElems; i++) {
            for (kk=kC[i][jj+1]; kk<kC[i][jj+2]; kk++) {
                row = iC[i][kk] + ii;
                lp->iA[k] = row;
                lp->A[k] = -C[i][kk];
                k++;
            }
            ii += Cm[i];
        }

        for (kk=kA[jj]; kk<kA[jj+1]; kk++) {
            row = iA[kk] + ii;
            lp->iA[k] = row;
            lp->A[k] = A[kk];
            k++;
        }

        j++;
        lp->kA[j] = k;
    }

    CALLOC( lp->b, lp->m, double );     
    for (ii=0; ii<NCellElems; ii++) {
        lp->b[ii] = 0;
    }
    for (i=0; i<NCellElems; i++) {
        for (kk=kC[i][0]; kk<kC[i][1]; kk++) {
            row = iC[i][kk] + ii;
            lp->b[row] = C[i][kk];
        }
        ii += Cm[i];
    }
    for (kk=0; kk<m; kk++) {
        row = kk + ii;
        lp->b[row] = b[kk];
    }

    MALLOC( lp->Q, lp->qnz, double );   /* values updated in nlupdate()*/
    MALLOC( lp->iQ,lp->qnz, int    );   
    MALLOC( lp->kQ,lp->n+1, int    );   
    j=0; k=0; lp->kQ[0] = 0; ii = -1;
    for (i=0; i<NCellElems; i++) {
        ii += Cm[i];
        for (kk=0; kk<Cm[i]-1; kk++) { 
            lp->iQ[k] = j; 
            k++;                

            lp->iQ[k] = ii;
            k++;

            j++;
            lp->kQ[j] = k;
        }

        ii -= Cm[i];
        for (kk=1; kk<=Cm[i]; kk++) { 
            lp->iQ[k] = ii+kk;
            k++;
        }
        ii += Cm[i];
        j++;
        lp->kQ[j] = k;
    }
    for (j=M+1; j<=lp->n; j++) {lp->kQ[j] = lp->kQ[j-1];}

    MALLOC( lp->l, lp->n, double );
    for (j=0; j<lp->n; j++) lp->l[j] = -HUGE_VAL;
    j = -1;
    for (i=0; i<NCellElems; i++) {
        j += Cm[i];
        lp->l[j] = 0.0;         /* variable t shouldn't hit 0 */
    }
    j++;
    if (l != NULL) {
        for (jj=0; jj<n; jj++, j++) {
            lp->l[j] = l[jj];
        }
    }

    MALLOC( lp->u, lp->n, double );
    for (j=0; j<lp->n; j++) lp->u[j] = HUGE_VAL;
    j = M;
    if (u != NULL) {
        for (jj=0; jj<n; jj++, j++) {
            lp->u[j] = u[jj];
        }
    }

    FREE( lp->y );

    MALLOC( lp->x, lp->n, double );
    jj=0;
    for (i=0; i<NCellElems; i++) {
                for (j=0; j<Cm[i]-1; j++) {
                        lp->x[jj] = 0;
                        jj++;
                }
                lp->x[jj] = 1;
                jj++;
        }
    for (j=M; j<lp->n; j++) lp->x[j] = 1;
    j = M;
    if (x0 != NULL) {
        for (jj=0; jj<n; jj++, j++) {
            lp->x[j] = x0[jj];
        }
    }

    MALLOC(lp->r, lp->m, double); 
    for (i=0; i<NCellElems;    i++) lp->r[i] = HUGE_VAL;
    for (   ; i<lp->m-(m-neq); i++) lp->r[i] = 0;
    for (   ; i<lp->m;         i++) lp->r[i] = HUGE_VAL;

	if (opnstr != NULL) {
		bad_opns = 0;
		if (!rd_opns()) {
			printf("loqo_options not set correctly\n");
			return 0;
		}
	}
	set_opns(lp);

    nlsetup(lp);

    status = solvelp(lp);

    if (x != NULL) for (j=0, jj=lp->n-n; j<n; j++,jj++) { x[j] = lp->x[jj]; }
    if (y_lin != NULL) for (i=0; i<m; i++) y_lin[i] = lp->y[M+NCellElems+i];
	if (y_cone != NULL) for (i=0; i<NCellElems; i++) y_cone[i] = 2*(lp->y[i]);

    inv_clo();
    closelp(lp);
    lp=NULL;

    return status;
}

double objval( double *x )
{
    int j, jj;
    double val = 0;
    for (j=M, jj=0; jj<n; j++, jj++) val += c[jj]*x[j];
    return val;
}

void objgrad( double *c, double *x )
{   /*  objective function is linear and so no need to update */
}

void hessian( double *Q, double *x, double *yy )
{
    int i, kk;
    double y, u, t, normu2;
    int j=0, k=0, ii = -1, jj;

    for (i=0; i<NCellElems; i++) {
        ii += Cm[i];
        y = lp->y[i];
        t = lp->x[ii];
        for (kk=0; kk<Cm[i]-1; kk++) { 
            u = lp->x[j];
            lp->Q[k] = 2*y/t;
            k++;                

            lp->Q[k] = -2*y*u/(t*t);
            k++;

            j++;
        }

        normu2 = 0;
        for (kk=1; kk<Cm[i]; kk++) { 
            u = lp->x[lp->iQ[k]];
            normu2 += u*u;
            lp->Q[k] = -2*y*u/(t*t);
            k++;
        }
        lp->Q[k] = y*2*normu2/(t*t*t);
        k++;

        j++;
    }

    /*  for animation only */

    if (x0 != NULL) for (j=0, jj=lp->n-n; j<n; j++,jj++) { x0[j] = lp->x[jj]; }
    /*    mexEvalString(cmdstr); */

     if (cmdstr != NULL) {
        mexCallMATLAB(0,NULL,1,&xk,cmdstr); 
     }

    /*  end of for animation only */
}

void conval( double *h, double *x )
{
    int i, j, ii, kk;
    double u, t, normu2;

    smx(lp->m, lp->n, lp->A, lp->kA, lp->iA, x, h);

    j = 0;
    for (i=0; i<NCellElems; i++) {
        normu2=0;
        for (kk=0; kk<Cm[i]-1; kk++) { 
            u = lp->x[j];
            normu2 += u*u;
            j++;
        }
        t = lp->x[j];
        j++;
        h[i] = t - normu2/t;
    }
}

void congrad( double *A, double *At, double *x )
{
    int i, j, k, ii, kk;
    double u, t, normu2;

    /******** here's the code for A *************/
    ii = -1;
    j = 0;
    k = 0;
    for (i=0; i<NCellElems; i++) {
        ii += Cm[i];
        t = lp->x[ii];
        normu2 = 0;
        for (kk=0; kk<Cm[i]-1; kk++) { 
            u = lp->x[j];
            normu2 += u*u;
            A[k] = -2*u/t;
            k += 2;
            j++;
        }
        A[k] = 1 + normu2/(t*t);
        k += 2;
        j++;
    }

    /******** here's the code for At *************/
    ii = -1;
    j = 0;
    k = 0;
    for (i=0; i<NCellElems; i++) {
        ii += Cm[i];
        t = lp->x[ii];
        normu2 = 0;
        for (kk=0; kk<Cm[i]-1; kk++) { 
            u = lp->x[j];
            normu2 += u*u;
            At[k] = -2*u/t;
            k++;
            j++;
        }
        At[k] = 1 + normu2/(t*t);
        k++;
        j++;
    }
}

int rd_opns() {
	char *s = opnstr;
	char *s1;
	keyword *kw;

	while (*s) {
		while (*s && *s==' ') s++;
		if (!*s) return 1;
		if ((*s >= 'a' && *s <= 'z') || (*s >= 'A' && *s <= 'Z')) {
			kw = (keyword *)binsearch(&s);
			if (bad_opns) return 0;
			while (*s && *s == ' ') s++;
			switch (kw->type) {
				case 0:
					kw->ivalue = 1; // depending on keyword, this may need to change
					break;
				case 1:
					if (*s == '=') {
						s++;
						while (*s && *s == ' ') s++;
						if (*s <= '9' && *s >= '0') {
							kw->ivalue = (int)strtol(s1 = s, &s, 10);
						} else {
							return 0;
						}
					} else {
						return 0;
					}
					break;
				case 2:
                    if (*s == '=') {
                        s++;
                        while (*s && *s == ' ') s++;
                        if (*s <= '9' && *s >= '0') {
                            kw->dvalue = strtod(s1 = s, &s);
                        } else {
                            return 0;
                        }
                    } else {
                        return 0;
                    }
                    break;
			}
		} else {
			return 0;
		}
	}
					
	return 1;
}

void *binsearch(char **sp) {
	keyword *kw = keywds;
	keyword *kw1;
	int kwsize = (int)sizeof(keyword);
	int n = 30; // number of keywords
	int n1;
	int c, c1, c2;
	char *s, *s1, *s2;

	s = *sp;
	while (n > 0) {
		kw1 = kw + (n1 = n >> 1);
		s2 = *(char **)kw1;
		for (s1 = s;; s1++) {
			c1 = tolower(*(unsigned char *)s1);
			if (!(c2 = *s2++)) {
				if (c1 <= ' ' || c1 == '=') {
					*sp = s1;
					return kw1;
				}
				break;
			}
			if (c1 != c2) break;
		}
		if (c1 == '=' || c1 < c2) n = n1;
		else {
			n -= n1 + 1;
			kw = kw1 + 1;
		}
	}
	bad_opns++;
	return 0;
}

void set_opns(LOQO *lp) {
	keyword *kw = keywds;
	int verbose, itnlim;
	int noreord, md, mlf;
	int primal, dual;
	int max, maximize, min, minimize;

	lp->bndpush = kw->dvalue; kw++;
	lp->convex = kw->ivalue; kw++;
	lp->dense = kw->ivalue; kw++;
	dual = kw->ivalue; kw++;
	lp->epsdiag = kw->dvalue; kw++;
	lp->epsnum = kw->dvalue; kw++;
	lp->epssol = kw->dvalue; kw++;
	lp->honor_bnds = kw->ivalue; kw++;
	lp->honor_bnds_init = kw->ivalue; kw++;
	lp->inftol = kw->dvalue; kw++;
	lp->inftol2 = kw->dvalue; kw++;
	itnlim = kw->ivalue; kw++;
	max = kw->ivalue; kw++;
	maximize = kw->ivalue; kw++;
	if (itnlim == 200) lp->itnlim = kw->ivalue; else lp->itnlim = itnlim; kw++;
	min = kw->ivalue; kw++;
	md = kw->ivalue; kw++;
	minimize = kw->ivalue; kw++;
	mlf = kw->ivalue; kw++;
	lp->mufactor = kw->dvalue; kw++;
	noreord = kw->ivalue; kw++;
	verbose = kw->ivalue; kw++;
	lp->pred_corr = kw->ivalue; kw++;
	primal = kw->ivalue; kw++;
	lp->sdp = kw->ivalue; kw++;
	lp->sf_req = kw->ivalue; kw++;
	lp->stablty = kw->dvalue; kw++;
	lp->steplen = kw->dvalue; kw++;
	if (kw->dvalue != -1.0) lp->timlim = kw->dvalue; kw++;
	if (verbose == 0) lp->verbose = kw->ivalue; else lp->verbose = verbose; kw++;
	lp->method=1;	
	if (noreord) lp->method = 0;
	else if (mlf) lp->method = 2;
		else if (md) lp->method = 1;

	lp->pdf=0;
	if (primal) lp->pdf = 1;
	if (dual) lp->pdf = 2;

	lp->max=1;
	if (max || maximize) lp->max = -1;
	if (min || minimize) lp->max = 1;
}

void mexFunction(
    int nlhs, mxArray *plhs[],
    int nrhs, const mxArray *prhs[]
)
{
    int i, j, k, jj;
    char how[40] = "Not Available";
    double xinit;
    unsigned int status;
    int m_done=0, n_done=0;
    static char *str[] = {
        "OPTIMAL SOLUTION FOUND\n",
        "SUBOPTIMAL SOLUTION FOUND\n",
        "ITERATION LIMIT\n",
        "PRIMAL INFEASIBLE\n",
        "DUAL INFEASIBLE\n",
        "PRIMAL and/or DUAL INFEASIBLE\n",
        "PRIMAL and/or DUAL INFEASIBLE (INCONSISTENT EQUATIONS)\n",
        "PRIMAL UNBOUNDED\n",
        "DUAL UNBOUNDED\n",
        "TIME LIMIT\n"
    };

    Cm = NULL; Cn = NULL; C = NULL; iC = NULL; kC = NULL; c = NULL;
    A = NULL; iA = NULL; kA = NULL; b = NULL; Q = NULL; iQ = NULL;
    kQ = NULL; l = NULL; u = NULL; x0 = NULL; x = NULL; y_lin = NULL; y_cone = NULL;
    neq = 0; m = 0; n = 0; 

    if (nrhs > 11 || nrhs < 1) {
            mexErrMsgTxt("Usage: [x,lambda_lin,lambda_cone,how] "
            	" = loqosoclp(C,H,c,A,b,l,u,x0,neqcstr,hookfn,loqo_options)");
            return;
    }
    switch (nrhs) {
	case 11:
			if (mxIsChar(prhs[10]) && mxGetM(prhs[10])==1) {
                 int buflen = mxGetN(prhs[10])+1;
                 MALLOC(opnstr, buflen, char);
				 mxGetString(prhs[10], opnstr, buflen);
            } else {
                 mexErrMsgTxt("Eleventh argument (loqo_options) must be "
								"a string.");
				 return;
			}
    case 10:
            if (mxIsChar(prhs[9]) && mxGetM(prhs[9])==1) {
                 int buflen = mxGetN(prhs[9])+1;
                 MALLOC(cmdstr, buflen, char);
                 mxGetString(prhs[9], cmdstr, buflen);
		 /* crappy matlab doesn't understand that empty strings are
		  * indeed strings.  Therefore we kludge with a space */
		 if (strcmp(cmdstr," ") == 0) cmdstr=NULL;
            } else {
                 mexErrMsgTxt("Tenth argument (hookfn) must be "
                              "a string.");
                 return;
            }
    case 9:
            if (mxGetM(prhs[8]) != 0 || mxGetN(prhs[8]) != 0) {
                    if (!mxIsNumeric(prhs[8]) || mxIsComplex(prhs[8]) 
                     ||  mxIsSparse(prhs[8])
                     || !(mxGetM(prhs[8])==1 && mxGetN(prhs[8])==1)) {
                         mexErrMsgTxt("Ninth argument (neqcstr) must be "
                                      "an integer scalar.");
                         return;
                    }
                    neq = *mxGetPr(prhs[8]);
            }
    case 8:
            if (mxGetM(prhs[7]) != 0 || mxGetN(prhs[7]) != 0) {
                if (!mxIsNumeric(prhs[7]) || mxIsComplex(prhs[7]) 
                 ||  mxIsSparse(prhs[7])
                 || !mxIsDouble(prhs[7]) 
                 ||  mxGetN(prhs[7])!=1 ) {
                     mexErrMsgTxt("Eighth argument (x0) must be "
                                  "a column vector.");
                     return;
                }
                x0 = mxGetPr(prhs[7]);
                n  = mxGetM(prhs[7]);
                n_done = 1;
                xk=(mxArray*)prhs[7]; /*space for intermediate result */
            } 
    case 7:
            if (mxGetM(prhs[6]) != 0 || mxGetN(prhs[6]) != 0) {
                    if (!mxIsNumeric(prhs[6]) || mxIsComplex(prhs[6]) 
                     ||  mxIsSparse(prhs[6])
                     || !mxIsDouble(prhs[6]) 
                     ||  mxGetN(prhs[6])!=1 ) {
                         mexErrMsgTxt("Seventh argument (u) must be "
                                      "a column vector.");
                         return;
                    }
                    if (n != 0 && n != mxGetM(prhs[6])) {
                         mexErrMsgTxt("Dimension error (arg 7 and later).");
                         return;
                    }
                    u = mxGetPr(prhs[6]);
                    n = mxGetM(prhs[6]);
                    if (!n_done) {
                            n_done = 1;
                    }
            }
    case 6:
            if (mxGetM(prhs[5]) != 0 || mxGetN(prhs[5]) != 0) {
                    if (!mxIsNumeric(prhs[5]) || mxIsComplex(prhs[5]) 
                     ||  mxIsSparse(prhs[5])
                                       || !mxIsDouble(prhs[5]) 
                     ||  mxGetN(prhs[5])!=1 ) {
                         mexErrMsgTxt("Sixth argument (l) must be "
                                      "a column vector.");
                         return;
                    }
                    if (n != 0 && n != mxGetM(prhs[5])) {
                         mexErrMsgTxt("Dimension error (arg 6 and later).");
                         return;
                    }
                    l = mxGetPr(prhs[5]);
                    n = mxGetM(prhs[5]);
                    if (!n_done) {
                            n_done = 1;
                    }
            }
    case 5:
            if (mxGetM(prhs[4]) != 0 || mxGetN(prhs[4]) != 0) {
                    if (!mxIsNumeric(prhs[4]) || mxIsComplex(prhs[4]) 
                     ||  mxIsSparse(prhs[4])
                     || !mxIsDouble(prhs[4]) 
                     ||  mxGetN(prhs[4])!=1 ) {
                         mexErrMsgTxt("Fifth argument (b) must be "
                                      "a column vector.");
                         return;
                    }
                    if (m != 0 && m != mxGetM(prhs[4])) {
                         mexErrMsgTxt("Dimension error (arg 5 and later).");
                         return;
                    }
                    b = mxGetPr(prhs[4]);
                    m = mxGetM(prhs[4]);
                    m_done = 1;
            }
    case 4:
            if (mxGetM(prhs[3]) != 0 || mxGetN(prhs[3]) != 0) {
                    if (!mxIsNumeric(prhs[3]) || mxIsComplex(prhs[3]) 
                     || !mxIsSparse(prhs[3]) ) {
                         mexErrMsgTxt("Fourth argument (A) must be "
                                      "a sparse matrix.");
                         return;
                    }
                    if (m != 0 && m != mxGetM(prhs[3])) {
                         mexErrMsgTxt("Dimension error (arg 4 and later).");
                         return;
                    }
                    if (n != 0 && n != mxGetN(prhs[3])) {
                         mexErrMsgTxt("Dimension error (arg 4 and later).");
                         return;
                    }
                    m = mxGetM(prhs[3]);
                    n = mxGetN(prhs[3]);
                    if (!m_done) {
                            m_done = 1;
                    }
                    if (!n_done) {
                            n_done = 1;
                    }
                    A = mxGetPr(prhs[3]);
                    iA = mxGetIr(prhs[3]);
                    kA = mxGetJc(prhs[3]);

            }
    case 3:
            if (mxGetM(prhs[2]) != 0 || mxGetN(prhs[2]) != 0) {
                    if (!mxIsNumeric(prhs[2]) || mxIsComplex(prhs[2]) 
                     ||  mxIsSparse(prhs[2])
                     || !mxIsDouble(prhs[2]) 
                     ||  mxGetN(prhs[2])!=1 ) {
                         mexErrMsgTxt("Second argument (c) must be "
                                      "a column vector.");
                         return;
                    }
                    if (n != 0 && n != mxGetM(prhs[2])) {
                         mexErrMsgTxt("Dimension error (arg 3 and later).");
                         return;
                    }
                    c = mxGetPr(prhs[2]);
                    n = mxGetM(prhs[2]);
                    if (!n_done) {
                            n_done = 1;
                    }
            }
    case 2:
            if (mxGetM(prhs[1]) != 0 || mxGetN(prhs[1]) != 0) {
                    if (!mxIsNumeric(prhs[1]) || mxIsComplex(prhs[1]) 
                     || !mxIsSparse(prhs[1]) ) {
                         mexErrMsgTxt("First argument (H) must be "
                                      "a sparse matrix.");
                         return;
                    }
                    if (n != 0 && n != mxGetM(prhs[1])) {
                         mexErrMsgTxt("Dimension error (arg 2 and later).");
                                     return;
                    }
                    if (n != 0 && n != mxGetN(prhs[1])) {
                         mexErrMsgTxt("Dimension error (arg 2 and later).");
                         return;
                    }
                    n  = mxGetN(prhs[1]);
                    if (!n_done) {
                            n_done = 1;
                    }
                    Q = mxGetPr(prhs[1]);
                    iQ = mxGetIr(prhs[1]);
                    kQ = mxGetJc(prhs[1]);
            }
    case 1:
            if (!mxIsCell(prhs[0]))
            {
                 mexErrMsgTxt("First argument (C) must be "
                              "a cell array.");
                 return;
            }
            NCellElems = mxGetNumberOfElements(prhs[0]); 
            MALLOC(Cm, NCellElems, int);
            MALLOC(Cn, NCellElems, int);
            MALLOC(C, NCellElems, double *);
            MALLOC(iC, NCellElems, int *);
            MALLOC(kC, NCellElems, int *);
            for (i=0; i<NCellElems; i++) {
                mxArray *a = mxGetCell(prhs[0],i); 
                Cm[i] = mxGetM(a);
                Cn[i] = mxGetN(a);
                C[i] = mxGetPr(a);
                iC[i] = mxGetIr(a);
                kC[i] = mxGetJc(a);
            }
            break;
    }

    if (nlhs > 4 || nlhs < 1) {
            mexErrMsgTxt("Usage: [x,lambda_lin,lambda_cone,how] "
                         "= loqosoclp(C,H,c,A,b,l,u,x0,neqcstr,hookfn,loqo_options)");
            return;
    }

    switch (nlhs) {
	case 4:
    case 3:
			plhs[2] = mxCreateDoubleMatrix(NCellElems, 1, mxREAL);
            y_cone = mxGetPr(plhs[2]);
    case 2:
            plhs[1] = mxCreateDoubleMatrix(m, 1, mxREAL);
            y_lin = mxGetPr(plhs[1]);
    case 1:
            plhs[0] = mxCreateDoubleMatrix(n, 1, mxREAL);
            x = mxGetPr(plhs[0]);
            break;
    }

    status = mloqo();

    switch (nlhs) {
    case 4:
            plhs[3] = mxCreateString(str[status]);
    }

        mexAtExit(cleanup);
}

void cleanup(void) {
    mexPrintf("MEX-file is terminating, destroying array\n");
}