mathnet-numerics

mathnet-numerics copied to clipboard

Since some iterative solvers may require an operation y <- A^t * x, SparseMatrix should provide a special TransposeMultiply method (without actually having to compute/allocate the transpose). Here's what I do:

/// <summary>
/// Multiplies a (m-by-n) matrix by a vector, y = alpha*A*x + beta*y. 
/// </summary>
/// <param name="x">Vector of length n (column count).</param>
/// <param name="y">Vector of length m (row count), containing the result.</param>
/// <param name="alpha">Scalar to multiply with matrix.</param>
/// <param name="beta">Scalar to multiply with vector y.</param>
/// <remarks>
/// Input values of vector y will be accumulated.
/// </remarks>
public void Multiply(double alpha, double[] x, double beta, double[] y)
{
    var ax = this.Values;
    var ap = this.RowPointers;
    var ai = this.ColumnIndices;

    double yi;

    int end, start = ap[0];

    for (int i = 0; i < nrows; i++)
    {
        end = ap[i + 1];

        yi = beta * y[i];
        for (int k = start; k < end; k++)
        {
            yi += alpha * ax[k] * x[ai[k]];
        }
        y[i] = yi;

        start = end;
    }
}

/// <summary>
/// Multiplies the transpose of a (m-by-n) matrix by a vector, y = alpha*A'*x + beta*y. 
/// </summary>
/// <param name="x">Vector of length m (column count of A').</param>
/// <param name="y">Vector of length n (row count of A'), containing the result.</param>
/// <param name="alpha">Scalar to multiply with matrix.</param>
/// <param name="beta">Scalar to multiply with vector y.</param>
/// <remarks>
/// Input values of vector y will be accumulated.
/// </remarks>
public void TransposeMultiply(double alpha, double[] x, double beta, double[] y)
{
    var ax = this.Values;
    var ap = this.RowPointers;
    var ai = this.ColumnIndices;

    // Scale y by beta
    for (int j = 0; j < ncols; j++)
    {
        y[j] = beta * y[j];
    }

    int end;
    double xi;

    for (int i = 0; i < nrows; i++)
    {
        xi = alpha * x[i];

        end = ap[i + 1];

        for (int k = ap[i]; k < end; k++)
        {
            y[ai[k]] += ax[k] * xi;
        }
    }
}

For dense matrices this would directly map to GEMM?

Yes (GEMV). This is the most general form, so you may decide if you want to use this one or other variants (without scaling factors, without update etc.).

Related: https://mathnetnumerics.codeplex.com/discussions/543114

Got fixed somewhere along the way...