C# ILArray.Clone方法代码示例

        /// <summary>
        /// multiply and fold array elements along first non singleton dimension
        /// </summary>
        /// <param name="inArray">N-dimensional double array</param>
        /// <returns>array having the first non singleton dimension 
		/// reduced to the length 1 with the result of the products of 
		/// corresponding elements of inArray in that dimension.
		/// The result will have the same number of dimensions as 
		/// inArray, but the first non singleton dimension having the 
		/// size 1.</returns>
		public static ILArray<double> prod(ILArray<double> inArray) {
            int[] newDims = inArray.Dimensions.ToIntArray();
            int nsDim = 0;
            while (nsDim < newDims.Length && newDims[nsDim] < 2)
                nsDim++;
            if (nsDim == newDims.Length)
                // scalar -> return copy
                return (ILArray<double>)inArray.Clone();
            return prod(inArray, nsDim);
        }

 /// <summary>
 /// find eigenvalues  and eigenvectors 
 /// </summary>
 /// <param name="A">input: square matrix, size [n x n]</param>
 /// <param name="V">output (optional): eigenvectors</param>   
 /// <param name="propsA">matrix properties, on input - if specified, 
 /// will be used to choose the proper method of solution. On exit will be 
 /// filled according to the properties of A.</param>
 /// <param name="balance">true: permute A in order to increase the 
 /// numerical stability, false: do not permute A.</param>
 /// <returns>eigenvalues as vector (if V is null) or as diagonoal 
 /// matrix (if V was requested, i.e. not equaled null).</returns>
 /// <remarks><para>The eigenvalues of A are found by use of the 
 /// Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
 /// <para>The arrays returned will be of complex inner type, 
 /// since no further constraints are set on the structure of 
 /// A (it may be nonsymmetric). Use 
 /// <see cref="ILNumerics.BuiltInFunctions.ILMath.eigSymm(ILArray&lt;double&gt;)"/> 
 /// or <see cref="ILNumerics.BuiltInFunctions.ILMath.eigSymm(ILArray&lt;double&gt;,ref ILArray&lt;double&gt;)"/> 
 /// functions for computing the real eigenvalues of symmetric 
 /// matrices explicitly.</para>
 /// <para>Depending on the parameter <paramref name="balance"/>, 
 /// A will be balanced first. This includes permutations and 
 /// scaling of A in order to improve the conditioning of the 
 /// eigenvalues.</para></remarks>
 /// <seealso cref="ILNumerics.BuiltInFunctions.ILMath.eig(ILArray&lt;double&gt;)"/>
 /// <seealso cref="ILNumerics.BuiltInFunctions.ILMath.eig(ILArray&lt;double&gt;,ref ILArray&lt;complex&gt;,ref MatrixProperties,bool)"/>
 /// <exception cref="ILNumerics.Exceptions.ILArgumentException">if a 
 /// is not square</exception>
 public static /*!HC:HCClsCmplx*/ ILArray<complex> eig(/*!HC:HCCls1*/ ILArray<double> A, ref /*!HC:HCClsCmplx*/ ILArray<complex> V, ref MatrixProperties propsA, bool balance) {
     if (A.IsEmpty) {
         V = /*!HC:HCClsCmplx*/ ILArray<complex> .empty(A.Dimensions); 
         return /*!HC:HCClsCmplx*/ ILArray<complex> .empty(A.Dimensions); 
     }
     /*!HC:HCClsCmplx*/ ILArray<complex> ret = null;  
     int n = A.Dimensions[0]; 
     bool createVR = (object.Equals(V,null))? false:true; 
     if (n != A.Dimensions[1]) 
         throw new ILArgumentException("eig: matrix A must be square!");
     propsA |= MatrixProperties.Square; 
     if (((propsA & MatrixProperties.Hermitian) != 0 || ILMath.ishermitian(A))) {
         propsA |= MatrixProperties.Hermitian; 
         /*!HC:HCCls1*/ ILArray<double> Vd = null; 
         if (createVR) 
             Vd = /*!HC:HCCls1*/ ILArray<double> .empty(0,0);
         /*!HC:HCCls1*/ ILArray<double> tmpRet = eigSymm(A,ref Vd);
         if (createVR)
             V = /*!HC:HCCls2Cmplx*/ ILMath.tocomplex (Vd); 
         ret = /*!HC:HCCls2Cmplx*/ ILMath.tocomplex (tmpRet);
     } else {
         // nonsymmetric case
         char bal = (balance)? 'B':'N', jobvr;  
         /*!HC:HCCls1*/ ILArray<double> tmpA = (/*!HC:HCCls1*/ ILArray<double> )A.Clone(); 
         /*!HC:HCArr1*/ double [] vr = null;
         /*!HC:HCArr1*/ double [] wr = ILMemoryPool.Pool.New</*!HC:HCArr1*/ double >(n); 
         /*!HC:HCArrWI*/ 
         double[] wi = ILMemoryPool.Pool.New<double>(n); 
         /*!HC:HCArrReal*/ double [] scale  = ILMemoryPool.Pool.New</*!HC:HCArrReal*/ double >(n);
         /*!HC:HCArrReal*/ double [] rconde = ILMemoryPool.Pool.New</*!HC:HCArrReal*/ double >(n); 
         /*!HC:HCArrReal*/ double [] rcondv = ILMemoryPool.Pool.New</*!HC:HCArrReal*/ double >(n); 
         /*!HC:HCArrReal*/ double abnrm = 0; 
         int ldvr, ilo = 0, ihi = 0, info = 0;
         if (createVR) {
             ldvr = n;
             vr = ILMemoryPool.Pool.New</*!HC:HCArr1*/ double >(n * n);
             jobvr = 'V'; 
         } else {
             ldvr = 1; 
             vr = new /*!HC:HCArr1*/ double [1]; 
             jobvr = 'N'; 
         }
         /*!HC:HC?geevx*/
         Lapack.dgeevx(bal,'N',jobvr,'N',n,tmpA.m_data,n,wr,wi,new /*!HC:HCArr1*/ double [1],1,vr,ldvr,ref ilo,ref ihi,scale,ref abnrm,rconde,rcondv,ref info);   
         if (info != 0) 
             throw new ILArgumentException("eig: error in Lapack '?geevx': (" + info + ")");
         // create eigenvalues 
         /*!HC:HCArrCmplx*/ complex [] retArr = ILMemoryPool.Pool.New</*!HC:HCArrCmplx*/ complex >(n); 
         for (int i = 0; i < n; i++) {
             /*!HC:HCSortEVal*/
             retArr[i].real = wr[i]; retArr[i].imag = wi[i]; 
         }
         ret = new /*!HC:HCClsCmplx*/ ILArray<complex> (retArr,n,1);
         if (createVR) {
             #region HCSortEVec
             complex [] VArr = ILMemoryPool.Pool.New< complex >(n * n);
             for (int c = 0; c < n; c++) {
                 if (wi[c] != 0 && wi[c+1] != 0 && c < n-1) {
                     ilo = n * c; ihi = ilo + n;
                     for (int r = 0; r < n; r++) {
                         VArr[ilo].real = vr[ilo];
                         VArr[ilo].imag = vr[ihi];  
                         VArr[ihi].real = vr[ilo];  
                         VArr[ihi].imag = -vr[ihi]; 
                         ilo++; ihi++; 
                     }
                     c++; 
                 } else {
                     ilo = n * c;
                     for (int r = 0; r < n; r++) {
                         VArr[ilo].real = vr[ilo];
//.........这里部分代码省略.........

        /// <summary>
        /// maximum
        /// </summary>
        /// <param name="A">input array, N-dimensional</param>
        /// <param name="I">return value. If this is an instance of an ILArray 
        /// (f.e. 'empty'), on return I will hold the indices into leadDim of  
        /// the maximum values. If, on entering the function, I is null, those indices 
        /// will not be computed and I will be ignored.</param>
        /// <param name="leadDim">index of dimension to operate along</param>
        /// <returns>ILArray of type double. If I was empty  having the dimension 'leadDim' 
        /// reduced to 1 and holding maximum values </returns>
        public static  ILArray<complex> max( ILArray<complex> A, ref  ILArray<double> I, int leadDim) {
			if (A.IsEmpty) {
                if (!object.Equals (I,null))
                    I =  ILArray<double> .empty(0,0); 
                return  ILArray<complex> .empty(A.Dimensions); 
            }
            ILDimension inDim = A.Dimensions; 
			int[] newDims = inDim.ToIntArray();
            if (leadDim == newDims.Length || inDim[leadDim] == 1)
				// scalar or sum over singleton -> return copy
                return ( ILArray<complex> )A.Clone();

			int newLength;
			newLength = inDim.NumberOfElements / newDims[leadDim];
			newDims[leadDim] = 1;
            complex [] retSystemArr;
            retSystemArr = new  complex [newLength];
            int leadDimLen = inDim[leadDim];
			int nrHigherDims = inDim.NumberOfElements / leadDimLen;
            #region HYCALPER GLOBAL_INIT

            complex result;
            complex curval;
            double [] indices = null;
            bool createIndices = false; 
            if (!Object.Equals(I,null)) {
                indices = new  double [retSystemArr.Length];
                createIndices = true; 
            } 
            #endregion HYCALPER GLOBAL_INIT
            double curabsval; double curabsmaxval; 

			// physical -> pointer arithmetic
            if (leadDim == 0) {
#region physical along 1st leading dimension
                unsafe {
                    fixed ( complex * pOutArr = retSystemArr)
                    fixed ( complex * pInArr = A.m_data)
                    fixed ( double * pIndices = indices) {
                        complex * lastElement; 
                        complex * tmpOut = pOutArr;
                        complex * tmpIn = pInArr;
                        if (createIndices) {
                            double * tmpInd = pIndices;
                            for (int h = nrHigherDims; h-- > 0; ) {
                                lastElement = tmpIn + leadDimLen;
                                curabsmaxval = double.MinValue; result = new complex(); 
                                while (tmpIn < lastElement) {
                                    curval = *tmpIn;
                                    curabsval = complex.Abs(curval);
                                            if (curabsval > curabsmaxval) {
                                                curabsmaxval = curabsval;
                                                result = curval;
                                            
                                        *tmpInd = ( double )(tmpIn - (lastElement - leadDimLen));
                                    }
                                    tmpIn++;
                                }
                                *(tmpOut++) = ( complex )result;
                                tmpInd++; 
                            }
                        } else {
                            double * tmpInd = pIndices;
                            for (int h = nrHigherDims; h-- > 0; ) {
                                lastElement = tmpIn + leadDimLen;
                                curabsmaxval = double.MinValue; result = new complex(); 
                                while (tmpIn < lastElement) {
                                    curval = *tmpIn++;
                                    curabsval = complex.Abs(curval);
                                            if (curabsval > curabsmaxval) {
                                                curabsmaxval = curabsval;
                                                result = curval;
                                            
                                    }
                                }
                                #region HYCALPER POSTLOOP
                                *(tmpOut++) = ( complex )result;
                                #endregion HYCALPER POSTLOOP
                            }
                        }
                    }
                }
#endregion physical along 1st leading dimension
            } else {
#region physical along abitrary dimension
				// sum along abitrary dimension 
                unsafe {
                    fixed ( complex * pOutArr = retSystemArr)
                    fixed ( complex * pInArr = A.m_data)
//.........这里部分代码省略.........

 /// <summary>
 /// Solve system of linear equations A*x = b, with A beeing a lower triangular matrix
 /// </summary>
 /// <param name="A">input matrix of Size [n x n], must be lower triangular. No check is made for that!</param>
 /// <param name="B">solution vector. Size [n x 1]</param>
 /// <param name="singularityDetect">output: this value gives the row of A, where a singularity has been detected (if any). If A is not singular, this will be a negative value.</param>
 /// <returns>Solution x solving A * x = b.</returns>
 /// <remarks><para>The solution will be determined via forward substitution</para>
 /// <para>Make sure, A and b are of correct size, since no checks are made for that!</para>
 /// <para>This function is used by ILMath.linsolve. There should be rare need for you to call this function directly.</para>
 /// <para>Elements of A above the main diagonal will not be accessed.</para>
 /// <para>If A has been found to be singular, the array returned will contain NaN values for corresponding elements!</para></remarks>
 internal static  ILArray<fcomplex> solveLowerTriangularSystem ( ILArray<fcomplex> A,  ILArray<fcomplex> B, ref int singularityDetect) {
     System.Diagnostics.Debug.Assert(B.Dimensions[1] >= 0); 
     System.Diagnostics.Debug.Assert(B.Dimensions[0] == A.Dimensions[1]); 
     System.Diagnostics.Debug.Assert(A.Dimensions[0] == A.Dimensions[1]); 
     singularityDetect = -1; 
     int n = A.Dimensions[0]; 
     int m = B.Dimensions[1]; 
     int spacingA0, spacingA1, info = 0; 
     char transA; 
     ILArray<fcomplex> ret = ( ILArray<fcomplex> )B.Clone(); 
     fcomplex [] retArr = ret.m_data;  
     if (isSuitableForLapack(A,out spacingA0,out spacingA1,out transA)) {
         // solve using Lapack
         unsafe {
             fixed ( fcomplex * ptrA = A.m_data) 
             fixed ( fcomplex * ptrB = ret.m_data) {
                 fcomplex * pA = ptrA + A.getBaseIndex(0); 
                 Lapack.ctrtrs ((transA == 'T')? 'U':'L',transA,'N',A.Dimensions[0],ret.Dimensions[1],(IntPtr)pA,spacingA1,(IntPtr)ptrB,n,ref info);
                 if (info < 0) 
                     throw new ILArgumentException ("linsolve: error inside Lapack function ?trtrs for argument: " + (-info));
                 if (info > 0) {
                     singularityDetect = info-1;
                     retArr = ret.m_data; 
                     for (spacingA0 = 0; spacingA0 < ret.Dimensions[1]; spacingA0++) {
                         info = spacingA0 * n + singularityDetect; 
                         for (spacingA1 = singularityDetect; spacingA1 < n; spacingA1++) {
                             retArr[info++] =  new fcomplex(float.NaN,float.NaN) ;     
                         }
                     }
                 } else {
                     singularityDetect = -1; 
                 }
             }
         }
         return ret; 
     }
     // must do it manually.... 
     System.Diagnostics.Debug.Assert(A.IsReference); 
     int[] Adim0 = A.m_indexOffset[0]; 
     int[] Adim1 = A.m_indexOffset[1];    
     fcomplex diagVal =  new fcomplex(0.0f,0.0f) ;
     fcomplex tmpVal;
     for (int b = 0; b < m; b++) {
         for (int r = 0; r < n; r++) {
             tmpVal =  new fcomplex(0.0f,0.0f) ; 
             diagVal = A.GetValue(r,r); 
             if (diagVal ==  new fcomplex(0.0f,0.0f) ) {
                 singularityDetect = r; 
                 for (; r < n; r++)
                     retArr[r + b * n] =  new fcomplex(float.PositiveInfinity,float.PositiveInfinity) ; 
             } else {
                 for (int c = 0; c < r; c++) {
                     tmpVal +=  A.m_data[Adim0[r] + Adim1[c]] * retArr[c + b * n]; 
                 }
                 retArr[r + b * n] = (retArr[r + b * n] - tmpVal) / diagVal; 
             }
         }
     }
     ret = new  ILArray<fcomplex> (retArr,n,m); 
     return ret; 
 }

 /// <summary>
 /// QR decomposition - raw Lapack output
 /// </summary>
 /// <param name="A">general input matrix A</param>
 /// <returns>orthonormal / unitary matrix Q and upper triangular 
 /// matrix R packed into single matrix. This is the output of the 
 /// lapack function ?geqrf.</returns>
 /// <remarks><para>Input matrix A will not be altered. </para>
 /// <para>The matrix returned is the direct output of the lapack 
 /// function [d,s,c,z]geqrf respectively. This mean that it contains 
 /// the decomposition factors Q and R, but they are cmbined into a 
 /// single matrix for performance reasons. If you need one of the factors, 
 /// you would use the overloaded function 
 /// <see cref="ILNumerics.BuiltInFunctions.ILMath.qr(ILArray&lt;double&gt;,ref ILArray&lt;double&gt;)"/> 
 /// instead, which returns those factors seperately.</para></remarks>
 public static  ILArray<complex> qr( ILArray<complex> A) {
     if (!A.IsMatrix) 
         throw new ILArgumentException("qr decomposition: A must be a matrix"); 
     int m = A.Dimensions[0], n = A.Dimensions[1]; 
     ILArray<complex> ret = ( ILArray<complex> )A.Clone();  
     complex [] tau = new  complex [(m<n)?m:n];  
     int info = 0; 
     Lapack.zgeqrf (m,n,ret.m_data,m,tau,ref info); 
     if (info < 0)
         throw new ILArgumentException("qr: an error occoured during decomposition"); 
     return ret; 
 }

// DO NOT EDIT INSIDE THIS REGION !! CHANGES WILL BE LOST !! 
        /// <summary>
        /// maximum
        /// </summary>
        /// <param name="A">input array, N-dimensional</param>
        /// <param name="I">return value. If this is an instance of an ILArray 
        /// (f.e. 'empty'), on return I will hold the indices into leadDim of  
        /// the maximum values. If, on entering the function, I is null, those indices 
        /// will not be computed and I will be ignored.</param>
        /// <param name="leadDim">index of dimension to operate along</param>
        /// <returns>ILArray of type double. If I was empty  having the dimension 'leadDim' 
        /// reduced to 1 and holding maximum values </returns>
        public static  ILLogicalArray min( ILArray<byte> A, ref  ILArray<double> I, int leadDim) {
			ILDimension inDim = A.Dimensions; 
			int[] newDims = inDim.ToIntArray();
            if (leadDim == newDims.Length || inDim[leadDim] == 1)
				// scalar or sum over singleton -> return copy
                return ( ILLogicalArray )A.Clone();

			int newLength;
			newLength = inDim.NumberOfElements / newDims[leadDim];
			newDims[leadDim] = 1;
            byte [] retSystemArr;
            retSystemArr = new  byte [newLength];
            int leadDimLen = inDim[leadDim];
			int nrHigherDims = inDim.NumberOfElements / leadDimLen;
            #region HYCALPER GLOBAL_INIT

            byte result;
            byte curval;
            double [] indices = null;
            bool createIndices = false; 
            if (!Object.Equals(I,null)) {
                indices = new  double [retSystemArr.Length];
                createIndices = true; 
            } 
            #endregion HYCALPER GLOBAL_INIT
            

			// physical -> pointer arithmetic
            if (leadDim == 0) {
#region physical along 1st leading dimension
                unsafe {
                    fixed ( byte * pOutArr = retSystemArr)
                    fixed ( byte * pInArr = A.m_data)
                    fixed ( double * pIndices = indices) {
                        byte * lastElement; 
                        byte * tmpOut = pOutArr;
                        byte * tmpIn = pInArr;
                        if (createIndices) {
                            double * tmpInd = pIndices;
                            for (int h = nrHigherDims; h-- > 0; ) {
                                lastElement = tmpIn + leadDimLen;
                                result = byte.MaxValue;
                                while (tmpIn < lastElement) {
                                    curval = *tmpIn;
                                    if (curval < result) {
                                                result = curval;
                                            
                                        *tmpInd = ( double )(tmpIn - (lastElement - leadDimLen));
                                    }
                                    tmpIn++; 
                                }
                                *(tmpOut++) = ( byte )result;
                                tmpInd++; 
                            }
                        } else {
                            double * tmpInd = pIndices;
                            for (int h = nrHigherDims; h-- > 0; ) {
                                lastElement = tmpIn + leadDimLen;
                                result = byte.MaxValue;
                                while (tmpIn < lastElement) {
                                    curval = *tmpIn++;
                                    if (curval < result) {
                                                result = curval;
                                            
                                    }
                                }
                                #region HYCALPER POSTLOOP
                                *(tmpOut++) = ( byte )result;
                                #endregion HYCALPER POSTLOOP
                            }
                        }
                    }
                }
#endregion physical along 1st leading dimension
            } else {
#region physical along abitrary dimension
				// sum along abitrary dimension 
                unsafe {
                    fixed ( byte * pOutArr = retSystemArr)
                    fixed ( byte * pInArr = A.m_data)
                    fixed ( double * pIndices = indices) {
                        byte * lastElementOut = newLength + pOutArr - 1;
                        int inLength = inDim.NumberOfElements -1;
                        byte * lastElementIn = pInArr + inLength; 
                        int inc = inDim.SequentialIndexDistance(leadDim);
                        byte * tmpOut = pOutArr;
                        int outLength = newLength - 1;
                        byte * leadEnd;
//.........这里部分代码省略.........

 /// <summary>
 /// Solve linear equation system
 /// </summary>
 /// <param name="A">Matrix A. Size [n x q]</param>
 /// <param name="B">'rigth hand side' B. Size [n x m]</param>
 /// <param name="props">Matrix properties. If defined, no checks are made for the structure of A. If the matrix A was found to be (close to or) singular, the 'MatrixProperties.Singular' flag in props will be set. This flag should be tested on return, in order to verify the reliability of the solution.</param>
 /// <returns>the solution x solving multiply(A,x) = B. Size [n x m]</returns>
 /// <remarks><para>depending on the <paramref name="props"/> parameter the equation system will be solved differently for special structures of A:
 /// <list type="bullet">
 /// <item>If A is square (q == n) and an <b>upper or lower triangular</b> matrix, the system will directly be solved via backward- or forward substitution. Therefore the Lapack function ?trtrs will be used, whenever the memory layout of A is suitable. This may be the case even for reference ILArray's! 
 /// <example><![CDATA[ILArray<double> A = ILMath.randn(4); // construct 4 x 4 matrix 
 /// A = A.T; // A is a reference array now! The transpose is fast and does not consume much memory
 /// // now construct a right side and solve the equations: 
 /// ILArray<double> B = new ILArray<double> (1.0,2.0,3.0).T; 
 /// ILMath.linsolve(A,B); // ... will be carried out via Lapack, even for all arrays involved beeing reference arrays! ]]></example></item>
 /// <item><para>if A is square and symmetric or hermitian, A will be decomposed into a triangular equation system using cholesky factorization and Lapack. The system is than solved using performant Lapack routines.</para>
 /// <para>if during the cholesky factorization A was found to be <b>not positive definite</b> - the corresponding flag in props will be cleaned and <c>null</c> will be returned.</para></item>
 /// <item>otherwise if A is square only, it will be decomposed into upper and lower triangular matrices using LU decomposition and Lapack. The triangular system is than solved using performant Lapack routines.</item>
 /// <item>otherwise, if A is of size [q x n] and q != n, the system is solved using QR decomposition. A may be rank deficient. The solution is computed by use of the Lapack routine '?gelsy' and may be a reference array.</item>
 /// </list></para>
 /// <para>Compatibility with Matlab<sup>(R)</sup>: If A is square, the algorithm used follows the same logic as Matlab up to Rel 14, with the exception of Hessenberg matrices wich are solved via LU factorization here. The un-squared case is handled differently. A direct Lapack driver function (?gelsy) is used here. Therefore the solutions might differ! However, the solution will of course fullfill the equation A * x = B without round off errrors. </para>
 /// <para>For specifiying the rank of A in the unsquare case (q != n), the eps member from <see cref="ILNumerics.Settings.ILSettings"/> class is used.</para></remarks>
 public static  ILArray<complex> linsolve( ILArray<complex> A,  ILArray<complex> B, ref MatrixProperties props) {
     if (object.Equals(A,null)) 
         throw new ILArgumentException("linsolve: input argument A must not be null!"); 
     if (object.Equals(B,null))
         throw new ILArgumentException("linsolve: input argument B must not be null!"); 
     if (A.IsEmpty || B.IsEmpty)
         return  ILArray<complex> .empty(A.Dimensions); 
     if (A.Dimensions[0] != B.Dimensions[0])
         throw new ILArgumentException("linsolve: number of rows for matrix A must match number of rows for RHS!");
     int info = 0, m = A.Dimensions[0]; 
     ILArray<complex> ret; 
     if (m == A.Dimensions[1]) {
         props |= MatrixProperties.Square; 
         if ((props & MatrixProperties.LowerTriangular) != 0) {
             ret = solveLowerTriangularSystem(A,B,ref info); 
             if (info > 0)
                 props |= MatrixProperties.Singular; 
             return ret;
         } 
         if ((props & MatrixProperties.UpperTriangular) != 0) {
             ret = solveUpperTriangularSystem(A,B, ref info); 
             if (info > 0)
                 props |= MatrixProperties.Singular; 
            return ret; 
         } 
         if ((props & MatrixProperties.Hermitian) != 0) {
             ILArray<complex> cholFact = A.copyUpperTriangle(m);
             Lapack.zpotrf ('U', m, cholFact.m_data, m, ref info); 
             if (info > 0) {
                 props ^= MatrixProperties.Hermitian; 
                 return null; 
             }  else {
                 // solve 
                 ret = ( ILArray<complex> )B.Clone(); 
                 Lapack.zpotrs ('U',m,B.Dimensions[1],cholFact.m_data,m,ret.m_data,m,ref info);   
                 return ret; 
             } 
         } else {
             // attempt complete (expensive) LU factorization 
             ILArray<complex> L = ( ILArray<complex> )A.Clone();  
             int [] pivInd = new int[m]; 
             Lapack.zgetrf (m, m, L.m_data, m, pivInd ,ref info); 
             if (info > 0)
                 props |= MatrixProperties.Singular; 
             ret = ( ILArray<complex> )B.Clone(); 
             Lapack.zgetrs ('N',m,B.Dimensions[1],L.m_data,m,pivInd,ret.m_data,m,ref info); 
             if (info < 0) 
                 throw new ILArgumentException("linsolve: failed to solve via lapack dgetrs");
             return ret; 
         }
     } else {
         // under- / overdetermined system
         int n = A.Dimensions[1], rank = 0, minMN = (m < n)? m:n, maxMN = (m > n)? m:n;
         int nrhs = B.Dimensions[1];
         if (B.Dimensions[0] != m) 
             throw new ILArgumentException("linsolve: right hand side matrix B must match input A!"); 
         ILArray<complex> tmpA = ( ILArray<complex> )A.Clone();
         if (m < n) {
             ret = new  ILArray<complex> (new  complex [n * nrhs],n,nrhs);
             ret["0:"+(m-1)+";:"] = B;
         } else {
             ret = ( ILArray<complex> )B.Clone();
         }
         int [] JPVT = new int [n]; 
         Lapack.zgelsy (m, n, B.Dimensions[1], tmpA.m_data, m, ret.m_data, 
                             maxMN, JPVT,  ILMath.MachineParameterDouble.eps ,
                             ref rank, ref info); 
         if (n < m) {
             ret = ret[ILMath.vector(0,n-1),null];
         }
         if (rank < minMN)
             props |= MatrixProperties.RankDeficient; 
         return ret; 
     }
 }

        /// <summary>
		/// Sum elements of A along dimension specified.
		/// </summary>
		/// <param name="A">N-dimensional array</param>
		/// <param name="leadDim">index of dimension to operate along</param>
		/// <returns>array, same size as A, but having the 'leadDim's dimension 
		/// reduced to the length 1 with the sum of all
		/// elements along that dimension.</returns>
        public static /*!HC:outCls1*/ ILArray<double> /*!HC:funcname*/ sum (/*!HC:inCls1*/ ILArray<double> A, int leadDim) {
            if (A.IsEmpty) 
                return /*!HC:outCls1*/ ILArray<double> .empty(A.Dimensions); 
            if (A.IsScalar)
                /*!HC:HCscalarOp*/
                return new /*!HC:outCls1*/ ILArray<double> (new /*!HC:inArr1*/ double []{A.GetValue(0)},1,1);
            if (leadDim >= A.Dimensions.NumberOfDimensions)
                throw new ILArgumentException("dimension parameter out of range!");
            ILDimension inDim = A.Dimensions; 
			int[] newDims = inDim.ToIntArray();
            /*!HC:singletonDimOp*/
            if (inDim[leadDim] == 1) return (/*!HC:outCls1*/ ILArray<double> )A.Clone();
			int newLength;
			/*!HC:outArr1*/ double [] retDblArr;
			// build ILDimension
			newLength = inDim.NumberOfElements / newDims[leadDim];
			newDims[leadDim] = 1;
			retDblArr = ILMemoryPool.Pool.New</*!HC:outArr1*/ double >(newLength);
            ILDimension newDimension = new ILDimension(newDims); 
            int incOut = newDimension.SequentialIndexDistance(leadDim); 
            int leadDimLen = inDim[leadDim]; 
			int nrHigherDims = inDim.NumberOfElements / leadDimLen;

			// physical -> pointer arithmetic
            if (leadDim == 0) {
#region physical along 1st leading dimension
                unsafe {
                    fixed (/*!HC:outArr1*/ double * pOutArr = retDblArr)
                    fixed (/*!HC:inArr1*/ double * pInArr = A.m_data) {
                        /*!HC:inArr1*/ double * lastElement;
                        /*!HC:outArr1*/ double * tmpOut = pOutArr;
                        /*!HC:inArr1*/ double * tmpIn = pInArr;
                        for (int h = nrHigherDims; h-- > 0; ) {
							lastElement = tmpIn + leadDimLen;
                            /*!HC:HCzero*/ 
                            *tmpOut = 0.0; 
                            while (tmpIn < lastElement) {
                                /*!HC:tmpOutStorage*/ *tmpOut += /*!HC:preEvalOp*/ (double) (*tmpIn++) /*!HC:postEvalOp*/ ;
                            }
						    /*!HC:operationResult*/ 
                            /**/
                            tmpOut++;
                        }
                    }
				}
#endregion
			} else {
#region physical along abitrary dimension
				// sum along abitrary dimension 
                unsafe {
                    fixed (/*!HC:outArr1*/ double * pOutArr = retDblArr)
                    fixed (/*!HC:inArr1*/ double * pInArr = A.m_data) {
                        /*!HC:outArr1*/ double * lastElementOut = newLength + pOutArr -1;
                        int inLength = inDim.NumberOfElements -1; 
                        /*!HC:inArr1*/ double * lastElementIn = pInArr + inLength; 
                        int inc = inDim.SequentialIndexDistance(leadDim); 
                        /*!HC:outArr1*/ double * tmpOut = pOutArr;
                        int outLength = newLength - 1;
                        /*!HC:inArr1*/ double * leadEnd; 
                        /*!HC:inArr1*/ double * tmpIn = pInArr;
                        for (int h = nrHigherDims; h--> 0; ) {
							leadEnd = tmpIn + leadDimLen * inc;
                            /*!HC:HCzero*/ 
                            *tmpOut = 0.0; 
                            while (tmpIn < leadEnd) {
                                /*!HC:tmpOutStorage*/ *tmpOut += /*!HC:preEvalOp*/ (double) (*tmpIn) /*!HC:postEvalOp*/ ;
                                tmpIn += inc; 
                            }
							/*!HC:operationResult*/ 
                            /**/
                            tmpOut  += inc;
                            if (tmpOut > lastElementOut)
                                tmpOut = pOutArr + ((tmpOut - pOutArr) - outLength);
                            if (tmpIn > lastElementIn)
                                tmpIn = pInArr + ((tmpIn - pInArr) - inLength); 
                        }
                    }
				}
#endregion
			}

			return new /*!HC:outCls1*/ ILArray<double> (retDblArr, newDims);;
        }

		/// <summary>
		/// Multiply elements of inArray along specified dimension.
		/// </summary>
		/// <param name="inArray">N-dimensional double array</param>
		/// <param name="leadDim">index of dimension to multiply elements along</param>
		/// <returns>array having the 'leadDim's dimension 
		/// reduced to the length of 1 with the result of the product of 
		/// corresponding elements of inArray of that dimension.</returns>
        public static ILArray<double> prod(ILArray<double> inArray, int leadDim) {
			ILDimension inDim = inArray.Dimensions; 
			int[] newDims = inDim.ToIntArray();
            if (leadDim == newDims.Length || inDim[leadDim] == 1)
				// scalar or sum over singleton -> return copy
                return (ILArray<double>)inArray.Clone();

            int newLength;
			double[] retDblArr;
			// build ILDimension
			newLength = inDim.NumberOfElements / newDims[leadDim];
			newDims[leadDim] = 1;
			retDblArr = ILMemoryPool.Pool.New<double>(newLength);
            int leadDimLen = inDim[leadDim];
			int nrHigherDims = inDim.NumberOfElements / leadDimLen;

			// physical -> pointer arithmetic
            if (leadDim == 0) {
#region physical along 1st leading dimension
                unsafe {
                    fixed (double* pOutArr = retDblArr,
                            pInArr = inArray.m_data) {
                        double* lastElement;
                        double* tmpOut = pOutArr;
                        double* tmpIn = pInArr;
                        for (int h = nrHigherDims; h-- > 0; ) {
							lastElement = tmpIn + leadDimLen;
                            *tmpOut = 1.0;
                            while (tmpIn < lastElement) {
                                *tmpOut *= *tmpIn++;
                            }
                            tmpOut++;
                        }
                    }
				}
#endregion
			} else {
#region physical along abitrary dimension
				// sum along abitrary dimension 
                unsafe {
                    fixed (double* pOutArr = retDblArr,
                            pInArr = inArray.m_data) {
                        double* lastElementOut = newLength + pOutArr -1;
                        int inLength = inDim.NumberOfElements -1; 
                        double* lastElementIn = pInArr + inLength; 
                        int inc = inDim.SequentialIndexDistance(leadDim); 
                        double* tmpOut = pOutArr;
                        int outLength = newLength - 1;
                        double* leadEnd; 
                        double* tmpIn = pInArr;
                        for (int h = nrHigherDims; h--> 0; ) {
							leadEnd = tmpIn + leadDimLen * inc;
                            *tmpOut = 1.0;
                            while (tmpIn < leadEnd) {
                                *tmpOut *= *tmpIn;
                                tmpIn += inc; 
                            }
                            tmpOut  += inc;
                            if (tmpOut > lastElementOut)
                                tmpOut = pOutArr + ((tmpOut - pOutArr) - outLength);
                            if (tmpIn > lastElementIn)
                                tmpIn = pInArr + ((tmpIn - pInArr) - inLength); 
                        }
                    }
				}
#endregion
			}

			return new ILArray<double>(retDblArr, newDims);;  
		}

 public void Test_SetRange_RangedPhysical() {
     int errorCode = 0; 
     try {
         ILArray<double> AOrig = ILMath.vector(1,24).Reshape(new ILDimension(4,3,2));  
         ILArray<double> ResOrig = new ILArray<double>(new double[]{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24},4,3,2);
         ILArray<double> A, Res, indb; 
         // test setrange for single element (sequential access) 
         for (int i = 0; i < 20; i++) {
             A = (ILArray<double>)AOrig.Clone();  
             Res = ((ILArray<double>)ResOrig.Clone()).Detach(); 
             A[i.ToString()] = 99;
             Res.m_data[i] = 99;
             if (!A.Equals(Res)) throw new Exception("A[\"" + i.ToString() + "\"] = 99 failed.");
             // test setrange via basearray
             A = (ILArray<double>)AOrig.Clone();  
             A[new int[]{i}] = 99;
             if (!A.Equals(Res)) throw new Exception("A[" + i.ToString() + "] = 99 failed.");
         }
         // test setrange for 2 single elements (sequential access) 
         errorCode = 1; 
         for (int i = 0; i < 20; i++) {
             A = (ILArray<double>)AOrig.Clone();  
             Res = ((ILArray<double>)ResOrig.Clone()).Detach(); 
             indb = new double[]{i,i+3}; 
             string inds = i.ToString() + "," + (i+3).ToString(); 
             Res.m_data[i] = 99;
             Res.m_data[i+3] = 99; 
             A[inds] = 99;
             if (!A.Equals(Res)) throw new Exception("A[\"" + i.ToString() + "," + (i+3).ToString() + "\"] = 99 failed.");
             // test setrange via basearray
             A = (ILArray<double>)AOrig.Clone();  
             A[indb] = 99;
             if (!A.Equals(Res)) throw new Exception("A[" + i.ToString() + "," + (i+3).ToString() + "] = 99 failed.");
         }
         errorCode = 2; 
         A = (ILArray<double>)AOrig.Clone();
         // reshaping range in second dimension 
         A["1;1:end"] = new double []{100,101,102,103,104};
         Res = ((ILArray<double>)ResOrig.Clone()).Detach(); 
         Res.m_data[5] = 100; 
         Res.m_data[9] = 101; 
         Res.m_data[13] = 102; 
         Res.m_data[17] = 103; 
         Res.m_data[21] = 104;
         if (!A.Equals(Res)) throw new Exception("A[1,1:end] = [100...104] failed.");
         // reshaping range in 2nd dimension - basearray indices
         indb = new double[]{1,2,3,4,5}; 
         A[1.0,indb] = new double[]{100,101,102,103,104};
         Res = ((ILArray<double>)ResOrig.Clone()).Detach(); 
         Res.m_data[5] = 100; 
         Res.m_data[9] = 101; 
         Res.m_data[13] = 102; 
         Res.m_data[17] = 103; 
         Res.m_data[21] = 104;
         if (!A.Equals(Res)) throw new Exception("A[1,[1,2,3,4,5]] = [100...104] failed.");
         // reshaping range in 2nd dimension - for reference storages 
         errorCode = 2; 
         A = ((ILArray<double>)AOrig.Clone()); 
         ILArray<double>.DetachReferences = ILDetachingBehavior.DetachSave;
         ILArray<double> ARef = (ILArray<double>)A.CreateReference();
         A.Dispose(); 
         // reshaping range in second dimension 
         ARef["1;1:end"] = new double[]{100,101,102,103,104};
         Res = ((ILArray<double>)ResOrig.Clone()).Detach(); 
         Res.m_data[5] = 100; 
         Res.m_data[9] = 101; 
         Res.m_data[13] = 102; 
         Res.m_data[17] = 103; 
         Res.m_data[21] = 104;
         if (!ARef.Equals(Res) || !ARef.IsReference) throw new Exception("A[1,1:end] = [100...104] failed.");
         // reshaping range in 2nd dimension - basearray indices
         indb = new double[]{1,2,3,4,5}; 
         ARef[1.0,indb] = new double[]{100,101,102,103,104};
         Res = ((ILArray<double>)ResOrig.Clone()).Detach(); 
         Res.m_data[5] = 100; 
         Res.m_data[9] = 101; 
         Res.m_data[13] = 102; 
         Res.m_data[17] = 103; 
         Res.m_data[21] = 104;
         if (!ARef.Equals(Res)) throw new Exception("A[1,[1,2,3,4,5]] = [100...104] failed.");
         // 
         errorCode = 3; 
         // test single element (sequential) via range for reference array
         for (int i = 0; i < 20; i++) {
             A = (ILArray<double>)AOrig.Clone(); 
             ARef = (ILArray<double>)A.CreateReference(); 
             A.Dispose(); 
             Res = ((ILArray<double>)ResOrig.Clone()).Detach(); 
             ARef[i.ToString()] = 99;
             Res.m_data[i] = 99;
             if (!ARef.Equals(Res) || !ARef.IsReference) throw new Exception("A[\"" + i.ToString() + "\"] = 99 failed.");
             // test setrange via basearray
             A = (ILArray<double>)AOrig.Clone(); 
             ARef = (ILArray<double>)A.CreateReference(); 
             A.Dispose(); 
             ARef[i] = 99;
             if (!ARef.Equals(Res) || !ARef.IsReference) throw new Exception("A[" + i.ToString() + "] = 99 failed.");
         }
         errorCode = 4; 
         
//.........这里部分代码省略.........

        public void Test_IndexAccessPhysicalSequential() {
            int errorCode = 1;
            try {
                // B = A[idx,0]
				ILArray<double> A = ILMath.reshape(ILMath.vector(0.0,23.0),2,3,4);
                A.MinimumRefDimensions = 2;
				ILArray<double> ind = new ILArray<double> (2,4,6,8); 
                ILArray<double> B = A[ind];  
				if (!B.Equals(ind))
					throw new Exception("Invalid value of result!"); 
				errorCode = 2;
                // B = A[idx,0]
				A = ILMath.reshape(ILMath.vector(0.0,23.0),2,3,4);
				ind = new ILArray<double> (2,4,6,8); 
				A.MinimumRefDimensions = 4; 
                B = A[ind];  
				if (!B.Equals(ind))
					throw new Exception("Invalid value of result!"); 
				errorCode = 3;
                ind = new ILArray<double> (new double[4]{2,4,6,8},2,2); 
                ILArray<double> Res = (ILArray<double>)ind.Clone(); 
                B = A[ind]; 
				if (!B.Equals(Res))
					throw new Exception("Invalid value of result!"); 
				errorCode = 4;
                ind = new ILArray<double> (new double[4]{2,4,6,8},2,2); 
                Res = (ILArray<double>)ind.T; 
                B = A[1,ind]; 
				if (!B.Equals(Res))
					throw new Exception("Invalid value of result!"); 
				Success();
                
            } catch (Exception e) {
				Error(errorCode,e.Message);
            }
        }

        /// <summary>
		/// Sum elements of A along dimension specified.
		/// </summary>
		/// <param name="A">N-dimensional array</param>
		/// <param name="leadDim">index of dimension to operate along</param>
		/// <returns>array, same size as A, but having the 'leadDim's dimension 
		/// reduced to the length 1 with the sum of all
		/// elements along that dimension.</returns>
        public static  ILArray<UInt16>  sum ( ILArray<UInt16> A, int leadDim) {
            if (leadDim >= A.Dimensions.NumberOfDimensions)
                throw new ILArgumentException("dimension parameter out of range!");
            if (A.IsEmpty)       
                return  ILArray<UInt16> .empty(A.Dimensions); 
            if (A.IsScalar) {
               
                return new  ILArray<UInt16> (new  UInt16 []{A.GetValue(0)},1,1);
            }
            ILDimension inDim = A.Dimensions; 
			int[] newDims = inDim.ToIntArray();
           
            if (inDim[leadDim] == 1) return ( ILArray<UInt16> )A.Clone();
			int newLength;
			 UInt16 [] retDblArr;
			// build ILDimension
			newLength = inDim.NumberOfElements / newDims[leadDim];
			newDims[leadDim] = 1;
			retDblArr = ILMemoryPool.Pool.New< UInt16 >(newLength);
            ILDimension newDimension = new ILDimension(newDims); 
            int incOut = newDimension.SequentialIndexDistance(leadDim); 
            int leadDimLen = inDim[leadDim];
            int posCounter; 
			int nrHigherDims = inDim.NumberOfElements / leadDimLen;
            if (A.IsReference) {
#region Reference storage
				// ========================  REFERENCE double Storage ===========
				if (A.IsMatrix) {
#region Matrix
					////////////////////////////   MATRIX   ///////////////////////
					unsafe {
						ILIndexOffset idxOffset = A.m_indexOffset;
						int secDim = (leadDim + 1) % 2;
						fixed (int* leadDimStart = idxOffset[leadDim],
								secDimStart = idxOffset[secDim]) {
							fixed ( UInt16 * pOutArr = retDblArr)
						    fixed ( UInt16 * pInArr = A.m_data) {
								 UInt16 * tmpOut = pOutArr;
								 UInt16 * lastElementOut = tmpOut + retDblArr.Length;
								 UInt16 * tmpIn = pInArr;
								int* secDimEnd = secDimStart + idxOffset[secDim].Length - 1;
								int* secDimIdx = secDimStart; 
								int* leadDimIdx = leadDimStart;
								int* leadDimEnd = leadDimStart + leadDimLen - 1;
								// start at first element
								while (secDimIdx <= secDimEnd) {
									tmpIn = pInArr + *secDimIdx++; 
									leadDimIdx = leadDimStart;
                                    *tmpOut = 0;
									while (leadDimIdx <= leadDimEnd) {
                                        UInt16 inVal = *(tmpIn + *leadDimIdx++); 
                                        /**/
										 *tmpOut +=  (UInt16) (inVal)  ;
									}
                                    
                                    /**/
                                    tmpOut++;
								}
							}
						}
					} 
#endregion
				} else {
					/////////////////////////////   ARBITRARY DIMENSIONS //////////
#region arbitrary size
					unsafe {
						ILIndexOffset idxOffset = A.m_indexOffset;
						int[] curPosition = new int[A.Dimensions.NumberOfDimensions];
						fixed (int* leadDimStart = idxOffset[leadDim]) {
							fixed ( UInt16 * pOutArr = retDblArr)
							fixed ( UInt16 * pInArr = A.m_data) {
								 UInt16 * tmpOut = pOutArr;
								 UInt16 * lastElementOut = tmpOut + retDblArr.Length - 1;
								 UInt16 * tmpIn = pInArr + A.m_indexOffset.Map(0);
								int* leadDimIdx = leadDimStart;
								int* leadDimEnd = leadDimStart + leadDimLen;
								int dimLen = curPosition.Length;
								int d, curD;
								// start at first element
                                posCounter = retDblArr.Length; 
								while (posCounter-->0) {
									leadDimIdx = leadDimStart;
                                    *tmpOut = 0;
									while (leadDimIdx < leadDimEnd){
                                        UInt16 inVal = *(tmpIn + *leadDimIdx++); 
                                        /**/
										 *tmpOut +=  (UInt16) (inVal)  ;
									 
                                    /**/
                                    }
                                    tmpOut += incOut;
                                    if (tmpOut > lastElementOut)
//.........这里部分代码省略.........

 protected void CreateMeshILArray(ILArray<double> x, ILArray<double> y, ILArray<double> z)
 {
     bounds = new Cuboid(x.MinValue, y.MinValue, z.MinValue, x.MaxValue, y.MaxValue, z.MaxValue);
     lengthU = x.Dimensions[0];
     lengthV = x.Dimensions[1];
     ILArray<double> xs, ys, zs;
     if (x.IsReference)
         xs = x.Clone() as ILArray<double>;
     else xs = x;
     if (y.IsReference)
         ys = y.Clone() as ILArray<double>;
     else ys = y;
     if (z.IsReference)
         zs = z.Clone() as ILArray<double>;
     else zs = z;
     //if (x.IsReference || y.IsReference || z.IsReference) throw new Exception("x, y and z must be solid arrays");
     double[] xa = xs.InternalArray4Experts;
     double[] ya = ys.InternalArray4Experts;
     double[] za = zs.InternalArray4Experts;
     Cuboid modelBounds = new Cuboid(new System.Windows.Media.Media3D.Point3D(-10, -10, -10), new System.Windows.Media.Media3D.Point3D(10, 10, 10));
     UpdateModelVertices(xa, ya, za, lengthU, lengthV);
     CreateVertsAndInds();
     colourMap = new ColourMap(ColourMapType.Jet, 256);
     colourMapIndices = FalseColourImage.IEnumerableToIndexArray(za, lengthU, lengthV, 256);
     SetColorFromIndices();
 }