docs: def solve(self,problem):
          """
          Solves optimization problem.
          
          Parameters
          ----------
          problem : Object
          """
      
          # Local vars
          norm2 = self.norm2
          norminf = self.norminf
          parameters = self.parameters
          
          # Parameters
          feastol = parameters['feastol']
          optol = parameters['optol']
          maxiter = parameters['maxiter']
          quiet = parameters['quiet']
          sigma = parameters['sigma']
          eps = parameters['eps']
          ls_maxiter = parameters['line_search_maxiter']
  
          # Problem
          problem = cast_problem(problem)
          self.problem = 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)
  
          # Constants
          self.A = problem.A
          self.AT = problem.A.T
          self.b = problem.b
          self.u = problem.u+feastol/10.
          self.l = problem.l-feastol/10.
          self.n = problem.get_num_primal_variables()
          self.m1 = problem.get_num_linear_equality_constraints()
          self.m2 = problem.get_num_nonlinear_equality_constraints()
          self.e = np.ones(self.n)
          self.I = eye(self.n,format='coo')
          self.Omm1 = coo_matrix((self.m1,self.m1))
          self.Omm2 = coo_matrix((self.m2,self.m2))
  
          # Initial primal
          if problem.x is None:
              self.x = (self.u + self.l)/2.
          else:
              self.x = np.maximum(np.minimum(problem.x,problem.u),problem.l)
  
          # Initial duals
          if problem.lam is None:
              self.lam = np.zeros(problem.get_num_linear_equality_constraints())
          else:
              self.lam = problem.lam.copy()
          if problem.nu is None:
              self.nu = np.zeros(problem.get_num_nonlinear_equality_constraints())
          else:
              self.nu = problem.nu.copy()
          self.mu = np.minimum(1./(self.u-self.x), 1.)
          self.pi = np.minimum(1./(self.x-self.l), 1.)
  
          # Init vector
          self.y = np.hstack((self.x,self.lam,self.nu,self.mu,self.pi))
  
          # Average violation of complementarity slackness
          self.eta_mu = (np.dot(self.mu,self.u-self.x)/self.x.size) if self.x.size else 0.
          self.eta_pi = (np.dot(self.pi,self.x-self.l)/self.x.size) if self.x.size else 0.
  
          # Objective scaling
          fdata = self.func(self.y)
          self.obj_sca = np.maximum(norminf(problem.gphi)/10.,1.)
          fdata = self.func(self.y)
  
          # Header
          if not quiet:
              print('\nSolver: inlp')
              print('------------')
                                     
          # Outer
          s = 0.
          self.k = 0
          while True:
  
              # Average violation of complementarity slackness
              self.eta_mu = (np.dot(self.mu,self.u-self.x)/self.x.size) if self.x.size else 0.
              self.eta_pi = (np.dot(self.pi,self.x-self.l)/self.x.size) if self.x.size else 0.
              
              # Init eval
              fdata = self.func(self.y)
  
              # 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('pres'),end=' ')
                  print('{0:^9s}'.format('dres'),end=' ')
                  print('{0:^9s}'.format('gmax'),end=' ')
                  print('{0:^8s}'.format('cu'),end=' ')
                  print('{0:^8s}'.format('cl'),end=' ')
                  print('{0:^8s}'.format('alpha'))
   
              # Inner
              while True:
                  
                  # Eval
                  fdata = self.func(self.y)
                  pres = norminf(np.hstack((fdata.rp1,fdata.rp2)))
                  dres = norminf(np.hstack((fdata.rd,fdata.ru,fdata.rl)))
                  gmax = norminf(fdata.GradF)                
                  compu = norminf(self.mu*(self.u-self.x))
                  compl = norminf(self.pi*(self.x-self.l))
                  phi = problem.phi
                  
                  # Show progress
                  if not quiet:
                      print('{0:^3d}'.format(self.k),end=' ')
                      print('{0:^9.2e}'.format(phi),end=' ')
                      print('{0:^9.2e}'.format(pres),end=' ')
                      print('{0:^9.2e}'.format(dres),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 self.k > 0 and pres < feastol and dres < optol and sigma*np.maximum(self.eta_mu,self.eta_pi) < optol:
                      self.set_status(self.STATUS_SOLVED)
                      self.set_error_msg('')
                      return
                      
                  # Done
                  if gmax < tau: 
                      break
  
                  # Done 
                  if pres < feastol and dres < optol and np.maximum(compu,compl) < optol:
                      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.rp1,fdata.rp2))
                  Hbar = coo_matrix((np.concatenate((problem.Hphi.data/self.obj_sca,
                                                     problem.H_combined.data,
                                                     D1.data,
                                                     D2.data)),
                                     (np.concatenate((problem.Hphi.row,
                                                      problem.H_combined.row,
                                                      D1.row,
                                                      D2.row)),
                                      np.concatenate((problem.Hphi.col,
                                                      problem.H_combined.col,
                                                      D1.col,
                                                      D2.col)))))
                  Jbar = bmat([[Hbar,None,None],
                               [-self.A,self.Omm1,None],
                               [-problem.J,None,self.Omm2]],
                              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(((self.u-self.x)[indices]/px[indices],np.inf)))
                  indices = px < 0
                  s2 = np.min(np.hstack(((self.l-self.x)[indices]/px[indices],np.inf)))
                  indices = pmu < 0
                  s3 = np.min(np.hstack((-self.mu[indices]/pmu[indices],np.inf)))
                  indices = ppi < 0
                  s4 = np.min(np.hstack((-self.pi[indices]/ppi[indices],np.inf)))
                  smax = (1.-eps)*np.min([s1,s2,s3,s4])
                  spmax = (1.-eps)*np.min([s1,s2])
                  sdmax = (1.-eps)*np.min([s3,s4])
                  
                  # Line search
                  try:
                      s, fdata = self.line_search(self.y, p, fdata.F, fdata.GradF, self.func, smax=smax, maxiter=ls_maxiter)
  
                      # Update point
                      self.y += s*p
                      self.x, self.lam, self.nu, self.mu, self.pi = self.extract_components(self.y)
                      
                  except OptSolverError_LineSearch:
                      sp = np.minimum(1., spmax)
                      sd = np.minimum(1., sdmax)
                      s = np.minimum(sp,sd)
  
                      # Update point
                      self.x += sp*px
                      self.lam += sd*pbar[self.x.size:self.x.size+self.lam.size]
                      self.nu += sd*pbar[self.x.size+self.lam.size:]
                      self.mu += sd*pmu
                      self.pi += sd*ppi
                      self.y = np.hstack((self.x,self.lam,self.nu,self.mu,self.pi))
  
                  # Update iters
                  self.k += 1                
  
                  # 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)
  
              # Update iters
              self.k += 1