SymbolicColumnStorage

CSparse.Tests/CSparse.Tests.csproj

@@ -41,9 +41,9 @@
   </ItemGroup>
 
   <ItemGroup>
-    <PackageReference Include="NUnit" Version="4.1.0" />
-    <PackageReference Include="NUnit3TestAdapter" Version="4.5.0" />
-    <PackageReference Include="Microsoft.NET.Test.Sdk" Version="17.9.0" />
+    <PackageReference Include="NUnit" Version="4.2.2" />
+    <PackageReference Include="NUnit3TestAdapter" Version="4.6.0" />
+    <PackageReference Include="Microsoft.NET.Test.Sdk" Version="17.11.1" />
   </ItemGroup>
 
   <ItemGroup>

CSparse.Tests/Ordering/TestAMD.cs

@@ -0,0 +1,89 @@
+using CSparse.Ordering;
+using CSparse.Storage;
+using NUnit.Framework;
+
+namespace CSparse.Tests.Ordering
+{
+    class TestAMD
+    {
+        [Test]
+        public void TestAMD1()
+        {
+            int[] ap = [0, 9, 15, 21, 27, 33, 39, 48, 57, 61, 70, 76, 82, 88, 94, 100, 106, 110, 119, 128, 137, 143, 152, 156, 160];
+            int[] ai = [
+            /* col  0 */  0, 5, 6, 12, 13, 17, 18, 19, 21,
+            /* col  1 */  1, 8, 9, 13, 14, 17,
+            /* col  2 */  2, 6, 11, 20, 21, 22,
+            /* col  3 */  3, 7, 10, 15, 18, 19,
+            /* col  4 */  4, 7, 9, 14, 15, 16,
+            /* col  5 */  0, 5, 6, 12, 13, 17,
+            /* col  6 */  0, 2, 5, 6, 11, 12, 19, 21, 23,
+            /* col  7 */  3, 4, 7, 9, 14, 15, 16, 17, 18,
+            /* col  8 */  1, 8, 9, 14,
+            /* col  9 */  1, 4, 7, 8, 9, 13, 14, 17, 18,
+            /* col 10 */  3, 10, 18, 19, 20, 21,
+            /* col 11 */  2, 6, 11, 12, 21, 23,
+            /* col 12 */  0, 5, 6, 11, 12, 23,
+            /* col 13 */  0, 1, 5, 9, 13, 17,
+            /* col 14 */  1, 4, 7, 8, 9, 14,
+            /* col 15 */  3, 4, 7, 15, 16, 18,
+            /* col 16 */  4, 7, 15, 16,
+            /* col 17 */  0, 1, 5, 7, 9, 13, 17, 18, 19,
+            /* col 18 */  0, 3, 7, 9, 10, 15, 17, 18, 19,
+            /* col 19 */  0, 3, 6, 10, 17, 18, 19, 20, 21,
+            /* col 20 */  2, 10, 19, 20, 21, 22,
+            /* col 21 */  0, 2, 6, 10, 11, 19, 20, 21, 22,
+            /* col 22 */  2, 20, 21, 22,
+            /* col 23 */  6, 11, 12, 23 ];
+
+            var S = new SymbolicColumnStorage(24, 24, ap, ai, false);
+
+            var p = AMD.Generate(S, ColumnOrdering.MinimumDegreeAtA);
+
+            int[] expected = [8, 16, 4, 14, 15, 1, 7, 9, 22, 23, 2, 11, 3, 12, 13, 17, 18, 0, 10, 19, 20, 6, 21, 5, 24];
+
+            Assert.That(Permutation.IsValid(p), Is.True);
+            Assert.That(p, Is.EqualTo(expected).AsCollection);
+        }
+
+        [Test]
+        public void TestAMD2()
+        {
+            int[] ap = [0, 9, 14, 20, 28, 33, 37, 44, 53, 58, 63, 63, 66, 69, 72, 75, 78, 82, 86, 91, 97, 101, 112, 112, 116];
+            int[] ai = [
+            /* col  0 */  0, 17, 18, 21, 5, 12, 5, 0, 13,
+            /* col  1 */  14, 1, 8, 13, 17,
+            /* col  2 */  2, 20, 11, 6, 11, 22,
+            /* col  3 */  3, 3, 10, 7, 18, 18, 15, 19,
+            /* col  4 */  7, 9, 15, 14, 16,
+            /* col  5 */  5, 13, 6, 17,
+            /* col  6 */  5, 0, 11, 6, 12, 6, 23,
+            /* col  7 */  3, 4, 9, 7, 14, 16, 15, 17, 18,
+            /* col  8 */  1, 9, 14, 14, 14,
+            /* col  9 */  7, 13, 8, 1, 17,
+            /* col 10 */
+            /* col 11 */  2, 12, 23,
+            /* col 12 */  5, 11, 12,
+            /* col 13 */  0, 13, 17,
+            /* col 14 */  1, 9, 14,
+            /* col 15 */  3, 15, 16,
+            /* col 16 */  16, 4, 4, 15,
+            /* col 17 */  13, 17, 19, 17,
+            /* col 18 */  15, 17, 19, 9, 10,
+            /* col 19 */  17, 19, 20, 0, 6, 10,
+            /* col 20 */  22, 10, 20, 21,
+            /* col 21 */  6, 2, 10, 19, 20, 11, 21, 22, 22, 22, 22,
+            /* col 22 */
+            /* col 23 */  12, 11, 12, 23 ];
+
+            var S = new SymbolicColumnStorage(24, 24, ap, ai, false);
+
+            var p = AMD.Generate(S, ColumnOrdering.MinimumDegreeAtA);
+
+            int[] expected = [10, 11, 23, 12, 2, 6, 8, 14, 15, 16, 4, 1, 9, 7, 18, 3, 5, 17, 0, 19, 20, 21, 13, 22, 24];
+
+            Assert.That(Permutation.IsValid(p), Is.True);
+            Assert.That(p, Is.EqualTo(expected).AsCollection);
+        }
+    }
+}

