/// <summary>
        /// Initializes a new instance of the <see cref="DenseQR"/> class. This object will compute the
        /// QR factorization when the constructor is called and cache it's factorization.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="method">The QR factorization method to use.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> row count is less then column count</exception>
        public static DenseQR Create(DenseMatrix matrix, QRMethod method = QRMethod.Full)
        {
            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(matrix);
            }

            var tau = new Complex32[Math.Min(matrix.RowCount, matrix.ColumnCount)];
            Matrix<Complex32> q;
            Matrix<Complex32> r;

            if (method == QRMethod.Full)
            {
                r = matrix.Clone();
                q = new DenseMatrix(matrix.RowCount);
                Control.LinearAlgebraProvider.QRFactor(((DenseMatrix) r).Values, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix) q).Values, tau);
            }
            else
            {
                q = matrix.Clone();
                r = new DenseMatrix(matrix.ColumnCount);
                Control.LinearAlgebraProvider.ThinQRFactor(((DenseMatrix) q).Values, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix) r).Values, tau);
            }

            return new DenseQR(q, r, method, tau);
        }

        /// <summary>
        /// Adds a scaled vector to another: <c>result = y + alpha*x</c>.
        /// </summary>
        /// <param name="y">The vector to update.</param>
        /// <param name="alpha">The value to scale <paramref name="x"/> by.</param>
        /// <param name="x">The vector to add to <paramref name="y"/>.</param>
        /// <param name="result">The result of the addition.</param>
        /// <remarks>This is similar to the AXPY BLAS routine.</remarks>
        public virtual void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result)
        {
            if (y == null)
            {
                throw new ArgumentNullException("y");
            }

            if (x == null)
            {
                throw new ArgumentNullException("x");
            }

            if (y.Length != x.Length)
            {
                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
            }

            if (alpha.IsZero())
            {
                CommonParallel.For(0, y.Length, index => result[index] = y[index]);
            }
            else if (alpha.IsOne())
            {
                CommonParallel.For(0, y.Length, index => result[index] = y[index] + x[index]);
            }
            else
            {
                CommonParallel.For(0, y.Length, index => result[index] = y[index] + (alpha * x[index]));
            }
        }

        /// <summary>
        /// Returns the smallest absolute value from the unsorted data array.
        /// Returns NaN if data is empty or any entry is NaN.
        /// </summary>
        /// <param name="data">Sample array, no sorting is assumed.</param>
        public static Complex32 MinimumMagnitudePhase(Complex32[] data)
        {
            if (data.Length == 0)
            {
                return new Complex32(float.NaN, float.NaN);
            }

            float minMagnitude = float.PositiveInfinity;
            Complex32 min = new Complex32(float.PositiveInfinity, float.PositiveInfinity);
            for (int i = 0; i < data.Length; i++)
            {
                float magnitude = data[i].Magnitude;
                if (float.IsNaN(magnitude))
                {
                    return new Complex32(float.NaN, float.NaN);
                }
                if (magnitude < minMagnitude || magnitude == minMagnitude && data[i].Phase < min.Phase)
                {
                    minMagnitude = magnitude;
                    min = data[i];
                }
            }

            return min;
        }

 /// <summary>
 /// Initializes a new instance of the <see cref="DenseVector"/> class with a given size
 /// and each element set to the given value;
 /// </summary>
 /// <param name="size">
 /// the size of the vector.
 /// </param>
 /// <param name="value">
 /// the value to set each element to.
 /// </param>
 /// <exception cref="ArgumentException">
 /// If <paramref name="size"/> is less than one.
 /// </exception>
 public DenseVector(int size, Complex32 value) : this(size)
 {
     for (var index = 0; index < Data.Length; index++)
     {
         Data[index] = value;
     }
 }

        public override double MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (rows <= 0)
            {
                throw new ArgumentException(Resources.ArgumentMustBePositive, "rows");
            }

            if (columns <= 0)
            {
                throw new ArgumentException(Resources.ArgumentMustBePositive, "columns");
            }

            if (matrix.Length < rows * columns)
            {
                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, rows * columns), "matrix");
            }

            var work = new float[rows];
            return SafeNativeMethods.c_matrix_norm((byte)norm, rows, columns, matrix, work);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserQR"/> class. This object will compute the
        /// QR factorization when the constructor is called and cache it's factorization.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="method">The QR factorization method to use.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        public static UserQR Create(Matrix<Complex32> matrix, QRMethod method = QRMethod.Full)
        {
            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(matrix);
            }

            Matrix<Complex32> q;
            Matrix<Complex32> r;

            var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount);
            var u = new Complex32[minmn][];

            if (method == QRMethod.Full)
            {
                r = matrix.Clone();
                q = Matrix<Complex32>.Build.SameAs(matrix, matrix.RowCount, matrix.RowCount);

                for (var i = 0; i < matrix.RowCount; i++)
                {
                    q.At(i, i, 1.0f);
                }

                for (var i = 0; i < minmn; i++)
                {
                    u[i] = GenerateColumn(r, i, i);
                    ComputeQR(u[i], r, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.MaxDegreeOfParallelism);
                }

                for (var i = minmn - 1; i >= 0; i--)
                {
                    ComputeQR(u[i], q, i, matrix.RowCount, i, matrix.RowCount, Control.MaxDegreeOfParallelism);
                }
            }
            else
            {
                q = matrix.Clone();

                for (var i = 0; i < minmn; i++)
                {
                    u[i] = GenerateColumn(q, i, i);
                    ComputeQR(u[i], q, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.MaxDegreeOfParallelism);
                }

                r = q.SubMatrix(0, matrix.ColumnCount, 0, matrix.ColumnCount);
                q.Clear();

                for (var i = 0; i < matrix.ColumnCount; i++)
                {
                    q.At(i, i, 1.0f);
                }

                for (var i = minmn - 1; i >= 0; i--)
                {
                    ComputeQR(u[i], q, i, matrix.RowCount, i, matrix.ColumnCount, Control.MaxDegreeOfParallelism);
                }
            }

            return new UserQR(q, r, method);
        }

        public override void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result)
        {
            if (y == null)
            {
                throw new ArgumentNullException("y");
            }

            if (x == null)
            {
                throw new ArgumentNullException("x");
            }

            if (y.Length != x.Length)
            {
                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
            }

            if (!ReferenceEquals(y, result))
            {
                Array.Copy(y, 0, result, 0, y.Length);
            }

            if (alpha == Complex32.Zero)
            {
                return;
            }

            SafeNativeMethods.c_axpy(y.Length, alpha, x, result);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="DenseQR"/> class. This object will compute the
        /// QR factorization when the constructor is called and cache it's factorization.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="method">The QR factorization method to use.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> row count is less then column count</exception>
        public DenseQR(DenseMatrix matrix, QRMethod method = QRMethod.Full)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(matrix);
            }

            Tau = new Complex32[Math.Min(matrix.RowCount, matrix.ColumnCount)];

            if (method == QRMethod.Full)
            {
                MatrixR = matrix.Clone();
                MatrixQ = new DenseMatrix(matrix.RowCount);
                Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Values, matrix.RowCount, matrix.ColumnCount,
                                                       ((DenseMatrix)MatrixQ).Values, Tau);
            }
            else
            {
                MatrixQ = matrix.Clone();
                MatrixR = new DenseMatrix(matrix.ColumnCount);
                Control.LinearAlgebraProvider.ThinQRFactor(((DenseMatrix)MatrixQ).Values, matrix.RowCount, matrix.ColumnCount,
                                                       ((DenseMatrix)MatrixR).Values, Tau);
            }
        }

 /// <summary>
 /// Initializes a new instance of the <see cref="DenseMatrix"/> class. This matrix is square with a given size.
 /// </summary>
 /// <param name="order">the size of the square matrix.</param>
 /// <exception cref="ArgumentException">
 /// If <paramref name="order"/> is less than one.
 /// </exception>
 public DenseMatrix(int order)
     : base(order)
 {
     _rowCount = order;
     _columnCount = order;
     Data = new Complex32[order * order];
 }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserQR"/> class. This object will compute the
        /// QR factorization when the constructor is called and cache it's factorization.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        public UserQR(Matrix<Complex32> matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
            }

            MatrixR = matrix.Clone();
            MatrixQ = matrix.CreateMatrix(matrix.RowCount, matrix.RowCount);

            for (var i = 0; i < matrix.RowCount; i++)
            {
                MatrixQ.At(i, i, 1.0f);
            }

            var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount);
            var u = new Complex32[minmn][];
            for (var i = 0; i < minmn; i++)
            {
                u[i] = GenerateColumn(MatrixR, i, i);
                ComputeQR(u[i], MatrixR, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
            }

            for (var i = minmn - 1; i >= 0; i--)
            {
                ComputeQR(u[i], MatrixQ, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
            }
        }

        /// <summary>
        /// Does a point wise add of two arrays <c>z = x + y</c>. This can be used
        /// to add vectors or matrices.
        /// </summary>
        /// <param name="x">The array x.</param>
        /// <param name="y">The array y.</param>
        /// <param name="result">The result of the addition.</param>
        /// <remarks>There is no equivalent BLAS routine, but many libraries
        /// provide optimized (parallel and/or vectorized) versions of this
        /// routine.</remarks>
        public virtual void AddArrays(Complex32[] x, Complex32[] y, Complex32[] result)
        {
            if (y == null)
            {
                throw new ArgumentNullException("y");
            }

            if (x == null)
            {
                throw new ArgumentNullException("x");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (y.Length != x.Length || y.Length != result.Length)
            {
                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
            }

            CommonParallel.For(0, y.Length, 4096, (a, b) =>
            {
                for (int i = a; i < b; i++)
                {
                    result[i] = x[i] + y[i];
                }
            });
        }

        /// <summary>
        /// Does a point wise add of two arrays <c>z = x + y</c>. This can be used 
        /// to add vectors or matrices.
        /// </summary>
        /// <param name="x">The array x.</param>
        /// <param name="y">The array y.</param>
        /// <param name="result">The result of the addition.</param>
        /// <remarks>There is no equivalent BLAS routine, but many libraries
        /// provide optimized (parallel and/or vectorized) versions of this
        /// routine.</remarks>
        public virtual void AddArrays(Complex32[] x, Complex32[] y, Complex32[] result)
        {
            if (y == null)
            {
                throw new ArgumentNullException("y");
            }

            if (x == null)
            {
                throw new ArgumentNullException("x");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (y.Length != x.Length || y.Length != result.Length)
            {
                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
            }

            if (Control.ParallelizeOperation(x.Length))
            {
                CommonParallel.For(0, y.Length, index => { result[index] = x[index] + y[index]; });
            }
            else
            {
                for (var index = 0; index < x.Length; index++)
                {
                    result[index] = x[index] + y[index];
                }
            }
        }

        public virtual Complex32 TheveninVoltage(Node referenceNode, Complex32 ?W = null)
        {
            Complex32 v = 0;
            Dipole compo1 = null;
            Node node1 = null;
            foreach (var item in Components)
            {
                if (item.Nodes[0] == referenceNode || item.Nodes[1] == referenceNode)
                {
                    compo1 = item;
                    break;
                }
            }
            node1 = referenceNode;
            do
            {
                if (compo1 is Branch)
                    v += ((Branch)compo1).TheveninVoltage(node1, W);                    
                else
                    v += compo1.voltage(node1, W);
                node1 = compo1.OtherNode(node1);
                compo1 = node1.OtherComponent(compo1);
            } while (InternalNodes.Contains(node1));

            return v;
        }

 public override Complex32 Current(Node referenceNode, Complex32? W = null)
 {
     Complex32 i = Voltage / Impedance(W);
     if (referenceNode == Nodes[0])
         return i;
     else
         return -i;
 }

 public override Complex32 voltage(Node ReferenceNode, Complex32 ?W)
 {
     if (ReferenceNode == Nodes[0])
         return Voltage;
     if (ReferenceNode == Nodes[1])
         return -Voltage;
     return Complex32.NaN;
 }

 public override Complex32 Impedance(Complex32? W = null)
 {
     if (W == null || W.Value.IsZero())
         return Complex32.PositiveInfinity;
     // 1/jWC
     // 1/SC
     return (W.Value * new Complex32((float)Value,0)).Reciprocal();
 }

 public virtual Complex32 Current(Node referenceNode, Complex32? W = null)
 {
     //el signo contrario al pensado, entra por neagtivo y sale por positivo
     if (referenceNode == Nodes[0])
         return -current;
     else
         return current;
 }

 /// <summary>
 /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class with all entries set to a particular value.
 /// </summary>
 /// <param name="rows">
 /// The number of rows.
 /// </param>
 /// <param name="columns">
 /// The number of columns.
 /// </param>
 /// <param name="value">The value which we assign to each element of the matrix.</param>
 public DiagonalMatrix(int rows, int columns, Complex32 value)
     : base(rows, columns)
 {
     Data = new Complex32[Math.Min(rows, columns)];
     for (var i = 0; i < Data.Length; i++)
     {
         Data[i] = value;
     }
 }

 /// <summary>
 /// Initializes a new instance of the <see cref="DenseMatrix"/> class with all entries set to a particular value.
 /// </summary>
 /// <param name="rows">
 /// The number of rows.
 /// </param>
 /// <param name="columns">
 /// The number of columns.
 /// </param>
 /// <param name="value">The value which we assign to each element of the matrix.</param>
 public DenseMatrix(int rows, int columns, Complex32 value)
     : base(rows, columns)
 {
     Data = new Complex32[rows * columns];
     for (var i = 0; i < Data.Length; i++)
     {
         Data[i] = value;
     }
 }

 public void CanAddArrays()
 {
     var result = new Complex32[_y.Length];
     Control.LinearAlgebraProvider.AddArrays(_x, _y, result);
     for (var i = 0; i < result.Length; i++)
     {
         Assert.AreEqual(_x[i] + _y[i], result[i]);
     }
 }