/*************************************************************************
        This function scales and copies NMain elements of GUnscaled into GScaled.
        Other NSlack components of GScaled are set to zero.
        *************************************************************************/
        private static void scalegradientandexpand(minbleicstate state,
            double[] gunscaled,
            ref double[] gscaled)
        {
            int i = 0;

            for(i=0; i<=state.nmain-1; i++)
            {
                gscaled[i] = gunscaled[i]*state.transforms[i];
            }
            for(i=0; i<=state.nslack-1; i++)
            {
                gscaled[state.nmain+i] = 0;
            }
        }

        /*************************************************************************
        This subroutine finalizes internal structures after emergency  termination
        from State.LSStart report (see comments on MinBLEICState for more information).

        INPUT PARAMETERS:
            State   -   structure after exit from LSStart report

          -- ALGLIB --
             Copyright 28.11.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void minbleicemergencytermination(minbleicstate state)
        {
            sactivesets.sasstopoptimization(state.sas);
        }

 /*************************************************************************
 Clears request fileds (to be sure that we don't forget to clear something)
 *************************************************************************/
 private static void clearrequestfields(minbleicstate state)
 {
     state.needf = false;
     state.needfg = false;
     state.xupdated = false;
     state.lsstart = false;
 }

        /*************************************************************************
        BLEIC results

        INPUT PARAMETERS:
            State   -   algorithm state

        OUTPUT PARAMETERS:
            X       -   array[0..N-1], solution
            Rep     -   optimization report. You should check Rep.TerminationType
                        in  order  to  distinguish  successful  termination  from
                        unsuccessful one:
                        * -8    internal integrity control  detected  infinite or
                                NAN   values   in   function/gradient.   Abnormal
                                termination signalled.
                        * -7   gradient verification failed.
                               See MinBLEICSetGradientCheck() for more information.
                        * -3   inconsistent constraints. Feasible point is
                               either nonexistent or too hard to find. Try to
                               restart optimizer with better initial approximation
                        *  1   relative function improvement is no more than EpsF.
                        *  2   scaled step is no more than EpsX.
                        *  4   scaled gradient norm is no more than EpsG.
                        *  5   MaxIts steps was taken
                        *  8   terminated by user who called minbleicrequesttermination().
                               X contains point which was "current accepted"  when
                               termination request was submitted.
                        More information about fields of this  structure  can  be
                        found in the comments on MinBLEICReport datatype.
           
          -- ALGLIB --
             Copyright 28.11.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void minbleicresults(minbleicstate state,
            ref double[] x,
            minbleicreport rep)
        {
            x = new double[0];

            minbleicresultsbuf(state, ref x, rep);
        }

        /*************************************************************************
        This subroutine restarts algorithm from new point.
        All optimization parameters (including constraints) are left unchanged.

        This  function  allows  to  solve multiple  optimization  problems  (which
        must have  same number of dimensions) without object reallocation penalty.

        INPUT PARAMETERS:
            State   -   structure previously allocated with MinBLEICCreate call.
            X       -   new starting point.

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

            n = state.nmain;
            
            //
            // First, check for errors in the inputs
            //
            alglib.ap.assert(alglib.ap.len(x)>=n, "MinBLEICRestartFrom: Length(X)<N");
            alglib.ap.assert(apserv.isfinitevector(x, n), "MinBLEICRestartFrom: X contains infinite or NaN values!");
            
            //
            // Set XC
            //
            for(i_=0; i_<=n-1;i_++)
            {
                state.xstart[i_] = x[i_];
            }
            
            //
            // prepare RComm facilities
            //
            state.rstate.ia = new int[6+1];
            state.rstate.ba = new bool[0+1];
            state.rstate.ra = new double[5+1];
            state.rstate.stage = -1;
            clearrequestfields(state);
            sactivesets.sasstopoptimization(state.sas);
        }

        /*************************************************************************
        This function turns on/off reporting.

        INPUT PARAMETERS:
            State   -   structure which stores algorithm state
            NeedXRep-   whether iteration reports are needed or not

        If NeedXRep is True, algorithm will call rep() callback function if  it is
        provided to MinBLEICOptimize().

          -- ALGLIB --
             Copyright 28.11.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void minbleicsetxrep(minbleicstate state,
            bool needxrep)
        {
            state.xrep = needxrep;
        }

        /*************************************************************************
        This function sets maximum step length

        IMPORTANT: this feature is hard to combine with preconditioning. You can't
        set upper limit on step length, when you solve optimization  problem  with
        linear (non-boundary) constraints AND preconditioner turned on.

        When  non-boundary  constraints  are  present,  you  have to either a) use
        preconditioner, or b) use upper limit on step length.  YOU CAN'T USE BOTH!
        In this case algorithm will terminate with appropriate error code.

        INPUT PARAMETERS:
            State   -   structure which stores algorithm state
            StpMax  -   maximum step length, >=0. Set StpMax to 0.0,  if you don't
                        want to limit step length.

        Use this subroutine when you optimize target function which contains exp()
        or  other  fast  growing  functions,  and optimization algorithm makes too
        large  steps  which  lead   to overflow. This function allows us to reject
        steps  that  are  too  large  (and  therefore  expose  us  to the possible
        overflow) without actually calculating function value at the x+stp*d.

          -- ALGLIB --
             Copyright 02.04.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void minbleicsetstpmax(minbleicstate state,
            double stpmax)
        {
            alglib.ap.assert(math.isfinite(stpmax), "MinBLEICSetStpMax: StpMax is not finite!");
            alglib.ap.assert((double)(stpmax)>=(double)(0), "MinBLEICSetStpMax: StpMax<0!");
            state.stpmax = stpmax;
        }

        /*************************************************************************
        Internal initialization subroutine
        *************************************************************************/
        private static void minbleicinitinternal(int n,
            double[] x,
            double diffstep,
            minbleicstate state)
        {
            int i = 0;
            double[,] c = new double[0,0];
            int[] ct = new int[0];

            state.nmain = n;
            state.optdim = 0;
            state.diffstep = diffstep;
            state.bndloriginal = new double[n];
            state.bndleffective = new double[n];
            state.hasbndl = new bool[n];
            state.bnduoriginal = new double[n];
            state.bndueffective = new double[n];
            state.hasbndu = new bool[n];
            state.xstart = new double[n];
            state.soriginal = new double[n];
            state.x = new double[n];
            state.g = new double[n];
            for(i=0; i<=n-1; i++)
            {
                state.bndloriginal[i] = Double.NegativeInfinity;
                state.hasbndl[i] = false;
                state.bnduoriginal[i] = Double.PositiveInfinity;
                state.hasbndu[i] = false;
                state.soriginal[i] = 1.0;
            }
            minbleicsetlc(state, c, ct, 0);
            minbleicsetinnercond(state, 0.0, 0.0, 0.0);
            minbleicsetoutercond(state, 1.0E-6, 1.0E-6);
            minbleicsetmaxits(state, 0);
            minbleicsetxrep(state, false);
            minbleicsetstpmax(state, 0.0);
            minbleicsetprecdefault(state);
            minbleicrestartfrom(state, x);
        }

 public override alglib.apobject make_copy()
 {
     minbleicstate _result = new minbleicstate();
     _result.nmain = nmain;
     _result.nslack = nslack;
     _result.epsg = epsg;
     _result.epsf = epsf;
     _result.epsx = epsx;
     _result.maxits = maxits;
     _result.xrep = xrep;
     _result.drep = drep;
     _result.stpmax = stpmax;
     _result.diffstep = diffstep;
     _result.sas = (sactivesets.sactiveset)sas.make_copy();
     _result.s = (double[])s.Clone();
     _result.prectype = prectype;
     _result.diagh = (double[])diagh.Clone();
     _result.x = (double[])x.Clone();
     _result.f = f;
     _result.g = (double[])g.Clone();
     _result.needf = needf;
     _result.needfg = needfg;
     _result.xupdated = xupdated;
     _result.lsstart = lsstart;
     _result.steepestdescentstep = steepestdescentstep;
     _result.boundedstep = boundedstep;
     _result.userterminationneeded = userterminationneeded;
     _result.teststep = teststep;
     _result.rstate = (rcommstate)rstate.make_copy();
     _result.ugc = (double[])ugc.Clone();
     _result.cgc = (double[])cgc.Clone();
     _result.xn = (double[])xn.Clone();
     _result.ugn = (double[])ugn.Clone();
     _result.cgn = (double[])cgn.Clone();
     _result.xp = (double[])xp.Clone();
     _result.fc = fc;
     _result.fn = fn;
     _result.fp = fp;
     _result.d = (double[])d.Clone();
     _result.cleic = (double[,])cleic.Clone();
     _result.nec = nec;
     _result.nic = nic;
     _result.lastgoodstep = lastgoodstep;
     _result.lastscaledgoodstep = lastscaledgoodstep;
     _result.maxscaledgrad = maxscaledgrad;
     _result.hasbndl = (bool[])hasbndl.Clone();
     _result.hasbndu = (bool[])hasbndu.Clone();
     _result.bndl = (double[])bndl.Clone();
     _result.bndu = (double[])bndu.Clone();
     _result.repinneriterationscount = repinneriterationscount;
     _result.repouteriterationscount = repouteriterationscount;
     _result.repnfev = repnfev;
     _result.repvaridx = repvaridx;
     _result.repterminationtype = repterminationtype;
     _result.repdebugeqerr = repdebugeqerr;
     _result.repdebugfs = repdebugfs;
     _result.repdebugff = repdebugff;
     _result.repdebugdx = repdebugdx;
     _result.repdebugfeasqpits = repdebugfeasqpits;
     _result.repdebugfeasgpaits = repdebugfeasgpaits;
     _result.xstart = (double[])xstart.Clone();
     _result.solver = (snnls.snnlssolver)solver.make_copy();
     _result.fbase = fbase;
     _result.fm2 = fm2;
     _result.fm1 = fm1;
     _result.fp1 = fp1;
     _result.fp2 = fp2;
     _result.xm1 = xm1;
     _result.xp1 = xp1;
     _result.gm1 = gm1;
     _result.gp1 = gp1;
     _result.cidx = cidx;
     _result.cval = cval;
     _result.tmpprec = (double[])tmpprec.Clone();
     _result.tmp0 = (double[])tmp0.Clone();
     _result.nfev = nfev;
     _result.mcstage = mcstage;
     _result.stp = stp;
     _result.curstpmax = curstpmax;
     _result.activationstep = activationstep;
     _result.work = (double[])work.Clone();
     _result.lstate = (linmin.linminstate)lstate.make_copy();
     _result.trimthreshold = trimthreshold;
     _result.nonmonotoniccnt = nonmonotoniccnt;
     _result.bufyk = (double[,])bufyk.Clone();
     _result.bufsk = (double[,])bufsk.Clone();
     _result.bufrho = (double[])bufrho.Clone();
     _result.buftheta = (double[])buftheta.Clone();
     _result.bufsize = bufsize;
     return _result;
 }

        /*************************************************************************
        This function projects gradient onto equality constrained subspace
        *************************************************************************/
        private static void makegradientprojection(minbleicstate state,
            ref double[] pg)
        {
            int i = 0;
            int nmain = 0;
            int nslack = 0;
            double v = 0;
            int i_ = 0;

            nmain = state.nmain;
            nslack = state.nslack;
            for(i=0; i<=nmain+nslack-1; i++)
            {
                if( state.activeconstraints[i] )
                {
                    pg[i] = 0;
                }
            }
            for(i=0; i<=state.cecnt-1; i++)
            {
                v = 0.0;
                for(i_=0; i_<=nmain+nslack-1;i_++)
                {
                    v += pg[i_]*state.cecurrent[i,i_];
                }
                for(i_=0; i_<=nmain+nslack-1;i_++)
                {
                    pg[i_] = pg[i_] - v*state.cecurrent[i,i_];
                }
            }
        }

        /*************************************************************************
        This function prepares equality constrained subproblem:

        1. X is used to activate constraints (if there are ones which are still
           inactive, but should be activated).
        2. constraints matrix CEOrt is copied to CECurrent and modified  according
           to the list of active bound constraints (corresponding elements are
           filled by zeros and reorthogonalized).
        3. XE - least squares solution of equality constraints - is recalculated
        4. X is copied to PX and projected onto equality constrained subspace
        5. inactive constraints are checked against PX - if there is at least one
           which should be activated, we activate it and move back to (2)
        6. as result, PX is feasible with respect to bound constraints - step (5)
           guarantees it. But PX can be infeasible with respect to equality ones,
           because step (2) is done without checks for consistency. As the final
           step, we check that PX is feasible. If not, we return False. True is
           returned otherwise.

        If this algorithm returned True, then:
        * X is not changed
        * PX contains projection of X onto constrained subspace
        * G is not changed
        * PG contains projection of G onto constrained subspace
        * PX is feasible with respect to all constraints
        * all constraints which are active at PX, are activated
        *************************************************************************/
        private static bool prepareconstraintmatrix(minbleicstate state,
            double[] x,
            double[] g,
            ref double[] px,
            ref double[] pg)
        {
            bool result = new bool();
            int i = 0;
            int nmain = 0;
            int nslack = 0;
            double v = 0;
            double ferr = 0;
            int i_ = 0;

            nmain = state.nmain;
            nslack = state.nslack;
            result = true;
            
            //
            // Step 1
            //
            additionalcheckforconstraints(state, x);
            
            //
            // Steps 2-5
            //
            do
            {
                
                //
                // Steps 2-3
                //
                rebuildcexe(state);
                
                //
                // Step 4
                //
                // Calculate PX, PG
                //
                for(i_=0; i_<=nmain+nslack-1;i_++)
                {
                    px[i_] = x[i_];
                }
                for(i_=0; i_<=nmain+nslack-1;i_++)
                {
                    px[i_] = px[i_] - state.xe[i_];
                }
                for(i_=0; i_<=nmain+nslack-1;i_++)
                {
                    pg[i_] = g[i_];
                }
                for(i=0; i<=nmain+nslack-1; i++)
                {
                    if( state.activeconstraints[i] )
                    {
                        px[i] = 0;
                        pg[i] = 0;
                    }
                }
                for(i=0; i<=state.cecnt-1; i++)
                {
                    v = 0.0;
                    for(i_=0; i_<=nmain+nslack-1;i_++)
                    {
                        v += px[i_]*state.cecurrent[i,i_];
                    }
                    for(i_=0; i_<=nmain+nslack-1;i_++)
                    {
                        px[i_] = px[i_] - v*state.cecurrent[i,i_];
                    }
                    v = 0.0;
                    for(i_=0; i_<=nmain+nslack-1;i_++)
                    {
                        v += pg[i_]*state.cecurrent[i,i_];
//.........这里部分代码省略.........

        /*************************************************************************
        This function rebuilds CECurrent and XE according to current set of
        active bound constraints.
        *************************************************************************/
        private static void rebuildcexe(minbleicstate state)
        {
            int i = 0;
            int j = 0;
            int k = 0;
            int nmain = 0;
            int nslack = 0;
            double v = 0;
            int i_ = 0;

            nmain = state.nmain;
            nslack = state.nslack;
            ablas.rmatrixcopy(state.cecnt, nmain+nslack+1, state.ceeffective, 0, 0, ref state.cecurrent, 0, 0);
            for(i=0; i<=state.cecnt-1; i++)
            {
                
                //
                // "Subtract" active bound constraints from I-th linear constraint
                //
                for(j=0; j<=nmain+nslack-1; j++)
                {
                    if( state.activeconstraints[j] )
                    {
                        state.cecurrent[i,nmain+nslack] = state.cecurrent[i,nmain+nslack]-state.cecurrent[i,j]*state.constrainedvalues[j];
                        state.cecurrent[i,j] = 0.0;
                    }
                }
                
                //
                // Reorthogonalize I-th constraint with respect to previous ones
                // NOTE: we also update right part, which is CECurrent[...,NMain+NSlack].
                //
                for(k=0; k<=i-1; k++)
                {
                    v = 0.0;
                    for(i_=0; i_<=nmain+nslack-1;i_++)
                    {
                        v += state.cecurrent[k,i_]*state.cecurrent[i,i_];
                    }
                    for(i_=0; i_<=nmain+nslack;i_++)
                    {
                        state.cecurrent[i,i_] = state.cecurrent[i,i_] - v*state.cecurrent[k,i_];
                    }
                }
                
                //
                // Calculate norm of I-th row of CECurrent. Fill by zeros, if it is
                // too small. Normalize otherwise.
                //
                // NOTE: we also scale last column of CECurrent (right part)
                //
                v = 0.0;
                for(i_=0; i_<=nmain+nslack-1;i_++)
                {
                    v += state.cecurrent[i,i_]*state.cecurrent[i,i_];
                }
                v = Math.Sqrt(v);
                if( (double)(v)>(double)(10000*math.machineepsilon) )
                {
                    v = 1/v;
                    for(i_=0; i_<=nmain+nslack;i_++)
                    {
                        state.cecurrent[i,i_] = v*state.cecurrent[i,i_];
                    }
                }
                else
                {
                    for(j=0; j<=nmain+nslack; j++)
                    {
                        state.cecurrent[i,j] = 0;
                    }
                }
            }
            for(j=0; j<=nmain+nslack-1; j++)
            {
                state.xe[j] = 0;
            }
            for(i=0; i<=nmain+nslack-1; i++)
            {
                if( state.activeconstraints[i] )
                {
                    state.xe[i] = state.xe[i]+state.constrainedvalues[i];
                }
            }
            for(i=0; i<=state.cecnt-1; i++)
            {
                v = state.cecurrent[i,nmain+nslack];
                for(i_=0; i_<=nmain+nslack-1;i_++)
                {
                    state.xe[i_] = state.xe[i_] + v*state.cecurrent[i,i_];
                }
            }
        }

        /*************************************************************************
        This function makes additional check for constraints which can be activated.

        We try activate constraints one by one, but it is possible that several
        constraints should be activated during one iteration. It this case only
        one of them (probably last) will be activated. This function will fix it -
        it will pass through constraints and activate those which are at the boundary
        or beyond it.

        It will return True, if at least one constraint was activated by this function.
        *************************************************************************/
        private static bool additionalcheckforconstraints(minbleicstate state,
            double[] x)
        {
            bool result = new bool();
            int i = 0;
            int nmain = 0;
            int nslack = 0;

            result = false;
            nmain = state.nmain;
            nslack = state.nslack;
            for(i=0; i<=nmain-1; i++)
            {
                if( !state.activeconstraints[i] )
                {
                    if( state.hasbndl[i] )
                    {
                        if( (double)(x[i])<=(double)(state.bndleffective[i]) )
                        {
                            state.activeconstraints[i] = true;
                            state.constrainedvalues[i] = state.bndleffective[i];
                            result = true;
                        }
                    }
                    if( state.hasbndu[i] )
                    {
                        if( (double)(x[i])>=(double)(state.bndueffective[i]) )
                        {
                            state.activeconstraints[i] = true;
                            state.constrainedvalues[i] = state.bndueffective[i];
                            result = true;
                        }
                    }
                }
            }
            for(i=0; i<=nslack-1; i++)
            {
                if( !state.activeconstraints[nmain+i] )
                {
                    if( (double)(x[nmain+i])<=(double)(0) )
                    {
                        state.activeconstraints[nmain+i] = true;
                        state.constrainedvalues[nmain+i] = 0;
                        result = true;
                    }
                }
            }
            return result;
        }

        /*************************************************************************
        This subroutine applies modifications to the target function given by
        its value F and gradient G at the projected point X which lies in the
        equality constrained subspace.

        Following modifications are applied:
        * modified barrier functions to handle inequality constraints
          (both F and G are modified)
        * projection of gradient into equality constrained subspace
          (only G is modified)
        * quadratic penalty for deviations from equality constrained subspace
          (both F and G are modified)

        It also calculates gradient norm (three different norms for three
        different types of gradient), feasibility and complementary slackness
        errors.

        INPUT PARAMETERS:
            State   -   optimizer state (we use its fields to get information
                        about constraints)
            X       -   point (projected into equality constrained subspace)
            R       -   residual from projection
            RNorm2  -   residual norm squared
            F       -   function value at X
            G       -   function gradient at X

        OUTPUT PARAMETERS:
            F       -   modified function value at X
            G       -   modified function gradient at X
            GNorm   -   2-norm of unmodified G
            MPGNorm -   2-norm of modified G
            MBA     -   minimum argument of barrier functions.
                        If X is strictly feasible, it is greater than zero.
                        If X lies on a boundary, it is zero.
                        It is negative for infeasible X.
            FIErr   -   2-norm of feasibility error with respect to
                        inequality/bound constraints
            CSErr   -   2-norm of complementarity slackness error
        *************************************************************************/
        private static void modifytargetfunction(minbleicstate state,
            double[] x,
            double[] r,
            double rnorm2,
            ref double f,
            ref double[] g,
            ref double gnorm,
            ref double mpgnorm)
        {
            double v = 0;
            int i = 0;
            int nmain = 0;
            int nslack = 0;
            bool hasconstraints = new bool();
            int i_ = 0;

            gnorm = 0;
            mpgnorm = 0;

            nmain = state.nmain;
            nslack = state.nslack;
            hasconstraints = false;
            
            //
            // GNorm
            //
            v = 0.0;
            for(i_=0; i_<=nmain+nslack-1;i_++)
            {
                v += g[i_]*g[i_];
            }
            gnorm = Math.Sqrt(v);
            
            //
            // Process equality constraints:
            // * modify F to handle penalty term for equality constraints
            // * project gradient on null space of equality constraints
            // * add penalty term for equality constraints to gradient
            //
            f = f+rnorm2;
            for(i=0; i<=nmain+nslack-1; i++)
            {
                if( state.activeconstraints[i] )
                {
                    g[i] = 0;
                }
            }
            for(i=0; i<=state.cecnt-1; i++)
            {
                v = 0.0;
                for(i_=0; i_<=nmain+nslack-1;i_++)
                {
                    v += g[i_]*state.cecurrent[i,i_];
                }
                for(i_=0; i_<=nmain+nslack-1;i_++)
                {
                    g[i_] = g[i_] - v*state.cecurrent[i,i_];
                }
            }
            for(i_=0; i_<=nmain+nslack-1;i_++)
            {
//.........这里部分代码省略.........

        /*************************************************************************
        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 13.10.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void minbleicsetprecdiag(minbleicstate state,
            double[] d)
        {
            int i = 0;

            alglib.ap.assert(alglib.ap.len(d)>=state.nmain, "MinBLEICSetPrecDiag: D is too short");
            for(i=0; i<=state.nmain-1; i++)
            {
                alglib.ap.assert(math.isfinite(d[i]), "MinBLEICSetPrecDiag: D contains infinite or NAN elements");
                alglib.ap.assert((double)(d[i])>(double)(0), "MinBLEICSetPrecDiag: D contains non-positive elements");
            }
            apserv.rvectorsetlengthatleast(ref state.diagh, state.nmain);
            state.prectype = 2;
            for(i=0; i<=state.nmain-1; i++)
            {
                state.diagh[i] = d[i];
            }
        }

        /*************************************************************************
                             BOUND CONSTRAINED OPTIMIZATION
               WITH ADDITIONAL LINEAR EQUALITY AND INEQUALITY CONSTRAINTS

        DESCRIPTION:
        The  subroutine  minimizes  function   F(x)  of N arguments subject to any
        combination of:
        * bound constraints
        * linear inequality constraints
        * linear equality constraints

        REQUIREMENTS:
        * user must provide function value and gradient
        * starting point X0 must be feasible or
          not too far away from the feasible set
        * grad(f) must be Lipschitz continuous on a level set:
          L = { x : f(x)<=f(x0) }
        * function must be defined everywhere on the feasible set F

        USAGE:

        Constrained optimization if far more complex than the unconstrained one.
        Here we give very brief outline of the BLEIC optimizer. We strongly recommend
        you to read examples in the ALGLIB Reference Manual and to read ALGLIB User Guide
        on optimization, which is available at http://www.alglib.net/optimization/

        1. User initializes algorithm state with MinBLEICCreate() call

        2. USer adds boundary and/or linear constraints by calling
           MinBLEICSetBC() and MinBLEICSetLC() functions.

        3. User sets stopping conditions with MinBLEICSetCond().

        4. User calls MinBLEICOptimize() function which takes algorithm  state and
           pointer (delegate, etc.) to callback function which calculates F/G.

        5. User calls MinBLEICResults() to get solution

        6. Optionally user may call MinBLEICRestartFrom() to solve another problem
           with same N but another starting point.
           MinBLEICRestartFrom() allows to reuse already initialized structure.


        INPUT PARAMETERS:
            N       -   problem dimension, N>0:
                        * if given, only leading N elements of X are used
                        * if not given, automatically determined from size ofX
            X       -   starting point, array[N]:
                        * it is better to set X to a feasible point
                        * but X can be infeasible, in which case algorithm will try
                          to find feasible point first, using X as initial
                          approximation.

        OUTPUT PARAMETERS:
            State   -   structure stores algorithm state

          -- ALGLIB --
             Copyright 28.11.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void minbleiccreate(int n,
            double[] x,
            minbleicstate state)
        {
            double[,] c = new double[0,0];
            int[] ct = new int[0];

            alglib.ap.assert(n>=1, "MinBLEICCreate: N<1");
            alglib.ap.assert(alglib.ap.len(x)>=n, "MinBLEICCreate: Length(X)<N");
            alglib.ap.assert(apserv.isfinitevector(x, n), "MinBLEICCreate: X contains infinite or NaN values!");
            minbleicinitinternal(n, x, 0.0, state);
        }

        /*************************************************************************
        Modification of the preconditioner: scale-based diagonal preconditioning.

        This preconditioning mode can be useful when you  don't  have  approximate
        diagonal of Hessian, but you know that your  variables  are  badly  scaled
        (for  example,  one  variable is in [1,10], and another in [1000,100000]),
        and most part of the ill-conditioning comes from different scales of vars.

        In this case simple  scale-based  preconditioner,  with H[i] = 1/(s[i]^2),
        can greatly improve convergence.

        IMPRTANT: you should set scale of your variables  with  MinBLEICSetScale()
        call  (before  or after MinBLEICSetPrecScale() call). Without knowledge of
        the scale of your variables scale-based preconditioner will be  just  unit
        matrix.

        INPUT PARAMETERS:
            State   -   structure which stores algorithm state

          -- ALGLIB --
             Copyright 13.10.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void minbleicsetprecscale(minbleicstate state)
        {
            state.prectype = 3;
        }

        /*************************************************************************
        The subroutine is finite difference variant of MinBLEICCreate().  It  uses
        finite differences in order to differentiate target function.

        Description below contains information which is specific to  this function
        only. We recommend to read comments on MinBLEICCreate() in  order  to  get
        more information about creation of BLEIC optimizer.

        INPUT PARAMETERS:
            N       -   problem dimension, N>0:
                        * if given, only leading N elements of X are used
                        * if not given, automatically determined from size of X
            X       -   starting point, array[0..N-1].
            DiffStep-   differentiation step, >0

        OUTPUT PARAMETERS:
            State   -   structure which stores algorithm state

        NOTES:
        1. algorithm uses 4-point central formula for differentiation.
        2. differentiation step along I-th axis is equal to DiffStep*S[I] where
           S[] is scaling vector which can be set by MinBLEICSetScale() call.
        3. we recommend you to use moderate values of  differentiation  step.  Too
           large step will result in too large truncation  errors, while too small
           step will result in too large numerical  errors.  1.0E-6  can  be  good
           value to start with.
        4. Numerical  differentiation  is   very   inefficient  -   one   gradient
           calculation needs 4*N function evaluations. This function will work for
           any N - either small (1...10), moderate (10...100) or  large  (100...).
           However, performance penalty will be too severe for any N's except  for
           small ones.
           We should also say that code which relies on numerical  differentiation
           is  less  robust and precise. CG needs exact gradient values. Imprecise
           gradient may slow  down  convergence, especially  on  highly  nonlinear
           problems.
           Thus  we  recommend to use this function for fast prototyping on small-
           dimensional problems only, and to implement analytical gradient as soon
           as possible.

          -- ALGLIB --
             Copyright 16.05.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void minbleiccreatef(int n,
            double[] x,
            double diffstep,
            minbleicstate state)
        {
            double[,] c = new double[0,0];
            int[] ct = new int[0];

            alglib.ap.assert(n>=1, "MinBLEICCreateF: N<1");
            alglib.ap.assert(alglib.ap.len(x)>=n, "MinBLEICCreateF: Length(X)<N");
            alglib.ap.assert(apserv.isfinitevector(x, n), "MinBLEICCreateF: X contains infinite or NaN values!");
            alglib.ap.assert(math.isfinite(diffstep), "MinBLEICCreateF: DiffStep is infinite or NaN!");
            alglib.ap.assert((double)(diffstep)>(double)(0), "MinBLEICCreateF: DiffStep is non-positive!");
            minbleicinitinternal(n, x, diffstep, state);
        }

        /*************************************************************************
        This function turns on/off line search reports.
        These reports are described in more details in developer-only  comments on
        MinBLEICState object.

        INPUT PARAMETERS:
            State   -   structure which stores algorithm state
            NeedDRep-   whether line search reports are needed or not

        This function is intended for private use only. Turning it on artificially
        may cause program failure.

          -- ALGLIB --
             Copyright 02.04.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void minbleicsetdrep(minbleicstate state,
            bool needdrep)
        {
            sta