CSparse.Tests/Ordering/TestDulmageMendelsohn.cs

@@ -1,6 +1,8 @@
 namespace CSparse.Tests.Ordering
 {
+    using CSparse.Double;
     using CSparse.Ordering;
+    using CSparse.Storage;
     using NUnit.Framework;
     using System;
 
@@ -33,5 +35,46 @@ public void TestGenerate2()
 
             Assert.That(dm.StructuralRank == n, Is.True);
         }
+
+        [Test]
+        public void TestGenerate3()
+        {
+            var A = SparseMatrix.OfRowMajor(8, 8,
+            [
+                11, 12,  0,  0,  0,  0,  0,  0,
+                 0, 22, 23,  0, 25, 26,  0,  0,
+                 0,  0, 33, 34,  0,  0, 37,  0,
+                 0,  0, 43, 44,  0,  0,  0, 48,
+                51,  0,  0,  0, 55, 56,  0,  0,
+                 0,  0,  0,  0,  0, 66, 67,  0,
+                 0,  0,  0,  0,  0, 76, 77,  0,
+                 0,  0,  0, 84,  0,  0, 87, 88
+            ]);
+
+            var S = SymbolicColumnStorage.Create(A);
+
+            var dm = DulmageMendelsohn.Generate(S, 1);
+
+            Assert.That(dm.StructuralRank, Is.EqualTo(8));
+            Assert.That(dm.Blocks, Is.EqualTo(3));
+
+            int[] expected = [0, 1, 4, 2, 3, 7, 5, 6];
+
+            Assert.That(dm.RowPermutation, Is.EqualTo(expected).AsCollection);
+            Assert.That(dm.ColumnPermutation, Is.EqualTo(expected).AsCollection);
+
+            expected = [0, 3, 6, 8];
+
+            Assert.That(dm.BlockRowPointers, Is.EqualTo(expected).AsCollection);
+            Assert.That(dm.BlockColumnPointers, Is.EqualTo(expected).AsCollection);
+
+            expected = [0, 0, 8, 8, 8];
+
+            Assert.That(dm.CoarseRowDecomposition, Is.EqualTo(expected).AsCollection);
+
+            expected = [0, 0, 0, 8, 8];
+
+            Assert.That(dm.CoarseColumnDecomposition, Is.EqualTo(expected).AsCollection);
+        }
     }
 }

CSparse/CSparse.csproj

@@ -11,17 +11,21 @@
     <Company />
     <Copyright>Copyright Christian Woltering © 2012-2024</Copyright>
     <Authors>Christian Woltering</Authors>
-    <AssemblyVersion>4.1.0.0</AssemblyVersion>
-    <FileVersion>4.1.0.0</FileVersion>
+    <AssemblyVersion>4.2.0.0</AssemblyVersion>
+    <FileVersion>4.2.0.0</FileVersion>
     <PackageTags>math sparse matrix lu cholesky qr decomposition factorization </PackageTags>
