CreateQFromSparseQR.cs

// See https://github.com/wo80/CSparse.NET/wiki/Factorization-members

namespace ConsoleApp
{
    using CSparse;
    using CSparse.Double;
    using CSparse.Double.Factorization;
    using System;
    using System.Diagnostics;
    using System.Globalization;
    using System.Threading;

    static class Test
    {
        public static void Run()
        {
            Thread.CurrentThread.CurrentCulture = CultureInfo.InvariantCulture;

            var A = (SparseMatrix)SparseMatrix.OfArray(new double[,]
            {
                { 7, 3, 1, 0 },
                { 0, 2, 0, 2 },
                { 0, 0, 1, 0 },
                { 1, 0, 0, 2 }
                //{ 0, 1, 1, 1 }
            });

            Run(A, true);
        }

        static void Run(SparseMatrix A, bool print = false)
        {
            // For underdetermined systems, SparseQR factorizes the transposed matrix
            // and the following code won't work.

            if (A.RowCount < A.ColumnCount) throw new Exception("Underdetermined system not supported.");

            var timer = Stopwatch.StartNew();

            PrintMatrix("A:", A);

            var qr = SparseQR.Create(A, ColumnOrdering.MinimumDegreeAtA);

            var t1 = timer.Elapsed;

            qr.GetFactors(out SparseMatrix H, out SparseMatrix R, out double[] beta, out int[] p, out int[] q);

            timer.Restart();

            var Q = CreateQ(H, beta);

            var t2 = timer.Elapsed;

            // Uncomment for benchmark.
            //Console.Write("N = {0,5}: A {1,4:0.0}% , H {2,4:0.0}% , R {3,4:0.0}% , Q {4,4:0.0}%, ", A.RowCount, D(A), D(H), D(R), D(Q));
            //Console.Write("fac time: {0,6:0}ms, Q time: {1,6:0}ms", t1.TotalMilliseconds, t2.TotalMilliseconds);
            //Console.WriteLine();

            if (print)
            {
                var result = Q.Multiply(R.PermuteColumns(Permutation.Invert(p)));

                result.DropZeros(1e-12);
                result.PermuteRows(Permutation.Invert(q));

                PrintMatrix("Q*R:", (SparseMatrix)result);
            }
        }

        static double D(SparseMatrix A)
        {
            return A.NonZerosCount * 100.0 / (A.RowCount * A.ColumnCount);
        }

