C# Matrix.At方法代码示例

        /// <summary>
        /// Initializes a new instance of the <see cref="UserEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public static UserEvd Create(Matrix<Complex32> matrix)
        {
            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            var eigenVectors = DenseMatrix.CreateIdentity(order);
            var blockDiagonal = matrix.CreateMatrix(order, order);
            var eigenValues = new LinearAlgebra.Complex.DenseVector(order);

            var isSymmetric = true;

            for (var i = 0; isSymmetric && i < order; i++)
            {
                for (var j = 0; isSymmetric && j < order; j++)
                {
                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
                }
            }

            if (isSymmetric)
            {
                var matrixCopy = matrix.ToArray();
                var tau = new Complex32[order];
                var d = new float[order];
                var e = new float[order];

                SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
                SymmetricDiagonalize(eigenVectors, d, e, order);
                SymmetricUntridiagonalize(eigenVectors, matrixCopy, tau, order);

                for (var i = 0; i < order; i++)
                {
                    eigenValues[i] = new Complex(d[i], e[i]);
                }
            }
            else
            {
                var matrixH = matrix.ToArray();
                NonsymmetricReduceToHessenberg(eigenVectors, matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(eigenVectors, eigenValues, matrixH, order);
            }

            for (var i = 0; i < eigenValues.Count; i++)
            {
                blockDiagonal.At(i, i, (Complex32) eigenValues[i]);
            }

            return new UserEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public UserEvd(Matrix<Complex> matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            MatrixEv = DenseMatrix.Identity(order);
            MatrixD = matrix.CreateMatrix(order, order);
            VectorEv = new DenseVector(order);
           
            IsSymmetric = true;

            for (var i = 0; IsSymmetric && i < order; i++)
            {
                for (var j = 0; IsSymmetric && j < order; j++)
                {
                    IsSymmetric &= matrix.At(i, j) == matrix.At(j, i).Conjugate();
                }
            }

            if (IsSymmetric)
            {
                var matrixCopy = matrix.ToArray();
                var tau = new Complex[order];
                var d = new double[order];
                var e = new double[order];

                SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
                SymmetricDiagonalize(d, e, order);
                SymmetricUntridiagonalize(matrixCopy, tau, order);

                for (var i = 0; i < order; i++)
                {
                    VectorEv[i] = new Complex(d[i], e[i]);
                }
            }
            else
            {
                var matrixH = matrix.ToArray();
                NonsymmetricReduceToHessenberg(matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(matrixH, order);
            }

            MatrixD.SetDiagonal(VectorEv);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserCholesky"/> class. This object will compute the
        /// Cholesky 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>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> is not positive definite.</exception>
        public UserCholesky(Matrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
            CholeskyFactor = matrix.Clone();
            for (var j = 0; j < CholeskyFactor.RowCount; j++)
            {
                var d = 0.0;
                for (var k = 0; k < j; k++)
                {
                    var s = 0.0;
                    for (var i = 0; i < k; i++)
                    {
                        s += CholeskyFactor.At(k, i) * CholeskyFactor.At(j, i);
                    }

                    s = (matrix.At(j, k) - s) / CholeskyFactor.At(k, k);
                    CholeskyFactor.At(j, k, s);
                    d += s * s;
                }

                d = matrix.At(j, j) - d;
                if (d <= 0.0)
                {
                    throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                }

                CholeskyFactor.At(j, j, Math.Sqrt(d));
                for (var k = j + 1; k < CholeskyFactor.RowCount; k++)
                {
                    CholeskyFactor.At(j, k, 0.0);
                }
            }
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserEvd"/> class. This object will compute the
        /// the eigenvalue decomposition when the constructor is called and cache it's decomposition.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If EVD algorithm failed to converge with matrix <paramref name="matrix"/>.</exception>
        public static UserEvd Create(Matrix<float> matrix)
        {
            if (matrix.RowCount != matrix.ColumnCount)
            {
                throw new ArgumentException(Resources.ArgumentMatrixSquare);
            }

            var order = matrix.RowCount;

            // Initialize matricies for eigenvalues and eigenvectors
            var eigenVectors = Matrix<float>.Build.SameAs(matrix, order, order);
            var blockDiagonal = Matrix<float>.Build.SameAs(matrix, order, order);
            var eigenValues = new LinearAlgebra.Complex.DenseVector(order);

            var isSymmetric = true;

            for (var i = 0; isSymmetric && i < order; i++)
            {
                for (var j = 0; isSymmetric && j < order; j++)
                {
                    isSymmetric &= matrix.At(i, j) == matrix.At(j, i);
                }
            }

            var d = new float[order];
            var e = new float[order];

            if (isSymmetric)
            {
                matrix.CopyTo(eigenVectors);
                d = eigenVectors.Row(order - 1).ToArray();

                SymmetricTridiagonalize(eigenVectors, d, e, order);
                SymmetricDiagonalize(eigenVectors, d, e, order);
            }
            else
            {
                var matrixH = matrix.ToArray();

                NonsymmetricReduceToHessenberg(eigenVectors, matrixH, order);
                NonsymmetricReduceHessenberToRealSchur(eigenVectors, matrixH, d, e, order);
            }

            for (var i = 0; i < order; i++)
            {
                blockDiagonal.At(i, i, d[i]);

                if (e[i] > 0)
                {
                    blockDiagonal.At(i, i + 1, e[i]);
                }
                else if (e[i] < 0)
                {
                    blockDiagonal.At(i, i - 1, e[i]);
                }
            }

            for (var i = 0; i < order; i++)
            {
                eigenValues[i] = new Complex(d[i], e[i]);
            }

            return new UserEvd(eigenVectors, eigenValues, blockDiagonal, isSymmetric);
        }

//.........这里部分代码省略.........
                        {
                            d[n] = x - (w/z);
                        }

                        e[n - 1] = 0.0f;
                        e[n] = 0.0f;
                        x = matrixH[n, n - 1];
                        s = Math.Abs(x) + Math.Abs(z);
                        p = x/s;
                        q = z/s;
                        r = (float) Math.Sqrt((p*p) + (q*q));
                        p = p/r;
                        q = q/r;

                        // Row modification
                        for (var j = n - 1; j < order; j++)
                        {
                            z = matrixH[n - 1, j];
                            matrixH[n - 1, j] = (q*z) + (p*matrixH[n, j]);
                            matrixH[n, j] = (q*matrixH[n, j]) - (p*z);
                        }

                        // Column modification
                        for (var i = 0; i <= n; i++)
                        {
                            z = matrixH[i, n - 1];
                            matrixH[i, n - 1] = (q*z) + (p*matrixH[i, n]);
                            matrixH[i, n] = (q*matrixH[i, n]) - (p*z);
                        }

                        // Accumulate transformations
                        for (var i = 0; i < order; i++)
                        {
                            z = eigenVectors.At(i, n - 1);
                            eigenVectors.At(i, n - 1, (q*z) + (p*eigenVectors.At(i, n)));
                            eigenVectors.At(i, n, (q*eigenVectors.At(i, n)) - (p*z));
                        }

                        // Complex pair
                    }
                    else
                    {
                        d[n - 1] = x + p;
                        d[n] = x + p;
                        e[n - 1] = z;
                        e[n] = -z;
                    }

                    n = n - 2;
                    iter = 0;

                    // No convergence yet
                }
                else
                {
                    // Form shift
                    x = matrixH[n, n];
                    y = 0.0f;
                    w = 0.0f;
                    if (l < n)
                    {
                        y = matrixH[n - 1, n - 1];
                        w = matrixH[n, n - 1]*matrixH[n - 1, n];
                    }

                    // Wilkinson's original ad hoc shift

        /// <summary>
        /// Symmetric Householder reduction to tridiagonal form.
        /// </summary>
        /// <param name="eigenVectors">The eigen vectors to work on.</param>
        /// <param name="d">Arrays for internal storage of real parts of eigenvalues</param>
        /// <param name="e">Arrays for internal storage of imaginary parts of eigenvalues</param>
        /// <param name="order">Order of initial matrix</param>
        /// <remarks>This is derived from the Algol procedures tred2 by 
        /// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 
        /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 
        /// Fortran subroutine in EISPACK.</remarks>
        static void SymmetricTridiagonalize(Matrix<float> eigenVectors, float[] d, float[] e, int order)
        {
            // Householder reduction to tridiagonal form.
            for (var i = order - 1; i > 0; i--)
            {
                // Scale to avoid under/overflow.
                var scale = 0.0f;
                var h = 0.0f;

                for (var k = 0; k < i; k++)
                {
                    scale = scale + Math.Abs(d[k]);
                }

                if (scale == 0.0f)
                {
                    e[i] = d[i - 1];
                    for (var j = 0; j < i; j++)
                    {
                        d[j] = eigenVectors.At(i - 1, j);
                        eigenVectors.At(i, j, 0.0f);
                        eigenVectors.At(j, i, 0.0f);
                    }
                }
                else
                {
                    // Generate Householder vector.
                    for (var k = 0; k < i; k++)
                    {
                        d[k] /= scale;
                        h += d[k]*d[k];
                    }

                    var f = d[i - 1];
                    var g = (float) Math.Sqrt(h);
                    if (f > 0)
                    {
                        g = -g;
                    }

                    e[i] = scale*g;
                    h = h - (f*g);
                    d[i - 1] = f - g;

                    for (var j = 0; j < i; j++)
                    {
                        e[j] = 0.0f;
                    }

                    // Apply similarity transformation to remaining columns.
                    for (var j = 0; j < i; j++)
                    {
                        f = d[j];
                        eigenVectors.At(j, i, f);
                        g = e[j] + (eigenVectors.At(j, j)*f);

                        for (var k = j + 1; k <= i - 1; k++)
                        {
                            g += eigenVectors.At(k, j)*d[k];
                            e[k] += eigenVectors.At(k, j)*f;
                        }

                        e[j] = g;
                    }

                    f = 0.0f;

                    for (var j = 0; j < i; j++)
                    {
                        e[j] /= h;
                        f += e[j]*d[j];
                    }

                    var hh = f/(h + h);

                    for (var j = 0; j < i; j++)
                    {
                        e[j] -= hh*d[j];
                    }

                    for (var j = 0; j < i; j++)
                    {
                        f = d[j];
                        g = e[j];

                        for (var k = j; k <= i - 1; k++)
                        {
                            eigenVectors.At(k, j, eigenVectors.At(k, j) - (f*e[k]) - (g*d[k]));
                        }
//.........这里部分代码省略.........

        /// <summary>
        /// Calculate Cholesky step
        /// </summary>
        /// <param name="data">Factor matrix</param>
        /// <param name="rowDim">Number of rows</param>
        /// <param name="firstCol">Column start</param>
        /// <param name="colLimit">Total columns</param>
        /// <param name="multipliers">Multipliers calculated previously</param>
        /// <param name="availableCores">Number of available processors</param>
        private static void DoCholeskyStep(Matrix<Complex32> data, int rowDim, int firstCol, int colLimit, Complex32[] multipliers, int availableCores)
        {
            var tmpColCount = colLimit - firstCol;

            if ((availableCores > 1) && (tmpColCount > 200))
            {
                var tmpSplit = firstCol + (tmpColCount / 3);
                var tmpCores = availableCores / 2;

                CommonParallel.Invoke(
                    () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
                    () => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
            }
            else
            {
                for (var j = firstCol; j < colLimit; j++)
                {
                    var tmpVal = multipliers[j];
                    for (var i = j; i < rowDim; i++)
                    {
                        data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal.Conjugate()));
                    }
                }
            }
        }

        /// <summary>
        /// Multiplies the transpose of this matrix with another matrix and places the results into the result matrix.
        /// </summary>
        /// <param name="other">The matrix to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        protected override void DoTransposeThisAndMultiply(Matrix<Complex> other, Matrix<Complex> result)
        {
            var diagonalOther = other as DiagonalMatrix;
            var diagonalResult = result as DiagonalMatrix;
            if (diagonalOther != null && diagonalResult != null)
            {
                var thisDataCopy = new Complex[diagonalResult._data.Length];
                var otherDataCopy = new Complex[diagonalResult._data.Length];
                Array.Copy(_data, thisDataCopy, (diagonalResult._data.Length > _data.Length) ? _data.Length : diagonalResult._data.Length);
                Array.Copy(diagonalOther._data, otherDataCopy, (diagonalResult._data.Length > diagonalOther._data.Length) ? diagonalOther._data.Length : diagonalResult._data.Length);
                Control.LinearAlgebraProvider.PointWiseMultiplyArrays(thisDataCopy, otherDataCopy, diagonalResult._data);
                return;
            }

            var denseOther = other.Storage as DenseColumnMajorMatrixStorage<Complex>;
            if (denseOther != null)
            {
                var dense = denseOther.Data;
                var diagonal = _data;
                var d = Math.Min(denseOther.RowCount, ColumnCount);
                if (d < ColumnCount)
                {
                    result.ClearSubMatrix(denseOther.RowCount, ColumnCount - denseOther.RowCount, 0, denseOther.ColumnCount);
                }
                int index = 0;
                for (int i = 0; i < denseOther.ColumnCount; i++)
                {
                    for (int j = 0; j < d; j++)
                    {
                        result.At(j, i, dense[index]*diagonal[j]);
                        index++;
                    }
                    index += (denseOther.RowCount - d);
                }
                return;
            }

            base.DoTransposeThisAndMultiply(other, result);
        }

        /// <summary>
        /// Negate each element of this matrix and place the results into the result matrix.
        /// </summary>
        /// <param name="result">The result of the negation.</param>
        protected override void DoNegate(Matrix<Complex> result)
        {
            var diagResult = result as DiagonalMatrix;
            if (diagResult != null)
            {
                Control.LinearAlgebraProvider.ScaleArray(-1, _data, diagResult._data);
                return;
            }

            result.Clear();
            for (var i = 0; i < _data.Length; i++)
            {
                result.At(i, i, -_data[i]);
            }
        }

        /// <summary>
        /// Adds another matrix to this matrix.
        /// </summary>
        /// <param name="other">The matrix to add to this matrix.</param>
        /// <param name="result">The matrix to store the result of the addition.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
        protected override void DoAdd(Matrix<float> other, Matrix<float> result)
        {
            var sparseOther = other as SparseMatrix;
            var sparseResult = result as SparseMatrix;
            if (sparseOther == null || sparseResult == null)
            {
                base.DoAdd(other, result);
                return;
            }

            if (ReferenceEquals(this, other))
            {
                if (!ReferenceEquals(this, result))
                {
                    CopyTo(result);
                }

                Control.LinearAlgebraProvider.ScaleArray(2.0f, _storage.Values, _storage.Values);
                return;
            }

            SparseMatrix left;

            if (ReferenceEquals(sparseOther, sparseResult))
            {
                left = this;
            }
            else if (ReferenceEquals(this, sparseResult))
            {
                left = sparseOther;
            }
            else
            {
                CopyTo(sparseResult);
                left = sparseOther;
            }

            var leftStorage = left._storage;
            for (var i = 0; i < leftStorage.RowCount; i++)
            {
                var endIndex = leftStorage.RowPointers[i + 1];
                for (var j = leftStorage.RowPointers[i]; j < endIndex; j++)
                {
                    var columnIndex = leftStorage.ColumnIndices[j];
                    var resVal = leftStorage.Values[j] + result.At(i, columnIndex);
                    result.At(i, columnIndex, resVal);
                }
            }
        }

        /// <summary>
        /// Puts the upper triangle of this matrix into the result matrix.
        /// </summary>
        /// <param name="result">Where to store the lower triangle.</param>
        private void UpperTriangleImpl(Matrix<float> result)
        {
            var rowPointers = _storage.RowPointers;
            var columnIndices = _storage.ColumnIndices;
            var values = _storage.Values;

            for (var row = 0; row < result.RowCount; row++)
            {
                var endIndex = rowPointers[row + 1];
                for (var j = rowPointers[row]; j < endIndex; j++)
                {
                    if (row <= columnIndices[j])
                    {
                        result.At(row, columnIndices[j], values[j]);
                    }
                }
            }
        }

        /// <summary>
        /// Multiplies this matrix with transpose of another matrix and places the results into the result matrix.
        /// </summary>
        /// <param name="other">The matrix to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        protected override void DoTransposeAndMultiply(Matrix<float> other, Matrix<float> result)
        {
            var otherSparse = other as SparseMatrix;
            var resultSparse = result as SparseMatrix;

            if (otherSparse == null || resultSparse == null)
            {
                base.DoTransposeAndMultiply(other, result);
                return;
            }

            resultSparse.Clear();

            var rowPointers = _storage.RowPointers;
            var values = _storage.Values;

            var otherStorage = otherSparse._storage;

            for (var j = 0; j < RowCount; j++)
            {
                var startIndexOther = otherStorage.RowPointers[j];
                var endIndexOther = otherStorage.RowPointers[j + 1];

                if (startIndexOther == endIndexOther)
                {
                    continue;
                }

                for (var i = 0; i < RowCount; i++)
                {
                    // Multiply row of matrix A on row of matrix B

                    var startIndexThis = rowPointers[i];
                    var endIndexThis = rowPointers[i + 1];

                    if (startIndexThis == endIndexThis)
                    {
                        continue;
                    }

                    var sum = 0f;
                    for (var index = startIndexOther; index < endIndexOther; index++)
                    {
                        var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]);
                        if (ind >= 0)
                        {
                            sum += otherStorage.Values[index]*values[ind];
                        }
                    }

                    resultSparse._storage.At(i, j, sum + result.At(i, j));
                }
            }
        }

        /// <summary>
        /// Subtracts another matrix from this matrix.
        /// </summary>
        /// <param name="other">The matrix to subtract to this matrix.</param>
        /// <param name="result">The matrix to store the result of subtraction.</param>
        /// <exception cref="ArgumentNullException">If the other matrix is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
        protected override void DoSubtract(Matrix<float> other, Matrix<float> result)
        {
            var sparseOther = other as SparseMatrix;
            var sparseResult = result as SparseMatrix;
            if (sparseOther == null || sparseResult == null)
            {
                base.DoSubtract(other, result);
                return;
            }

            if (ReferenceEquals(this, other))
            {
                result.Clear();
                return;
            }

            var otherStorage = sparseOther._storage;

            if (ReferenceEquals(this, sparseResult))
            {
                for (var i = 0; i < otherStorage.RowCount; i++)
                {
                    var endIndex = otherStorage.RowPointers[i + 1];
                    for (var j = otherStorage.RowPointers[i]; j < endIndex; j++)
                    {
                        var columnIndex = otherStorage.ColumnIndices[j];
                        var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j];
                        result.At(i, columnIndex, resVal);
                    }
                }
            }
            else
            {
                if (!ReferenceEquals(sparseOther, sparseResult))
                {
                    sparseOther.CopyTo(sparseResult);
                }

                sparseResult.Negate(sparseResult);

                var rowPointers = _storage.RowPointers;
                var columnIndices = _storage.ColumnIndices;
                var values = _storage.Values;

                for (var i = 0; i < RowCount; i++)
                {
                    var endIndex = rowPointers[i + 1];
                    for (var j = rowPointers[i]; j < endIndex; j++)
                    {
                        var columnIndex = columnIndices[j];
                        var resVal = sparseResult.At(i, columnIndex) + values[j];
                        result.At(i, columnIndex, resVal);
                    }
                }
            }
        }

        /// <summary>
        /// Pointwise multiplies this matrix with another matrix and stores the result into the result matrix.
        /// </summary>
        /// <param name="other">The matrix to pointwise multiply with this one.</param>
        /// <param name="result">The matrix to store the result of the pointwise multiplication.</param>
        protected override void DoPointwiseMultiply(Matrix<float> other, Matrix<float> result)
        {
            result.Clear();

            var rowPointers = _storage.RowPointers;
            var columnIndices = _storage.ColumnIndices;
            var values = _storage.Values;

            for (var i = 0; i < RowCount; i++)
            {
                var endIndex = rowPointers[i + 1];
                for (var j = rowPointers[i]; j < endIndex; j++)
                {
                    var resVal = values[j]*other.At(i, columnIndices[j]);
                    if (resVal != 0f)
                    {
                        result.At(i, columnIndices[j], resVal);
                    }
                }
            }
        }

        /// <summary>
        /// Pointwise divide this matrix by another matrix and stores the result into the result matrix.
        /// </summary>
        /// <param name="divisor">The matrix to pointwise divide this one by.</param>
        /// <param name="result">The matrix to store the result of the pointwise division.</param>
        protected override void DoPointwiseDivide(Matrix<float> divisor, Matrix<float> result)
        {
            result.Clear();

            var rowPointers = _storage.RowPointers;
            var columnIndices = _storage.ColumnIndices;
            var values = _storage.Values;

            for (var i = 0; i < RowCount; i++)
            {
                var endIndex = rowPointers[i + 1];
                for (var j = rowPointers[i]; j < endIndex; j++)
                {
                    if (values[j] != 0f)
                    {
                        result.At(i, columnIndices[j], values[j]/divisor.At(i, columnIndices[j]));
                    }
                }
            }
        }