C# MatrixFixed.Datablock方法代码示例

        public unsafe QR(MatrixFixed M)
        {
            qrdc_out_ = new MatrixFixed(M.Columns, M.Rows);
            qraux_ = new Vector(M.Columns);
            jpvt_ = new int[M.Rows];
            Q_ = null;
            R_ = null;

            // Fill transposed O/P matrix
            int c = M.Columns;
            int r = M.Rows;
            for (int i = 0; i < r; ++i)
                for (int j = 0; j < c; ++j)
                    qrdc_out_[j,i] = M[i,j];

            int do_pivot = 0; // Enable[!=0]/disable[==0] pivoting.
            for (int i = 0; i < jpvt_.Length; i++) jpvt_[i] = 0;

            Vector work = new Vector(M.Rows);

            fixed (float* data = qrdc_out_.Datablock())
            {
                fixed (float* data2 = qraux_.Datablock())
                {
                    fixed (int* data3 = jpvt_)
                    {
                        fixed (float* data4 = work.Datablock())
                        {
                            Netlib.dqrdc_(data,       // On output, UT is R, below diag is mangled Q
                                          &r, &r, &c,
                                          data2,      // Further information required to demangle Q
                                          data3,
                                          data4,
                                          &do_pivot);
                        }
                    }
                }
            }
        }

        //subtraction
        public static MatrixFixed operator -(MatrixFixed lmtx, MatrixFixed rmtx)
        {
            int left_rows = lmtx.Rows;
            int left_columns = lmtx.Columns;
            int right_rows = rmtx.Rows;
            int right_columns = rmtx.Columns;

            if (left_rows != right_rows || left_columns != right_columns)
                throw new MatrixFixedException("Attempt to subtract matrixes of different sizes");
            //return null;

            MatrixFixed result = new MatrixFixed(left_rows, left_columns);

            float[,] result_matrix = result.Datablock();
            float[,] left_matrix = lmtx.Datablock();
            float[,] right_matrix = rmtx.Datablock();

            for (int i = 0; i < left_rows; i++)
                for (int j = 0; j < left_columns; j++)
                    result_matrix[i, j] = left_matrix[i, j] - right_matrix[i, j];

            return result;
        }

        //devision by const
        public static MatrixFixed operator /(MatrixFixed mtx, float val)
        {
            int rows = mtx.Rows;
            int cols = mtx.Columns;

            MatrixFixed result = new MatrixFixed(rows, cols);

            float[,] result_matrix = result.Datablock();
            float[,] matrix_data = mtx.Datablock();

            if (val == 0)
                throw new MatrixFixedException("Attempt to devide the matrix by zero");
            //return null;

            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < cols; j++)
                {
                    result_matrix[i, j] = matrix_data[i, j] / val;
                }
            }

            return result;
        }

        /// <summary>
        /// Cholesky decomposition.
        /// Make cholesky decomposition of M optionally computing
        /// the reciprocal condition number.  If mode is estimate_condition, the
        /// condition number and an approximate nullspace are estimated, at a cost
        /// of a factor of (1 + 18/n).  Here's a table of 1 + 18/n:
        ///<pre>
        /// n:              3      5     10     50    100    500   1000
        /// slowdown:     7.0f    4.6    2.8    1.4   1.18   1.04   1.02
        /// </summary>
        /// <param name="M"></param>
        /// <param name="mode"></param>
        public unsafe void init(MatrixFixed M, Operation mode)
        {
            A_ = new MatrixFixed(M);

            int n = M.Columns;
            //assert(n == (int)(M.Rows()));
            num_dims_rank_def_ = -1;
            int num_dims_rank_def_temp = num_dims_rank_def_;

            // BJT: This warning is pointless - it often doesn't detect non symmetry and
            // if you know what you're doing you don't want to be slowed down
            // by a cerr
            /*
               if (Math.Abs(M[0,n-1] - M[n-1,0]) > 1e-8) 
               {
                   Debug.WriteLine("cholesky: WARNING: unsymmetric: " + M);
               }
            */

            if (mode != Operation.estimate_condition) 
            {
                // Quick factorization
                fixed (float* data = A_.Datablock())
                {
                    Netlib.dpofa_(data, &n, &n, &num_dims_rank_def_temp);                    
                }
                //if ((mode == Operation.verbose) && (num_dims_rank_def_temp != 0))
                //    Debug.WriteLine("cholesky:: " + Convert.ToString(num_dims_rank_def_temp) + " dimensions of non-posdeffness");
            } 
            else 
            {
                Vector nullvector = new Vector(n);
                float rcond_temp = rcond_;
                fixed (float* data = A_.Datablock())
                {
                    fixed (float* data2 = nullvector.Datablock())
                    {
                        Netlib.dpoco_(data, &n, &n, &rcond_temp, data2, &num_dims_rank_def_temp);
                    }
                }
                rcond_ = rcond_temp;
            }
            num_dims_rank_def_ = num_dims_rank_def_temp;
        }

        /// <summary>
        /// Compute inverse.\  Not efficient.
        /// </summary>
        /// <returns></returns>
        public unsafe MatrixFixed Inverse()
        {
            int n = A_.Columns;
            MatrixFixed I = new MatrixFixed(A_);
            float[] det = new float[2];
            int job = 01;
            fixed (float* data = I.Datablock())
            {
                fixed (float* data2 = det)
                {
                    Netlib.dpodi_(data, &n, &n, data2, &job);
                }
            }

            // Copy lower triangle into upper
            for (int i = 0; i < n; ++i)
                for (int j = i+1; j < n; ++j)
                    I[i,j] = I[j,i];

            return I;
        }

 /// <summary>
 /// Compute determinant.
 /// </summary>
 /// <returns></returns>
 public unsafe float Determinant()
 {
     int n = A_.Columns;
     MatrixFixed I = new MatrixFixed(A_);
     float[] det = new float[2];
     int job = 10;
     fixed (float* data = I.Datablock())
     {
         fixed (float* data2 = det)
         {
             Netlib.dpodi_(data, &n, &n, data2, &job);
         }
     }
     return det[0] * (float)Math.Pow(10.0, det[1]);
 }