Can I construct sparse matrix from multiple diagonals? about csparse.net HOT 4 CLOSED

/// <summary>
/// Create a sparse matrix from given diagonals.
/// </summary>
/// <typeparam name="T"></typeparam>
/// <param name="A">The input diagonals stored column-wise.</param>
/// <param name="diags">The diagonal offsets.</param>
/// <param name="rowCount">The target matrix row count.</param>
/// <param name="columnCount">The target matrix column count.</param>
/// <returns></returns>
/// <exception cref="ArgumentException"></exception>
static CompressedColumnStorage<T> spdiags<T>(DenseColumnMajorStorage<T> A, int[] diags,
    int rowCount, int columnCount) where T : struct, IEquatable<T>, IFormattable
{
    int k = diags.Length;

    if (A.ColumnCount != k)
    {
        throw new ArgumentException("Columns of A must correspond to diagonals.");
    }

    int dim = Math.Min(rowCount, columnCount);

    var coo = new CoordinateStorage<T>(rowCount, columnCount, k * dim);

    // Stores current column of A.
    var column = new T[A.RowCount];

    for (int i = 0; i < k; i++)
    {
        A.Column(i, column);

        int d = diags[i];

        int ik = d < 0 ? -d : 0;
        int jk = d > 0 ?  d : 0;

        int len = dim - Math.Abs(d);

        int offset = d > 0 ? d : 0;

        for (int j = 0; j < len; j++)
        {
            coo.At(ik++, jk++, column[offset + j]);
        }
    }

    return CompressedColumnStorage<T>.OfIndexed(coo, true);
}

var A = DenseMatrix.OfColumnMajor(4, 3, new double[]
{
    -1, -1, -1, -1,
     2,  2,  2,  2,
    -1, -1, -1, -1
});

var B = spdiags(A, new int[] { -1, 0, 1 }, 4, 4);

public static CompressedColumnStorage<T> OfDiagonals(DenseColumnMajorStorage<T> A, int[] diags,
    int rowCount, int columnCount)
{
    int k = diags.Length;

    if (A.ColumnCount != k)
    {
        throw new ArgumentException("Columns of A must correspond to diagonals.");
    }

    // Upper limit for storage size.
    int size = k * Math.Min(rowCount, columnCount);

    // Current non-zeros count.
    int nz = 0;

    var result = Create(rowCount, columnCount, size);

    var ap = result.ColumnPointers;
    var ai = result.RowIndices;
    var ax = result.Values;

    // Fill each column of the result matrix.
    for (int col = 0; col < columnCount; col++)
    {
        ap[col] = nz;

        // Add diagonals at specified offsets.
        for (int j = 0; j < k; j++)
        {
            int row = col - diags[j];

            if (row >= 0 && row < rowCount)
            {
                ai[nz] = row;
                ax[nz] = A.At(col, j);

                nz++;
            }
        }
    }

    ap[columnCount] = nz;

    Helper.SortIndices(result);

    return result;
}

// Copied from CompressedColumnStorage.cs
static CompressedColumnStorage<T> Create(int rowCount, int columnCount, int valueCount)
{
    if (typeof(T) == typeof(double))
    {
        return new CSparse.Double.SparseMatrix(rowCount, columnCount, valueCount)
            as CompressedColumnStorage<T>;
    }

    if (typeof(T) == typeof(Complex))
    {
        return new CSparse.Complex.SparseMatrix(rowCount, columnCount, valueCount)
            as CompressedColumnStorage<T>;
    }

    throw new NotSupportedException();
}