/*************************************************************************
        This  subroutine returns L1 penalty for violation of active general linear
        constraints (violation of boundary or inactive linear constraints  is  not
        added to penalty).

        Penalty term is equal to:
            
            Penalty = SUM( Abs((C_i*x-R_i)/Alpha_i) )
            
        Here:
        * summation is performed for I=0...NEC+NIC-1, ActiveSet[N+I]>0
          (only for rows of CLEIC which are in active set)
        * C_i is I-th row of CLEIC
        * R_i is corresponding right part
        * S is a scale matrix
        * Alpha_i = ||S*C_i|| - is a scaling coefficient which "normalizes"
          I-th summation term according to its scale.

        INPUT PARAMETERS:
            S       -   active set object
            X       -   array[N], candidate point

          -- ALGLIB --
             Copyright 21.12.2012 by Bochkanov Sergey
        *************************************************************************/
        public static double sasactivelcpenalty1(sactiveset state,
            double[] x)
        {
            double result = 0;
            int i = 0;
            int j = 0;
            int n = 0;
            int nec = 0;
            int nic = 0;
            double v = 0;
            double alpha = 0;
            double p = 0;

            alglib.ap.assert(state.algostate==1, "SASActiveLCPenalty1: is not in optimization mode");
            sasrebuildbasis(state);
            n = state.n;
            nec = state.nec;
            nic = state.nic;
            
            //
            // Calculate penalty term.
            //
            result = 0;
            for(i=0; i<=nec+nic-1; i++)
            {
                if( state.activeset[n+i]>0 )
                {
                    alpha = 0;
                    p = -state.cleic[i,n];
                    for(j=0; j<=n-1; j++)
                    {
                        v = state.cleic[i,j];
                        p = p+v*x[j];
                        alpha = alpha+math.sqr(v*state.s[j]);
                    }
                    alpha = Math.Sqrt(alpha);
                    if( (double)(alpha)!=(double)(0) )
                    {
                        result = result+Math.Abs(p/alpha);
                    }
                }
            }
            return result;
        }

        /*************************************************************************
        This  subroutine  performs  correction of some (possibly infeasible) point
        with respect to a) current active set, b) all boundary  constraints,  both
        active and inactive:

        0) we calculate L1 penalty term for violation of active linear constraints
           (one which is returned by SASActiveLCPenalty1() function).
        1) first, it performs projection (orthogonal with respect to scale  matrix
           S) of X into current active set: X -> X1.
        2) next, we perform projection with respect to  ALL  boundary  constraints
           which are violated at X1: X1 -> X2.
        3) X is replaced by X2.

        The idea is that this function can preserve and enforce feasibility during
        optimization, and additional penalty parameter can be used to prevent algo
        from leaving feasible set because of rounding errors.

        INPUT PARAMETERS:
            S       -   active set object
            X       -   array[N], candidate point
            
        OUTPUT PARAMETERS:
            X       -   "improved" candidate point:
                        a) feasible with respect to all boundary constraints
                        b) feasibility with respect to active set is retained at
                           good level.
            Penalty -   penalty term, which can be added to function value if user
                        wants to penalize violation of constraints (recommended).
                        
        NOTE: this function is not intended to find exact  projection  (i.e.  best
              approximation) of X into feasible set. It just improves situation  a
              bit.
              Regular  use  of   this function will help you to retain feasibility
              - if you already have something to start  with  and  constrain  your
              steps is such way that the only source of infeasibility are roundoff
              errors.

          -- ALGLIB --
             Copyright 21.12.2012 by Bochkanov Sergey
        *************************************************************************/
        public static void sascorrection(sactiveset state,
            double[] x,
            ref double penalty)
        {
            int i = 0;
            int j = 0;
            int n = 0;
            double v = 0;
            int i_ = 0;

            penalty = 0;

            alglib.ap.assert(state.algostate==1, "SASCorrection: is not in optimization mode");
            sasrebuildbasis(state);
            n = state.n;
            apserv.rvectorsetlengthatleast(ref state.corrtmp, n);
            
            //
            // Calculate penalty term.
            //
            penalty = sasactivelcpenalty1(state, x);
            
            //
            // Perform projection 1.
            //
            // This projecton is given by:
            //
            //     x_proj = x - S*S*As'*(As*x-b)
            //
            // where x is original x before projection, S is a scale matrix,
            // As is a matrix of equality constraints (active set) which were
            // orthogonalized with respect to inner product given by S (i.e. we
            // have As*S*S'*As'=I), b is a right part of the orthogonalized
            // constraints.
            //
            // NOTE: you can verify that x_proj is strictly feasible w.r.t.
            //       active set by multiplying it by As - you will get
            //       As*x_proj = As*x - As*x + b = b.
            //
            //       This formula for projection can be obtained by solving
            //       following minimization problem.
            //
            //           min ||inv(S)*(x_proj-x)||^2 s.t. As*x_proj=b
            //       
            //
            for(i_=0; i_<=n-1;i_++)
            {
                state.corrtmp[i_] = x[i_];
            }
            for(i=0; i<=state.basissize-1; i++)
            {
                v = -state.sbasis[i,n];
                for(j=0; j<=n-1; j++)
                {
                    v = v+state.sbasis[i,j]*state.corrtmp[j];
                }
                for(j=0; j<=n-1; j++)
                {
                    state.corrtmp[j] = state.corrtmp[j]-v*state.sbasis[i,j]*math.sqr(state.s[j]);
                }
//.........这里部分代码省略.........

 public override alglib.apobject make_copy()
 {
     sactiveset _result = new sactiveset();
     _result.n = n;
     _result.algostate = algostate;
     _result.xc = (double[])xc.Clone();
     _result.hasxc = hasxc;
     _result.s = (double[])s.Clone();
     _result.h = (double[])h.Clone();
     _result.activeset = (int[])activeset.Clone();
     _result.basisisready = basisisready;
     _result.sbasis = (double[,])sbasis.Clone();
     _result.pbasis = (double[,])pbasis.Clone();
     _result.ibasis = (double[,])ibasis.Clone();
     _result.basissize = basissize;
     _result.constraintschanged = constraintschanged;
     _result.hasbndl = (bool[])hasbndl.Clone();
     _result.hasbndu = (bool[])hasbndu.Clone();
     _result.bndl = (double[])bndl.Clone();
     _result.bndu = (double[])bndu.Clone();
     _result.cleic = (double[,])cleic.Clone();
     _result.nec = nec;
     _result.nic = nic;
     _result.mtx = (double[])mtx.Clone();
     _result.mtas = (int[])mtas.Clone();
     _result.cdtmp = (double[])cdtmp.Clone();
     _result.corrtmp = (double[])corrtmp.Clone();
     _result.unitdiagonal = (double[])unitdiagonal.Clone();
     _result.solver = (snnls.snnlssolver)solver.make_copy();
     _result.scntmp = (double[])scntmp.Clone();
     _result.tmp0 = (double[])tmp0.Clone();
     _result.tmpfeas = (double[])tmpfeas.Clone();
     _result.tmpm0 = (double[,])tmpm0.Clone();
     _result.rctmps = (double[])rctmps.Clone();
     _result.rctmpg = (double[])rctmpg.Clone();
     _result.rctmprightpart = (double[])rctmprightpart.Clone();
     _result.rctmpdense0 = (double[,])rctmpdense0.Clone();
     _result.rctmpdense1 = (double[,])rctmpdense1.Clone();
     _result.rctmpisequality = (bool[])rctmpisequality.Clone();
     _result.rctmpconstraintidx = (int[])rctmpconstraintidx.Clone();
     _result.rctmplambdas = (double[])rctmplambdas.Clone();
     _result.tmpbasis = (double[,])tmpbasis.Clone();
     return _result;
 }

        /*************************************************************************
        This function recalculates constraints - activates  and  deactivates  them
        according to gradient value at current point.

        Algorithm  assumes  that  we  want  to make Quasi-Newton step from current
        point with diagonal Quasi-Newton matrix H. Constraints are  activated  and
        deactivated in such way that we won't violate any constraint by step.

        After call to  this  function  active set is ready to  try  preconditioned
        steepest descent step (SASDescentDirection-SASExploreDirection-SASMoveTo).

        Only already "active" and "candidate" elements of ActiveSet are  examined;
        constraints which are not active are not examined.

        INPUT PARAMETERS:
            State       -   active set object
            GC          -   array[N], gradient at XC
            
        OUTPUT PARAMETERS:
            State       -   active set object, with new set of constraint

          -- ALGLIB --
             Copyright 26.09.2012 by Bochkanov Sergey
        *************************************************************************/
        public static void sasreactivateconstraintsprec(sactiveset state,
            double[] gc)
        {
            alglib.ap.assert(state.algostate==1, "SASReactivateConstraintsPrec: must be in optimization mode");
            reactivateconstraints(state, gc, state.h);
        }

        /*************************************************************************
        This  subroutine  calculates  preconditioned  descent direction subject to
        current active set.

        INPUT PARAMETERS:
            State   -   active set object
            G       -   array[N], gradient
            H       -   array[N], Hessian matrix
            HA      -   active constraints orthogonalized in such way
                        that HA*inv(H)*HA'= I.
            Normalize-  whether we need normalized descent or not
            D       -   possibly preallocated buffer; automatically resized.
            
        OUTPUT PARAMETERS:
            D       -   descent direction projected onto current active set.
                        Components of D which correspond to active boundary
                        constraints are forced to be exactly zero.
                        In case D is non-zero and Normalize is True, it is
                        normalized to have unit norm.
                        
        NOTE: if we have N active constraints, D is explicitly set to zero.

          -- ALGLIB --
             Copyright 21.12.2012 by Bochkanov Sergey
        *************************************************************************/
        private static void constraineddescent(sactiveset state,
            double[] g,
            double[] h,
            double[,] ha,
            bool normalize,
            ref double[] d)
        {
            int i = 0;
            int j = 0;
            int n = 0;
            double v = 0;
            int nactive = 0;
            int i_ = 0;

            alglib.ap.assert(state.algostate==1, "SAS: internal error in ConstrainedDescent() - not in optimization mode");
            alglib.ap.assert(state.basisisready, "SAS: internal error in ConstrainedDescent() - no basis");
            n = state.n;
            apserv.rvectorsetlengthatleast(ref d, n);
            
            //
            // Calculate preconditioned constrained descent direction:
            //
            //     d := -inv(H)*( g - HA'*(HA*inv(H)*g) )
            //
            // Formula above always gives direction which is orthogonal to rows of HA.
            // You can verify it by multiplication of both sides by HA[i] (I-th row),
            // taking into account that HA*inv(H)*HA'= I (by definition of HA - it is
            // orthogonal basis with inner product given by inv(H)).
            //
            nactive = 0;
            for(i=0; i<=n-1; i++)
            {
                if( state.activeset[i]>0 )
                {
                    d[i] = 0;
                    nactive = nactive+1;
                }
                else
                {
                    d[i] = g[i];
                }
            }
            for(i=0; i<=state.basissize-1; i++)
            {
                v = 0.0;
                for(j=0; j<=n-1; j++)
                {
                    v = v+ha[i,j]*d[j]/h[j];
                }
                for(i_=0; i_<=n-1;i_++)
                {
                    d[i_] = d[i_] - v*ha[i,i_];
                }
                nactive = nactive+1;
            }
            v = 0.0;
            for(i=0; i<=n-1; i++)
            {
                if( state.activeset[i]>0 )
                {
                    d[i] = 0;
                }
                else
                {
                    d[i] = -(d[i]/h[i]);
                    v = v+math.sqr(d[i]);
                }
            }
            v = Math.Sqrt(v);
            if( nactive>=n )
            {
                v = 0;
                for(i=0; i<=n-1; i++)
                {
                    d[i] = 0;
//.........这里部分代码省略.........

        /*************************************************************************
        Another variation of SASSetLC(), which accepts  linear  constraints  using
        another representation.

        Linear constraints are inactive by default (after initial creation).

        INPUT PARAMETERS:
            State   -   SAS structure
            CLEIC   -   linear constraints, array[NEC+NIC,N+1].
                        Each row of C represents one constraint:
                        * first N elements correspond to coefficients,
                        * last element corresponds to the right part.
                        First NEC rows store equality constraints, next NIC -  are
                        inequality ones.
                        All elements of C (including right part) must be finite.
            NEC     -   number of equality constraints, NEC>=0
            NIC     -   number of inequality constraints, NIC>=0

        NOTE 1: linear (non-bound) constraints are satisfied only approximately:
        * there always exists some minor violation (about Epsilon in magnitude)
          due to rounding errors
        * numerical differentiation, if used, may  lead  to  function  evaluations
          outside  of the feasible  area,   because   algorithm  does  NOT  change
          numerical differentiation formula according to linear constraints.
        If you want constraints to be  satisfied  exactly, try to reformulate your
        problem  in  such  manner  that  all constraints will become boundary ones
        (this kind of constraints is always satisfied exactly, both in  the  final
        solution and in all intermediate points).

          -- ALGLIB --
             Copyright 28.11.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void sassetlcx(sactiveset state,
            double[,] cleic,
            int nec,
            int nic)
        {
            int n = 0;
            int i = 0;
            int j = 0;

            alglib.ap.assert(state.algostate==0, "SASSetLCX: you may change constraints only in modification mode");
            n = state.n;
            
            //
            // First, check for errors in the inputs
            //
            alglib.ap.assert(nec>=0, "SASSetLCX: NEC<0");
            alglib.ap.assert(nic>=0, "SASSetLCX: NIC<0");
            alglib.ap.assert(alglib.ap.cols(cleic)>=n+1 || nec+nic==0, "SASSetLCX: Cols(CLEIC)<N+1");
            alglib.ap.assert(alglib.ap.rows(cleic)>=nec+nic, "SASSetLCX: Rows(CLEIC)<NEC+NIC");
            alglib.ap.assert(apserv.apservisfinitematrix(cleic, nec+nic, n+1), "SASSetLCX: CLEIC contains infinite or NaN values!");
            
            //
            // Store constraints
            //
            apserv.rmatrixsetlengthatleast(ref state.cleic, nec+nic, n+1);
            state.nec = nec;
            state.nic = nic;
            for(i=0; i<=nec+nic-1; i++)
            {
                for(j=0; j<=n; j++)
                {
                    state.cleic[i,j] = cleic[i,j];
                }
            }
            
            //
            // Mark state as changed
            //
            state.constraintschanged = true;
        }

        /*************************************************************************
        This subroutine turns on optimization mode:
        1. feasibility in X is enforced  (in case X=S.XC and constraints  have not
           changed, algorithm just uses X without any modifications at all)
        2. constraints are marked as "candidate" or "inactive"

        INPUT PARAMETERS:
            S   -   active set object
            X   -   initial point (candidate), array[N]. It is expected that X
                    contains only finite values (we do not check it).
            
        OUTPUT PARAMETERS:
            S   -   state is changed
            X   -   initial point can be changed to enforce feasibility
            
        RESULT:
            True in case feasible point was found (mode was changed to "optimization")
            False in case no feasible point was found (mode was not changed)

          -- ALGLIB --
             Copyright 21.12.2012 by Bochkanov Sergey
        *************************************************************************/
        public static bool sasstartoptimization(sactiveset state,
            double[] x)
        {
            bool result = new bool();
            int n = 0;
            int nec = 0;
            int nic = 0;
            int i = 0;
            int j = 0;
            double v = 0;
            int i_ = 0;

            alglib.ap.assert(state.algostate==0, "SASStartOptimization: already in optimization mode");
            result = false;
            n = state.n;
            nec = state.nec;
            nic = state.nic;
            
            //
            // Enforce feasibility and calculate set of "candidate"/"active" constraints.
            // Always active equality constraints are marked as "active", all other constraints
            // are marked as "candidate".
            //
            apserv.ivectorsetlengthatleast(ref state.activeset, n+nec+nic);
            for(i=0; i<=n-1; i++)
            {
                if( state.hasbndl[i] && state.hasbndu[i] )
                {
                    if( (double)(state.bndl[i])>(double)(state.bndu[i]) )
                    {
                        return result;
                    }
                }
            }
            for(i_=0; i_<=n-1;i_++)
            {
                state.xc[i_] = x[i_];
            }
            if( state.nec+state.nic>0 )
            {
                
                //
                // General linear constraints are present; general code is used.
                //
                apserv.rvectorsetlengthatleast(ref state.tmp0, n);
                apserv.rvectorsetlengthatleast(ref state.tmpfeas, n+state.nic);
                apserv.rmatrixsetlengthatleast(ref state.tmpm0, state.nec+state.nic, n+state.nic+1);
                for(i=0; i<=state.nec+state.nic-1; i++)
                {
                    for(i_=0; i_<=n-1;i_++)
                    {
                        state.tmpm0[i,i_] = state.cleic[i,i_];
                    }
                    for(j=n; j<=n+state.nic-1; j++)
                    {
                        state.tmpm0[i,j] = 0;
                    }
                    if( i>=state.nec )
                    {
                        state.tmpm0[i,n+i-state.nec] = 1.0;
                    }
                    state.tmpm0[i,n+state.nic] = state.cleic[i,n];
                }
                for(i_=0; i_<=n-1;i_++)
                {
                    state.tmpfeas[i_] = state.xc[i_];
                }
                for(i=0; i<=state.nic-1; i++)
                {
                    v = 0.0;
                    for(i_=0; i_<=n-1;i_++)
                    {
                        v += state.cleic[i+state.nec,i_]*state.xc[i_];
                    }
                    state.tmpfeas[i+n] = Math.Max(state.cleic[i+state.nec,n]-v, 0.0);
                }
                if( !optserv.findfeasiblepoint(ref state.tmpfeas, state.bndl, state.hasbndl, state.bndu, state.hasbndu, n, state.nic, state.tmpm0, state.nec+state.nic, 1.0E-6, ref i, ref j) )
                {
//.........这里部分代码省略.........

        /*************************************************************************
        This function sets linear constraints for SAS object.

        Linear constraints are inactive by default (after initial creation).

        INPUT PARAMETERS:
            State   -   SAS structure
            C       -   linear constraints, array[K,N+1].
                        Each row of C represents one constraint, either equality
                        or inequality (see below):
                        * first N elements correspond to coefficients,
                        * last element corresponds to the right part.
                        All elements of C (including right part) must be finite.
            CT      -   type of constraints, array[K]:
                        * if CT[i]>0, then I-th constraint is C[i,*]*x >= C[i,n+1]
                        * if CT[i]=0, then I-th constraint is C[i,*]*x  = C[i,n+1]
                        * if CT[i]<0, then I-th constraint is C[i,*]*x <= C[i,n+1]
            K       -   number of equality/inequality constraints, K>=0

        NOTE 1: linear (non-bound) constraints are satisfied only approximately:
        * there always exists some minor violation (about Epsilon in magnitude)
          due to rounding errors
        * numerical differentiation, if used, may  lead  to  function  evaluations
          outside  of the feasible  area,   because   algorithm  does  NOT  change
          numerical differentiation formula according to linear constraints.
        If you want constraints to be  satisfied  exactly, try to reformulate your
        problem  in  such  manner  that  all constraints will become boundary ones
        (this kind of constraints is always satisfied exactly, both in  the  final
        solution and in all intermediate points).

          -- ALGLIB --
             Copyright 28.11.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void sassetlc(sactiveset state,
            double[,] c,
            int[] ct,
            int k)
        {
            int n = 0;
            int i = 0;
            int i_ = 0;

            alglib.ap.assert(state.algostate==0, "SASSetLC: you may change constraints only in modification mode");
            n = state.n;
            
            //
            // First, check for errors in the inputs
            //
            alglib.ap.assert(k>=0, "SASSetLC: K<0");
            alglib.ap.assert(alglib.ap.cols(c)>=n+1 || k==0, "SASSetLC: Cols(C)<N+1");
            alglib.ap.assert(alglib.ap.rows(c)>=k, "SASSetLC: Rows(C)<K");
            alglib.ap.assert(alglib.ap.len(ct)>=k, "SASSetLC: Length(CT)<K");
            alglib.ap.assert(apserv.apservisfinitematrix(c, k, n+1), "SASSetLC: C contains infinite or NaN values!");
            
            //
            // Handle zero K
            //
            if( k==0 )
            {
                state.nec = 0;
                state.nic = 0;
                state.constraintschanged = true;
                return;
            }
            
            //
            // Equality constraints are stored first, in the upper
            // NEC rows of State.CLEIC matrix. Inequality constraints
            // are stored in the next NIC rows.
            //
            // NOTE: we convert inequality constraints to the form
            // A*x<=b before copying them.
            //
            apserv.rmatrixsetlengthatleast(ref state.cleic, k, n+1);
            state.nec = 0;
            state.nic = 0;
            for(i=0; i<=k-1; i++)
            {
                if( ct[i]==0 )
                {
                    for(i_=0; i_<=n;i_++)
                    {
                        state.cleic[state.nec,i_] = c[i,i_];
                    }
                    state.nec = state.nec+1;
                }
            }
            for(i=0; i<=k-1; i++)
            {
                if( ct[i]!=0 )
                {
                    if( ct[i]>0 )
                    {
                        for(i_=0; i_<=n;i_++)
                        {
                            state.cleic[state.nec+state.nic,i_] = -c[i,i_];
                        }
                    }
                    else
                    {
//.........这里部分代码省略.........

        /*************************************************************************
        Modification  of  the  preconditioner:  diagonal of approximate Hessian is
        used.

        INPUT PARAMETERS:
            State   -   structure which stores algorithm state
            D       -   diagonal of the approximate Hessian, array[0..N-1],
                        (if larger, only leading N elements are used).

        NOTE 1: D[i] should be positive. Exception will be thrown otherwise.

        NOTE 2: you should pass diagonal of approximate Hessian - NOT ITS INVERSE.

          -- ALGLIB --
             Copyright 21.12.2012 by Bochkanov Sergey
        *************************************************************************/
        public static void sassetprecdiag(sactiveset state,
            double[] d)
        {
            int i = 0;

            alglib.ap.assert(state.algostate==0, "SASSetPrecDiag: you may change preconditioner only in modification mode");
            alglib.ap.assert(alglib.ap.len(d)>=state.n, "SASSetPrecDiag: D is too short");
            for(i=0; i<=state.n-1; i++)
            {
                alglib.ap.assert(math.isfinite(d[i]), "SASSetPrecDiag: D contains infinite or NAN elements");
                alglib.ap.assert((double)(d[i])>(double)(0), "SASSetPrecDiag: D contains non-positive elements");
            }
            for(i=0; i<=state.n-1; i++)
            {
                state.h[i] = d[i];
            }
        }

        /*************************************************************************
        This function sets scaling coefficients for SAS object.

        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

        During orthogonalization phase, scale is used to calculate drop tolerances
        (whether vector is significantly non-zero or not).

        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 21.12.2012 by Bochkanov Sergey
        *************************************************************************/
        public static void sassetscale(sactiveset state,
            double[] s)
        {
            int i = 0;

            alglib.ap.assert(state.algostate==0, "SASSetScale: you may change scale only in modification mode");
            alglib.ap.assert(alglib.ap.len(s)>=state.n, "SASSetScale: Length(S)<N");
            for(i=0; i<=state.n-1; i++)
            {
                alglib.ap.assert(math.isfinite(s[i]), "SASSetScale: S contains infinite or NAN elements");
                alglib.ap.assert((double)(s[i])!=(double)(0), "SASSetScale: S contains zero elements");
            }
            for(i=0; i<=state.n-1; i++)
            {
                state.s[i] = Math.Abs(s[i]);
            }
        }

        /*************************************************************************
        This   subroutine   is   used  to initialize active set. By default, empty
        N-variable model with no constraints is  generated.  Previously  allocated
        buffer variables are reused as much as possible.

        Two use cases for this object are described below.

        CASE 1 - STEEPEST DESCENT:

            SASInit()
            repeat:
                SASReactivateConstraints()
                SASDescentDirection()
                SASExploreDirection()
                SASMoveTo()
            until convergence

        CASE 1 - PRECONDITIONED STEEPEST DESCENT:

            SASInit()
            repeat:
                SASReactivateConstraintsPrec()
                SASDescentDirectionPrec()
                SASExploreDirection()
                SASMoveTo()
            until convergence

          -- ALGLIB --
             Copyright 21.12.2012 by Bochkanov Sergey
        *************************************************************************/
        public static void sasinit(int n,
            sactiveset s)
        {
            int i = 0;

            s.n = n;
            s.algostate = 0;
            
            //
            // Constraints
            //
            s.constraintschanged = true;
            s.nec = 0;
            s.nic = 0;
            apserv.rvectorsetlengthatleast(ref s.bndl, n);
            apserv.bvectorsetlengthatleast(ref s.hasbndl, n);
            apserv.rvectorsetlengthatleast(ref s.bndu, n);
            apserv.bvectorsetlengthatleast(ref s.hasbndu, n);
            for(i=0; i<=n-1; i++)
            {
                s.bndl[i] = Double.NegativeInfinity;
                s.bndu[i] = Double.PositiveInfinity;
                s.hasbndl[i] = false;
                s.hasbndu[i] = false;
            }
            
            //
            // current point, scale
            //
            s.hasxc = false;
            apserv.rvectorsetlengthatleast(ref s.xc, n);
            apserv.rvectorsetlengthatleast(ref s.s, n);
            apserv.rvectorsetlengthatleast(ref s.h, n);
            for(i=0; i<=n-1; i++)
            {
                s.xc[i] = 0.0;
                s.s[i] = 1.0;
                s.h[i] = 1.0;
            }
            
            //
            // Other
            //
            apserv.rvectorsetlengthatleast(ref s.unitdiagonal, n);
            for(i=0; i<=n-1; i++)
            {
                s.unitdiagonal[i] = 1.0;
            }
        }

        /*************************************************************************
        This subroutine calculates scaled norm of  vector  after  projection  onto
        subspace of active constraints. Most often this function is used  to  test
        stopping conditions.

        INPUT PARAMETERS:
            S       -   active set object
            D       -   vector whose norm is calculated
            
        RESULT:
            Vector norm (after projection and scaling)
            
        NOTE: projection is performed first, scaling is performed after projection
                        
        NOTE: if we have N active constraints, zero value (exact zero) is returned

          -- ALGLIB --
             Copyright 21.12.2012 by Bochkanov Sergey
        *************************************************************************/
        public static double sasscaledconstrainednorm(sactiveset state,
            double[] d)
        {
            double result = 0;
            int i = 0;
            int n = 0;
            double v = 0;
            int nactive = 0;
            int i_ = 0;

            alglib.ap.assert(state.algostate==1, "SASMoveTo: is not in optimization mode");
            n = state.n;
            apserv.rvectorsetlengthatleast(ref state.scntmp, n);
            
            //
            // Prepare basis (if needed)
            //
            sasrebuildbasis(state);
            
            //
            // Calculate descent direction
            //
            nactive = 0;
            for(i=0; i<=n-1; i++)
            {
                if( state.activeset[i]>0 )
                {
                    state.scntmp[i] = 0;
                    nactive = nactive+1;
                }
                else
                {
                    state.scntmp[i] = d[i];
                }
            }
            if( nactive+state.basissize>=n )
            {
                
                //
                // Quick exit if number of active constraints is N or larger
                //
                result = 0.0;
                return result;
            }
            for(i=0; i<=state.basissize-1; i++)
            {
                v = 0.0;
                for(i_=0; i_<=n-1;i_++)
                {
                    v += state.ibasis[i,i_]*state.scntmp[i_];
                }
                for(i_=0; i_<=n-1;i_++)
                {
                    state.scntmp[i_] = state.scntmp[i_] - v*state.ibasis[i,i_];
                }
            }
            v = 0.0;
            for(i=0; i<=n-1; i++)
            {
                v = v+math.sqr(state.s[i]*state.scntmp[i]);
            }
            result = Math.Sqrt(v);
            return result;
        }

        /*************************************************************************
        This function explores search direction and calculates bound for  step  as
        well as information for activation of constraints.

        INPUT PARAMETERS:
            State       -   SAS structure which stores current point and all other
                            active set related information
            D           -   descent direction to explore

        OUTPUT PARAMETERS:
            StpMax      -   upper  limit  on  step  length imposed by yet inactive
                            constraints. Can be  zero  in  case  some  constraints
                            can be activated by zero step.  Equal  to  some  large
                            value in case step is unlimited.
            CIdx        -   -1 for unlimited step, in [0,N+NEC+NIC) in case of
                            limited step.
            VVal        -   value which is assigned to X[CIdx] during activation.
                            For CIdx<0 or CIdx>=N some dummy value is assigned to
                            this parameter.
        *************************************************************************/
        public static void sasexploredirection(sactiveset state,
            double[] d,
            ref double stpmax,
            ref int cidx,
            ref double vval)
        {
            int n = 0;
            int nec = 0;
            int nic = 0;
            int i = 0;
            double prevmax = 0;
            double vc = 0;
            double vd = 0;
            int i_ = 0;

            stpmax = 0;
            cidx = 0;
            vval = 0;

            alglib.ap.assert(state.algostate==1, "SASExploreDirection: is not in optimization mode");
            n = state.n;
            nec = state.nec;
            nic = state.nic;
            cidx = -1;
            vval = 0;
            stpmax = 1.0E50;
            for(i=0; i<=n-1; i++)
            {
                if( state.activeset[i]<=0 )
                {
                    alglib.ap.assert(!state.hasbndl[i] || (double)(state.xc[i])>=(double)(state.bndl[i]), "SASExploreDirection: internal error - infeasible X");
                    alglib.ap.assert(!state.hasbndu[i] || (double)(state.xc[i])<=(double)(state.bndu[i]), "SASExploreDirection: internal error - infeasible X");
                    if( state.hasbndl[i] && (double)(d[i])<(double)(0) )
                    {
                        prevmax = stpmax;
                        stpmax = apserv.safeminposrv(state.xc[i]-state.bndl[i], -d[i], stpmax);
                        if( (double)(stpmax)<(double)(prevmax) )
                        {
                            cidx = i;
                            vval = state.bndl[i];
                        }
                    }
                    if( state.hasbndu[i] && (double)(d[i])>(double)(0) )
                    {
                        prevmax = stpmax;
                        stpmax = apserv.safeminposrv(state.bndu[i]-state.xc[i], d[i], stpmax);
                        if( (double)(stpmax)<(double)(prevmax) )
                        {
                            cidx = i;
                            vval = state.bndu[i];
                        }
                    }
                }
            }
            for(i=nec; i<=nec+nic-1; i++)
            {
                if( state.activeset[n+i]<=0 )
                {
                    vc = 0.0;
                    for(i_=0; i_<=n-1;i_++)
                    {
                        vc += state.cleic[i,i_]*state.xc[i_];
                    }
                    vc = vc-state.cleic[i,n];
                    vd = 0.0;
                    for(i_=0; i_<=n-1;i_++)
                    {
                        vd += state.cleic[i,i_]*d[i_];
                    }
                    if( (double)(vd)<=(double)(0) )
                    {
                        continue;
                    }
                    if( (double)(vc)<(double)(0) )
                    {
                        
                        //
                        // XC is strictly feasible with respect to I-th constraint,
                        // we can perform non-zero step because there is non-zero distance
                        // between XC and bound.
//.........这里部分代码省略.........

        /*************************************************************************
        This subroutine turns off optimization mode.

        INPUT PARAMETERS:
            S   -   active set object
            
        OUTPUT PARAMETERS:
            S   -   state is changed

        NOTE: this function can be called many times for optimizer which was
              already stopped.

          -- ALGLIB --
             Copyright 21.12.2012 by Bochkanov Sergey
        *************************************************************************/
        public static void sasstopoptimization(sactiveset state)
        {
            state.algostate = 0;
        }

        /*************************************************************************
        This subroutine moves current point to XN, which can be:
        a) point in the direction previously explored  with  SASExploreDirection()
           function (in this case NeedAct/CIdx/CVal are used)
        b) point in arbitrary direction, not necessarily previously  checked  with
           SASExploreDirection() function.

        Step may activate one constraint. It is assumed than XN  is  approximately
        feasible (small error as  large  as several  ulps  is  possible).   Strict
        feasibility  with  respect  to  bound  constraints  is  enforced    during
        activation, feasibility with respect to general linear constraints is  not
        enforced.

        This function activates boundary constraints, such that both is True:
        1) XC[I] is not at the boundary
        2) XN[I] is at the boundary or beyond it

        INPUT PARAMETERS:
            S       -   active set object
            XN      -   new point.
            NeedAct -   True in case one constraint needs activation
            CIdx    -   index of constraint, in [0,N+NEC+NIC).
                        Ignored if NeedAct is false.
                        This value is calculated by SASExploreDirection().
            CVal    -   for CIdx in [0,N) this field stores value which is
                        assigned to XC[CIdx] during activation. CVal is ignored in
                        other cases.
                        This value is calculated by SASExploreDirection().
            
        OUTPUT PARAMETERS:
            S       -   current point and list of active constraints are changed.

        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

        NOTE: in general case State.XC<>XN because activation of  constraints  may
              slightly change current point (to enforce feasibility).

          -- ALGLIB --
             Copyright 21.12.2012 by Bochkanov Sergey
        *************************************************************************/
        public static int sasmoveto(sactiveset state,
            double[] xn,
            bool needact,
            int cidx,
            double cval)
        {
            int result = 0;
            int n = 0;
            int nec = 0;
            int nic = 0;
            int i = 0;
            bool wasactivation = new bool();

            alglib.ap.assert(state.algostate==1, "SASMoveTo: is not in optimization mode");
            n = state.n;
            nec = state.nec;
            nic = state.nic;
            
            //
            // Save previous state, update current point
            //
            apserv.rvectorsetlengthatleast(ref state.mtx, n);
            apserv.ivectorsetlengthatleast(ref state.mtas, n+nec+nic);
            for(i=0; i<=n-1; i++)
            {
                state.mtx[i] = state.xc[i];
                state.xc[i] = xn[i];
            }
            for(i=0; i<=n+nec+nic-1; i++)
            {
                state.mtas[i] = state.activeset[i];
            }
            
            //
            // Activate constraints
            //
            wasactivation = false;
            if( needact )
            {
                
                //
                // Activation
                //
                alglib.ap.assert(cidx>=0 && cidx<n+nec+nic, "SASMoveTo: incorrect CIdx");