/*************************************************************************
        Having feasible current point XC and possibly infeasible candidate   point
        XN,  this  function  performs  longest  step  from  XC to XN which retains
        feasibility. In case XN is found to be infeasible, at least one constraint
        is activated.

        For example, if we have:
          XC=0.5
          XN=1.2
          x>=0, x<=1
        then this function will move us to X=1.0 and activate constraint "x<=1".

        INPUT PARAMETERS:
            State   -   MinQP state.
            XC      -   current point, must be feasible with respect to
                        all constraints
            XN      -   candidate point, can be infeasible with respect to some
                        constraints. Must be located in the subspace of current
                        active set, i.e. it is feasible with respect to already
                        active constraints.
            Buf     -   temporary buffer, automatically resized if needed

        OUTPUT PARAMETERS:
            State   -   this function changes following fields of State:
                        * State.ActiveSet
                        * State.ActiveC     -   active linear constraints
            XC      -   new position

        RESULT:
            >0, in case at least one inactive non-candidate constraint was activated
            =0, in case only "candidate" constraints were activated
            <0, in case no constraints were activated by the step


          -- ALGLIB --
             Copyright 29.02.2012 by Bochkanov Sergey
        *************************************************************************/
        private static int minqpboundedstepandactivation(minqpstate state,
            double[] xn,
            ref double[] buf)
        {
            int result = 0;
            int n = 0;
            double stpmax = 0;
            int cidx = 0;
            double cval = 0;
            bool needact = new bool();
            double v = 0;
            int i_ = 0;

            n = state.n;
            apserv.rvectorsetlengthatleast(ref buf, n);
            for(i_=0; i_<=n-1;i_++)
            {
                buf[i_] = xn[i_];
            }
            for(i_=0; i_<=n-1;i_++)
            {
                buf[i_] = buf[i_] - state.sas.xc[i_];
            }
            sactivesets.sasexploredirection(state.sas, buf, ref stpmax, ref cidx, ref cval);
            needact = (double)(stpmax)<=(double)(1);
            v = Math.Min(stpmax, 1.0);
            for(i_=0; i_<=n-1;i_++)
            {
                buf[i_] = v*buf[i_];
            }
            for(i_=0; i_<=n-1;i_++)
            {
                buf[i_] = buf[i_] + state.sas.xc[i_];
            }
            result = sactivesets.sasmoveto(state.sas, buf, needact, cidx, cval);
            return result;
        }

        /*************************************************************************
        Fast version of MinQPSetQuadraticTerm(), which doesn't check its arguments.

        It accepts additional parameter - shift S, which allows to "shift"  matrix
        A by adding s*I to A. S must be positive (although it is not checked).

        For internal use only.

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void minqpsetquadratictermfast(minqpstate state,
            double[,] a,
            bool isupper,
            double s)
        {
            int i = 0;
            int j = 0;
            int n = 0;
            double v = 0;
            int j0 = 0;
            int j1 = 0;

            n = state.n;
            state.akind = 0;
            cqmodels.cqmseta(state.a, a, isupper, 1.0);
            if( (double)(s)>(double)(0) )
            {
                apserv.rvectorsetlengthatleast(ref state.tmp0, n);
                for(i=0; i<=n-1; i++)
                {
                    state.tmp0[i] = a[i,i]+s;
                }
                cqmodels.cqmrewritedensediagonal(state.a, state.tmp0);
            }
            
            //
            // Estimate norm of A
            // (it will be used later in the quadratic penalty function)
            //
            state.absamax = 0;
            state.absasum = 0;
            state.absasum2 = 0;
            for(i=0; i<=n-1; i++)
            {
                if( isupper )
                {
                    j0 = i;
                    j1 = n-1;
                }
                else
                {
                    j0 = 0;
                    j1 = i;
                }
                for(j=j0; j<=j1; j++)
                {
                    v = Math.Abs(a[i,j]);
                    state.absamax = Math.Max(state.absamax, v);
                    state.absasum = state.absasum+v;
                    state.absasum2 = state.absasum2+v*v;
                }
            }
        }

        /*************************************************************************
        Fast version of MinQPSetStartingPoint(), which doesn't check its arguments.
        For internal use only.

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void minqpsetstartingpointfast(minqpstate state,
            double[] x)
        {
            int n = 0;
            int i_ = 0;

            n = state.n;
            for(i_=0; i_<=n-1;i_++)
            {
                state.startx[i_] = x[i_];
            }
            state.havex = true;
        }

        /*************************************************************************
        This function solves quadratic programming problem.

        Prior to calling this function you should choose solver by means of one of
        the following functions:

        * MinQPSetAlgoQuickQP() - for QuickQP solver
        * MinQPSetAlgoBLEIC() - for BLEIC-QP solver

        These functions also allow you to control stopping criteria of the solver.
        If you did not set solver,  MinQP  subpackage  will  automatically  select
        solver for your problem and will run it with default stopping criteria.

        However, it is better to set explicitly solver and its stopping criteria.

        INPUT PARAMETERS:
            State   -   algorithm state

        You should use MinQPResults() function to access results after calls
        to this function.

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey.
             Special thanks to Elvira Illarionova  for  important  suggestions  on
             the linearly constrained QP algorithm.
        *************************************************************************/
        public static void minqpoptimize(minqpstate state)
        {
            int n = 0;
            int i = 0;
            int nbc = 0;
            int currentsolver = 0;

            n = state.n;
            state.repterminationtype = -5;
            state.repinneriterationscount = 0;
            state.repouteriterationscount = 0;
            state.repncholesky = 0;
            state.repnmv = 0;
            
            //
            // check correctness of constraints
            //
            for(i=0; i<=n-1; i++)
            {
                if( state.havebndl[i] && state.havebndu[i] )
                {
                    if( (double)(state.bndl[i])>(double)(state.bndu[i]) )
                    {
                        state.repterminationtype = -3;
                        return;
                    }
                }
            }
            
            //
            // count number of bound and linear constraints
            //
            nbc = 0;
            for(i=0; i<=n-1; i++)
            {
                if( state.havebndl[i] )
                {
                    nbc = nbc+1;
                }
                if( state.havebndu[i] )
                {
                    nbc = nbc+1;
                }
            }
            
            //
            // Initial point:
            // * if we have starting point in StartX, we just have to bound it
            // * if we do not have StartX, deduce initial point from boundary constraints
            //
            if( state.havex )
            {
                for(i=0; i<=n-1; i++)
                {
                    state.xs[i] = state.startx[i];
                    if( state.havebndl[i] && (double)(state.xs[i])<(double)(state.bndl[i]) )
                    {
                        state.xs[i] = state.bndl[i];
                    }
                    if( state.havebndu[i] && (double)(state.xs[i])>(double)(state.bndu[i]) )
                    {
                        state.xs[i] = state.bndu[i];
                    }
                }
            }
            else
            {
                for(i=0; i<=n-1; i++)
                {
                    if( state.havebndl[i] && state.havebndu[i] )
                    {
                        state.xs[i] = 0.5*(state.bndl[i]+state.bndu[i]);
                        continue;
                    }
//.........这里部分代码省略.........

        /*************************************************************************
        This function tells solver to use BLEIC-based algorithm and sets  stopping
        criteria for the algorithm.

        ALGORITHM FEATURES:

        * supports dense and sparse QP problems
        * supports boundary and general linear equality/inequality constraints
        * can solve all types of problems  (convex,  semidefinite,  nonconvex)  as
          long as they are bounded from below under constraints.
          Say, it is possible to solve "min{-x^2} subject to -1<=x<=+1".
          Of course, global  minimum  is found only  for  positive  definite   and
          semidefinite  problems.  As  for indefinite ones - only local minimum is
          found.

        ALGORITHM OUTLINE:

        * BLEIC-QP solver is just a driver function for MinBLEIC solver; it solves
          quadratic  programming   problem   as   general   linearly   constrained
          optimization problem, which is solved by means of BLEIC solver  (part of
          ALGLIB, active set method).
          
        ALGORITHM LIMITATIONS:

        * unlike QuickQP solver, this algorithm does not perform Newton steps  and
          does not use Level 3 BLAS. Being general-purpose active set  method,  it
          can activate constraints only one-by-one. Thus, its performance is lower
          than that of QuickQP.
        * its precision is also a bit  inferior  to  that  of   QuickQP.  BLEIC-QP
          performs only LBFGS steps (no Newton steps), which are good at detecting
          neighborhood of the solution, buy need many iterations to find  solution
          with more than 6 digits of precision.

        INPUT PARAMETERS:
            State   -   structure which stores algorithm state
            EpsG    -   >=0
                        The  subroutine  finishes  its  work   if   the  condition
                        |v|<EpsG is satisfied, where:
                        * |.| means Euclidian norm
                        * v - scaled constrained gradient vector, v[i]=g[i]*s[i]
                        * g - gradient
                        * s - scaling coefficients set by MinQPSetScale()
            EpsF    -   >=0
                        The  subroutine  finishes its work if exploratory steepest
                        descent  step  on  k+1-th iteration  satisfies   following
                        condition:  |F(k+1)-F(k)|<=EpsF*max{|F(k)|,|F(k+1)|,1}
            EpsX    -   >=0
                        The  subroutine  finishes its work if exploratory steepest
                        descent  step  on  k+1-th iteration  satisfies   following
                        condition:  
                        * |.| means Euclidian norm
                        * v - scaled step vector, v[i]=dx[i]/s[i]
                        * dx - step vector, dx=X(k+1)-X(k)
                        * s - scaling coefficients set by MinQPSetScale()
            MaxIts  -   maximum number of iterations. If MaxIts=0, the  number  of
                        iterations is unlimited. NOTE: this  algorithm uses  LBFGS
                        iterations,  which  are  relatively  cheap,  but   improve
                        function value only a bit. So you will need many iterations
                        to converge - from 0.1*N to 10*N, depending  on  problem's
                        condition number.

        IT IS VERY IMPORTANT TO CALL MinQPSetScale() WHEN YOU USE THIS  ALGORITHM
        BECAUSE ITS STOPPING CRITERIA ARE SCALE-DEPENDENT!

        Passing EpsG=0, EpsF=0 and EpsX=0 and MaxIts=0 (simultaneously) will lead
        to automatic stopping criterion selection (presently it is  small    step
        length, but it may change in the future versions of ALGLIB).

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void minqpsetalgobleic(minqpstate state,
            double epsg,
            double epsf,
            double epsx,
            int maxits)
        {
            alglib.ap.assert(math.isfinite(epsg), "MinQPSetAlgoBLEIC: EpsG is not finite number");
            alglib.ap.assert((double)(epsg)>=(double)(0), "MinQPSetAlgoBLEIC: negative EpsG");
            alglib.ap.assert(math.isfinite(epsf), "MinQPSetAlgoBLEIC: EpsF is not finite number");
            alglib.ap.assert((double)(epsf)>=(double)(0), "MinQPSetAlgoBLEIC: negative EpsF");
            alglib.ap.assert(math.isfinite(epsx), "MinQPSetAlgoBLEIC: EpsX is not finite number");
            alglib.ap.assert((double)(epsx)>=(double)(0), "MinQPSetAlgoBLEIC: negative EpsX");
            alglib.ap.assert(maxits>=0, "MinQPSetAlgoBLEIC: negative MaxIts!");
            state.algokind = 2;
            if( (((double)(epsg)==(double)(0) && (double)(epsf)==(double)(0)) && (double)(epsx)==(double)(0)) && maxits==0 )
            {
                epsx = 1.0E-6;
            }
            state.qpbleicsettingsuser.epsg = epsg;
            state.qpbleicsettingsuser.epsf = epsf;
            state.qpbleicsettingsuser.epsx = epsx;
            state.qpbleicsettingsuser.maxits = maxits;
        }

        /*************************************************************************
        This function sets boundary constraints for QP solver

        Boundary constraints are inactive by default (after initial creation).
        After  being  set,  they  are  preserved  until explicitly turned off with
        another SetBC() call.

        INPUT PARAMETERS:
            State   -   structure stores algorithm state
            BndL    -   lower bounds, array[N].
                        If some (all) variables are unbounded, you may specify
                        very small number or -INF (latter is recommended because
                        it will allow solver to use better algorithm).
            BndU    -   upper bounds, array[N].
                        If some (all) variables are unbounded, you may specify
                        very large number or +INF (latter is recommended because
                        it will allow solver to use better algorithm).
                        
        NOTE: it is possible to specify BndL[i]=BndU[i]. In this case I-th
        variable will be "frozen" at X[i]=BndL[i]=BndU[i].

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void minqpsetbc(minqpstate state,
            double[] bndl,
            double[] bndu)
        {
            int i = 0;
            int n = 0;

            n = state.n;
            alglib.ap.assert(alglib.ap.len(bndl)>=n, "MinQPSetBC: Length(BndL)<N");
            alglib.ap.assert(alglib.ap.len(bndu)>=n, "MinQPSetBC: Length(BndU)<N");
            for(i=0; i<=n-1; i++)
            {
                alglib.ap.assert(math.isfinite(bndl[i]) || Double.IsNegativeInfinity(bndl[i]), "MinQPSetBC: BndL contains NAN or +INF");
                alglib.ap.assert(math.isfinite(bndu[i]) || Double.IsPositiveInfinity(bndu[i]), "MinQPSetBC: BndU contains NAN or -INF");
                state.bndl[i] = bndl[i];
                state.havebndl[i] = math.isfinite(bndl[i]);
                state.bndu[i] = bndu[i];
                state.havebndu[i] = math.isfinite(bndu[i]);
            }
        }

        /*************************************************************************
        Fast version of MinQPSetQuadraticTerm(), which doesn't check its arguments.

        It accepts additional parameter - shift S, which allows to "shift"  matrix
        A by adding s*I to A. S must be positive (although it is not checked).

        For internal use only.

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void minqpsetquadratictermfast(minqpstate state,
            double[,] a,
            bool isupper,
            double s)
        {
            int k = 0;
            int n = 0;
            int i_ = 0;

            
            //
            // We store off-diagonal part of A in the lower triangle of DenseA.
            // Diagonal elements of A are stored in the DiagA.
            // Diagonal of DenseA and uppper triangle are used as temporaries.
            //
            // Why such complex storage? Because it:
            // 1. allows us to easily recover from exceptions (lower triangle
            //    is unmodified during execution as well as DiagA, and on entry
            //    we will always find unmodified matrix)
            // 2. allows us to make Cholesky decomposition in the upper triangle
            //    of DenseA or to do other SPD-related operations.
            //
            n = state.n;
            state.akind = 0;
            apserv.rmatrixsetlengthatleast(ref state.densea, n, n);
            apserv.rvectorsetlengthatleast(ref state.diaga, n);
            if( isupper )
            {
                for(k=0; k<=n-2; k++)
                {
                    state.diaga[k] = a[k,k]+s;
                    for(i_=k+1; i_<=n-1;i_++)
                    {
                        state.densea[i_,k] = a[k,i_];
                    }
                }
                state.diaga[n-1] = a[n-1,n-1]+s;
            }
            else
            {
                state.diaga[0] = a[0,0]+s;
                for(k=1; k<=n-1; k++)
                {
                    for(i_=0; i_<=k-1;i_++)
                    {
                        state.densea[k,i_] = a[k,i_];
                    }
                    state.diaga[k] = a[k,k]+s;
                }
            }
        }

        /*************************************************************************
                            CONSTRAINED QUADRATIC PROGRAMMING

        The subroutine creates QP optimizer. After initial creation,  it  contains
        default optimization problem with zero quadratic and linear terms  and  no
        constraints. You should set quadratic/linear terms with calls to functions
        provided by MinQP subpackage.

        INPUT PARAMETERS:
            N       -   problem size
            
        OUTPUT PARAMETERS:
            State   -   optimizer with zero quadratic/linear terms
                        and no constraints

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void minqpcreate(int n,
            minqpstate state)
        {
            int i = 0;

            ap.assert(n>=1, "MinQPCreate: N<1");
            
            //
            // initialize QP solver
            //
            state.n = n;
            state.akind = -1;
            state.repterminationtype = 0;
            state.b = new double[n];
            state.bndl = new double[n];
            state.bndu = new double[n];
            state.workbndl = new double[n];
            state.workbndu = new double[n];
            state.havebndl = new bool[n];
            state.havebndu = new bool[n];
            state.startx = new double[n];
            state.xorigin = new double[n];
            state.xc = new double[n];
            state.gc = new double[n];
            for(i=0; i<=n-1; i++)
            {
                state.b[i] = 0.0;
                state.workbndl[i] = Double.NegativeInfinity;
                state.workbndu[i] = Double.PositiveInfinity;
                state.havebndl[i] = false;
                state.havebndu[i] = false;
                state.startx[i] = 0.0;
                state.xorigin[i] = 0.0;
            }
            state.havex = false;
            minqpsetalgocholesky(state);
        }

        /*************************************************************************
        This function solves quadratic programming problem.
        You should call it after setting solver options with MinQPSet...() calls.

        INPUT PARAMETERS:
            State   -   algorithm state

        You should use MinQPResults() function to access results after calls
        to this function.

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void minqpoptimize(minqpstate state)
        {
            int n = 0;
            int i = 0;
            int j = 0;
            int k = 0;
            int nbc = 0;
            int nlc = 0;
            int nactive = 0;
            int nfree = 0;
            double f = 0;
            double fprev = 0;
            double v = 0;
            bool b = new bool();
            int i_ = 0;

            n = state.n;
            state.repterminationtype = -5;
            state.repinneriterationscount = 0;
            state.repouteriterationscount = 0;
            state.repncholesky = 0;
            state.repnmv = 0;
            
            //
            // check correctness of constraints
            //
            for(i=0; i<=n-1; i++)
            {
                if( state.havebndl[i] & state.havebndu[i] )
                {
                    if( (double)(state.bndl[i])>(double)(state.bndu[i]) )
                    {
                        state.repterminationtype = -3;
                        return;
                    }
                }
            }
            
            //
            // count number of bound and linear constraints
            //
            nbc = 0;
            nlc = 0;
            for(i=0; i<=n-1; i++)
            {
                if( state.havebndl[i] )
                {
                    nbc = nbc+1;
                }
                if( state.havebndu[i] )
                {
                    nbc = nbc+1;
                }
            }
            
            //
            // Our formulation of quadratic problem includes origin point,
            // i.e. we have F(x-x_origin) which is minimized subject to
            // constraints on x, instead of having simply F(x).
            //
            // Here we make transition from non-zero origin to zero one.
            // In order to make such transition we have to:
            // 1. subtract x_origin from x_start
            // 2. modify constraints
            // 3. solve problem
            // 4. add x_origin to solution
            //
            // There is alternate solution - to modify quadratic function
            // by expansion of multipliers containing (x-x_origin), but
            // we prefer to modify constraints, because it is a) more precise
            // and b) easier to to.
            //
            // Parts (1)-(2) are done here. After this block is over,
            // we have:
            // * XC, which stores shifted XStart (if we don't have XStart,
            //   value of XC will be ignored later)
            // * WorkBndL, WorkBndU, which store modified boundary constraints.
            //
            for(i=0; i<=n-1; i++)
            {
                state.xc[i] = state.startx[i]-state.xorigin[i];
                if( state.havebndl[i] )
                {
                    state.workbndl[i] = state.bndl[i]-state.xorigin[i];
                }
                if( state.havebndu[i] )
                {
//.........这里部分代码省略.........

        /*************************************************************************
        This  function  calculates x'*A*x for given X.

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey
        *************************************************************************/
        private static double minqpxtax(minqpstate state,
            double[] x)
        {
            double result = 0;
            int n = 0;
            int i = 0;
            int j = 0;

            n = state.n;
            ap.assert(state.akind==-1 | state.akind==0, "MinQPXTAX: internal error");
            result = 0;
            
            //
            // zero A
            //
            if( state.akind==-1 )
            {
                result = 0.0;
                return result;
            }
            
            //
            // dense A
            //
            if( state.akind==0 )
            {
                result = 0;
                for(i=0; i<=n-1; i++)
                {
                    for(j=0; j<=i-1; j++)
                    {
                        result = result+state.densea[i,j]*x[i]*x[j];
                    }
                    result = result+0.5*state.diaga[i]*math.sqr(x[i]);
                }
                return result;
            }
            return result;
        }

        /*************************************************************************
        This  function  calculates gradient of quadratic function at XC and stores
        it in the GC.

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey
        *************************************************************************/
        private static void minqpgrad(minqpstate state)
        {
            int n = 0;
            int i = 0;
            double v = 0;
            int i_ = 0;

            n = state.n;
            ap.assert(state.akind==-1 | state.akind==0, "MinQPGrad: internal error");
            
            //
            // zero A
            //
            if( state.akind==-1 )
            {
                for(i_=0; i_<=n-1;i_++)
                {
                    state.gc[i_] = state.b[i_];
                }
                return;
            }
            
            //
            // dense A
            //
            if( state.akind==0 )
            {
                for(i_=0; i_<=n-1;i_++)
                {
                    state.gc[i_] = state.b[i_];
                }
                state.gc[0] = state.gc[0]+state.diaga[0]*state.xc[0];
                for(i=1; i<=n-1; i++)
                {
                    v = 0.0;
                    for(i_=0; i_<=i-1;i_++)
                    {
                        v += state.densea[i,i_]*state.xc[i_];
                    }
                    state.gc[i] = state.gc[i]+v+state.diaga[i]*state.xc[i];
                    v = state.xc[i];
                    for(i_=0; i_<=i-1;i_++)
                    {
                        state.gc[i_] = state.gc[i_] + v*state.densea[i,i_];
                    }
                }
                return;
            }
        }

        /*************************************************************************
        Optimum of A subject to:
        a) active boundary constraints (given by ActiveSet[] and corresponding
           elements of XC)
        b) active linear constraints (given by C, R, LagrangeC)

        INPUT PARAMETERS:
            A       -   main quadratic term of the model;
                        although structure may  store  linear  and  rank-K  terms,
                        these terms are ignored and rewritten  by  this  function.
            ANorm   -   estimate of ||A|| (2-norm is used)
            B       -   array[N], linear term of the model
            XN      -   possibly preallocated buffer
            Tmp     -   temporary buffer (automatically resized)
            Tmp1    -   temporary buffer (automatically resized)

        OUTPUT PARAMETERS:
            A       -   modified quadratic model (this function changes rank-K
                        term and linear term of the model)
            LagrangeC-  current estimate of the Lagrange coefficients
            XN      -   solution

        RESULT:
            True on success, False on failure (non-SPD model)

          -- ALGLIB --
             Copyright 20.06.2012 by Bochkanov Sergey
        *************************************************************************/
        private static bool minqpconstrainedoptimum(minqpstate state,
            cqmodels.convexquadraticmodel a,
            double anorm,
            double[] b,
            ref double[] xn,
            ref double[] tmp,
            ref bool[] tmpb,
            ref double[] lagrangec)
        {
            bool result = new bool();
            int itidx = 0;
            int i = 0;
            double v = 0;
            double feaserrold = 0;
            double feaserrnew = 0;
            double theta = 0;
            int n = 0;
            int i_ = 0;

            n = state.n;
            
            //
            // Rebuild basis accroding to current active set.
            // We call SASRebuildBasis() to make sure that fields of SAS
            // store up to date values.
            //
            sactivesets.sasrebuildbasis(state.sas);
            
            //
            // Allocate temporaries.
            //
            apserv.rvectorsetlengthatleast(ref tmp, Math.Max(n, state.sas.basissize));
            apserv.bvectorsetlengthatleast(ref tmpb, n);
            apserv.rvectorsetlengthatleast(ref lagrangec, state.sas.basissize);
            
            //
            // Prepare model
            //
            for(i=0; i<=state.sas.basissize-1; i++)
            {
                tmp[i] = state.sas.pbasis[i,n];
            }
            theta = 100.0*anorm;
            for(i=0; i<=n-1; i++)
            {
                if( state.sas.activeset[i]>0 )
                {
                    tmpb[i] = true;
                }
                else
                {
                    tmpb[i] = false;
                }
            }
            cqmodels.cqmsetactiveset(a, state.sas.xc, tmpb);
            cqmodels.cqmsetq(a, state.sas.pbasis, tmp, state.sas.basissize, theta);
            
            //
            // Iterate until optimal values of Lagrange multipliers are found
            //
            for(i=0; i<=state.sas.basissize-1; i++)
            {
                lagrangec[i] = 0;
            }
            feaserrnew = math.maxrealnumber;
            result = true;
            for(itidx=1; itidx<=maxlagrangeits; itidx++)
            {
                
                //
                // Generate right part B using linear term and current
                // estimate of the Lagrange multipliers.
//.........这里部分代码省略.........

        /*************************************************************************
        This function sets scaling coefficients.

        ALGLIB optimizers use scaling matrices to test stopping  conditions  (step
        size and gradient are scaled before comparison with tolerances).  Scale of
        the I-th variable is a translation invariant measure of:
        a) "how large" the variable is
        b) how large the step should be to make significant changes in the function

        BLEIC-based QP solver uses scale for two purposes:
        * to evaluate stopping conditions
        * for preconditioning of the underlying BLEIC solver

        INPUT PARAMETERS:
            State   -   structure stores algorithm state
            S       -   array[N], non-zero scaling coefficients
                        S[i] may be negative, sign doesn't matter.

          -- ALGLIB --
             Copyright 14.01.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void minqpsetscale(minqpstate state,
            double[] s)
        {
            int i = 0;

            alglib.ap.assert(alglib.ap.len(s)>=state.n, "MinQPSetScale: Length(S)<N");
            for(i=0; i<=state.n-1; i++)
            {
                alglib.ap.assert(math.isfinite(s[i]), "MinQPSetScale: S contains infinite or NAN elements");
                alglib.ap.assert((double)(s[i])!=(double)(0), "MinQPSetScale: S contains zero elements");
                state.s[i] = Math.Abs(s[i]);
            }
        }

 public override alglib.apobject make_copy()
 {
     minqpstate _result = new minqpstate();
     _result.n = n;
     _result.qqpsettingsuser = (qqpsolver.qqpsettings)qqpsettingsuser.make_copy();
     _result.qqpsettingscurrent = (qqpsolver.qqpsettings)qqpsettingscurrent.make_copy();
     _result.qpbleicsettingsuser = (qpbleicsolver.qpbleicsettings)qpbleicsettingsuser.make_copy();
     _result.qpbleicsettingscurrent = (qpbleicsolver.qpbleicsettings)qpbleicsettingscurrent.make_copy();
     _result.algokind = algokind;
     _result.akind = akind;
     _result.a = (cqmodels.convexquadraticmodel)a.make_copy();
     _result.sparsea = (sparse.sparsematrix)sparsea.make_copy();
     _result.sparseaupper = sparseaupper;
     _result.absamax = absamax;
     _result.absasum = absasum;
     _result.absasum2 = absasum2;
     _result.b = (double[])b.Clone();
     _result.bndl = (double[])bndl.Clone();
     _result.bndu = (double[])bndu.Clone();
     _result.s = (double[])s.Clone();
     _result.havebndl = (bool[])havebndl.Clone();
     _result.havebndu = (bool[])havebndu.Clone();
     _result.xorigin = (double[])xorigin.Clone();
     _result.startx = (double[])startx.Clone();
     _result.havex = havex;
     _result.cleic = (double[,])cleic.Clone();
     _result.nec = nec;
     _result.nic = nic;
     _result.xs = (double[])xs.Clone();
     _result.repinneriterationscount = repinneriterationscount;
     _result.repouteriterationscount = repouteriterationscount;
     _result.repncholesky = repncholesky;
     _result.repnmv = repnmv;
     _result.repterminationtype = repterminationtype;
     _result.tmp0 = (double[])tmp0.Clone();
     _result.qpbleicfirstcall = qpbleicfirstcall;
     _result.qpbleicbuf = (qpbleicsolver.qpbleicbuffers)qpbleicbuf.make_copy();
     _result.qqpbuf = (qqpsolver.qqpbuffers)qqpbuf.make_copy();
     _result.qpcholeskybuf = (qpcholeskysolver.qpcholeskybuffers)qpcholeskybuf.make_copy();
     _result.estimator = (normestimator.normestimatorstate)estimator.make_copy();
     return _result;
 }

        /*************************************************************************
        This function tells solver to use Cholesky-based algorithm. This algorithm
        was deprecated in ALGLIB 3.9.0 because its performance is inferior to that
        of BLEIC-QP or  QuickQP  on  high-dimensional  problems.  Furthermore,  it
        supports only dense convex QP problems.

        This solver is no longer active by default.

        We recommend you to switch to BLEIC-QP or QuickQP solver.

        INPUT PARAMETERS:
            State   -   structure which stores algorithm state

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void minqpsetalgocholesky(minqpstate state)
        {
            state.algokind = 1;
        }

        /*************************************************************************
                            CONSTRAINED QUADRATIC PROGRAMMING

        The subroutine creates QP optimizer. After initial creation,  it  contains
        default optimization problem with zero quadratic and linear terms  and  no
        constraints. You should set quadratic/linear terms with calls to functions
        provided by MinQP subpackage.

        You should also choose appropriate QP solver and set it  and  its stopping
        criteria by means of MinQPSetAlgo??????() function. Then, you should start
        solution process by means of MinQPOptimize() call. Solution itself can  be
        obtained with MinQPResults() function.

        INPUT PARAMETERS:
            N       -   problem size
            
        OUTPUT PARAMETERS:
            State   -   optimizer with zero quadratic/linear terms
                        and no constraints

          -- ALGLIB --
             Copyright 11.01.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void minqpcreate(int n,
            minqpstate state)
        {
            int i = 0;

            alglib.ap.assert(n>=1, "MinQPCreate: N<1");
            
            //
            // initialize QP solver
            //
            state.n = n;
            state.nec = 0;
            state.nic = 0;
            state.repterminationtype = 0;
            state.absamax = 1;
            state.absasum = 1;
            state.absasum2 = 1;
            state.akind = 0;
            state.sparseaupper = false;
            cqmodels.cqminit(n, state.a);
            state.b = new double[n];
            state.bndl = new double[n];
            state.bndu = new double[n];
            state.havebndl = new bool[n];
            state.havebndu = new bool[n];
            state.s = new double[n];
            state.startx = new double[n];
            state.xorigin = new double[n];
            state.xs = new double[n];
            for(i=0; i<=n-1; i++)
            {
                state.bndl[i] = Double.NegativeInfinity;
                state.bndu[i] = Double.PositiveInfinity;
                state.havebndl[i] = false;
                state.havebndu[i] = false;
                state.b[i] = 0.0;
                state.startx[i] = 0.0;
                state.xorigin[i] = 0.0;
                state.s[i] = 1.0;
            }
            state.havex = false;
            minqpsetalgocholesky(state);
            normestimator.normestimatorcreate(n, n, 5, 5, state.estimator);
            qqpsolver.qqploaddefaults(n, state.qqpsettingsuser);
            qpbleicsolver.qpbleicloaddefaults(n, state.qpbleicsettingsuser);
            state.qpbleicfirstcall = true;
        }

//.........这里部分代码省略.........
          to perform several constrained Newton steps. If  Cholesky  decomposition
          failed (matrix is indefinite even under constraints),  we  perform  more
          CG iterations until we converge to such set of constraints  that  system
          matrix becomes  positive  definite.  Constrained  Newton  steps  greatly
          increase convergence speed and precision.
        * algorithm interleaves CG and Newton iterations which  allows  to  handle
          indefinite matrices (CG phase) and quickly converge after final  set  of
          constraints is found (Newton phase). Combination of CG and Newton phases
          is called "outer iteration".
        * it is possible to turn off Newton  phase  (beneficial  for  semidefinite
          problems - Cholesky decomposition will fail too often)
          
        ALGORITHM LIMITATIONS:

        * algorithm does not support general  linear  constraints;  only  boundary
          ones are supported
        * Cholesky decomposition for sparse problems  is  performed  with  Skyline
          Cholesky solver, which is intended for low-profile matrices. No profile-
          reducing reordering of variables is performed in this version of ALGLIB.
        * problems with near-zero negative eigenvalues (or exacty zero  ones)  may
          experience about 2-3x performance penalty. The reason is  that  Cholesky
          decomposition can not be performed until we identify directions of