        static SparseMatrix CreateQ(SparseMatrix H, double[] beta)
        {
            int rows = H.RowCount;
            int cols = H.ColumnCount;

            var hp = H.ColumnPointers;

            int begin = hp[0], end, nnz = 0;

            // Get max nnz of Householder columns.
            for (int i = 0; i < cols; i++)
            {
                end = hp[i + 1];
                nnz = Math.Max(end - begin, nnz);
                begin = end;
            }

            // Dyadic product will have at most nnz^2 entries.
            nnz = nnz * nnz;

            if (!CheckStorageRequirements(rows, cols, nnz))
            {
                Console.WriteLine("Exit.");
            }

            // DyadicProduct assumes row indices to be sorted.
            Helper.SortIndices(H);

            // Allocate all storage outside the loop.
            var E = SparseMatrix.CreateIdentity(rows);
            var D = new SparseMatrix(rows, rows, nnz);
            var P = new SparseMatrix(rows, rows, nnz);
            var Q = (SparseMatrix)E.Clone();

            var Qtemp = new SparseMatrix(rows, rows, nnz);

            // Disable auto-trim for matrix addition and multiplication.
            SparseMatrix.AutoTrimStorage = false;

            for (int i = 0; i < cols; i++)
            {
                DyadicProduct(i, H, D);

                nnz = E.NonZerosCount + D.NonZerosCount;

                if (nnz > Qtemp.Values.Length)
                {
                    Array.Resize(ref Qtemp.RowIndices, 2 * nnz);
                    Array.Resize(ref Qtemp.Values, 2 * nnz);
                }

                E.Add(1.0, -beta[i], D, Qtemp);

                Q.Multiply(Qtemp, P);

                P.CopyTo(Q);
            }

            return Q;

        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="i">The current column index.</param>
        /// <param name="H">The Householder matrix.</param>
        /// <param name="Q">The Q matrix of the QR factorization.</param>
        /// <remarks>
        /// IMPORTANT:
        /// - the row indices of the Householder matrix <paramref name="H"/> have to be sorted.
        /// - the output matrix <paramref name="Q"/> has to provide enough space to store the non-zero entries.
        /// </remarks>
        static void DyadicProduct(int i, SparseMatrix H, SparseMatrix Q)
        {
            const double EPS = 1e-12;

            // The Householder matrix
            var hp = H.ColumnPointers;
            var hi = H.RowIndices;
            var hx = H.Values;

            // The Q matrix for current column of H.
            var qp = Q.ColumnPointers;
            var qi = Q.RowIndices;
            var qx = Q.Values;

            int columns = Q.ColumnCount;

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

            // Reset the column pointers.
            for (int j = 0; j < columns; j++)
            {
                qp[j] = 0;
            }

            int start = hp[i];
            int end = hp[i + 1];

            // Let v be the i-th column of the Householder matrix (lower triangle).
            // Set the bottom-right submatrix of Q to v*v'.

            // Traverse v'.
            for (int j = start; j < end; j++)
            {
                // Current value of v'.
                double a = hx[j];

                if (Math.Abs(a - EPS) > 0)
                {
                    // Traverse v.
                    for (int k = start; k < end; k++)
                    {
                        qi[nz] = hi[k];
                        qx[nz] = a * hx[k];

                        nz++;
                    }
                }

                // Store the nz count of current column.
                qp[hi[j] + 1] = nz;
            }

            qp[columns] = nz;

            // There might be gaps in the column indices of Q.
            nz = qp[i];
            for (int j = i + 1; j < columns; j++)
            {
                if (qp[j] > 0)
                {
                    nz = qp[j];
                }
                else
                {
                    qp[j] = nz;
                }
            }

        }

        static bool CheckStorageRequirements(int rows, int cols, int nnz)
        {
            const int NZ_THRESHOLD = 100 * 100; // Matrix dimensions > 100 x 100

            if (nnz > NZ_THRESHOLD && nnz > rows * cols / 2)
            {
                Console.WriteLine("WARNING: matrix Q will be dense.");
            }

            const int MAX_BYTES = 200 * 1024 * 1024; // 200 MB

            return ((cols + 1) * 4 + nnz * 4 + nnz * 8) < MAX_BYTES;

        }

        static void PrintMatrix(string caption, SparseMatrix A)
        {
            Console.WriteLine(caption);

            for (int y = 0; y < A.RowCount; y++)
            {
                for (int x = 0; x < A.ColumnCount; x++)
                {
                    Console.Write("{0,5:0.###}", A.At(y, x));
                }
                Console.WriteLine();
            }

            Console.WriteLine();
        }
    }
    
    static class ExtensionMethods
    {
        public static void CopyTo(this SparseMatrix source, SparseMatrix target)
        {
            int m = source.RowCount;
            int n = source.ColumnCount;
            int nnz = source.NonZerosCount;

            if (target.RowCount != m || target.ColumnCount != n)
            {
                throw new ArgumentException("Invalid dimensions.", nameof(target));
            }

            if (target.NonZerosCount < nnz)
            {
                target.RowIndices = (int[])source.RowIndices.Clone();
                target.Values = (double[])source.Values.Clone();
            }
            else
            {
                Array.Copy(source.RowIndices, 0, target.RowIndices, 0, nnz);
                Array.Copy(source.Values, 0, target.Values, 0, nnz);
            }

            source.ColumnPointers.CopyTo(target.ColumnPointers, 0);
        }
    }
}

But keep in mind, that the above code is just for demonstration purpose. I haven't used it in production and wouldn't recommend doing so.