public static void Main()
 {
     var a = new complex(1,2);
     var b = new complex(2,3);
     var c = a+b;
     if(c == new complex(3,5)) c.print("test complex addition passed");
     else c.print("test complex addition failed");
 }

        public static void fftc1d(ref complex[] a) {
            int n;


            n = ap.len(a);
            fft.fftc1d(ref a, n);

            return;
        }

 public void AddTest()
 {
     IComplex compl = new complex();
       complex a = new complex(2, 3);
       complex b = new complex(3, 4);
       complex actual = compl.Add(a, b);
       complex expected = a + b;
       Assert.AreEqual(actual.im, expected.im);
       Assert.AreEqual(actual.re, expected.re);
 }

 public void SubtractTest()
 {
     IComplex compl = new complex();
       complex a = new complex(20, 30);
       complex b = new complex(30, 40);
       complex actual = compl.Subtract(a, b);
       complex expected = a - b;
       Assert.AreEqual(actual.im, expected.im);
       Assert.AreEqual(actual.re, expected.re);
 }

 public void MultiplyTest()
 {
     IComplex compl = new complex();
       complex a = new complex(20, 30);
       complex b = new complex(30, 40);
       complex actual = compl.Multiply(a, b);
       complex expected = a * b;
       Assert.AreEqual(actual.im, expected.im);
       Assert.AreEqual(actual.re, expected.re);
 }

 public void DivideTest()
 {
     IComplex compl = new complex();
       complex a = new complex(20, 30);
       complex b = new complex(30, 40);
       complex actual = compl.Divide(a, b);
       complex expected = a / b;
       Assert.AreEqual(actual.im, expected.im);
       Assert.AreEqual(actual.re, expected.re);
 }

        private void Test_conj() {
            // test for common matrix 4x3
            ILArray<complex> A = new complex[,] {
                {new complex(1,2), new complex(3,4), new complex(5,6)},
                {new complex(7,8), new complex(9,10), new complex(11,12)},
                {new complex(13,14), new complex(15,16), new complex(17,18)},
                {new complex(19,20), new complex(21,22), new complex(23,24)}};
            ILArray<fcomplex> Af = ILMath.tofcomplex(A); 
            Test_isConj(A,ILMath.conj(A)); 
            Test_isfConj(Af,ILMath.conj(Af));
            // test on reference matrix
            Test_isConj(A,ILMath.conj(A.R)); 
            Test_isfConj(Af,ILMath.conj(Af.R));
            // test on 2x2x3 matrix (reference)
            A = A.Reshape(new ILDimension(2,2,3)).R;
            if (!A.IsReference) 
                throw new Exception("Unable to test conj() for reference 2x2x3!"); 
            Af = Af.Reshape(new ILDimension(2,2,3)).R; 
            if (!Af.IsReference) 
                throw new Exception("Unable to test conj() for reference 2x2x3!"); 
            Test_isConj(A,ILMath.conj(A)); 
            Test_isfConj(Af,ILMath.conj(Af));
            // test those on dense array 
            Test_isConj(A,ILMath.conj(A.C)); 
            Test_isfConj(Af,ILMath.conj(Af.C));
            // test for empty
            A = ILArray<complex>.empty();
            Af = ILArray<fcomplex>.empty();
            Test_isConj(A,ILMath.conj(A)); 
            Test_isfConj(Af,ILMath.conj(Af)); 
            // test on vector
            A = new complex[] { new complex(1,2), new complex(3,4), new complex(5,6)}; 
            Af = ILMath.tofcomplex(A); 
            Test_isConj(A,ILMath.conj(A)); 
            Test_isfConj(Af,ILMath.conj(Af)); 
            // transpose 
            Test_isConj(A.T,ILMath.conj(A.T)); 
            Test_isfConj(Af.T,ILMath.conj(Af.T)); 
            // test on scalar 
            A = A[2]; Af = Af[2]; 
            Test_isConj(A,ILMath.conj(A)); 
            Test_isfConj(Af,ILMath.conj(Af)); 

        }

    /*************************************************************************
    Dense solver. Same as RMatrixMixedSolve(), but for complex matrices.

    Algorithm features:
    * automatic detection of degenerate cases
    * condition number estimation
    * iterative refinement
    * O(N^2) complexity

    INPUT PARAMETERS
        A       -   array[0..N-1,0..N-1], system matrix
        LUA     -   array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
        P       -   array[0..N-1], pivots array, CMatrixLU result
        N       -   size of A
        B       -   array[0..N-1], right part

    OUTPUT PARAMETERS
        Info    -   same as in RMatrixSolveM
        Rep     -   same as in RMatrixSolveM
        X       -   same as in RMatrixSolveM

      -- ALGLIB --
         Copyright 27.01.2010 by Bochkanov Sergey
    *************************************************************************/
    public static void cmatrixmixedsolve(complex[,] a, complex[,] lua, int[] p, int n, complex[] b, out int info, out densesolverreport rep, out complex[] x)
    {
        info = 0;
        rep = new densesolverreport();
        x = new complex[0];
        densesolver.cmatrixmixedsolve(a, lua, p, n, b, ref info, rep.innerobj, ref x);
        return;
    }

 public static void smp_cmatrixsolve(complex[,] a, int n, complex[] b, out int info, out densesolverreport rep, out complex[] x)
 {
     info = 0;
     rep = new densesolverreport();
     x = new complex[0];
     densesolver._pexec_cmatrixsolve(a, n, b, ref info, rep.innerobj, ref x);
     return;
 }

        /*************************************************************************
        Basic LU solver for ScaleA*PLU*x = y.

        This subroutine assumes that:
        * L is well-scaled, and it is U which needs scaling by ScaleA.
        * A=PLU is well-conditioned, so no zero divisions or overflow may occur

          -- ALGLIB --
             Copyright 27.01.2010 by Bochkanov Sergey
        *************************************************************************/
        private static void cbasiclusolve(complex[,] lua,
            int[] p,
            double scalea,
            int n,
            ref complex[] xb,
            ref complex[] tmp)
        {
            int i = 0;
            complex v = 0;
            int i_ = 0;

            for(i=0; i<=n-1; i++)
            {
                if( p[i]!=i )
                {
                    v = xb[i];
                    xb[i] = xb[p[i]];
                    xb[p[i]] = v;
                }
            }
            for(i=1; i<=n-1; i++)
            {
                v = 0.0;
                for(i_=0; i_<=i-1;i_++)
                {
                    v += lua[i,i_]*xb[i_];
                }
                xb[i] = xb[i]-v;
            }
            xb[n-1] = xb[n-1]/(scalea*lua[n-1,n-1]);
            for(i=n-2; i>=0; i--)
            {
                for(i_=i+1; i_<=n-1;i_++)
                {
                    tmp[i_] = scalea*lua[i,i_];
                }
                v = 0.0;
                for(i_=i+1; i_<=n-1;i_++)
                {
                    v += tmp[i_]*xb[i_];
                }
                xb[i] = (xb[i]-v)/(scalea*lua[i,i]);
            }
        }

        /*************************************************************************
        Internal LU solver

          -- ALGLIB --
             Copyright 27.01.2010 by Bochkanov Sergey
        *************************************************************************/
        private static void cmatrixlusolveinternal(complex[,] lua,
            int[] p,
            double scalea,
            int n,
            complex[,] a,
            bool havea,
            complex[,] b,
            int m,
            ref int info,
            densesolverreport rep,
            ref complex[,] x)
        {
            int i = 0;
            int j = 0;
            int k = 0;
            int rfs = 0;
            int nrfs = 0;
            complex[] xc = new complex[0];
            complex[] y = new complex[0];
            complex[] bc = new complex[0];
            complex[] xa = new complex[0];
            complex[] xb = new complex[0];
            complex[] tx = new complex[0];
            double[] tmpbuf = new double[0];
            complex v = 0;
            double verr = 0;
            double mxb = 0;
            double scaleright = 0;
            bool smallerr = new bool();
            bool terminatenexttime = new bool();
            int i_ = 0;

            info = 0;
            x = new complex[0,0];

            alglib.ap.assert((double)(scalea)>(double)(0));
            
            //
            // prepare: check inputs, allocate space...
            //
            if( n<=0 || m<=0 )
            {
                info = -1;
                return;
            }
            for(i=0; i<=n-1; i++)
            {
                if( p[i]>n-1 || p[i]<i )
                {
                    info = -1;
                    return;
                }
            }
            x = new complex[n, m];
            y = new complex[n];
            xc = new complex[n];
            bc = new complex[n];
            tx = new complex[n];
            xa = new complex[n+1];
            xb = new complex[n+1];
            tmpbuf = new double[2*n+2];
            
            //
            // estimate condition number, test for near singularity
            //
            rep.r1 = rcond.cmatrixlurcond1(lua, n);
            rep.rinf = rcond.cmatrixlurcondinf(lua, n);
            if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) || (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
            {
                for(i=0; i<=n-1; i++)
                {
                    for(j=0; j<=m-1; j++)
                    {
                        x[i,j] = 0;
                    }
                }
                rep.r1 = 0;
                rep.rinf = 0;
                info = -3;
                return;
            }
            info = 1;
            
            //
            // solve
            //
            for(k=0; k<=m-1; k++)
            {
                
                //
                // copy B to contiguous storage
                //
                for(i_=0; i_<=n-1;i_++)
                {
//.........这里部分代码省略.........

        /*************************************************************************
        Dense solver. Same as RMatrixLUSolve(), but for  HPD matrices  represented
        by their Cholesky decomposition.

        Algorithm features:
        * automatic detection of degenerate cases
        * O(N^2) complexity
        * condition number estimation
        * matrix is represented by its upper or lower triangle

        No iterative refinement is provided because such partial representation of
        matrix does not allow efficient calculation of extra-precise  matrix-vector
        products for large matrices. Use RMatrixSolve or RMatrixMixedSolve  if  you
        need iterative refinement.

        INPUT PARAMETERS
            CHA     -   array[0..N-1,0..N-1], Cholesky decomposition,
                        SPDMatrixCholesky result
            N       -   size of A
            IsUpper -   what half of CHA is provided
            B       -   array[0..N-1], right part

        OUTPUT PARAMETERS
            Info    -   same as in RMatrixSolve
            Rep     -   same as in RMatrixSolve
            X       -   same as in RMatrixSolve

          -- ALGLIB --
             Copyright 27.01.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void hpdmatrixcholeskysolve(complex[,] cha,
            int n,
            bool isupper,
            complex[] b,
            ref int info,
            densesolverreport rep,
            ref complex[] x)
        {
            complex[,] bm = new complex[0,0];
            complex[,] xm = new complex[0,0];
            int i_ = 0;

            info = 0;
            x = new complex[0];

            if( n<=0 )
            {
                info = -1;
                return;
            }
            bm = new complex[n, 1];
            for(i_=0; i_<=n-1;i_++)
            {
                bm[i_,0] = b[i_];
            }
            hpdmatrixcholeskysolvem(cha, n, isupper, bm, 1, ref info, rep, ref xm);
            x = new complex[n];
            for(i_=0; i_<=n-1;i_++)
            {
                x[i_] = xm[i_,0];
            }
        }

 /*************************************************************************
 Single-threaded stub. HPC ALGLIB replaces it by multithreaded code.
 *************************************************************************/
 public static void _pexec_hpdmatrixsolve(complex[,] a,
     int n,
     bool isupper,
     complex[] b,
     ref int info,
     densesolverreport rep,
     ref complex[] x)
 {
     hpdmatrixsolve(a,n,isupper,b,ref info,rep,ref x);
 }

        /*************************************************************************
        1-dimensional real inverse FFT.

        Algorithm has O(N*logN) complexity for any N (composite or prime).

        INPUT PARAMETERS
            F   -   array[0..floor(N/2)] - frequencies from forward real FFT
            N   -   problem size

        OUTPUT PARAMETERS
            A   -   inverse DFT of a input array, array[0..N-1]

        NOTE:
            F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just  one
        half of frequencies array is needed - elements from 0 to floor(N/2).  F[0]
        is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd,  then
        F[floor(N/2)] has no special properties.

        Relying on properties noted above, FFTR1DInv subroutine uses only elements
        from 0th to floor(N/2)-th. It ignores imaginary part of F[0],  and in case
        N is even it ignores imaginary part of F[floor(N/2)] too.

        When you call this function using full arguments list - "FFTR1DInv(F,N,A)"
        - you can pass either either frequencies array with N elements or  reduced
        array with roughly N/2 elements - subroutine will  successfully  transform
        both.

        If you call this function using reduced arguments list -  "FFTR1DInv(F,A)"
        - you must pass FULL array with N elements (although higher  N/2 are still
        not used) because array size is used to automatically determine FFT length


          -- ALGLIB --
             Copyright 01.06.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void fftr1dinv(complex[] f,
            int n,
            ref double[] a)
        {
            int i = 0;
            double[] h = new double[0];
            complex[] fh = new complex[0];

            a = new double[0];

            ap.assert(n>0, "FFTR1DInv: incorrect N!");
            ap.assert(ap.len(f)>=(int)Math.Floor((double)n/(double)2)+1, "FFTR1DInv: Length(F)<Floor(N/2)+1!");
            ap.assert(math.isfinite(f[0].x), "FFTR1DInv: F contains infinite or NAN values!");
            for(i=1; i<=(int)Math.Floor((double)n/(double)2)-1; i++)
            {
                ap.assert(math.isfinite(f[i].x) & math.isfinite(f[i].y), "FFTR1DInv: F contains infinite or NAN values!");
            }
            ap.assert(math.isfinite(f[(int)Math.Floor((double)n/(double)2)].x), "FFTR1DInv: F contains infinite or NAN values!");
            if( n%2!=0 )
            {
                ap.assert(math.isfinite(f[(int)Math.Floor((double)n/(double)2)].y), "FFTR1DInv: F contains infinite or NAN values!");
            }
            
            //
            // Special case: N=1, FFT is just identity transform.
            // After this block we assume that N is strictly greater than 1.
            //
            if( n==1 )
            {
                a = new double[1];
                a[0] = f[0].x;
                return;
            }
            
            //
            // inverse real FFT is reduced to the inverse real FHT,
            // which is reduced to the forward real FHT,
            // which is reduced to the forward real FFT.
            //
            // Don't worry, it is really compact and efficient reduction :)
            //
            h = new double[n];
            a = new double[n];
            h[0] = f[0].x;
            for(i=1; i<=(int)Math.Floor((double)n/(double)2)-1; i++)
            {
                h[i] = f[i].x-f[i].y;
                h[n-i] = f[i].x+f[i].y;
            }
            if( n%2==0 )
            {
                h[(int)Math.Floor((double)n/(double)2)] = f[(int)Math.Floor((double)n/(double)2)].x;
            }
            else
            {
                h[(int)Math.Floor((double)n/(double)2)] = f[(int)Math.Floor((double)n/(double)2)].x-f[(int)Math.Floor((double)n/(double)2)].y;
                h[(int)Math.Floor((double)n/(double)2)+1] = f[(int)Math.Floor((double)n/(double)2)].x+f[(int)Math.Floor((double)n/(double)2)].y;
            }
            fftr1d(h, n, ref fh);
            for(i=0; i<=n-1; i++)
            {
                a[i] = (fh[i].x-fh[i].y)/n;
            }
        }

    /*************************************************************************
    1-dimensional complex inverse FFT.

    Array size N may be arbitrary number (composite or prime).  Algorithm  has
    O(N*logN) complexity for any N (composite or prime).

    See FFTC1D() description for more information about algorithm performance.

    INPUT PARAMETERS
        A   -   array[0..N-1] - complex array to be transformed
        N   -   problem size

    OUTPUT PARAMETERS
        A   -   inverse DFT of a input array, array[0..N-1]
                A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), k = 0..N-1)


      -- ALGLIB --
         Copyright 29.05.2009 by Bochkanov Sergey
    *************************************************************************/
    public static void fftc1dinv(ref complex[] a, int n)
    {

        fft.fftc1dinv(ref a, n);
        return;
    }

        /*************************************************************************
        Dense solver. Same as RMatrixSolveM(), but for Hermitian positive definite
        matrices.

        Algorithm features:
        * automatic detection of degenerate cases
        * condition number estimation
        * O(N^3+M*N^2) complexity
        * matrix is represented by its upper or lower triangle

        No iterative refinement is provided because such partial representation of
        matrix does not allow efficient calculation of extra-precise  matrix-vector
        products for large matrices. Use RMatrixSolve or RMatrixMixedSolve  if  you
        need iterative refinement.

        COMMERCIAL EDITION OF ALGLIB:

          ! Commercial version of ALGLIB includes two  important  improvements  of
          ! this function, which can be used from C++ and C#:
          ! * Intel MKL support (lightweight Intel MKL is shipped with ALGLIB)
          ! * multicore support
          !
          ! Intel MKL gives approximately constant  (with  respect  to  number  of
          ! worker threads) acceleration factor which depends on CPU  being  used,
          ! problem  size  and  "baseline"  ALGLIB  edition  which  is  used   for
          ! comparison.
          !
          ! Say, on SSE2-capable CPU with N=1024, HPC ALGLIB will be:
          ! * about 2-3x faster than ALGLIB for C++ without MKL
          ! * about 7-10x faster than "pure C#" edition of ALGLIB
          ! Difference in performance will be more striking  on  newer  CPU's with
          ! support for newer SIMD instructions. Generally,  MKL  accelerates  any
          ! problem whose size is at least 128, with best  efficiency achieved for
          ! N's larger than 512.
          !
          ! Commercial edition of ALGLIB also supports multithreaded  acceleration
          ! of this function. We should note that Cholesky decomposition is harder
          ! to parallelize than, say, matrix-matrix product - this  algorithm  has
          ! several synchronization points which  can  not  be  avoided.  However,
          ! parallelism starts to be profitable starting from N=500.
          !
          ! In order to use multicore features you have to:
          ! * use commercial version of ALGLIB
          ! * call  this  function  with  "smp_"  prefix,  which  indicates  that
          !   multicore code will be used (for multicore support)
          !
          ! We recommend you to read 'Working with commercial version' section  of
          ! ALGLIB Reference Manual in order to find out how to  use  performance-
          ! related features provided by commercial edition of ALGLIB.

        INPUT PARAMETERS
            A       -   array[0..N-1,0..N-1], system matrix
            N       -   size of A
            IsUpper -   what half of A is provided
            B       -   array[0..N-1,0..M-1], right part
            M       -   right part size

        OUTPUT PARAMETERS
            Info    -   same as in RMatrixSolve.
                        Returns -3 for non-HPD matrices.
            Rep     -   same as in RMatrixSolve
            X       -   same as in RMatrixSolve

          -- ALGLIB --
             Copyright 27.01.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void hpdmatrixsolvem(complex[,] a,
            int n,
            bool isupper,
            complex[,] b,
            int m,
            ref int info,
            densesolverreport rep,
            ref complex[,] x)
        {
            complex[,] da = new complex[0,0];
            double sqrtscalea = 0;
            int i = 0;
            int j = 0;
            int j1 = 0;
            int j2 = 0;
            int i_ = 0;

            info = 0;
            x = new complex[0,0];

            
            //
            // prepare: check inputs, allocate space...
            //
            if( n<=0 || m<=0 )
            {
                info = -1;
                return;
            }
            da = new complex[n, n];
            
            //
            // 1. scale matrix, max(|A[i,j]|)
            // 2. factorize scaled matrix
//.........这里部分代码省略.........

 public static void smp_hpdmatrixsolvem(complex[,] a, int n, bool isupper, complex[,] b, int m, out int info, out densesolverreport rep, out complex[,] x)
 {
     info = 0;
     rep = new densesolverreport();
     x = new complex[0,0];
     densesolver._pexec_hpdmatrixsolvem(a, n, isupper, b, m, ref info, rep.innerobj, ref x);
     return;
 }

        /*************************************************************************
        Dense solver. Same as RMatrixLUSolveM(), but for HPD matrices  represented
        by their Cholesky decomposition.

        Algorithm features:
        * automatic detection of degenerate cases
        * O(M*N^2) complexity
        * condition number estimation
        * matrix is represented by its upper or lower triangle

        No iterative refinement is provided because such partial representation of
        matrix does not allow efficient calculation of extra-precise  matrix-vector
        products for large matrices. Use RMatrixSolve or RMatrixMixedSolve  if  you
        need iterative refinement.

        INPUT PARAMETERS
            CHA     -   array[0..N-1,0..N-1], Cholesky decomposition,
                        HPDMatrixCholesky result
            N       -   size of CHA
            IsUpper -   what half of CHA is provided
            B       -   array[0..N-1,0..M-1], right part
            M       -   right part size

        OUTPUT PARAMETERS
            Info    -   same as in RMatrixSolve
            Rep     -   same as in RMatrixSolve
            X       -   same as in RMatrixSolve

          -- ALGLIB --
             Copyright 27.01.2010 by Bochkanov Sergey
        *************************************************************************/
        public static void hpdmatrixcholeskysolvem(complex[,] cha,
            int n,
            bool isupper,
            complex[,] b,
            int m,
            ref int info,
            densesolverreport rep,
            ref complex[,] x)
        {
            complex[,] emptya = new complex[0,0];
            double sqrtscalea = 0;
            int i = 0;
            int j = 0;
            int j1 = 0;
            int j2 = 0;

            info = 0;
            x = new complex[0,0];

            
            //
            // prepare: check inputs, allocate space...
            //
            if( n<=0 || m<=0 )
            {
                info = -1;
                return;
            }
            
            //
            // 1. scale matrix, max(|U[i,j]|)
            // 2. factorize scaled matrix
            // 3. solve
            //
            sqrtscalea = 0;
            for(i=0; i<=n-1; i++)
            {
                if( isupper )
                {
                    j1 = i;
                    j2 = n-1;
                }
                else
                {
                    j1 = 0;
                    j2 = i;
                }
                for(j=j1; j<=j2; j++)
                {
                    sqrtscalea = Math.Max(sqrtscalea, math.abscomplex(cha[i,j]));
                }
            }
            if( (double)(sqrtscalea)==(double)(0) )
            {
                sqrtscalea = 1;
            }
            sqrtscalea = 1/sqrtscalea;
            hpdmatrixcholeskysolveinternal(cha, sqrtscalea, n, isupper, emptya, false, b, m, ref info, rep, ref x);
        }

        /*************************************************************************
        Internal Cholesky solver

          -- ALGLIB --
             Copyright 27.01.2010 by Bochkanov Sergey
        *************************************************************************/
        private static void hpdmatrixcholeskysolveinternal(complex[,] cha,
            double sqrtscalea,
            int n,
            bool isupper,
            complex[,] a,
            bool havea,
            complex[,] b,
            int m,
            ref int info,
            densesolverreport rep,
            ref complex[,] x)
        {
            int i = 0;
            int j = 0;
            int k = 0;
            complex[] xc = new complex[0];
            complex[] y = new complex[0];
            complex[] bc = new complex[0];
            complex[] xa = new complex[0];
            complex[] xb = new complex[0];
            complex[] tx = new complex[0];
            double v = 0;
            double mxb = 0;
            double scaleright = 0;
            int i_ = 0;

            info = 0;
            x = new complex[0,0];

            alglib.ap.assert((double)(sqrtscalea)>(double)(0));
            
            //
            // prepare: check inputs, allocate space...
            //
            if( n<=0 || m<=0 )
            {
                info = -1;
                return;
            }
            x = new complex[n, m];
            y = new complex[n];
            xc = new complex[n];
            bc = new complex[n];
            tx = new complex[n+1];
            xa = new complex[n+1];
            xb = new complex[n+1];
            
            //
            // estimate condition number, test for near singularity
            //
            rep.r1 = rcond.hpdmatrixcholeskyrcond(cha, n, isupper);
            rep.rinf = rep.r1;
            if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) )
            {
                for(i=0; i<=n-1; i++)
                {
                    for(j=0; j<=m-1; j++)
                    {
                        x[i,j] = 0;
                    }
                }
                rep.r1 = 0;
                rep.rinf = 0;
                info = -3;
                return;
            }
            info = 1;
            
            //
            // solve
            //
            for(k=0; k<=m-1; k++)
            {
                
                //
                // copy B to contiguous storage
                //
                for(i_=0; i_<=n-1;i_++)
                {
                    bc[i_] = b[i_,k];
                }
                
                //
                // Scale right part:
                // * MX stores max(|Bi|)
                // * ScaleRight stores actual scaling applied to B when solving systems
                //   it is chosen to make |scaleRight*b| close to 1.
                //
                mxb = 0;
                for(i=0; i<=n-1; i++)
                {
                    mxb = Math.Max(mxb, math.abscomplex(bc[i]));
                }
                if( (double)(mxb)==(double)(0) )
//.........这里部分代码省略.........

 /*************************************************************************
 Single-threaded stub. HPC ALGLIB replaces it by multithreaded code.
 *************************************************************************/
 public static void _pexec_cmatrixsolvem(complex[,] a,
     int n,
     complex[,] b,
     int m,
     bool rfs,
     ref int info,
     densesolverreport rep,
     ref complex[,] x)
 {
     cmatrixsolvem(a,n,b,m,rfs,ref info,rep,ref x);
 }