small update to the doc

doc/clean.sh

@@ -0,0 +1,8 @@
+sed -n '
+1h
+1!H
+$ {
+  g
+  s/\\begin{fulllineitems}\n*\\end{fulllineitems}//g
+  p
+}' build/latex/ABCDSolver.tex

doc/clean_full.sh

@@ -0,0 +1,16 @@
+#!/usr/bin/env bash
+
+sed '/^a test$/{
+       $!{ N        # append the next line when not on the last line
+         s/^a test\nPlease do not$/not a test\nBe/
+                    # now test for a successful substitution, otherwise
+                    #+  unpaired "a test" lines would be mis-handled
+         t sub-yes  # branch_on_substitute (goto label :sub-yes)
+         :sub-not   # a label (not essential; here to self document)
+                    # if no substituion, print only the first line
+         P          # pattern_first_line_print
+         D          # pattern_ltrunc(line+nl)_top/cycle
+         :sub-yes   # a label (the goto target of the 't' branch)
+                    # fall through to final auto-pattern_print (2 lines)
+       }    
+     }' alpha.txt  

doc/source/algorithm.rst

@@ -2,25 +2,26 @@
 The Algorithm
 =============
 
-The package `ABCD Solver` is a distributed hybrid (iterative/direct)
-solver for sparse linear systems :math:`Ax = b` where :math:`A` is a
-double precision matrix :math:`m \times n`, :math:`x` an
-:math:`n`-vector and :math:`b` an :math:`m`-vector.
-`ABCD Solver` uses two methods to solve the linear system:
+The package **Augmented Block Cimmino Distributed Solver** (`ABCD
+Solver`) is a distributed hybrid (iterative/direct) solver for sparse
+linear systems :math:`Ax = b` where :math:`A` is a double precision
+matrix :math:`m \times n`, :math:`x` an :math:`n`-vector and :math:`b`
+an :math:`m`-vector.  `ABCD Solver` uses two methods to solve the
+linear system:
 
 - *Regular Block Cimmino*: A block-projection technique that iterates
    to solve the linear system. During the iterations it solves a set
    of small problems (augmented systems built using the partitions of
    the original system).
-- *Augmented Block Cimmino*: A pseudo-direct technique that augments
-   the original system and through a succession of direct solves finds
-   the solution.
+- *Augmented Block Cimmino*: A pseudo-direct solver that augments the
+   original system and constructs the solution directly out of
+   independent solves using the augmented subsystems.
 
 The regular block Cimmino
 -------------------------
 
 The block Cimmino method is an iterative method that uses block-row
-projections. To solve the :math:`Ax = b`, where :math:`A` is an
+projections. To solve :math:`Ax = b`, where :math:`A` is an
 :math:`m\times n` sparse matrix, :math:`x` is an :math:`n`-vector and
 :math:`b` is an :math:`m`-vector, we subdivide the system into strips of
 rows as follows:
@@ -45,7 +46,7 @@ rows as follows:
 Let :math:`P_{\mathcal{R}(A_i^T)}` be the projector onto the range of
 :math:`A_i^T` and :math:`{A_i}^+` be the Moore-Penrose pseudo-inverse of the
 partition :math:`A_i`. The block Cimmino algorithm then computes a solution
-iteratively from an initial estimate :math:`x^{(0)}` according to:
+iteratively from an initial estimate :math:`x^{(0)}`, viz:
 
 .. math::
    :nowrap:
@@ -57,10 +58,10 @@ iteratively from an initial estimate :math:`x^{(0)}` according to:
         \end{array}
     \end{eqnarray*}
 
-where we see the independence of the set of :math:`p` equations, which is why
-the method is so attractive in a parallel environment.
+where we can see the independence of the set of :math:`p` equations,
+which is why the method is so attractive in a parallel environment.
 
-With the above notations, the iteration equations are thus:
+With the above notations, the iteration equation is thus:
 
 .. math::
    :nowrap:
@@ -73,16 +74,18 @@ With the above notations, the iteration equations are thus:
           \label{something}
     \end{eqnarray*}
 
