context: namespace CSparse.Storage
  {
      using CSparse.Properties;
      using System;
      using System.Collections.Generic;
      using System.Numerics;
  
      /// <summary>
      /// Compressed sparse column storage.
      /// </summary>
      /// <typeparam name="T"></typeparam>
      [Serializable]
      public abstract class CompressedColumnStorage<T> : Matrix<T>
          where T : struct, IEquatable<T>, IFormattable
      {
          /// <summary>
          /// Gets or sets a value indicating whether the storage should be
          /// automatically resized to non-zeros count. Defaults to true.
          /// </summary>
          /// <remarks>
          /// Affects only sparse matrix addition and multiplication.
          /// </remarks>
          public static bool AutoTrimStorage { get; set; } = true;
  
          /// <summary>
          /// Column pointers with last entry equal number of non-zeros (size = ColumnCount + 1)
          /// </summary>
          public int[] ColumnPointers;
  
          /// <summary>
          /// Row indices (size >= NonZerosCount)
          /// </summary>
          public int[] RowIndices;
  
          /// <summary>
          /// Numerical values (size >= NonZerosCount)
          /// </summary>
          public T[] Values;
  
          /// <summary>
          /// Gets the number of non-zero entries.
          /// </summary>
          public int NonZerosCount
          {
              get { return ColumnPointers[columns]; }
          }
  
          /// <summary>
          /// Initializes a new instance of the <see cref="CompressedColumnStorage{T}"/> class.
          /// </summary>
          /// <param name="rowCount">The number of rows.</param>
          /// <param name="columnCount">The number of columns.</param>
          /// <remarks>By default, no arrays are allocated. The user will have to assign the storage arrays.</remarks>
          protected CompressedColumnStorage(int rowCount, int columnCount)
              : base(rowCount, columnCount)
          {
              // No array initialization here
          }
  
          /// <summary>
          /// Initializes a new instance of the <see cref="CompressedColumnStorage{T}"/> class.
          /// </summary>
          /// <param name="rowCount">The number of rows.</param>
          /// <param name="columnCount">The number of columns.</param>
          /// <param name="valueCount">The number of non-zero values.</param>
          protected CompressedColumnStorage(int rowCount, int columnCount, int valueCount)
              : base(rowCount, columnCount)
          {
              if (valueCount < 0)
              {
                  throw new ArgumentOutOfRangeException(nameof(valueCount), Resources.ValueNonNegative);
              }
  
              this.ColumnPointers = new int[columnCount + 1];
              this.RowIndices = new int[valueCount];
              this.Values = new T[valueCount];
          }
  
          /// <summary>
          /// Initializes a new instance of the <see cref="CompressedColumnStorage{T}"/> class. Based on other CCS arrays
          /// </summary>
          /// <param name="rowCount">The number of rows.</param>
          /// <param name="columnCount">The number of columns.</param>
          /// <param name="values"></param>
          /// <param name="rowIndices"></param>
          /// <param name="columnPointers"></param>
          /// <remarks>The provided arrays will be re-used (not cloned).</remarks>
          public CompressedColumnStorage(int rowCount, int columnCount, T[] values, int[] rowIndices, int[] columnPointers)
              : base(rowCount, columnCount)
          {
              if (columnPointers.Length != columnCount + 1)
              {
                  throw new ArgumentException("columnPointers.Length must equal (columnCount + 1)");
              }
  
              if (values.Length != rowIndices.Length)
              {
                  throw new ArgumentException("rowIndices.Length must equal values.Length");
              }
  
              this.ColumnPointers = columnPointers;
              this.RowIndices = rowIndices;
              this.Values = values;
          }
  
          #region Public static functions
  
          /// <summary>
          /// Create a new sparse matrix as a copy of the given other matrix.
          /// </summary>
          public static CompressedColumnStorage<T> OfMatrix(Matrix<T> matrix)
          {
              var c = Converter.FromEnumerable_<T>(matrix.EnumerateIndexed(), matrix.RowCount, matrix.ColumnCount);
  
              return Converter.ToCompressedColumnStorage_(c);
          }
  
          /// <summary>
          /// Create a new sparse matrix as a copy of the given two-dimensional array.
          /// </summary>
          public static CompressedColumnStorage<T> OfArray(T[,] array)
          {
              var c = Converter.FromDenseArray_(array);
  
              return Converter.ToCompressedColumnStorage_(c);
          }
  
          /// <summary>
          /// Create a new sparse matrix as a copy of the given two-dimensional array.
          /// </summary>
          public static CompressedColumnStorage<T> OfJaggedArray(T[][] array)
          {
              var c = Converter.FromJaggedArray_(array);
  
              return Converter.ToCompressedColumnStorage_(c);
          }
  
          /// <summary>
          /// Create a new sparse matrix as a copy of the given coordinate storage.
          /// </summary>
          public static CompressedColumnStorage<T> OfIndexed(CoordinateStorage<T> coordinateStorage, bool inplace = false)
          {
              return Converter.ToCompressedColumnStorage_(coordinateStorage, true, inplace);
          }
  
          /// <summary>
          /// Create a new sparse matrix as a copy of the given indexed enumerable.
          /// </summary>
          public static CompressedColumnStorage<T> OfIndexed(int rows, int columns, IEnumerable<Tuple<int, int, T>> enumerable)
          {
              var c = Converter.FromEnumerable_<T>(enumerable, rows, columns);
  
              return Converter.ToCompressedColumnStorage_(c);
          }
  
          /// <summary>
          /// Create a new sparse matrix as a copy of the given array (row-major).
          /// </summary>
          public static CompressedColumnStorage<T> OfRowMajor(int rows, int columns, T[] rowMajor)
          {
              var c = Converter.FromRowMajorArray_<T>(rowMajor, rows, columns);
  
              return Converter.ToCompressedColumnStorage_(c);
          }
  
          /// <summary>
          /// Create a new sparse matrix as a copy of the given array (column-major).
          /// </summary>
          public static CompressedColumnStorage<T> OfColumnMajor(int rows, int columns, T[] columnMajor)
          {
              var c = Converter.FromColumnMajorArray_<T>(columnMajor, rows, columns);
  
              return Converter.ToCompressedColumnStorage_(c);
          }
  
          /// <summary>
          /// Create a new square sparse matrix with the diagonal as a copy of the given array.
          /// </summary>
          public static CompressedColumnStorage<T> OfDiagonalArray(T[] diagonal)
          {
              int order = diagonal.Length;
  
              var A = Create(order, order, order);
  
              var ap = A.ColumnPointers;
              var ai = A.RowIndices;
              var ax = A.Values;
  
              for (int i = 0; i < order; i++)
              {
                  ap[i] = i;
                  ai[i] = i;
                  ax[i] = diagonal[i];
              }
  
              ap[order] = order;
  
              return A;
          }
  
          /// <summary>
          /// Create a new square sparse matrix and initialize each diagonal value to the same provided value.
          /// </summary>
          public static CompressedColumnStorage<T> CreateDiagonal(int order, T value)
          {
              var A = Create(order, order, order);
  
              var ap = A.ColumnPointers;
              var ai = A.RowIndices;
              var ax = A.Values;
  
              for (int i = 0; i < order; i++)
              {
                  ap[i] = i;
                  ai[i] = i;
                  ax[i] = value;
              }
  
              ap[order] = order;
  
              return A;
          }
  
          /// <summary>
          /// Create a new square sparse identity matrix where each diagonal value is set to One.
          /// </summary>
          public static CompressedColumnStorage<T> CreateIdentity(int order)
          {
              return CreateDiagonal(order, One);
          }
  
          #endregion
  
          /// <inheritdoc />
          public override T At(int row, int column)
          {
              int index = ColumnPointers[column];
              int length = ColumnPointers[column + 1] - index;
              int pos = Array.BinarySearch(RowIndices, index, length, row);
              return pos >= 0 ? Values[pos] : Zero;
          }
  
          /// <inheritdoc />
          public override void Clear()
          {
              Array.Clear(ColumnPointers, 0, ColumnPointers.Length);
              Array.Clear(Values, 0, Values.Length);
          }
  
          /// <inheritdoc />
          public override T[] Row(int rowIndex)
          {
              var target = new T[ColumnCount];
  
              Row(rowIndex, target);
  
              return target;
          }
  
          /// <inheritdoc />
          public override void Row(int rowIndex, T[] target)
          {
              if (target.Length != ColumnCount)
              {
                  throw new Exception();
              }
  
              int i, n = ColumnCount;
  
              var ap = ColumnPointers;
              var ai = RowIndices;
              var ax = Values;
  
              for (int k = 0; k < n; k++)
              {
                  // Check if columns contain row index.
                  i = Array.BinarySearch(ai, ap[k], ap[k + 1] - ap[k], rowIndex);
  
                  if (i >= 0)
                  {
                      target[k] = ax[i];
                  }
              }
          }
  
          /// <inheritdoc />
          public override T[] Column(int columnIndex)
          {
              var target = new T[RowCount];
  
              Column(columnIndex, target);
  
              return target;
          }
  
          /// <inheritdoc />
          public override void Column(int columnIndex, T[] target)
          {
              if (target.Length != RowCount)
              {
                  throw new Exception();
              }
  
              var ap = ColumnPointers;
              var ai = RowIndices;
              var ax = Values;
  
              int colEnd = ap[columnIndex + 1];
  
              for (int k = ap[columnIndex]; k < colEnd; k++)
              {
                  target[ai[k]] = ax[k];
              }
          }
  
          #region Linear Algebra (Matrix)
  
          /// <summary>
          /// Returns the transpose of this matrix.
          /// </summary>
          public CompressedColumnStorage<T> Transpose()
          {
              return this.Transpose(false);
          }
  
          /// <summary>
          /// Transpose this matrix and store the result in given matrix.
          /// </summary>
          /// <param name="result">Storage for the tranposed matrix.</param>
          public void Transpose(CompressedColumnStorage<T> result)
          {
              this.Transpose(result, false);
          }
  
          /// <summary>
          /// Returns the transpose of this matrix.
          /// </summary>
          /// <param name="storage">A value indicating, whether the transpose should be done on storage level (without complex conjugation).</param>
          public CompressedColumnStorage<T> Transpose(bool storage)
          {
              var result = CompressedColumnStorage<T>.Create(columns, rows, this.NonZerosCount);
              this.Transpose(result, storage);
              return result;
          }
  
          /// <summary>
          /// Transpose this matrix and store the result in given matrix.
          /// </summary>
          /// <param name="result">Storage for the tranposed matrix.</param>
          /// <param name="storage">A value indicating, whether the transpose should be done on storage level (without complex conjugation).</param>
          public virtual void Transpose(CompressedColumnStorage<T> result, bool storage)
          {
              int i, j, p;
  
              var cx = result.Values;
              var cp = result.ColumnPointers;
              var ci = result.RowIndices;
  
              int[] w = new int[rows];
  
              for (p = 0; p < ColumnPointers[columns]; p++)
              {
                  // Row counts.
                  w[RowIndices[p]]++;
              }
  
              // Row pointers.
              Helper.CumulativeSum(cp, w, rows);
  
  
gt:             for (i = 0; i < columns; i++)