C# Vector.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);
                        }
                    }
                }
            }
        }

        private unsafe void init(MatrixFixed M, float zero_out_tol)
        {
            m_ = M.Rows;
            n_ = M.Columns;
            U_ = new MatrixFixed(m_, n_);
            W_ = new DiagMatrix(n_);
            Winverse_ = new DiagMatrix(n_);
            V_ = new MatrixFixed(n_, n_);

            //assert(m_ > 0);  
            //assert(n_ > 0);
		
            int n = M.Rows;    
            int p = M.Columns;
            int mm = Netlib.min(n+1,p);

            // Copy source matrix into fortran storage
            // SVD is slow, don't worry about the cost of this transpose.
            Vector X = Vector.fortran_copy(M);

            // Make workspace vectors
            Vector work = new Vector(n);
            work.Fill(0);
            Vector uspace = new Vector(n*p);
            uspace.Fill(0);
            Vector vspace = new Vector(p*p);
            vspace.Fill(0);
            Vector wspace = new Vector(mm);
            wspace.Fill(0); // complex fortran routine actually _wants_ complex W!
            Vector espace = new Vector(p);
            espace.Fill(0);
    
            // Call Linpack SVD
            int info = 0;
            int job = 21;

            fixed (float* data = X.Datablock())
            {
                fixed (float* data2 = wspace.Datablock())
                {
                    fixed (float* data3 = espace.Datablock())
                    {
                        fixed (float* data4 = uspace.Datablock())
                        {
                            fixed (float* data5 = vspace.Datablock())
                            {
                                fixed (float* data6 = work.Datablock())
                                {
                                    Netlib.dsvdc_(data, &n, &n, &p,
                                             data2,
                                             data3,
                                             data4, &n,
                                             data5, &p,
                                             data6,
                                             &job, &info);
                                }
                            }
                        }
                    }
                }
            }

            // Error return?
            if (info != 0) 
            {
                // If info is non-zero, it contains the number of singular values
                // for this the SVD algorithm failed to converge. The condition is
                // not bogus. Even if the returned singular values are sensible,
                // the singular vectors can be utterly wrong.

                // It is possible the failure was due to NaNs or infinities in the
                // matrix. Check for that now.
                M.assert_finite();

                // If we get here it might be because
                // 1. The scalar type has such
                // extreme precision that too few iterations were performed to
                // converge to within machine precision (that is the svdc criterion).
                // One solution to that is to increase the maximum number of
                // iterations in the netlib code.
                //
                // 2. The LINPACK dsvdc_ code expects correct IEEE rounding behaviour,
                // which some platforms (notably x86 processors)
                // have trouble doing. For example, gcc can output
                // code in -O2 and static-linked code that causes this problem.
                // One solution to this is to persuade gcc to output slightly different code
                // by adding and -fPIC option to the command line for v3p\netlib\dsvdc.c. If
                // that doesn't work try adding -ffloat-store, which should fix the problem
                // at the expense of being significantly slower for big problems. Note that
                // if this is the cause, vxl/vnl/tests/test_svd should have failed.
                //
                // You may be able to diagnose the problem here by printing a warning message.
                Debug.WriteLine("__FILE__ : suspicious return value (" + Convert.ToString(info) + ") from SVDC" +
                                "__FILE__ : M is " + Convert.ToString(M.Rows) + "x" + Convert.ToString(M.Columns));

                valid_ = false;
            }
            else
                valid_ = true;

//.........这里部分代码省略.........

        /// <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>
        /// Solve least squares problem M x = b.
        /// </summary>
        /// <param name="b"></param>
        /// <returns></returns>
        public unsafe Vector Solve(Vector b)
        {
            //assert(b.size() == A_.Columns());

            int n = A_.Columns;
            Vector ret = new Vector(b);
            fixed (float* data = A_.Datablock())
            {
                fixed (float* data2 = ret.Datablock())
                {
                    Netlib.dposl_(data, &n, &n, data2);
                }
            }
            return ret;
        }

        public static Vector operator* (MatrixFixed m, Vector v) 
        {
            Vector result = new Vector(m.Rows);        // Temporary
            MatrixFixed mm = m;                        // Drop const for get()
            float[] result_data = result.Datablock();
            float[] v_data = v.Datablock();
            float[,] mm_data = mm.Datablock();
            int vsize = v.size();

            for (int i = 0; i < m.Rows; i++) 
            {                                                     // For each index
                result_data[i] = 0;                               // Initialize element value
                for (int k = 0; k < vsize; k++)                // Loop over column values
                    result_data[i] += (mm_data[i,k] * v_data[k]); // Multiply
            }
            return result;
        }

        /// <summary>
        /// Return residual vector d of M x = b -> d = Q'b.
        /// </summary>
        /// <param name="b"></param>
        /// <returns></returns>
        public unsafe Vector QtB(Vector b)
        {
            int n = qrdc_out_.Columns;
            int p = qrdc_out_.Rows;
            float[] b_data = b.Datablock();
            Vector QtB = new Vector(n);

            // see comment above
            int JOB = 1000;

            int info = 0;

            fixed (float* data = qrdc_out_.Datablock())
            {
                fixed (float* data2 = qraux_.Datablock())
                {
                    fixed (float* data3 = b_data)
                    {
                        fixed (float* data4 = QtB.Datablock())
                        {
                            Netlib.dqrsl_(data, &n, &n, &p, data2, data3,
                                   (float*)0,         // A: Qb
                                   data4,              // B: Q'b
                                   (float*)0,         // C: x
                                   (float*)0,         // D: residual
                                   (float*)0,         // E: Ax
                                   &JOB,
                                   &info);
                        }
                    }
                }
            }
     
            if (info > 0)
                Debug.WriteLine(" __FILE__ : VNL::QR<T>::QtB() -- matrix is rank-def by " + Convert.ToString(info));
  
            return QtB;
        }

        /// <summary>
        /// JOB: ABCDE decimal
        /// A     B     C     D              E
        /// ---   ---   ---   ---            ---
        /// Qb    Q'b   x     norm(A*x - b)  A*x
        /// 
        /// Solve equation M x = b for x using the computed decomposition.
        /// </summary>
        /// <param name="b"></param>
        /// <returns></returns>
        public unsafe Vector Solve(Vector b)
        {
            int n = qrdc_out_.Columns;
            int p = qrdc_out_.Rows;
            float[] b_data = b.Datablock();
            Vector QtB = new Vector(n);
            Vector x = new Vector(p);

            // see comment above
            int JOB = 100;

            int info = 0;

            fixed (float* data = qrdc_out_.Datablock())
            {
                fixed (float* data2 = qraux_.Datablock())
                {
                    fixed (float* data3 = b_data)
                    {
                        fixed (float* data4 = QtB.Datablock())
                        {
                            fixed (float* data5 = x.Datablock())
                            {

                                Netlib.dqrsl_(data, &n, &n, &p, data2, data3,
                                       (float*)0, data4, data5,
                                       (float*)0, // residual*
                                       (float*)0, // Ax*
                                       &JOB,
                                       &info);
                            }
                        }
                    }
                }
            }

            if (info > 0)
                Debug.WriteLine("__FILE__ : VNL::QR<T>::Solve() : matrix is rank-deficient by " + Convert.ToString(info));
  
            return x;
        }