-The iteration matrix for block Cimmino, :math:`H = I - Q`, is then a
-sum of projectors :math:`H = \omega
+The iteration matrix for block Cimmino method is :math:`H = I - Q`,
+which corresponds to a sum of projectors :math:`H = \omega
 \sum_{i=1}^p{\mathcal{P}_{\mathcal{R}(A_i^T)}}`. It is thus symmetric
 and positive definite and so we can solve
 
 .. math::
     H x ~=~ \xi,
 
 where :math:`\xi = \omega \sum_{i=1}^p{A_i^+ b_i}`
-using conjugate gradient or block conjugate gradient methods.  As :math:`\omega` appears on both sides of the equation, we can set it to one.
+using conjugate gradient or block conjugate gradient methods.  As the
+relaxation scalar :math:`\omega` appears on both sides of the
+equation, we can set it to one.
 
 At each step of the conjugate gradient algorithm we must solve for the
 :math:`p` projections viz.
@@ -94,7 +97,7 @@ At each step of the conjugate gradient algorithm we must solve for the
         A_i u_i ~=~ r_i, ~~~~  (r_i = {b_i - A_i x^{(k)}}),~~~ i = 1, .... p.
     \end{equation}
 
-In our approach we choose to solve these equations using the augmented system
+In our implementation we choose to solve these equations using the augmented system approach
 
 .. math::
    :nowrap:
@@ -105,21 +108,25 @@ In our approach we choose to solve these equations using the augmented system
         &=&  \left ( \begin{array}{l} 0 \\ r_i \end{array} \right )
     \end{eqnarray*}
 
-that we will solve, at each iteration, using a direct method and gives :math:`u_i = A_i^+ r_i` the projection we need for the partition :math:`A_i`.
-We use the multifrontal parallel solver :math:`MUMPS` to do this.
+that we perform, at each iteration, using a direct method, giving
+:math:`u_i = A_i^+ r_i` the projection needed for each partition
+:math:`A_i`.  We use the multifrontal parallel solver :math:`MUMPS` to
+do this.
 
 Running our solver in the regular mode will go through the following steps:
 
 - Partition the system into strips of rows (:math:`A_i` and :math:`b_i` for :math:`i = 1, \dots p`)
 - Create the augmented systems
 - Analyze and factorize the augmented systems using the direct solver :math:`MUMPS`
-- Run a block conjugate gradient with an implicit matrix :math:`H`, and at each iteration compute the matrix-vector product as a sum of projections. These projects being a set of solves using the direct solver.
+- Run a block conjugate gradient with an implicit iteration matrix
+  :math:`H`, which resumes into independent augmented systems direct
+  solves at each iteration.
 
 
 The augmented block Cimmino
 ---------------------------
 
-To understand the algorith, suppose that we have a matrix :math:`A` with three partitions, described as follow:
+To understand the algorithm, suppose that we have a matrix :math:`A` with three partitions, described as follows:
 
 .. math::
    :nowrap:
@@ -139,10 +146,11 @@ Where :math:`A_{i,j}` is the sub-part of :math:`A_i`, the :math:`i`-th partition
 
 The goal of the augmented block Cimmino algorithm is to make these
 three partitions mutually orthogonal to each other, meaning that the
-product of each couple of partitions is zero. We consider two
-different ways to augment the matrix to obtain these zero matrix products.
+inner product of each couple of partitions is zero. We consider two
+different ways to augment the matrix to obtain these zero matrix inner
+products.
 
-  * The first way to augment the matrix to make all the partitions mutually orthogonal to each other is by putting the product :math:`C_{ij} = A_{ij}A_{ji}^T` on the right of the partition :math:`A_i` and put :math:`-I` on the right of :math:`A_j` viz.
+  * The first way to augment the matrix to make all the partitions mutually orthogonal to each other is obtained by putting the product :math:`C_{ij} = A_{ij}A_{ji}^T` on the right of the partition :math:`A_i` and adding :math:`-I` on the right of :math:`A_j` viz.
 
   .. math::
     \bar{A} =
@@ -154,7 +162,9 @@ different ways to augment the matrix to obtain these zero matrix products.
     \end{array}\right].
 
     
-  * We can repeat the submatrices :math:`A_{ij}` and :math:`A_{ji}`, reversing the signs of one of them to obtain the augmented matrix :math:`\bar{A}` as in the following
+  * The second way is to repeat the submatrices :math:`A_{ij}` and
+    :math:`A_{ji}`, reversing the signs of one of them to obtain the
+    augmented matrix :math:`\bar{A}` as in the following
 
   .. math::
     \bar{A} =
@@ -168,6 +178,6 @@ different ways to augment the matrix to obtain these zero matrix products.
   This way :math:`\bar{A}_i\bar{A}_j^T` is zero for any pair :math:`i/j`, hence the new matrix has mutually orthogonal partitions.
 
 Notice that we augment the matrix upper-down and shift the
-augmentation at each step. This way, we do not create any new
+augmentation at each step. This is done such that we do not create any new
 interconnections between the new partitions. A simple check shows that
 :math:`\bar{A}_i \bar{A}_j^T` is zero for any pair :math:`i/j`.

doc/source/api.rst

@@ -1,11 +1,10 @@
 **********************************************
-API (work in progress)
+API
 **********************************************
 
-..
-   .. doxygenclass:: abcd
-       :project: abcd
-       :members: 
+.. doxygenclass:: abcd
+    :project: abcd
+    :members: 
 
-   .. doxygenfile:: defaults.h
-      :project: abcd
+.. doxygenfile:: defaults.h
+   :project: abcd

doc/source/installation.rst

@@ -2,9 +2,9 @@
 Installation
 ============
 
-The ABCD Solver depends on a few libraries: ``MUMPS 5.0``, ``Sparselib++ (custom)``, ``PaToH``, ``lapack`` and ``Boost::MPI``.
+The ABCD Solver depends on a few external libraries: ``MUMPS 5.0``, ``Sparselib++ (custom)``, ``PaToH``, ``lapack`` and ``Boost::MPI``.
 
-The installation can be done by writing the following commands:
+The installation can be done by typing the following commands in your terminal
 
 .. code-block:: bash
 

doc/source/xml/abcd_8h.xml

@@ -1004,9 +1004,9 @@
 <codeline lineno="384"><highlight class="normal"><sp/><sp/><sp/><sp/>std::vector&lt;double&gt;<sp/>drow_;</highlight></codeline>
 <codeline lineno="385"><highlight class="normal"><sp/><sp/><sp/><sp/>std::vector&lt;double&gt;<sp/>dcol_;</highlight></codeline>
 <codeline lineno="386"><highlight class="normal"></highlight></codeline>
-<codeline lineno="387"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/***************************************************************************</highlight></codeline>
+<codeline lineno="387"><highlight class="normal"><sp/><sp/><sp/><sp/></highlight><highlight class="comment">/**************************************************************************</highlight></codeline>
 <codeline lineno="388"><highlight class="comment"><sp/><sp/><sp/><sp/><sp/>*<sp/>The<sp/>matrix<sp/>object<sp/>itself</highlight></codeline>
-<codeline lineno="389"><highlight class="comment"><sp/><sp/><sp/><sp/>***************************************************************************/</highlight><highlight class="normal"></highlight></codeline>
+<codeline lineno="389"><highlight class="comment"><sp/><sp/><sp/><sp/>**************************************************************************/</highlight><highlight class="normal"></highlight></codeline>
 <codeline lineno="390"><highlight class="normal"><sp/><sp/><sp/><sp/><ref refid="struct_m_u_m_p_s" kindref="compound">MUMPS</ref><sp/>mumps_S;</highlight></codeline>
 <codeline lineno="391"><highlight class="normal"><sp/><sp/><sp/><sp/>Coord_Mat_double<sp/>S;</highlight></codeline>
 <codeline lineno="392"><highlight class="normal"><sp/><sp/><sp/><sp/>std::vector&lt;int&gt;<sp/>S_rows;</highlight></codeline>

include/abcd.h

@@ -258,13 +258,13 @@ class abcd
     abcd();
     ~abcd();
 
+
+private:
     int gqr(MV_ColMat_double &P, MV_ColMat_double &AP, MV_ColMat_double &R, int s, bool use_a);
     int gqr(MV_ColMat_double &p, MV_ColMat_double &ap, MV_ColMat_double &r, CompCol_Mat_double g, int s, bool use_a);
     void gmgs2(MV_ColMat_double &p, MV_ColMat_double &ap, MV_ColMat_double &r, int s, bool use_a);
     void gmgs2(MV_ColMat_double &p, MV_ColMat_double &ap, MV_ColMat_double &r, CompCol_Mat_double g, int s, bool use_a);
-
 
-private:
     // Types to be used localy
     double nrmA;
     double nrmB;
@@ -384,9 +384,9 @@ class abcd
     std::vector<double> drow_;
     std::vector<double> dcol_;
 
-    /***************************************************************************
+    /**************************************************************************
      * The matrix object itself
-    ***************************************************************************/
+    **************************************************************************/
     MUMPS mumps_S;
     Coord_Mat_double S;
     std::vector<int> S_rows;

src/preprocess/cij_augment.cpp

@@ -63,7 +63,7 @@ void abcd::cijAugmentMatrix(std::vector<CompCol_Mat_double> &M)
             CompCol_Mat_double C_ij;
             {
                 CompRow_Mat_double A_ij(sub_matrix(M[i], intersect));
-                CompRow_Mat_double A_ji(sub_matrix(M[j], intersect)));
+                CompRow_Mat_double A_ji(sub_matrix(M[j], intersect));
                 CompRow_Mat_double A_jiT = csr_transpose(A_ji);
                 C_ij = spmm(A_ij, A_jiT);
             }