-    <Version>4.1.0</Version>
+    <Version>4.2.0</Version>
     <AssemblyName>CSparse</AssemblyName>
     <RootNamespace>CSparse</RootNamespace>
     <PackageLicenseExpression>LGPL-2.1-only</PackageLicenseExpression>
     <PackageProjectUrl>https://github.com/wo80/CSparse.NET</PackageProjectUrl>
     <RepositoryUrl>https://github.com/wo80/CSparse.NET.git</RepositoryUrl>
     <RepositoryType>git</RepositoryType>
     <PackageReleaseNotes>
+		Version 4.2.0
+
+		* Make SymbolicColumnStorage class public and update StronglyConnectedComponents and DulmageMendelsohn decomposition accordingly.
+
 		Version 4.1.0
 
 		* Add overload for creating a sparse matrix from an enumerable of ValueTuple.

CSparse/Ordering/AMD.cs

@@ -25,6 +25,22 @@ public static class AMD
         /// </remarks>
         public static int[] Generate<T>(CompressedColumnStorage<T> A, ColumnOrdering order)
             where T : struct, IEquatable<T>, IFormattable
+        {
+            return Generate(SymbolicColumnStorage.Create(A), order);
+        }
+
+        /// <summary>
+        /// Generate minimum degree ordering of A+A' (if A is symmetric) or A'A.
+        /// </summary>
+        /// <param name="A">Column-compressed matrix</param>
+        /// <param name="order">Column ordering method</param>
+        /// <returns>amd(A+A') if A is symmetric, or amd(A'A) otherwise, null on 
+        /// error or for natural ordering</returns>
+        /// <remarks>
+        /// See Chapter 7.1 (Fill-reducing orderings: Minimum degree ordering) in 
+        /// "Direct Methods for Sparse Linear Systems" by Tim Davis.
+        /// </remarks>
+        public static int[] Generate(SymbolicColumnStorage A, ColumnOrdering order)
         {
             int[] Cp, Ci, P, W, nv, next, head, elen, degree, w, hhead;
 
@@ -42,7 +58,7 @@ public static int[] Generate<T>(CompressedColumnStorage<T> A, ColumnOrdering ord
                 return Permutation.Create(n);
             }
 
-            var C = ConstructMatrix(SymbolicColumnStorage.Create(A), order);
+            var C = ConstructMatrix(A, order);
 
             Cp = C.ColumnPointers;
             cnz = Cp[n];
@@ -139,7 +155,7 @@ public static int[] Generate<T>(CompressedColumnStorage<T> A, ColumnOrdering ord
                             Ci[p] = -(j + 2); // first entry is now CS_FLIP(j)
                         }
                     }
