C# lincgstate类代码示例

 /*************************************************************************
 Clears request fileds (to be sure that we don't forgot to clear something)
 *************************************************************************/
 private static void updateitersdata(lincgstate state)
 {
     state.repiterationscount = 0;
     state.repnmv = 0;
     state.repterminationtype = 0;
 }

        /*************************************************************************
        Procedure for restart function LinCGIteration

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgrestart(lincgstate state)
        {
            state.rstate.ia = new int[0+1];
            state.rstate.ra = new double[2+1];
            state.rstate.stage = -1;
            clearrfields(state);
        }

        /*************************************************************************
        This function sets frequency of residual recalculations.

        Algorithm updates residual r_k using iterative formula,  but  recalculates
        it from scratch after each 10 iterations. It is done to avoid accumulation
        of numerical errors and to stop algorithm when r_k starts to grow.

        Such low update frequence (1/10) gives very  little  overhead,  but  makes
        algorithm a bit more robust against numerical errors. However, you may
        change it 

        INPUT PARAMETERS:
            Freq    -   desired update frequency, Freq>=0.
                        Zero value means that no updates will be done.

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgsetrupdatefreq(lincgstate state,
            int freq)
        {
            alglib.ap.assert(!state.running, "LinCGSetRUpdateFreq: you can not change update frequency when LinCGIteration() is running");
            alglib.ap.assert(freq>=0, "LinCGSetRUpdateFreq: non-positive Freq");
            state.itsbeforerupdate = freq;
        }

        /*************************************************************************
        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 MinCGOptimize().

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgsetxrep(lincgstate state,
            bool needxrep)
        {
            state.xrep = needxrep;
        }

        /*************************************************************************
        CG-solver: results.

        This function must be called after LinCGSolve

        INPUT PARAMETERS:
            State   -   algorithm state

        OUTPUT PARAMETERS:
            X       -   array[N], solution
            Rep     -   optimization report:
                        * Rep.TerminationType completetion code:
                            * -5    input matrix is either not positive definite,
                                    too large or too small                            
                            * -4    overflow/underflow during solution
                                    (ill conditioned problem)
                            *  1    ||residual||<=EpsF*||b||
                            *  5    MaxIts steps was taken
                            *  7    rounding errors prevent further progress,
                                    best point found is returned
                        * Rep.IterationsCount contains iterations count
                        * NMV countains number of matrix-vector calculations

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgresults(lincgstate state,
            ref double[] x,
            lincgreport rep)
        {
            int i_ = 0;

            x = new double[0];

            alglib.ap.assert(!state.running, "LinCGResult: you can not get result, because function LinCGIteration has been launched!");
            if( alglib.ap.len(x)<state.n )
            {
                x = new double[state.n];
            }
            for(i_=0; i_<=state.n-1;i_++)
            {
                x[i_] = state.rx[i_];
            }
            rep.iterationscount = state.repiterationscount;
            rep.nmv = state.repnmv;
            rep.terminationtype = state.repterminationtype;
            rep.r2 = state.r2;
        }

        /*************************************************************************
        This function sets restart frequency. By default, algorithm  is  restarted
        after N subsequent iterations.

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgsetrestartfreq(lincgstate state,
            int srf)
        {
            alglib.ap.assert(!state.running, "LinCGSetRestartFreq: you can not change restart frequency when LinCGIteration() is running");
            alglib.ap.assert(srf>0, "LinCGSetRestartFreq: non-positive SRF");
            state.itsbeforerestart = srf;
        }

        /*************************************************************************
        Reverse communication version of linear CG.

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static bool lincgiteration(lincgstate state)
        {
            bool result = new bool();
            int i = 0;
            double uvar = 0;
            double bnorm = 0;
            double v = 0;
            int i_ = 0;

            
            //
            // Reverse communication preparations
            // I know it looks ugly, but it works the same way
            // anywhere from C++ to Python.
            //
            // This code initializes locals by:
            // * random values determined during code
            //   generation - on first subroutine call
            // * values from previous call - on subsequent calls
            //
            if( state.rstate.stage>=0 )
            {
                i = state.rstate.ia[0];
                uvar = state.rstate.ra[0];
                bnorm = state.rstate.ra[1];
                v = state.rstate.ra[2];
            }
            else
            {
                i = -983;
                uvar = -989;
                bnorm = -834;
                v = 900;
            }
            if( state.rstate.stage==0 )
            {
                goto lbl_0;
            }
            if( state.rstate.stage==1 )
            {
                goto lbl_1;
            }
            if( state.rstate.stage==2 )
            {
                goto lbl_2;
            }
            if( state.rstate.stage==3 )
            {
                goto lbl_3;
            }
            if( state.rstate.stage==4 )
            {
                goto lbl_4;
            }
            if( state.rstate.stage==5 )
            {
                goto lbl_5;
            }
            if( state.rstate.stage==6 )
            {
                goto lbl_6;
            }
            if( state.rstate.stage==7 )
            {
                goto lbl_7;
            }
            
            //
            // Routine body
            //
            alglib.ap.assert(alglib.ap.len(state.b)>0, "LinCGIteration: B is not initialized (you must initialize B by LinCGSetB() call");
            state.running = true;
            state.repnmv = 0;
            clearrfields(state);
            updateitersdata(state);
            
            //
            // Start 0-th iteration
            //
            for(i_=0; i_<=state.n-1;i_++)
            {
                state.rx[i_] = state.startx[i_];
            }
            for(i_=0; i_<=state.n-1;i_++)
            {
                state.x[i_] = state.rx[i_];
            }
            state.repnmv = state.repnmv+1;
            clearrfields(state);
            state.needvmv = true;
            state.rstate.stage = 0;
            goto lbl_rcomm;
        lbl_0:
            state.needvmv = false;
//.........这里部分代码省略.........

        /*************************************************************************
        Procedure for solution of A*x=b with sparse A.

        INPUT PARAMETERS:
            State   -   algorithm state
            A       -   sparse matrix in the CRS format (you MUST contvert  it  to 
                        CRS format by calling SparseConvertToCRS() function).
            IsUpper -   whether upper or lower triangle of A is used:
                        * IsUpper=True  => only upper triangle is used and lower
                                           triangle is not referenced at all 
                        * IsUpper=False => only lower triangle is used and upper
                                           triangle is not referenced at all
            B       -   right part, array[N]

        RESULT:
            This function returns no result.
            You can get solution by calling LinCGResults()
            
        NOTE: this function uses lightweight preconditioning -  multiplication  by
              inverse of diag(A). If you want, you can turn preconditioning off by
              calling LinCGSetPrecUnit(). However, preconditioning cost is low and
              preconditioner  is  very  important  for  solution  of  badly scaled
              problems.

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgsolvesparse(lincgstate state,
            sparse.sparsematrix a,
            bool isupper,
            double[] b)
        {
            int n = 0;
            int i = 0;
            double v = 0;
            double vmv = 0;
            int i_ = 0;

            n = state.n;
            alglib.ap.assert(alglib.ap.len(b)>=state.n, "LinCGSetB: Length(B)<N");
            alglib.ap.assert(apserv.isfinitevector(b, state.n), "LinCGSetB: B contains infinite or NaN values!");
            
            //
            // Allocate temporaries
            //
            apserv.rvectorsetlengthatleast(ref state.tmpd, n);
            
            //
            // Compute diagonal scaling matrix D
            //
            if( state.prectype==0 )
            {
                
                //
                // Default preconditioner - inverse of matrix diagonal
                //
                for(i=0; i<=n-1; i++)
                {
                    v = sparse.sparsegetdiagonal(a, i);
                    if( (double)(v)>(double)(0) )
                    {
                        state.tmpd[i] = 1/Math.Sqrt(v);
                    }
                    else
                    {
                        state.tmpd[i] = 1;
                    }
                }
            }
            else
            {
                
                //
                // No diagonal scaling
                //
                for(i=0; i<=n-1; i++)
                {
                    state.tmpd[i] = 1;
                }
            }
            
            //
            // Solve
            //
            lincgrestart(state);
            lincgsetb(state, b);
            while( lincgiteration(state) )
            {
                
                //
                // Process different requests from optimizer
                //
                if( state.needmv )
                {
                    sparse.sparsesmv(a, isupper, state.x, ref state.mv);
                }
                if( state.needvmv )
                {
                    sparse.sparsesmv(a, isupper, state.x, ref state.mv);
                    vmv = 0.0;
//.........这里部分代码省略.........

        /*************************************************************************
        This  function  changes  preconditioning  settings  of  LinCGSolveSparse()
        function.  LinCGSolveSparse() will use diagonal of the  system  matrix  as
        preconditioner. This preconditioning mode is active by default.

        INPUT PARAMETERS:
            State   -   structure which stores algorithm state

          -- ALGLIB --
             Copyright 19.11.2012 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgsetprecdiag(lincgstate state)
        {
            alglib.ap.assert(!state.running, "LinCGSetPrecDiag: you can not change preconditioner, because function LinCGIteration is running!");
            state.prectype = 0;
        }

        /*************************************************************************
        This function sets stopping criteria.

        INPUT PARAMETERS:
            EpsF    -   algorithm will be stopped if norm of residual is less than 
                        EpsF*||b||.
            MaxIts  -   algorithm will be stopped if number of iterations is  more 
                        than MaxIts.

        OUTPUT PARAMETERS:
            State   -   structure which stores algorithm state

        NOTES:
        If  both  EpsF  and  MaxIts  are  zero then small EpsF will be set to small 
        value.

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgsetcond(lincgstate state,
            double epsf,
            int maxits)
        {
            alglib.ap.assert(!state.running, "LinCGSetCond: you can not change stopping criteria when LinCGIteration() is running");
            alglib.ap.assert(math.isfinite(epsf) && (double)(epsf)>=(double)(0), "LinCGSetCond: EpsF is negative or contains infinite or NaN values");
            alglib.ap.assert(maxits>=0, "LinCGSetCond: MaxIts is negative");
            if( (double)(epsf)==(double)(0) && maxits==0 )
            {
                state.epsf = defaultprecision;
                state.maxits = maxits;
            }
            else
            {
                state.epsf = epsf;
                state.maxits = maxits;
            }
        }

        /*************************************************************************
        This function sets right part. By default, right part is zero.

        INPUT PARAMETERS:
            B       -   right part, array[N].

        OUTPUT PARAMETERS:
            State   -   structure which stores algorithm state

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgsetb(lincgstate state,
            double[] b)
        {
            int i_ = 0;

            alglib.ap.assert(!state.running, "LinCGSetB: you can not set B, because function LinCGIteration is running!");
            alglib.ap.assert(alglib.ap.len(b)>=state.n, "LinCGSetB: Length(B)<N");
            alglib.ap.assert(apserv.isfinitevector(b, state.n), "LinCGSetB: B contains infinite or NaN values!");
            for(i_=0; i_<=state.n-1;i_++)
            {
                state.b[i_] = b[i_];
            }
        }

        /*************************************************************************
        This function sets starting point.
        By default, zero starting point is used.

        INPUT PARAMETERS:
            X       -   starting point, array[N]

        OUTPUT PARAMETERS:
            State   -   structure which stores algorithm state

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgsetstartingpoint(lincgstate state,
            double[] x)
        {
            int i_ = 0;

            alglib.ap.assert(!state.running, "LinCGSetStartingPoint: you can not change starting point because LinCGIteration() function is running");
            alglib.ap.assert(state.n<=alglib.ap.len(x), "LinCGSetStartingPoint: Length(X)<N");
            alglib.ap.assert(apserv.isfinitevector(x, state.n), "LinCGSetStartingPoint: X contains infinite or NaN values!");
            for(i_=0; i_<=state.n-1;i_++)
            {
                state.startx[i_] = x[i_];
            }
        }

        /*************************************************************************
        This function initializes linear CG Solver. This solver is used  to  solve
        symmetric positive definite problems. If you want  to  solve  nonsymmetric
        (or non-positive definite) problem you may use LinLSQR solver provided  by
        ALGLIB.

        USAGE:
        1. User initializes algorithm state with LinCGCreate() call
        2. User tunes solver parameters with  LinCGSetCond() and other functions
        3. Optionally, user sets starting point with LinCGSetStartingPoint()
        4. User  calls LinCGSolveSparse() function which takes algorithm state and
           SparseMatrix object.
        5. User calls LinCGResults() to get solution
        6. Optionally, user may call LinCGSolveSparse()  again  to  solve  another
           problem  with different matrix and/or right part without reinitializing
           LinCGState structure.
          
        INPUT PARAMETERS:
            N       -   problem dimension, N>0

        OUTPUT PARAMETERS:
            State   -   structure which stores algorithm state

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgcreate(int n,
            lincgstate state)
        {
            int i = 0;

            alglib.ap.assert(n>0, "LinCGCreate: N<=0");
            state.n = n;
            state.prectype = 0;
            state.itsbeforerestart = n;
            state.itsbeforerupdate = 10;
            state.epsf = defaultprecision;
            state.maxits = 0;
            state.xrep = false;
            state.running = false;
            
            //
            // * allocate arrays
            // * set RX to NAN (just for the case user calls Results() without 
            //   calling SolveSparse()
            // * set starting point to zero
            // * we do NOT initialize B here because we assume that user should
            //   initializate it using LinCGSetB() function. In case he forgets
            //   to do so, exception will be thrown in the LinCGIteration().
            //
            state.rx = new double[state.n];
            state.startx = new double[state.n];
            state.b = new double[state.n];
            for(i=0; i<=state.n-1; i++)
            {
                state.rx[i] = Double.NaN;
                state.startx[i] = 0.0;
                state.b[i] = 0;
            }
            state.cx = new double[state.n];
            state.p = new double[state.n];
            state.r = new double[state.n];
            state.cr = new double[state.n];
            state.z = new double[state.n];
            state.cz = new double[state.n];
            state.x = new double[state.n];
            state.mv = new double[state.n];
            state.pv = new double[state.n];
            updateitersdata(state);
            state.rstate.ia = new int[0+1];
            state.rstate.ra = new double[2+1];
            state.rstate.stage = -1;
        }

 public override alglib.apobject make_copy()
 {
     lincgstate _result = new lincgstate();
     _result.rx = (double[])rx.Clone();
     _result.b = (double[])b.Clone();
     _result.n = n;
     _result.prectype = prectype;
     _result.cx = (double[])cx.Clone();
     _result.cr = (double[])cr.Clone();
     _result.cz = (double[])cz.Clone();
     _result.p = (double[])p.Clone();
     _result.r = (double[])r.Clone();
     _result.z = (double[])z.Clone();
     _result.alpha = alpha;
     _result.beta = beta;
     _result.r2 = r2;
     _result.meritfunction = meritfunction;
     _result.x = (double[])x.Clone();
     _result.mv = (double[])mv.Clone();
     _result.pv = (double[])pv.Clone();
     _result.vmv = vmv;
     _result.startx = (double[])startx.Clone();
     _result.epsf = epsf;
     _result.maxits = maxits;
     _result.itsbeforerestart = itsbeforerestart;
     _result.itsbeforerupdate = itsbeforerupdate;
     _result.xrep = xrep;
     _result.xupdated = xupdated;
     _result.needmv = needmv;
     _result.needmtv = needmtv;
     _result.needmv2 = needmv2;
     _result.needvmv = needvmv;
     _result.needprec = needprec;
     _result.repiterationscount = repiterationscount;
     _result.repnmv = repnmv;
     _result.repterminationtype = repterminationtype;
     _result.running = running;
     _result.tmpd = (double[])tmpd.Clone();
     _result.rstate = (rcommstate)rstate.make_copy();
     return _result;
 }

        /*************************************************************************
        Procedure for solution of A*x=b with sparse A.

        INPUT PARAMETERS:
            State   -   algorithm state
            A       -   sparse matrix in the CRS format (you MUST contvert  it  to 
                        CRS format by calling SparseConvertToCRS() function).
            IsUpper -   whether upper or lower triangle of A is used:
                        * IsUpper=True  => only upper triangle is used and lower
                                           triangle is not referenced at all 
                        * IsUpper=False => only lower triangle is used and upper
                                           triangle is not referenced at all
            B       -   right part, array[N]

        RESULT:
            This function returns no result.
            You can get solution by calling LinCGResults()

          -- ALGLIB --
             Copyright 14.11.2011 by Bochkanov Sergey
        *************************************************************************/
        public static void lincgsolvesparse(lincgstate state,
            sparse.sparsematrix a,
            bool isupper,
            double[] b)
        {
            double vmv = 0;
            int i_ = 0;

            alglib.ap.assert(alglib.ap.len(b)>=state.n, "LinCGSetB: Length(B)<N");
            alglib.ap.assert(apserv.isfinitevector(b, state.n), "LinCGSetB: B contains infinite or NaN values!");
            lincgrestart(state);
            lincgsetb(state, b);
            while( lincgiteration(state) )
            {
                if( state.needmv )
                {
                    sparse.sparsesmv(a, isupper, state.x, ref state.mv);
                }
                if( state.needvmv )
                {
                    sparse.sparsesmv(a, isupper, state.x, ref state.mv);
                    vmv = 0.0;
                    for(i_=0; i_<=state.n-1;i_++)
                    {
                        vmv += state.x[i_]*state.mv[i_];
                    }
                    state.vmv = vmv;
                }
                if( state.needprec )
                {
                    for(i_=0; i_<=state.n-1;i_++)
                    {
                        state.pv[i_] = state.x[i_];
                    }
                }
            }
        }