C# GeneralMatrix.Subtract方法代码示例

        protected override GeneralMatrix CalculateNextHessianApproximation(GeneralMatrix previousH, 
			double[]prevX, double[]curX, double[]prevGrad, double[]curGrad)
        {
            GeneralMatrix currentH = new GeneralMatrix(_nDim,_nDim);
            GeneralMatrix cX = new GeneralMatrix(curX,_nDim);
            GeneralMatrix pX = new GeneralMatrix(prevX,_nDim);
            GeneralMatrix cG = new GeneralMatrix(curGrad,_nDim);
            GeneralMatrix pG = new GeneralMatrix(prevGrad,_nDim);

            GeneralMatrix dX = cX.Subtract(pX);
            GeneralMatrix dG = cG.Subtract(pG);

            double aK1 = 1/(dX.Transpose().Multiply(dG).GetElement(0,0));
            GeneralMatrix aK2 = dX.Multiply(dX.Transpose());

            GeneralMatrix aK = aK2.Multiply(aK1);

            double bK1 = -1/(dG.Transpose().Multiply(previousH).Multiply(dG).GetElement(0,0));
            GeneralMatrix bK2 = previousH.Multiply(dG).Multiply(dG.Transpose()).Multiply(previousH.Transpose());

            GeneralMatrix bK =bK2.Multiply(bK1);

            currentH = previousH.Add(aK).Add(bK);

            return currentH;
        }

        protected override GeneralMatrix CalculateNextHessianApproximation(GeneralMatrix pH, 
			double[]prevX, double[]curX, double[]prevGrad, double[]curGrad)
        {
            GeneralMatrix cH = new GeneralMatrix(_nDim,_nDim);
            GeneralMatrix cX = new GeneralMatrix(curX,_nDim);
            GeneralMatrix pX = new GeneralMatrix(prevX,_nDim);
            GeneralMatrix cG = new GeneralMatrix(curGrad,_nDim);
            GeneralMatrix pG = new GeneralMatrix(prevGrad,_nDim);

            GeneralMatrix sigma = cX.Subtract(pX);
            GeneralMatrix gamma = cG.Subtract(pG);

            double sigmaTGamma = sigma.Transpose().Multiply(gamma).GetElement(0,0);

            GeneralMatrix hGammaSigmaT = pH.Multiply(gamma.Multiply(sigma.Transpose()));
            GeneralMatrix sigmaGammaTH = sigma.Multiply(gamma.Transpose().Multiply(pH));
            double gammaTHGamma = (gamma.Transpose().Multiply(pH.Multiply(gamma))).GetElement(0,0);
            GeneralMatrix sigmaSigmaT  = sigma.Multiply(sigma.Transpose());

            GeneralMatrix term1 = (hGammaSigmaT.Add(sigmaGammaTH)).Multiply(1/sigmaTGamma);
            GeneralMatrix term2 = (sigmaSigmaT.Multiply(1/sigmaTGamma)).Multiply(1+gammaTHGamma/sigmaTGamma);

            return pH.Subtract(term1).Add(term2);
        }

        public static void Main(System.String[] argv)
        {
            GeneralMatrix A, B, C, Z, O, I, R, S, X, SUB, M, T, SQ, DEF, SOL;
            int errorCount = 0;
            int warningCount = 0;
            double tmp;
            double[] columnwise = new double[]{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
            double[] rowwise = new double[]{1.0, 4.0, 7.0, 10.0, 2.0, 5.0, 8.0, 11.0, 3.0, 6.0, 9.0, 12.0};
            double[][] avals = {new double[]{1.0, 4.0, 7.0, 10.0}, new double[]{2.0, 5.0, 8.0, 11.0}, new double[]{3.0, 6.0, 9.0, 12.0}};
            double[][] rankdef = avals;
            double[][] tvals = {new double[]{1.0, 2.0, 3.0}, new double[]{4.0, 5.0, 6.0}, new double[]{7.0, 8.0, 9.0}, new double[]{10.0, 11.0, 12.0}};
            double[][] subavals = {new double[]{5.0, 8.0, 11.0}, new double[]{6.0, 9.0, 12.0}};
            double[][] rvals = {new double[]{1.0, 4.0, 7.0}, new double[]{2.0, 5.0, 8.0, 11.0}, new double[]{3.0, 6.0, 9.0, 12.0}};
            double[][] pvals = {new double[]{1.0, 1.0, 1.0}, new double[]{1.0, 2.0, 3.0}, new double[]{1.0, 3.0, 6.0}};
            double[][] ivals = {new double[]{1.0, 0.0, 0.0, 0.0}, new double[]{0.0, 1.0, 0.0, 0.0}, new double[]{0.0, 0.0, 1.0, 0.0}};
            double[][] evals = {new double[]{0.0, 1.0, 0.0, 0.0}, new double[]{1.0, 0.0, 2e-7, 0.0}, new double[]{0.0, - 2e-7, 0.0, 1.0}, new double[]{0.0, 0.0, 1.0, 0.0}};
            double[][] square = {new double[]{166.0, 188.0, 210.0}, new double[]{188.0, 214.0, 240.0}, new double[]{210.0, 240.0, 270.0}};
            double[][] sqSolution = {new double[]{13.0}, new double[]{15.0}};
            double[][] condmat = {new double[]{1.0, 3.0}, new double[]{7.0, 9.0}};
            int rows = 3, cols = 4;
            int invalidld = 5; /* should trigger bad shape for construction with val */
            int raggedr = 0; /* (raggedr,raggedc) should be out of bounds in ragged array */
            int raggedc = 4;
            int validld = 3; /* leading dimension of intended test Matrices */
            int nonconformld = 4; /* leading dimension which is valid, but nonconforming */
            int ib = 1, ie = 2, jb = 1, je = 3; /* index ranges for sub GeneralMatrix */
            int[] rowindexset = new int[]{1, 2};
            int[] badrowindexset = new int[]{1, 3};
            int[] columnindexset = new int[]{1, 2, 3};
            int[] badcolumnindexset = new int[]{1, 2, 4};
            double columnsummax = 33.0;
            double rowsummax = 30.0;
            double sumofdiagonals = 15;
            double sumofsquares = 650;

            /// <summary>Constructors and constructor-like methods:
            /// double[], int
            /// double[][]
            /// int, int
            /// int, int, double
            /// int, int, double[][]
            /// Create(double[][])
            /// Random(int,int)
            /// Identity(int)
            ///
            /// </summary>

            print("\nTesting constructors and constructor-like methods...\n");
            try
            {
                /// <summary>check that exception is thrown in packed constructor with invalid length *</summary>
                A = new GeneralMatrix(columnwise, invalidld);
                errorCount = try_failure(errorCount, "Catch invalid length in packed constructor... ", "exception not thrown for invalid input");
            }
            catch (System.ArgumentException e)
            {
                try_success("Catch invalid length in packed constructor... ", e.Message);
            }
            try
            {
                /// <summary>check that exception is thrown in default constructor
                /// if input array is 'ragged' *
                /// </summary>
                A = new GeneralMatrix(rvals);
                tmp = A.GetElement(raggedr, raggedc);
            }
            catch (System.ArgumentException e)
            {
                try_success("Catch ragged input to default constructor... ", e.Message);
            }
            catch (System.IndexOutOfRangeException e)
            {
                errorCount = try_failure(errorCount, "Catch ragged input to constructor... ", "exception not thrown in construction...ArrayIndexOutOfBoundsException thrown later");
                System.Console.Out.WriteLine(e.Message);
            }
            try
            {
                /// <summary>check that exception is thrown in Create
                /// if input array is 'ragged' *
                /// </summary>
                A = GeneralMatrix.Create(rvals);
                tmp = A.GetElement(raggedr, raggedc);
            }
            catch (System.ArgumentException e)
            {
                try_success("Catch ragged input to Create... ", e.Message);
                System.Console.Out.WriteLine(e.Message);
            }
            catch (System.IndexOutOfRangeException e)
            {
                errorCount = try_failure(errorCount, "Catch ragged input to Create... ", "exception not thrown in construction...ArrayIndexOutOfBoundsException thrown later");
                System.Console.Out.WriteLine(e.Message);
            }

            A = new GeneralMatrix(columnwise, validld);
            B = new GeneralMatrix(avals);
            tmp = B.GetElement(0, 0);
            avals[0][0] = 0.0;
            C = B.Subtract(A);
            avals[0][0] = tmp;
//.........这里部分代码省略.........

 /// <summary>Check norm of difference of Matrices. *</summary>
 private static void check(GeneralMatrix X, GeneralMatrix Y)
 {
     double eps = System.Math.Pow(2.0, - 52.0);
     if (X.Norm1() == 0.0 & Y.Norm1() < 10 * eps)
         return ;
     if (Y.Norm1() == 0.0 & X.Norm1() < 10 * eps)
         return ;
     if (X.Subtract(Y).Norm1() > 1000 * eps * System.Math.Max(X.Norm1(), Y.Norm1()))
     {
         throw new System.SystemException("The norm of (X-Y) is too large: " + X.Subtract(Y).Norm1().ToString());
     }
 }

        public void Substraction()
        {
            double[] columnwise = new double[]{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
            double[][] avals = {new double[]{1.0, 4.0, 7.0, 10.0}, new double[]{2.0, 5.0, 8.0, 11.0}, new double[]{3.0, 6.0, 9.0, 12.0}};

            int validld = 3; /* leading dimension of intended test Matrices */

            //one dimensional array of doubles packed by columns ala Fortran
            //is passed to the constructor
            //the integer value tell the constructor how to
            //breakup one domensional value into mulidimensional column
            GeneralMatrix A = new GeneralMatrix(columnwise, validld);
            GeneralMatrix B = new GeneralMatrix(avals);
            double tmp = B.GetElement(0, 0);
            avals[0][0] = 0.0;
            GeneralMatrix C = B.Subtract(A);
            avals[0][0] = tmp;
            B = GeneralMatrix.Create(avals);
            tmp = B.GetElement(0, 0);
            avals[0][0] = 0.0;
            Assert.IsTrue(tmp==B.GetElement(0, 0));
        }

        /// <summary>Check norm of difference of Matrices.
        /// </summary>	
        public static bool Check(GeneralMatrix X, GeneralMatrix Y)
        {
            bool result=false;

            double eps = System.Math.Pow(2.0, - 52.0);
            if (X.Norm1() == 0.0 & Y.Norm1() < 10 * eps)
                result = true ;
            else if (Y.Norm1() == 0.0 & X.Norm1() < 10 * eps)
                result = true ;
            else if (X.Subtract(Y).Norm1() > 1000 * eps * System.Math.Max(X.Norm1(), Y.Norm1()))
            {
                throw new System.SystemException("The norm of (X-Y) is too large: " + X.Subtract(Y).Norm1().ToString());
            }
            else result = true;

            return result;
        }

 //General
 /// <summary>Check norm of difference of Matrices. *</summary>
 private static bool check(GeneralMatrix x, GeneralMatrix y)
 {
     double eps = Math.Pow(2.0, -52.0);
     if (x.Norm1() == 0.0 & y.Norm1() < 10 * eps)
         return true;
     if (y.Norm1() == 0.0 & x.Norm1() < 10 * eps)
         return true;
     return (x.Subtract(y).Norm1() > 1000*eps*Math.Max(x.Norm1(), y.Norm1()));
 }

 public void InitializeArrays()
 {
     B = new GeneralMatrix(avals);
     M = new GeneralMatrix(2, 3, 0.0);
     S = new GeneralMatrix(columnwise, nonconformld);
     R = GeneralMatrix.Random(rows, cols);
     Z = new GeneralMatrix(rows, cols);
     A = R.Copy();
     C = A.Subtract(B);
     O = new GeneralMatrix(rows,cols,1.0);
     twos = new GeneralMatrix(rows,cols,2.0);
 }

 public void Substract()
 {
     A=R;
     Assert.AreEqual(0.0,A.Subtract(R).Norm1());
 }