-                    for (q = 0, p = 0; p < cnz; ) // scan all of memory
+                    for (q = 0, p = 0; p < cnz;) // scan all of memory
                     {
                         if ((j = FLIP(Ci[p++])) >= 0) // found object j
                         {
@@ -303,7 +319,7 @@ public static int[] Generate<T>(CompressedColumnStorage<T> A, ColumnOrdering ord
                         eln = elen[i];
                         for (p = Cp[i] + 1; p <= Cp[i] + ln - 1; p++) w[Ci[p]] = mark;
                         jlast = i;
-                        for (j = next[i]; j != -1; ) // compare i with all j
+                        for (j = next[i]; j != -1;) // compare i with all j
                         {
                             ok = (W[j] == ln) && (elen[j] == eln);
                             for (p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++)
@@ -411,13 +427,13 @@ private static bool KeepOffDiag(int i, int j)
         private static SymbolicColumnStorage ConstructMatrix(SymbolicColumnStorage A, ColumnOrdering order)
         {
             SymbolicColumnStorage result = null;
-            
+
             // Compute A'
             var AT = A.Transpose();
 
             int m = A.RowCount;
             int n = A.ColumnCount;
-            
+
             if (order == ColumnOrdering.MinimumDegreeAtPlusA)
             {
                 if (n != m)
@@ -444,7 +460,7 @@ private static SymbolicColumnStorage ConstructMatrix(SymbolicColumnStorage A, Co
                 {
                     // Column j of AT starts here.
                     p = colptr[j];
-                    
+
                     // New column j starts here.
                     colptr[j] = p2;
 

CSparse/Ordering/DulmageMendelsohn.cs

@@ -22,7 +22,7 @@ public sealed class DulmageMendelsohn
         private int nb;    // number of blocks in fine dmperm decomposition
 
         /// <summary>
-        /// Create a new Decomposition instance. 
+        /// Create a new decomposition instance. 
         /// </summary>
         private DulmageMendelsohn(int m, int n)
         {
@@ -97,22 +97,29 @@ public int Singletons
         public int[] CoarseColumnDecomposition => cc;
 
         /// <summary>
-        /// Compute coarse and then fine Dulmage-Mendelsohn decomposition. seed
-        /// optionally selects a randomized algorithm.
+        /// Compute coarse and fine Dulmage-Mendelsohn decomposition.
         /// </summary>
         /// <param name="matrix">column-compressed matrix</param>
-        /// <param name="seed">0: natural, -1: reverse, random order otherwise</param>
+        /// <param name="seed"> The seed optionally selects a randomized algorithm (0 = default (natural), -1 = reverse, random order otherwise).</param>
         /// <returns>Dulmage-Mendelsohn analysis</returns>
         public static DulmageMendelsohn Generate<T>(CompressedColumnStorage<T> matrix, int seed = 0)
              where T : struct, IEquatable<T>, IFormattable
+        {
+            return Generate(SymbolicColumnStorage.Create(matrix), seed);
+        }
+
+        /// <summary>
+        /// Compute coarse and fine Dulmage-Mendelsohn decomposition.
+        /// </summary>
+        /// <param name="A">The matrix represented by <see cref="SymbolicColumnStorage"/>.</param>
+        /// <param name="seed"> The seed optionally selects a randomized algorithm (0 = default (natural), -1 = reverse, random order otherwise).</param>
+        /// <returns>Dulmage-Mendelsohn analysis</returns>
+        public static DulmageMendelsohn Generate(SymbolicColumnStorage A, int seed = 0)
         {
             int i, j, k, cnz, nc, nb1, nb2;
             int[] Cp, ps, rs;
             bool ok;
 
-            // We are not interested in the actual matrix values.
-            var A = SymbolicColumnStorage.Create(matrix);
-
             // Maximum matching
             int m = A.RowCount;
             int n = A.ColumnCount;
@@ -147,7 +154,7 @@ public static DulmageMendelsohn Generate<T>(CompressedColumnStorage<T> matrix, i
             // Fine decomposition
             int[] pinv = Permutation.Invert(p); // pinv=p'
 
-            var C = SymbolicColumnStorage.Create(matrix);
+            var C = A.Clone();
             A.Permute(pinv, q, C); // C=A(p,q) (it will hold A(R2,C2))
 
             Cp = C.ColumnPointers;

CSparse/Ordering/StronglyConnectedComponents.cs

@@ -48,19 +48,31 @@ private StronglyConnectedComponents(int m, int n)
         public static StronglyConnectedComponents Generate<T>(CompressedColumnStorage<T> matrix)
              where T : struct, IEquatable<T>, IFormattable
         {
-            return Generate(SymbolicColumnStorage.Create(matrix, false), matrix.ColumnCount);
+            return Generate(SymbolicColumnStorage.Create(matrix, false));
         }
-        
+
         /// <summary>
-        /// Find strongly connected components of A.
+        /// Compute strongly connected components of A.
         /// </summary>
-        /// <param name="A"></param>
-        /// <param name="n"></param>
-        /// <returns></returns>
-        internal static StronglyConnectedComponents Generate(SymbolicColumnStorage A, int n)
+        /// <param name="A">The matrix represented by <see cref="SymbolicColumnStorage"/>.</param>
+        /// <returns>Strongly connected components</returns>
+        public static StronglyConnectedComponents Generate(SymbolicColumnStorage A)
+        {
+            return Generate(A, A.ColumnCount);
+        }
+
+        /// <summary>
+        /// Compute strongly connected components of A.
+        /// </summary>
+        /// <param name="A">The matrix represented by <see cref="SymbolicColumnStorage"/>.</param>
+        /// <param name="size">The size of the matrix.</param>
+        /// <returns>Strongly connected components</returns>
+        public static StronglyConnectedComponents Generate(SymbolicColumnStorage A, int size)
         {
             // matrix A temporarily modified, then restored
 
+            int n = size;
+
             int i, k, b, nb = 0, top;
             int[] xi, p, r, Ap, ATp;
 

CSparse/Storage/CompressedColumnStorage.cs

@@ -40,10 +40,7 @@ public abstract class CompressedColumnStorage<T> : Matrix<T>
         /// <summary>
         /// Gets the number of non-zero entries.
         /// </summary>
-        public int NonZerosCount
-        {
-            get { return ColumnPointers[columns]; }
-        }
+        public int NonZerosCount => ColumnPointers[columns];
 
         /// <summary>
         /// Initializes a new instance of the <see cref="CompressedColumnStorage{T}"/> class.