def solve(self,problem):
        """
        Solves optimization problem.

        Parameters
        ----------
        problem : Object
        """
    
        # Local vars
        norm2 = self.norm2
        norminf = self.norminf
        parameters = self.parameters
        
        # Parameters
        tol = parameters['tol']
        maxiter = parameters['maxiter']
        quiet = parameters['quiet']
        sigma = parameters['sigma']
        eps = parameters['eps']
        eps_cold = parameters['eps_cold']

        # Problem
        if not isinstance(problem,QuadProblem):
            problem = cast_problem(problem)
            quad_problem = QuadProblem(None,None,None,None,None,None,problem=problem)
        else:
            quad_problem = problem
        self.problem = problem
        self.quad_problem = quad_problem

        # Linsolver
        self.linsolver = new_linsolver(parameters['linsolver'],'symmetric')

        # Reset
        self.reset()

        # Checks
        if not np.all(problem.l <= problem.u):
            raise OptSolverError_NoInterior(self)
    
        # Data
        self.H = quad_problem.H
        self.g = quad_problem.g
        self.A = quad_problem.A
        self.AT = quad_problem.A.T
        self.b = quad_problem.b
        self.l = quad_problem.l-tol/10.
        self.u = quad_problem.u+tol/10.
        self.n = quad_problem.H.shape[0]
        self.m = quad_problem.A.shape[0]
        self.e = np.ones(self.n)
        self.I = eye(self.n,format='coo')
        self.Onm = coo_matrix((self.n,self.m))
        self.Omm = coo_matrix((self.m,self.m))

        # Initial primal
        if quad_problem.x is None:
            self.x = (self.u + self.l)/2.
        else:
            self.x = np.maximum(np.minimum(quad_problem.x,problem.u),problem.l)

        # Initial duals
        if quad_problem.lam is None:
            self.lam = np.zeros(self.m)
        else:
            self.lam = quad_problem.lam.copy()
        if quad_problem.mu is None:
            self.mu = np.ones(self.x.size)*eps_cold
        else:
            self.mu = np.maximum(quad_problem.mu,eps)
        if quad_problem.pi is None:
            self.pi = np.ones(self.x.size)*eps_cold
        else:
            self.pi = np.maximum(quad_problem.pi,eps)

        # Check interior
        try:
            assert(np.all(self.l < self.x)) 
            assert(np.all(self.x < self.u))
            assert(np.all(self.mu > 0))
            assert(np.all(self.pi > 0))
        except AssertionError:
            raise OptSolverError_Infeasibility(self)

        # Init vector
        self.y = np.hstack((self.x,self.lam,self.mu,self.pi))

        # Complementarity measures
        self.eta_mu = np.dot(self.mu,self.u-self.x)/self.x.size
        self.eta_pi = np.dot(self.pi,self.x-self.l)/self.x.size

        # Objective scaling
        fdata = self.func(self.y)
        self.obj_sca = np.maximum(norminf(self.g+self.H*self.x)/10.,1.)
        self.H = self.H/self.obj_sca
        self.g = self.g/self.obj_sca
        fdata = self.func(self.y)

        # Header
        if not quiet:
            print('\nSolver: IQP')
            print('-----------')
                                   
        # Outer
        s = 0.
        self.k = 0
        while True:

            # Complementarity measures
            self.eta_mu = np.dot(self.mu,self.u-self.x)/self.x.size
            self.eta_pi = np.dot(self.pi,self.x-self.l)/self.x.size
            
            # Init eval
            fdata = self.func(self.y)
            fmax = norminf(fdata.f)
            gmax = norminf(fdata.GradF)
            
            # Done
            if fmax < tol and sigma*np.maximum(self.eta_mu,self.eta_pi) < tol:
                self.set_status(self.STATUS_SOLVED)
                self.set_error_msg('')
                return

            # Target
            tau = sigma*norminf(fdata.GradF)
           
            # Header
            if not quiet:
                if self.k > 0:
                    print('')
                print('{0:^3s}'.format('iter'), end=' ')
                print('{0:^9s}'.format('phi'), end=' ')
                print('{0:^9s}'.format('fmax'), end=' ')
                print('{0:^9s}'.format('gmax'), end=' ')
                print('{0:^8s}'.format('cu'), end=' ')
                print('{0:^8s}'.format('cl'), end=' ')
                print('{0:^8s}'.format('s'))
 
            # Inner
            while True:
                
                # Eval
                fdata = self.func(self.y)
                fmax = norminf(fdata.f)
                gmax = norminf(fdata.GradF)
                compu = norminf(self.mu*(self.u-self.x))
                compl = norminf(self.pi*(self.x-self.l))
                phi = (0.5*np.dot(self.x,self.H*self.x)+np.dot(self.g,self.x))*self.obj_sca
                
                # Show progress
                if not quiet:
                    print('{0:^3d}'.format(self.k), end=' ')
                    print('{0:^9.2e}'.format(phi), end=' ')
                    print('{0:^9.2e}'.format(fmax), end=' ')
                    print('{0:^9.2e}'.format(gmax), end=' ')
                    print('{0:^8.1e}'.format(compu), end=' ')
                    print('{0:^8.1e}'.format(compl), end=' ')
                    print('{0:^8.1e}'.format(s))
                
                # Done
                if gmax < tau:
                    break

                # Done
                if fmax < tol and np.maximum(compu,compl) < tol:
                    break

                # Maxiters
                if self.k >= maxiter:
                    raise OptSolverError_MaxIters(self)
                    
                # Search direction
                ux = self.u-self.x
                xl = self.x-self.l
                D1 = spdiags(self.mu/ux,0,self.n,self.n,format='coo')
                D2 = spdiags(self.pi/xl,0,self.n,self.n,format='coo')
                fbar = np.hstack((-fdata.rd+fdata.ru/ux-fdata.rl/xl,fdata.rp))
                if self.A.shape[0] > 0:
                    Jbar = bmat([[tril(self.H)+D1+D2,None],
                                 [-self.A,self.Omm]],format='coo')
                else:
                    Jbar = bmat([[tril(self.H)+D1+D2]],
                                format='coo')
                try:
                    if not self.linsolver.is_analyzed():
                        self.linsolver.analyze(Jbar)
                    pbar = self.linsolver.factorize_and_solve(Jbar,fbar)
                except RuntimeError:
                    raise OptSolverError_BadLinSystem(self)
                px = pbar[:self.n]
                pmu = (-fdata.ru + self.mu*px)/ux
                ppi = (-fdata.rl - self.pi*px)/xl
                p = np.hstack((pbar,pmu,ppi))

                # Steplength bounds
                indices = px > 0
                s1 = np.min(np.hstack(((1.-eps)*(self.u-self.x)[indices]/px[indices],np.inf)))
                indices = px < 0
                s2 = np.min(np.hstack(((eps-1.)*(self.x-self.l)[indices]/px[indices],np.inf)))
                indices = pmu < 0
                s3 = np.min(np.hstack(((eps-1.)*self.mu[indices]/pmu[indices],np.inf)))
                indices = ppi < 0
                s4 = np.min(np.hstack(((eps-1.)*self.pi[indices]/ppi[indices],np.inf)))
                smax = np.min([s1,s2,s3,s4])
                
                # Line search
                s,fdata = self.line_search(self.y,p,fdata.F,fdata.GradF,self.func,smax)

                # Update x
                self.y += s*p
                self.k += 1
                self.x,self.lam,self.mu,self.pi = self.extract_components(self.y)

                # Check
                try:
                    assert(np.all(self.x < self.u))
                    assert(np.all(self.x > self.l))
                    assert(np.all(self.mu > 0))
                    assert(np.all(self.pi > 0))
                except AssertionError:
                    raise OptSolverError_Infeasibility(self)
Solves optimization problem.

        Parameters
        ----------
        problem : Object