diff --git a/docs/pdf/Makefile b/docs/pdf/Makefile
index 9c0301eb..2c928afa 100644
--- a/docs/pdf/Makefile
+++ b/docs/pdf/Makefile
@@ -85,8 +85,9 @@
 
 TOPFILE   = userguide.tex
 SECFILE   = title.tex abstract.tex overview.tex conventions.tex distribution.tex \
-	building.tex gettingstarted.tex  advanced.tex errors.tex \
-	listofroutines.tex bibliography.tex userinterface.tex background.tex
+		building.tex background.tex gettingstarted.tex userinterface.tex \
+            errors.tex bibliography.tex license.tex
+#		advanced.tex listofroutines.tex 
 FIGDIR    = figures
 
 XPDFFLAGS = 
diff --git a/docs/pdf/abstract.tex b/docs/pdf/abstract.tex
index 1589cb4a..ab4e28ba 100644
--- a/docs/pdf/abstract.tex
+++ b/docs/pdf/abstract.tex
@@ -5,7 +5,7 @@ It implements various versions of one-level additive and of multi-level additive
 and hybrid Schwarz algorithms. In the multi-level case, a purely algebraic approach
 is applied to generate coarse-level corrections, so that no geometric background is needed
 concerning the matrix to be preconditioned. The matrix is required to be square, real
-or complex, with a symmetric sparsity pattern 
+or complex, with a symmetric sparsity pattern. 
 
 MLD2P4 has been designed to provide scalable and easy-to-use preconditioners in the
 context of the PSBLAS (Parallel Sparse Basic Linear Algebra Subprograms)
@@ -14,12 +14,9 @@ available in this framework. MLD2P4 enables the user to easily specify different
 of a generic algebraic multilevel Schwarz preconditioner, thus allowing to search
 for the ``best'' preconditioner for the problem at hand. The package has been designed 
 employing object-oriented techniques, using Fortran 95 and MPI, with interfaces to
-additional external libraries such as UMFPACK, SuperLU and SuperLU\_Dist, that
+additional third party libraries such as UMFPACK, SuperLU and SuperLU\_Dist, that
 can be exploited in building multi-level preconditioners.
 
-The software is freely  distributable,  under the terms of the license
-in Appendix~\ref{sec:distribution}. 
-
 This guide provides a brief description of the functionalities and
 the user interface of MLD2P4.
 \end{abstract}
diff --git a/docs/pdf/background.tex b/docs/pdf/background.tex
index 5e23a6db..a6abf84a 100644
--- a/docs/pdf/background.tex
+++ b/docs/pdf/background.tex
@@ -1,6 +1,6 @@
 \section{Multi-level Domain Decomposition Background\label{sec:background}}
-\markboth{\underline{MLD2P4 User's and Reference Guide}}
-         {\underline{\ref{sec:background} Background}}
+\markboth{\textsc{MLD2P4 User's and Reference Guide}}
+         {\textsc{\ref{sec:background} Multi-level Domain Decomposition Background}}
 
 \emph{Domain Decomposition} (DD) preconditioners, coupled with Krylov iterative
 solvers, are widely used in the parallel solution of large and sparse linear systems.
@@ -17,7 +17,7 @@ solution of the original problem from the local solutions
 \emph{Additive Schwarz} preconditioners are DD preconditioners using overlapping
 submatrices, i.e.\ with some common rows, to couple the local information
 related to the submatrices (see, e.g., \cite{dd2_96}).
-The main motivations for choosing Additive Schwarz preconditioners are their
+The main motivation for choosing Additive Schwarz preconditioners is their
 intrinsic parallelism. A drawback of these
 preconditioners is that the number of iterations of the preconditioned solvers
 generally grows with the number of submatrices. This may be a serious limitation
@@ -29,16 +29,16 @@ in a coarse space, which globally couples the information related to the single
 submatrices. 
 
 \emph{Two-level Schwarz} preconditioners are obtained
-by combining basic (one-level) Schwarz preconditioners with coarse-level
-corrections. In this context, the one-level preconditioner is often
-called smoother. Different two-level preconditioners are obtained by varying the
-choice of the smoother, of the coarse-level correction and the
+by combining basic (one-level) Schwarz preconditioners with a coarse-level
+correction. In this context, the one-level preconditioner is often
+called `smoother'. Different two-level preconditioners are obtained by varying the
+choice of the smoother and of the coarse-level correction, and the
 way they are combined \cite{dd2_96}. The same reasoning can be applied starting
 from the coarse-level system, i.e.\ a coarse-space correction can be built
 from this system, thus obtaining \emph{multi-level} preconditioners.
 
 It is worth noting that optimal preconditioners do not necessarily correspond
-to minimum execution times. Indeed, to obtain effective multilevel preconditioners
+to minimum execution times. Indeed, to obtain effective multi-level preconditioners
 a tradeoff between optimality of convergence and the cost of building and applying
 the coarse-space corrections must be achieved. The choice of the number of levels,
 i.e.\ of the coarse-space corrections, also affects the effectiveness of the
@@ -74,7 +74,7 @@ we first provide a sketch of the AS preconditioners.
 
 Given the linear system \Ref{system1},
 where $A=(a_{ij}) \in \Re^{n \times n}$ is a
-nonsingular sparse matrix with a symmetric non-zero pattern,
+nonsingular sparse matrix with a symmetric nonzero pattern,
 let $G=(W,E)$ be the adjacency graph of $A$, where $W=\{1, 2, \ldots, n\}$
 and $E=\{(i,j) : a_{ij} \neq 0\}$ are the vertex set and the edge set of $G$,
 respectively. Two vertices are called adjacent if there is an edge connecting
@@ -84,7 +84,7 @@ Given a 0-overlap (or non-overlapping) partition of $W$,
 i.e.\ a set of $m$ disjoint nonempty sets $W_i^0 \subset W$ such that
 $\cup_{i=1}^m W_i^0 = W$, a $\delta$-overlap
 partition of $W$ is obtained by considering the sets
-$W_i^\delta \supset W_i^{\delta-1}$, obtained by including the vertices that
+$W_i^\delta \supset W_i^{\delta-1}$ obtained by including the vertices that
 are adjacent to any vertex in $W_i^{\delta-1}$.
 
 Let $n_i^\delta$ be the size of $W_i^\delta$ and $R_i^{\delta} \in 
@@ -107,12 +107,16 @@ to a vector $v \in \Re^n$ within a Krylov solver requires the following
 three steps:
 \begin{enumerate}
 	\item restriction of $v$ as $v_i = R_i^{\delta} v$, $i=1,\ldots,m$;
-	\item (approximate) solution of the linear systems $A_i^\delta w_i = v_i$,
+	\item solution of the linear systems $A_i^\delta w_i = v_i$,
 	      $i=1,\ldots,m$;
 	\item prolongation and sum of the $w_i$'s, i.e. $w = \sum_{i=1}^m (R_i^{\delta})^T w_i$.
 \end{enumerate}
+Note that the linear systems at step 2 are usually solved approximately,
+e.g.\ using incomplete LU factorizations such as ILU($p$), MILU($p$) and
+ILU($p,t$) \cite[Chapter 10]{Saad_book}.
+
 A variant of the classical AS preconditioner that outperforms it
-in terms of both convergence rate and of computation and communication
+in terms of convergence rate and of computation and communication
 time on parallel distributed-memory computers is the so-called \emph{Restricted AS
 (RAS)} preconditioner~\cite{CAI_SARKIS,EFSTATHIOU}. It
 is obtained by zeroing the components of $w_i$ corresponding to the
@@ -207,16 +211,21 @@ On the other hand, the use of few processors often leads to local sumatrices tha
 are too expensive to be processed on single processors, because of memory and/or
 computing requirements. Therefore, it seems natural to use a recursive approach,
 in which the coarse-level correction is re-applied starting from the current
-coarse-level system. The corresponding preconditioners are called \emph{multi-level}.
-One more reason for the multi-level approach is that it may significantly
-reduce the computational cost of preconditioning with respect to the two-level case
-(see \cite[Chapter 3]{dd2_96}). Additive and hybrid multilevel preconditioners
-are obtained as direct extensions of the two-level counterparts. 
-The algorithm for applying a multi-level version of the two-level hybrid 
-post-smoothed preconditioner is reported in Figure~\ref{fig:mlhpost_alg}.
-Other combinations of the smoothers and coarse-level corrections are possible, leading to variants
-of the previous algorithms. For a detailed descrition of them, the reader is
+coarse-level system. The corresponding preconditioners, called \emph{multi-level}
+preconditioners, can significantly reduce the computational cost of preconditioning
+with respect to the two-level case (see \cite[Chapter 3]{dd2_96}). 
+Additive and hybrid multilevel preconditioners
+are obtained as direct extensions of the two-level counterparts.
+For a detailed descrition of them, the reader is
 referred to \cite[Chapter 3]{dd2_96}.
+The algorithm for the application of a multi-level hybrid 
+post-smoothed preconditioner $M$ to a vector $v$, i.e.\ for the
+computation of $w=M^{-1}v$ is reported, for
+example, in Figure~\ref{fig:mlhpost_alg}. Here the number of levels
+is denoted by $nlev$ and the levels are numbered in increasing order starting
+from the finest one, i.e.\ the finest level is level 1; the coarse matrix
+and the corresponding basic preconditioner at each level $l$ are denoted by $A_l$ and
+$M_l$, respectively, with $A_1=A$.
 % 
 \begin{figure}[t]
 \begin{center}
@@ -224,40 +233,41 @@ referred to \cite[Chapter 3]{dd2_96}.
 \begin{minipage}{.85\textwidth} {\small
 \begin{tabbing}
 \quad \=\quad \=\quad \=\quad \\[-1mm]
-! assigned the finest matrix\\
-$A_1 \leftarrow A$;\\[1mm]
-! defined the number of levels $nlev$ \\[1mm]
-! defined $nlev-1$ prolongators\\
-$R_l^T, l=2, \ldots, nlev$;\\[1mm]
-! defined $nlev-1$ coarser matrices\\
-$A_l \leftarrow R_lA_{l-1}R_l^T, \; l=2, \ldots, nlev$;\\[1mm]
-! defined the $nlev-1$ basic Schwarz preconditioners\\
-$M_l$, basic preconditioner for $A_l \; l=1, \ldots, nlev-1$;\\[1mm]
-! assigned a vector $v$\\
-$v_1 \leftarrow v$; \\[2mm]
+%
+%! assign the finest matrix\\
+%$A_1 \leftarrow A$;\\[1mm]
+%! define the number of levels $nlev$ \\[1mm]
+%! define $nlev-1$ prolongators\\
+%$R_l^T, l=2, \ldots, nlev$;\\[1mm]
+%! define $nlev-1$ coarser matrices\\
+%$A_l \leftarrow R_lA_{l-1}R_l^T, \; l=2, \ldots, nlev$;\\[1mm]
+%! define the $nlev-1$ basic Schwarz preconditioners\\
+%$M_l$, basic preconditioner for $A_l \; l=1, \ldots, nlev-1$;\\[1mm]
+%$! assign a vector $v$\\
+%
+$v_1 = v$; \\[2mm]
 \textbf{for $l=2, nlev$ do}\\[1mm]
 \> ! transfer $v_{l-1}$ to the next coarser level\\
-\>  $v_l \leftarrow R_lv_{l-1}$; \\[1mm]
+\>  $v_l = R_lv_{l-1}$ \\[1mm]
 \textbf{endfor} \\[2mm]
 ! apply the coarsest-level correction\\[1mm]
-$y_{nlev} \leftarrow A_{nlev}^{-1}*v_{nlev}$;\\[2mm]
+$y_{nlev} = A_{nlev}^{-1} v_{nlev}$\\[2mm]
 \textbf{for $l=nlev -1 , 1, -1$ do}\\[1mm]
 \> ! transfer $y_{l+1}$ to the next finer level\\
-\> $y_l \leftarrow R_{l+1}^T*y_{l+1}$;\\[1mm]
+\> $y_l = R_{l+1}^T y_{l+1}$;\\[1mm]
 \> ! compute the residual at the current level\\
-\> $r_l \leftarrow v_l-A_l^{-1}*y_l$;\\[1mm]
-\> ! apply the basic Schwarz preconditioner to $r_l$\\
-\> $r_l \leftarrow M_l^{-1}*r_l$\\[1mm]
+\> $r_l = v_l-A_l^{-1} y_l$;\\[1mm]
+\> ! apply the basic Schwarz preconditioner to the residual\\
+\> $r_l = M_l^{-1} r_l$\\[1mm]
 \> ! update $y_l$\\
-\> $y_l \leftarrow y_l+r_l$\\
+\> $y_l = y_l+r_l$\\
 \textbf{endfor} \\[1mm]
-! preconditioned vector
-$w \leftarrow y_1$;
+$w = y_1$;
 \end{tabbing}
 }
 \end{minipage}
 }
-\caption{Multi-level hybrid post-smoothed preconditioner.\label{fig:mlhpost_alg}}
+\caption{Application of the multi-level hybrid post-smoothed preconditioner.\label{fig:mlhpost_alg}}
 \end{center}
 \end{figure}
 %
@@ -265,7 +275,7 @@ $w \leftarrow y_1$;
 
 \subsection{Smoothed Aggregation\label{sec:aggregation}}
 
-To define the restriction operator $R_C$, which is used to compute
+In order to define the restriction operator $R_C$, which is used to compute
 the coarse-level matrix $A_C$, MLD2P4 uses the \emph{smoothed aggregation}
 algorithm described in \cite{BREZINA_VANEK,VANEK_MANDEL_BREZINA}.
 The basic idea of this algorithm is to build a coarse set of vertices
@@ -283,16 +293,15 @@ Three main steps can be identified in the smoothed aggregation procedure:
 %\textbf{NOTA: Controllare cosa fa trilinos dopo il primo passo.}
  
 To perform the coarsening step, we have implemented the aggregation algorithm sketched
-in \cite{apnum_07}. According to \cite{BREZINA_VANEK}, a modification of this algorithm
-has been actually considered,
+in \cite{apnum_07}. According to \cite{VANEK_MANDEL_BREZINA}, a modification of
+this algorithm has been actually considered,
 in which each aggregate $N_r$ is made of vertices of $W$ that are \emph{strongly coupled}
 to a certain root vertex $r \in W$, i.e.\
 \[  N_r = \left\{s \in W: |a_{rs}| > \theta \sqrt{|a_{rr}a_{ss}|} \right\}
     \cup \left\{ r \right\} ,
 \]
-for a given $\theta \in [0,1]$. \textbf{L'ALGORITMO USA IL MAGGIORE STRETTO? ALTRIMENTI
-CON THETA=0 AGGREGHIAMO ANCHE I VERTICI CON COEFFICIENTE CORRISPONDENTE NULLO}
-Since the previous algorithm has a sequential nature, a \emph{decoupled} version of
+for a given $\theta \in [0,1]$.
+Since this algorithm has a sequential nature, a \emph{decoupled} version of
 it has been chosen, where each processor $i$ independently applies the algorithm to
 the set of vertices $W_i^0$ assigned to it in the initial data distribution. This
 version is embarrassingly parallel, since it does not require any data communication.
diff --git a/docs/pdf/bibliography.tex b/docs/pdf/bibliography.tex
index 9599bdb9..f05597d5 100644
--- a/docs/pdf/bibliography.tex
+++ b/docs/pdf/bibliography.tex
@@ -1,7 +1,8 @@
-\section{Bibliography\label{sec:bib}}
-\markboth{\underline{MLD2P4 User's and Reference Guide}}
-         {\underline{\ref{sec:bib} Bibliography}}
-\let\refname\relax
+%\section{Bibliography\label{sec:bib}}
+\markboth{\textsc{MLD2P4 User's and Reference Guide}}
+         {\textsc{References}}
+
+%\let\refname\relax
 
 \begin{thebibliography}{99}
 %
@@ -146,6 +147,11 @@ ACM Transactions on Mathematical Software, 29, 2, 2003, 110--140.
 %{\em Fortran 95/2003 explained.}
 %{Oxford University Press}, 2004.
 %
+\bibitem{Saad_book}
+Y.~Saad,
+\emph{Iterative methods for sparse linear systems}, 2nd edition,
+SIAM, 2003
+
 \bibitem{dd2_96}
 B.~Smith, P.~Bjorstad, W.~Gropp,
 {\em Domain Decomposition: Parallel Multilevel Methods for Elliptic
diff --git a/docs/pdf/building.tex b/docs/pdf/building.tex
index 6b28bf2f..7fc35840 100644
--- a/docs/pdf/building.tex
+++ b/docs/pdf/building.tex
@@ -1,9 +1,13 @@
-\section{Configuring and Building MLD2P4\label{sec:configuring}}
-\markboth{\underline{MLD2P4 User's and Reference Guide}}
-         {\underline{\ref{sec:configuring} Configuring and Building}}
-    - uso di GNU autoconf e automake \\
-    - software di base necessario (MPI, BLACS, BLAS, PSBLAS, UMFPACK ? - specificare versioni)\\
-    - software opzionale (SuperLU, SuperLUdist - specificare versioni e opzioni di configure)\\
-    - sistemi operativi e compilatori su cui MLD2P4 e' stato costruito con successo \\
+\section{Configuring and Building MLD2P4\label{sec:building}}
+\markboth{\textsc{MLD2P4 User's and Reference Guide}}
+         {\textsc{\ref{sec:building} Configuring and Building MLD2P4}}
+    - uso di GNU autoconf e automake, con riferimento alla compilazione delle varie versioni
+      (reale/complessa, singola/doppia precisione) \\
+    - software di base (MPI, BLACS, BLAS, PSBLAS - specificare versioni)\\
+    - software `third party" (UMFPACK, SuperLU, SuperLUdist - specificare versioni e opzioni di 
+    configure, dire che UMFPACK e SuperLU sono richiesti, rispettivamente, dalle versioni in 
+    singola ed in doppia precisione.)\\
+    - sistemi operativi e compilatori su cui MLD2P4 e' stato implementato con successo \\
     - sono previste opzioni di configurazione per il debugging o per il profiling? \\
-    - albero delle directory generato al momento dell'installazione\\
+    - albero delle directory generato al momento dell'installazione, con brevissima
+    descrizione del contenuto delle directory - aggiungere la directory \texttt{examples}\\
diff --git a/docs/pdf/conventions.tex b/docs/pdf/conventions.tex
index 2f59b148..32058331 100644
--- a/docs/pdf/conventions.tex
+++ b/docs/pdf/conventions.tex
@@ -1,8 +1,8 @@
 \section{Notational Conventions\label{sec:conventions}}
-\markboth{\underline{MLD2P4 User's and Reference Guide}}
-         {\underline{\ref{sec:conventions} Notational conventions}}
+\markboth{\textsc{MLD2P4 User's and Reference Guide}}
+         {\textsc{\ref{sec:conventions} Notational Conventions}}
 
     - caratteri tipografici usati nella guida (vedi guida ML recente e guida Aztec) \\
-    - convenzioni sui nomi di routine (differenza nei nomi tra high-level e medium-level),
-      strutture dati, moduli, costanti, etc. (vedi guida psblas) \\
+    - convenzioni sui nomi di routine (anche differenza nei nomi tra high-level e 	 
+      medium-level?), strutture dati, moduli, costanti, etc. (vedi guida psblas) \\
     - versione reale e complessa, singola e doppia precisione\\
\ No newline at end of file
diff --git a/docs/pdf/distribution.tex b/docs/pdf/distribution.tex
index f6cd18b8..8bc01b41 100644
--- a/docs/pdf/distribution.tex
+++ b/docs/pdf/distribution.tex
@@ -1,44 +1,7 @@
-\section{License\label{sec:distribution}}
-\markboth{\underline{MLD2P4 User's and Reference Guide}}
-         {\underline{\ref{sec:distribution} License}}
+\section{Code Distribution\label{sec:distribution}}
+\markboth{\textsc{MLD2P4 User's and Reference Guide}}
+         {\textsc{\ref{sec:distribution} Code Distribution}}
 
-The MLD2P4 is freely distributable under the following copyright
-terms: {\small
-\begin{verbatim} 
-                         MLD2P4  version 1.0
-MultiLevel Domain Decomposition Parallel Preconditioners Package
-           based on PSBLAS (Parallel Sparse BLAS version 2.3)
-
-(C) Copyright 2008
-
-                    Salvatore Filippone  University of Rome Tor Vergata       
-                    Alfredo Buttari      University of Rome Tor Vergata
-                    Pasqua D'Ambra       ICAR-CNR, Naples
-                    Daniela di Serafino  Second University of Naples
-
-
-Redistribution and use in source and binary forms, with or without
-modification, are permitted provided that the following conditions
-are met:
-  1. Redistributions of source code must retain the above copyright
-     notice, this list of conditions and the following disclaimer.
-  2. Redistributions in binary form must reproduce the above copyright
-     notice, this list of conditions, and the following disclaimer in the
-     documentation and/or other materials provided with the distribution.
-  3. The name of the MLD2P4 group or the names of its contributors may
-     not be used to endorse or promote products derived from this
-     software without specific written permission.
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
-TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
-BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-POSSIBILITY OF SUCH DAMAGE.
-\end{verbatim}
-}
\ No newline at end of file
+\noindent
+- sito web sul quale e' disponibile il software\\
+- modalita' di uso con cenno alla licenza (i dettagli sulla licenza vanno in appendice)
\ No newline at end of file
diff --git a/docs/pdf/errors.tex b/docs/pdf/errors.tex
index a427ca97..28fb8f40 100644
--- a/docs/pdf/errors.tex
+++ b/docs/pdf/errors.tex
@@ -1,6 +1,6 @@
 \section{Error Handling}\label{sec:errors}
-\markboth{\underline{MLD2P4 User's and Reference Guide}}
-         {\underline{\ref{sec:errors} Error handling}}
+\markboth{\textsc{MLD2P4 User's and Reference Guide}}
+         {\textsc{\ref{sec:errors} Error handling}}
 
 Error handling
     - Breve descrizione con rinvio alla guida di PSBLAS
diff --git a/docs/pdf/gettingstarted.tex b/docs/pdf/gettingstarted.tex
index 03a57217..4a919564 100644
--- a/docs/pdf/gettingstarted.tex
+++ b/docs/pdf/gettingstarted.tex
@@ -1,6 +1,6 @@
 \section{Getting Started\label{sec:started}}
-\markboth{\underline{MLD2P4 User's and Reference Guide}}
-         {\underline{\ref{sec:started} Getting started}}
+\markboth{\textsc{MLD2P4 User's and Reference Guide}}
+         {\textsc{\ref{sec:started} Getting Started}}
 
 We describe the basics for building and applying MLD2P4 one-level and multi-level
 Schwarz preconditioners with the Krylov solvers included in PSBLAS \cite{PSBLASGUIDE}.
@@ -17,15 +17,16 @@ The following steps are required:
 	a preconditioner type chosen by the user}. This is performed by the routine
 	\verb|mld_precinit|, which also sets defaults for each preconditioner
 	type selected by the user. The defaults associated to each preconditioner
-	type are listed in Table~\ref{tab:precinit}, where the strings used by
+	type are given in Table~\ref{tab:precinit}, where the strings used by
 	\verb|mld_precinit| to identify the preconditioner types are also given.
+	Note that these strings are valid also if uppercase letters are substituted by
+	corresponding lowercase ones.
 \item \emph{Modify the selected preconditioner type, by properly setting
   preconditioner parameters.} This is performed by the routine \verb|mld_precset|.
   This routine must be called only if the user wants to modify the default values
   of the parameters associated to the selected preconditioner type, to obtain a variant
-  of the preconditioner.
-  Examples of use of \verb|mld_precset| is given in Section~\ref{sec:examples};
-  a complete list of all the
+  of the preconditioner. Examples of use of \verb|mld_precset| are given in
+  Section~\ref{sec:examples}; a complete list of all the
   preconditioner parameters and their allowed and default values is provided in 
   Section~\ref{sec:userinterface}, Tables~\ref{tab:p_type}-\ref{tab:p_coarse}. 
 \item \emph{Build the preconditioner for a given matrix.} This is performed by
@@ -39,39 +40,43 @@ The following steps are required:
   be performed when the preconditioner is no more used.
 \end{enumerate}
 A detailed description of the above routines is given in Section~\ref{sec:userinterface}.
-
-Note that the Fortran 95 module \verb|mld_prec_mod| must be used in the program
-calling the MLD2P4 routines; this requires also  the use of the
-\verb|psb_base_mod| for the sparse matrix and communication descriptor
-data types, as well as for the kind parameters for vectors, and the
-use of the module \verb|psb_krylov_mod| for interfacing with the
-Krylov solvers. Note that the include path for MLD2P4 must override
-those for the base PSBLAS, e.g. they must come first in the sequence
-passed to the compiler, as the MLD2P4 version of the Krylov interfaces
-must override that of PSBLAS. This will change in the future when the
-support for the \verb|class| statement becomes widespread in Fortran
-compilers. 
 Examples showing the basic use of MLD2P4 are reported in Section~\ref{sec:examples}.
 
-\noindent
-\textbf{Remark.} The coarsest-level solver used by the default two-level
+Note that the Fortran 95 module \verb|mld_prec_mod|, containing the definition of the 
+preconditioner data type and the interfaces to the routines of MLD2P4,
+must be used in any program calling such routines.
+The modules \verb|psb_base_mod|, for the sparse matrix and communication descriptor
+data types, and \verb|psb_krylov_mod|, for interfacing with the
+Krylov solvers, must be also used (see Section~\ref{sec:examples}).
+
+\ \\
+\textbf{Remark 1.} The coarsest-level solver used by the default two-level
 preconditioner has been chosen by taking into account that, on parallel
 machines, it often leads to the smallest execution time when applied to
 linear systems coming from finite-difference discretizations of basic
 elliptic PDE problems, considered as standard tests for multi-level Schwarz
 preconditioners \cite{aaecc_07,apnum_07}. However, this solver does
-not necessarily  to the smallest number of iterations of the
+not necessarily correspond to the smallest number of iterations of the
 preconditioned Krylov method, which is usually obtained by applying a
-direct solver, e.g.\ based on the LU factorization, on a matrix
-replicated at the coarsest level (see Section~\ref{sec:userinterface}
-for coarsest-level solvers available in MLD2P4). 
+a direct solver to the coarsest-level system, e.g.\ based on the LU
+factorization (see Section~\ref{sec:userinterface}
+for the coarsest-level solvers available in MLD2P4). 
+
+\ \\
+\textbf{Remark 2.} The include path for MLD2P4 must override
+those for PSBLAS, e.g.\ the latter must come first in the sequence
+passed to the compiler, as the MLD2P4 version of the Krylov solver
+interfaces must override that of PSBLAS. This will change in the future
+when the support for the \verb|class| statement becomes widespread in Fortran
+compilers. 
+
 
 \begin{table}[th]
-{
 \begin{center}
-\begin{tabular}{|l|l|p{6.7cm}|}
+%{\small
+\begin{tabular}{|l|l|p{6.4cm}|}
 \hline
-\emph{Type}       & \emph{String} & \emph{Default preconditioner} \\ \hline
+\textsc{type}       & \textsc{string} & \textsc{default preconditioner} \\ \hline
 No preconditioner &\verb|'NOPREC'|& (Considered only to use the PSBLAS
                                     Krylov solvers with no preconditioner.) \\
 Diagonal          & \verb|'DIAG'| & --- \\
@@ -80,18 +85,23 @@ Additive Schwarz  & \verb|'AS'|   & Restricted Additive Schwarz (RAS),
                                     with overlap 1 and ILU(0) on the local blocks. \\ 
 Multilevel        &\verb|'ML'|    & Multi-level hybrid preconditioner (additive on the
                                     same level and multiplicative through the levels),
-                                    with post-smoothing only. Number of levels: 2;
-                                    post-smoother: RAS with overlap 1 and with ILU(0)
-                                    on the local blocks; coarsest matrix: distributed 
-                                    among the processors; (approximate) coarse-level
-                                    solver: 4 sweeps of the block-Jacobi solver, 
-                                    with the UMFPACK LU factorization
-                                    on the blocks (double precision versions) or
-                                    \textbf{XXXXXXXXX} (single precision versions)\\
+                                    with post-smoothing only.
+                                    Number of levels: 2.
+	                                  Post-smoother: RAS with overlap 1 and ILU(0)
+                                    on the local blocks.
+                                    Aggregation: smoothed aggregation with
+                                    threshold $\theta = 0$.
+                                    Coarsest matrix: distributed among the processors.
+                                    Coarse-level solver: 
+                                    4 sweeps of the block-Jacobi solver, 
+                                    with LU factorization of the blocks
+                                    (UMFPACK for the double precision versions and
+                                    SuperLU for the single precision ones)         \\
 \hline
 \end{tabular}
+%}
 \end{center}
-}
+
 \caption{Preconditioner types, corresponding strings and default choices.
 \label{tab:precinit}}
 \end{table}
@@ -103,31 +113,32 @@ multi-level preconditioner available in the real double precision version
 of MLD2P4 (see Table~\ref{tab:precinit}). This preconditioner is chosen
 by simply specifying \verb|'ML'| as second argument of \verb|mld_precinit|
 (a call to \verb|mld_precset| is not needed) and is applied with the BiCGSTAB
-solver provided by PSBLAS. The setup and application of the default multi-level
-preconditioners for the real single precision and the complex, single and double
-precision, versions are obtained with straightforward modifications of the example.
+solver provided by PSBLAS. As previously observed, the modules \verb|psb_base_mod|,
+\verb|mld_prec_mod| and \verb|psb_krylov_mod| must be used by the example program.
  
-
 The part of the code concerning the
 reading and assembling of the sparse matrix and the right-hand side vector, performed
 through the PSBLAS routines for sparse matrix and vector management, is not reported
-here for brevity; the statements concerning the deallocation of the PSBLAS data structure
-are neglected too.
-The complete code can be found in the example program file \verb|example_ml.f90|
-in the directory \textbf{XXXXXX (COMPLETARE. DIRE CHE I FILE IN REALTA' SONO DUE, UNO CON
-LA GENERAZIONE DELLA MATRICE ED UNO CON LA LETTURA).} Note that the modules \verb|psb_base_mod|
-and \verb|psb_util_mod| at the beginning of the code are required by PSBLAS.
-\textbf{O psb\_base\_mod} E' RICHIESTO ANCHE DA MLD2P4?)
-For details on the use of the PSBLAS routines, see the PSBLAS User's Guide \cite{PSBLASGUIDE}.
+here for brevity; the statements concerning the deallocation of the PSBLAS
+data structure are neglected too.
+The complete code can be found in the example program file \verb|mld_dexample_ml.f90|,
+in the directory \verb|examples/fileread| of the MLD2P4 tree (see
+Section~\ref{sec:building}).
+For details on the use of the PSBLAS routines, see the PSBLAS User's
+Guide \cite{PSBLASGUIDE}.
 
-\textbf{LE FIGURE SONO DECENTRATE, NONOSTANTE IL CENTER. CI VUOLE UNA MINIPAGE?}
+The setup and application of the default multi-level
+preconditioners for the real single precision and the complex, single and double
+precision, versions are obtained with straightforward modifications of the previous
+example (see Section~\ref{sec:userinterface} for details). If these versions are installed,
+the corresponding Fortran 95 codes are available in \verb|examples/fileread/| with these.
 
 \begin{figure}[tbp]
 \begin{center}
+\begin{minipage}{.90\textwidth} 
 {\small
 \begin{verbatim}
   use psb_base_mod
-  use psb_util_mod 
   use mld_prec_mod
   use psb_krylov_mod
 ... ...
@@ -138,6 +149,8 @@ For details on the use of the PSBLAS routines, see the PSBLAS User's Guide \cite
   type(psb_desc_type)   :: desc_A
 ! preconditioner
   type(mld_dprec_type)  :: P
+! right-hand side and solution vectors
+  real(kind(1.d0))      :: b(:), x(:)
 ... ...
 !
 ! initialize the parallel environment
@@ -145,20 +158,18 @@ For details on the use of the PSBLAS routines, see the PSBLAS User's Guide \cite
   call psb_info(ictxt,iam,np)
 ... ...
 !
-! read and assemble the matrix A and the right-hand
-! side b using PSBLAS routines for sparse matrix /
-! vector management
+! read and assemble the matrix A and the right-hand side b 
+! using PSBLAS routines for sparse matrix / vector management 
 ... ...
 !
-! initialize the default multi-level preconditioner,
-! i.e. two-level hybrid Schwarz, using RAS (with
-! overlap 1 and ILU(0) on the blocks) as post-smoother
-! and 4 block-Jacobi sweeps (with UMFPACK LU on the
-! blocks) as distributed coarse-level solver
+! initialize the default multi-level preconditioner, i.e. hybrid
+! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks) 
+! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU
+! on the blocks) as distributed coarse-level solver
   call mld_precinit(P,'ML',info)
 !
 ! build the preconditioner
-  call psb_precbld(A,P,desc_A,info)
+  call mld_precbld(A,desc_A,P,info)
 !
 ! set the solver parameters and the initial guess
   ... ...
@@ -178,85 +189,94 @@ For details on the use of the PSBLAS routines, see the PSBLAS User's Guide \cite
   stop
 \end{verbatim}
 }
+\end{minipage}
 \caption{Setup and application of the default multi-level Schwarz preconditioner.
 \label{fig:ex_default}}
 \end{center}
 \end{figure}
 
-Different versions of multilevel preconditioner can be obtained by changing
+Different versions of multi-level preconditioners can be obtained by changing
 the default values of the preconditioner parameters. The code reported in
 Figure~\ref{fig:ex_3lh} shows how to set a three-level hybrid Schwarz
 preconditioner, which uses block Jacobi with ILU(0) on the
-local blocks as post-smoother, a coarsest matrix replicated on the processors,
-and the LU factorization from UMFPACK~\cite{UMFPACK}, version 4.4, as coarse-level solver.
+local blocks as post-smoother, has a coarsest matrix replicated on the processors,
+and solves the coarsest-level system with the LU factorization from UMFPACK~\cite{UMFPACK}.
 The number of levels is specified by using \verb|mld_precinit|; the other
 preconditioner parameters are set by calling \verb|mld_precset|. Note that
 the type of multilevel framework (i.e.\ multiplicative among the levels
 with post-smoothing only) is not specified since it is the default 
-set by \verb|mld_precinit|. Figure~\ref{fig:ex_3la} shows how to
+set by \verb|mld_precinit|. 
+
+Figure~\ref{fig:ex_3la} shows how to
 set a three-level additive Schwarz preconditioner,
-which applies RAS, with overlap 1 and ILU(0) on the blocks,
-as pre- and post-smoother, and five block-Jacobi sweeps, with
+which uses RAS, with overlap 1 and ILU(0) on the blocks, 
+as pre- and post-smoother, and applies five block-Jacobi sweeps, with
 the UMFPACK LU factorization on the blocks, as distributed coarsest-level
 solver. Again, \verb|mld_precset| is used only to set
 non-default values of the parameters (see Tables~\ref{tab:p_type}-\ref{tab:p_coarse}).
 In both cases, the construction and the application of the preconditioner
 are carried out as for the default multi-level preconditioner.
 The code fragments shown in in Figures~\ref{fig:ex_3lh}-\ref{fig:ex_3la} are
-included in the example program file \verb|example_ml.f90|.
-\textbf{LO STESSO PROGRAMMA CONTIENE I TRE ESEMPI, CON UN SWITCH TRA L'UNO E L'ALTRO
-O FACCIAMO 3 PROGRAMMI DISTINTI? NON RICORDO CHE COSA ABBIAMO DECISO.
-PASQUA: ABBIAMO DETTO CHE ERA PREFERIBILE UN UNICO PROGRAMMA CON SWITCH.}
+included in the example program file \verb|mld_dexample_ml.f90| too.
 
 Finally, Figure~\ref{fig:ex_1l} shows the setup of a one-level
-additive Schwarz preconditioner, i.e. RAS with overlap 2. The corresponding code,
-including also the application of the preconditioner is in the example
-program file \verb|example_1lev.f90|.
+additive Schwarz preconditioner, i.e.\ RAS with overlap 2. The corresponding
+example program is available in \verb|mld_dexample_1lev.f90|.
 
-\begin{figure}[tbp]
+For all the previous preconditioners, example programs where the sparse matrix and
+the right-hand side are generated by discretizing a Poisson equation with Dirichlet
+boundary conditions are also available in the directory \verb|examples/pdegen|.
+
+
+\begin{figure}[tbh]
 \begin{center}
+\begin{minipage}{.90\textwidth} 
 {\small
 \begin{verbatim}
 ... ...
-! set a three-level hybrid Schwarz preconditioner,
-! which uses block Jacobi (with ILU(0) on the blocks)
-! as post-smoother, a coarsest matrix replicated on the
-! processors, and the LU factorization from UMFPACK
-! as coarse-level solver
+! set a three-level hybrid Schwarz preconditioner, which uses 
+! block Jacobi (with ILU(0) on the blocks) as post-smoother,
+! a coarsest matrix replicated on the processors, and the 
+! LU factorization from UMFPACK as coarse-level solver
   call mld_precinit(P,'ML',info,nlev=3)
   call_mld_precset(P,mld_smoother_type_,'BJAC',info)
-  call mld_precset(P,mld_coarse_mat,'REPL')
-  call mld_precset(P,mld_coarse_solve,'UMF')
+  call mld_precset(P,mld_coarse_mat_,'REPL',info)
+  call mld_precset(P,mld_coarse_solve_,'UMF',info)
 ... ...
 \end{verbatim}
 }
+\end{minipage}
+
 \caption{Setup of a hybrid three-level Schwarz preconditioner.\label{fig:ex_3lh}}
 \end{center}
 \end{figure}
 
-\begin{figure}[htb]
+\begin{figure}[tbh]
 \begin{center}
+\begin{minipage}{.90\textwidth} 
 {\small
 \begin{verbatim}
 ... ...
-! set a three-level additive Schwarz preconditioner,
-! which uses RAS (with overlap 1 and ILU(0) on the blocks)
-! as pre- and post-smoother, and 5 block-Jacobi sweeps
-! (with UMFPACK LU on the blocks) as distributed
-! coarsest-level solver
+! set a three-level additive Schwarz preconditioner, which uses 
+! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and 
+! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU
+! on the blocks) as distributed coarsest-level solver
   call mld_precinit(P,'ML',info,nlev=3)
   call mld_precset(P,mld_ml_type_,'ADD',info)
   call_mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
-  call mld_precset(P,mld_coarse_sweeps_,5)
+  call mld_precset(P,mld_coarse_sweeps_,5,info)
 ... ...
 \end{verbatim}
 }
+\end{minipage}
+
 \caption{Setup of an additive three-level Schwarz preconditioner.\label{fig:ex_3la}}
 \end{center}
 \end{figure}
 
-\begin{figure}[htb]
+\begin{figure}[tbh]
 \begin{center}
+\begin{minipage}{.90\textwidth} 
 {\small
 \begin{verbatim}
 ... ...
@@ -266,12 +286,13 @@ program file \verb|example_1lev.f90|.
 ... ...
 \end{verbatim}
 }
+\end{minipage}
 \caption{Setup of a one-level Schwarz preconditioner.\label{fig:ex_1l}}
 \end{center}
 \end{figure}
 
 \ \\
-\textbf{Remark.} Any PSBLAS-based program using the basic preconditioners
+\textbf{Remark 3.} Any PSBLAS-based program using the basic preconditioners
 implemented in PSBLAS 2.0, i.e.\ the diagonal and block-Jacobi ones,
 can use the diagonal and block-Jacobi preconditioners
 implemented in MLD2P4 without any change in the code.
diff --git a/docs/pdf/license.tex b/docs/pdf/license.tex
new file mode 100644
index 00000000..8f305b24
--- /dev/null
+++ b/docs/pdf/license.tex
@@ -0,0 +1,44 @@
+\section{License\label{sec:distribution}}
+\markboth{\textsc{MLD2P4 User's and Reference Guide}}
+         {\textsc{\ref{sec:distribution} License}}
+
+The MLD2P4 is freely distributable under the following copyright
+terms: {\small
+\begin{verbatim} 
+                         MLD2P4  version 1.0
+MultiLevel Domain Decomposition Parallel Preconditioners Package
+           based on PSBLAS (Parallel Sparse BLAS version 2.3)
+
+(C) Copyright 2008
+
+                    Salvatore Filippone  University of Rome Tor Vergata       
+                    Alfredo Buttari      University of Rome Tor Vergata
+                    Pasqua D'Ambra       ICAR-CNR, Naples
+                    Daniela di Serafino  Second University of Naples
+
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+  1. Redistributions of source code must retain the above copyright
+     notice, this list of conditions and the following disclaimer.
+  2. Redistributions in binary form must reproduce the above copyright
+     notice, this list of conditions, and the following disclaimer in the
+     documentation and/or other materials provided with the distribution.
+  3. The name of the MLD2P4 group or the names of its contributors may
+     not be used to endorse or promote products derived from this
+     software without specific written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
+BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.
+\end{verbatim}
+}
\ No newline at end of file
diff --git a/docs/pdf/overview.tex b/docs/pdf/overview.tex
index 76de1adf..1a675ad4 100644
--- a/docs/pdf/overview.tex
+++ b/docs/pdf/overview.tex
@@ -1,8 +1,7 @@
 \section{General Overview\label{sec:overview}}
-
-\markboth{\underline{MLD2P4 User's and Reference Guide}}
-         {\underline{\ref{sec:overview} General Overview}}
-
+\markboth{\textsc{MLD2P4 User's and Reference Guide}}
+         {\textsc{\ref{sec:overview} General Overview}}
+            
 The \textsc{Multi-Level Domain Decomposition Parallel Preconditioners Package based on
 PSBLAS (MLD2P4}) provides \emph{multi-level Schwarz preconditioners}~\cite{dd2_96},
 to be used in the iterative solutions of sparse linear systems:
@@ -31,9 +30,9 @@ as algebraic coarsening strategy~\cite{BREZINA_VANEK,VANEK_MANDEL_BREZINA}.
 The package is written in \emph{Fortran~95}, following an
 \emph{object-oriented approach} through the exploitation of features
 such as abstract data type creation, functional 
-overloading and dynamic memory management.% , while providing a smooth
+overloading and dynamic memory management.
+% , while providing a smooth
 % path towards the integration in legacy application codes. 
-% \textbf{NON MI PIACE QUESTO PERIODO, E' TROPPO LUNGO. RIUSCITE A SCRIVERLO MEGLIO?}
 The parallel implementation is based
 on a Single Program Multiple Data (SPMD) paradigm for distributed-memory architectures. 
 Single and double precision implementations of MLD2P4 are available for both the
@@ -42,7 +41,7 @@ real and the complex case, that can be used through a single interface.
 
 MLD2P4 has been designed to implement scalable and easy-to-use multilevel preconditioners
 in the context of the \emph{PSBLAS (Parallel Sparse BLAS)
-  computational framework}~\cite{psblas_00}. 
+computational framework}~\cite{psblas_00}. 
 PSBLAS is a library originally developed to address the parallel implementation of
 iterative solvers for sparse linear system, by providing basic linear algebra
 operators and data management facilities for distributed sparse matrices; it
@@ -71,7 +70,7 @@ by expert users to build new versions of multi-level Schwarz preconditioners.
 We provide here a description of the upper-layer routines, but not of the
 medium-layer ones.
 
-This guide is organized as follows: \textbf{ORGANIZZAZIONE DELLA GUIDA}
+This guide is organized as follows. \textbf{ORGANIZZAZIONE DELLA GUIDA}
 
 %%% Local Variables: 
 %%% mode: latex
diff --git a/docs/pdf/userguide.tex b/docs/pdf/userguide.tex
index 4725534f..50a834da 100644
--- a/docs/pdf/userguide.tex
+++ b/docs/pdf/userguide.tex
@@ -1,4 +1,4 @@
-\documentclass[a4paper]{article}
+\documentclass[a4paper,twoside,11pt]{article}
 \usepackage{pstricks}
 \usepackage{fancybox}
 \usepackage{amsfonts}
@@ -112,6 +112,7 @@
 
 \include{overview}
 \include{conventions}
+\include{distribution}
 \include{building}
 \include{background}
 \include{gettingstarted}
@@ -121,12 +122,11 @@
 %\include{listofroutines}
 \cleardoublepage
 \appendix
-\include{distribution}
-
+\include{license}
 \cleardoublepage
-
 \include{bibliography}
 
+
 \end{document}
 %%% Local Variables: 
 %%% mode: latex
diff --git a/docs/pdf/userinterface.tex b/docs/pdf/userinterface.tex
index 97af115e..02ed7a26 100644
--- a/docs/pdf/userinterface.tex
+++ b/docs/pdf/userinterface.tex
@@ -1,395 +1,420 @@
-\section{User Interface\label{sec:userinterface}}
-\markboth{\underline{MLD2P4 User's and Reference Guide}}
-         {\underline{\ref{sec:userinterface} User Interface}}
-
-The basic user interface of MLD2P4 consists of six routines. The four routines \verb|mld_precinit|,
-\verb|mld_precset|, \verb|mld_precbld| and \verb|mld_precaply| encapsulate all the functionalities for the setup and application of any one-level and multi-level
-preconditioner implemented in the package.
-The routine \verb|mld_precfree| deallocates the preconditioner data structure, while
-\verb|mld_precdescr| prints a description of the preconditioner setup by the user.
-
-For each routine, the same user interface is overloaded with
-respect to the real/complex case and the single/double precision;
-arguments with appropriate data types must be passed to the routine,
-i.e.
-\begin{itemize}
-\item the sparse matrix data structure, containing the matrix to be
-  preconditioned, must be of type \verb|mld_|\emph{x}\verb|spmat_type|
-	with \emph{x} = \verb|s| for real single precision, \emph{x} = \verb|d|
-	for real double precision, \emph{x} = \verb|c| for complex single precision,
-	\emph{x} = \verb|z| for complex double precision;
-\item the preconditioner data structure must be of type
-  \verb|mld_|\emph{x}\verb|prec_type|, with \emph{x} =    
-  \verb|s|, \verb|d|, \verb|c|, \verb|z|, according to the sparse
-  matrix data structure;
-\item the arrays containing the vectors $v$ and $w$ involved in
-  the preconditioner application $w=M^{-1}v$ must be of type   
-  \emph{type}\verb|(|\emph{kind\_parameter}\verb|)|, with \emph{type} =
-  \verb|real|, \verb|complex| and \emph{kind\_parameter} = \verb|kind(1.e0)|,
-  \verb|kind(1.d0)|, according to the sparse matrix and preconditioner
-  data structure; note that the PSBLAS module provides the constants \verb|psb_spk_|
-  = \verb|kind(1.e0)| and \verb|psb_dpk_| = \verb|kind(1.d0)|;
-\item real parameters defining the preconditioner must be declared
-  according to the precision of the previous data structures
-  (see Section~\ref{sec:precset}).
-\end{itemize}
-A description of each routine is given in the remainder of this section.
-
-
-\subsection{Subroutine mld\_precinit\label{sec:precinit}}
-
-\begin{center}
-\verb|mld_precinit(p,ptype,info)| \\
-\verb|mld_precinit(p,ptype,info,nlev)| \\
-\end{center}
-
-\noindent
-This routine allocates and initializes the preconditioner data structure,
-according to the preconditioner type chosen by the user.
-
-\subsubsection*{Arguments}
-
-\begin{tabular}{p{1.2cm}p{11.5cm}}
-\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(inout)|.\\
-              & The preconditioner data structure. Note that \emph{x}
-                must be chosen according to the real/complex, single/double
-                precision version of MLD2P4 under use.\\
-\verb|ptype|  & \verb|character(len=*), intent(in)|.\\
-              & The type of preconditioner. Its values are specified in Table~\ref{tab:precinit}.\\
-\verb|info|   & \verb|integer, intent(out)|.\\
-              & Error code. See Section~\ref{sec:errors} for details.\\
-\verb|nlev|   & \verb|integer, optional, intent(in)|.\\
-              & The number of levels of the multilevel preconditioner.
-                If \verb|nlev| is not present and \verb|ptype|='ML'/'ml', 
-                then \verb|nlev|=2 is assumed. Otherwise, \verb|nlev| is ignored.
-\end{tabular}
-
-
-\subsection{Subroutine mld\_precset\label{sec:precset}}
-
-\begin{center}
-\verb|mld_precset(p,what,val,info)|\\
-\end{center}
-
-\noindent
-This routine sets the parameters defining the preconditioner. More
-precisely, the parameter identified by \verb|what| is assigned the value
-contained in \verb|val|.
-
-\subsubsection*{Arguments}
-
-\begin{tabular}{p{1.2cm}p{11.5cm}}
-\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(inout)|.\\
-              & The preconditioner data structure. Note that \emph{x} must
-                be chosen according to the real/complex, single/double precision
-                 version of MLD2P4 under use.\\
-\verb|what|   & \verb|integer, intent(in)|. \\
-              & The number identifying the parameter to be set.
-                A mnemonic constant has been associated to each of these
-                numbers, as reported in Tables~\ref{tab:p_type}-\ref{tab:p_coarse}.\\
-\verb|val |   & \verb|integer| \emph{or} \verb|character(len=*)| \emph{or}
-                \verb|real(psb_spk_)| \emph{or} \verb|real(psb_dpk_)|,
-                \verb|intent(in)|.\\
-              & The value of the parameter to be set. The list of allowed
-                values and the corresponding data types is given in
-                Table~\ref{tab:params}.\\
-\verb|info|   & \verb|integer, intent(out)|.\\
-              & Error code. See Section~\ref{sec:errors} for details.\\
-%\verb|ilev|   & \verb|integer, optional, intent(in)|.\\
-%              & For the multilevel preconditioner, the level at which the
-%                preconditioner parameter has to be set.
-%                The levels are numbered in increasing
-%                order starting from the finest one, i.e.\ level 1 is the finest level.
-%                If \verb|ilev| is not present, the parameter identified by \verb|what|
-%                is set at all the appropriate levels (see Table~\ref{tab:params}).
-\end{tabular}
-
-\ \\
-A variety of (one-level and multi-level) preconditioner can be obtained
-by a suitable setting of the preconditioner parameters. These parameters
-can be logically divided into four groups, i.e.\ parameters defining
-\begin{enumerate}
-	\item the type of multi-level preconditioner;
-	\item the one-level preconditioner to be used as smoother;
-	\item the aggregation algorithm;
-	\item the coarse-space correction at the coarsest level.
-\end{enumerate}
-
-A list of the parameters that can be set, along with their allowed and
-default values, is given in Tables~\ref{tab:p_type}-\ref{tab:p_coarse}. 
-%Note that the routine allows to set different features of the
-%preconditioner at each level through the use of \verb|ilev|.
-%This should be done by users with experience in the field of
-%multi-level preconditioners. Non-expert users are recommended
-%to call \verb| mld_precset| without specifying \verb|ilev|.
-\textbf{CORREGGERE LA ROUTINE E LA DOC INTERNA - ilev NON E' PIU'
-ACCESSIBILE ALL'UTENTE.}
-
-
-\begin{sidewaystable}
-\begin{center}
-\begin{tabular}{|l|l|p{2cm}|l|p{7cm}|}
-\hline
-\verb|what|              & \emph{data type}        &  \verb|val|      &  \emph{default}  &
-\emph{comments} \\ \hline
-%\multicolumn{5}{|c|}{\emph{type of the multi-level preconditioner}}\\ \hline
-\verb|mld_ml_type_|      & \verb|character(len=*)|
-                         & 'ADD' \ \ \ 'MULT'   
-                         & 'MULT'
-                         & basic multi-level framework: additive or multiplicative
-                           among the levels always additive inside a level)         \\  
-\verb|mld_smoother_type_|& \verb|character(len=*)|
-                         & 'DIAG' \ \ \ 'BJAC' \ \ \ 'AS'
-                         & 'AS'
-                         & basic one-level preconditioner (i.e.\ smoother) of the
-                           multi-level preconditioner \\
-\verb|mld_smoother_pos_| & \verb|character(len=*)|
-                         & 'PRE' \ \ \ 'POST' \ \ \ 'TWOSIDE'
-                         & 'POST'
-                         & ``position'' of the smoother: pre-smoother, post-smoother, 
-                           pre-/post-smoother \\
-\hline
-\end{tabular}
-\end{center}
-\caption{Parameters defining the type of multi-level preconditioner.
-\label{tab:p_type}}                       
-\end{sidewaystable}
-                   
-\begin{sidewaystable}
-\begin{center}
-\begin{tabular}{|l|l|p{2cm}|l|p{7cm}|}
-\hline
-\verb|what|              & \emph{data type}        &  \verb|val|      &  \emph{default}  &
-\emph{comments} \\ \hline
-%\multicolumn{5}{|c|}{\emph{basic one-level preconditioner (smoother)}} \\ \hline
-\verb|mld_sub_ovr|       & \verb|integer|
-                         & any number $\ge 0$
-                         & 1
-                         & number of overlap in the basic Schwarz preconditioner   \\
-\verb|mld_sub_restr_|    & \verb|character(len=*)|
-                         & 'HALO' \ \ \ 'NONE'
-                         & 'HALO'
-                         & type of restriction operator used in basic Schwarz preconditioner: 'HALO' for taking into account contributions from the overlap \\
-\verb|mld_sub_prol_|     & \verb|character(len=*)|
-                         & 'SUM' \ \ \ 'NONE'
-                         & 'NONE'
-                         & type of prolongator operator used in basic Schwarz preconditioner: 'NONE' for neglecting contributions from the overlap    \\
-\verb|mld_sub_solve_|    & \verb|character(len=*)|
-                         & 'ILU' \ \ \ 'MILU' \ \ \ 'ILUT' \ \ \ 'UMF' \ \ \ 'SLU'
-                         & 'UMF'
-                         & available local solver: 'ILU' for incomplete LU, 'MILU' for modified incomplete LU, 'ILUT' 
-for incomplete LU with threshold, 'UMF' for complete LU using UMFPACK~\cite{UMFPACK} version 4.4, 'SLU' for complete LU using SuperLU~\cite{SUPERLU}, version 3.0  \\    
-\verb|mld_sub_fillin_|   & \verb|integer|
-                         & any number $\ge 0$
-                         & 0
-                         & fill-in level for 'ILU', 'MILU' and 'ILUT' of local blocks\\
-\verb|mld_sub_thresh_|   & \verb|real| 
-                         & any number $\ge 0.$
-                         & 0.
-                         & drop tolerance for 'ILUT' \\
-\verb|mld_sub_ren_|      & \verb|character(len=*)|
-                         & 'RENUM\_NONE', 'RENUM\_GLOBAL' %, 'RENUM\_GPS' 
-                         &
-                         & reordering algorithm for the local blocks \\
-\hline
-
-\end{tabular}
-\end{center}
-\caption{Parameters defining the basic one-level preconditioner (smoother).
-\label{tab:p_smoother}}  
-\end{sidewaystable}                     
-                   
-\begin{sidewaystable}
-\begin{center}
-\begin{tabular}{|l|l|p{2cm}|l|p{7cm}|}
-\hline
-\verb|what|              & \emph{data type}        &  \verb|val|      &  \emph{default}  &
-\emph{comments} \\ \hline
-%\multicolumn{5}{|c|}{\emph{aggregation algorithm}} \\ \hline
-\verb|mld_aggr_alg_|     & \verb|character(len=*)|
-                         & 'DEC'
-                         & 'DEC'
-                         & define the aggregation scheme. Now, only decoupled aggregation is available\\
-\verb|mld_aggr_kind_|    & \verb|character(len=*)|
-                         & 'SMOOTH', 'RAW'
-                         & 'SMOOTH'
-                         & define the type of aggregation technique (smoothed or nonsmoothed).    \\
-\verb|mld_aggr_thresh_|  & \verb|real|
-                         & any number $\in [0, 1]$
-                         & 0.
-                         & dropping threshold in aggregation    \\
-\verb|mld_aggr_eig_|     & \verb|character(len=*)|
-                         & 'A\_NORMI'
-                         &
-                         & define the algorithm to evaluate the maximum eigenvalue of $D^{-1}A$ for smoothed
-aggregation. Currently only the infinity norm of the matrix A is available\\
-\hline
-\end{tabular}
-\end{center}
-\caption{Parameters defining the aggregation algorithm.
-\label{tab:p_aggregation}} 
-\end{sidewaystable}
-                     
-\begin{sidewaystable}
-\begin{center}
-\begin{tabular}{|l|l|p{2cm}|l|p{7cm}|}
-\hline
-\verb|what|              & \emph{data type}        &  \verb|val|      &  \emph{default}  &
-\emph{comments} \\ \hline
-%\multicolumn{5}{|c|}{\emph{coarse-space correction at the coarsest level}}\\ \hline
-\verb|mld_coarse_mat_|   & \verb|character(len=*)|
-                         & 'DISTR', 'REPL'
-                         & 'DISTR'
-                         & Coarse matrix: distributed or replicated    \\
-\verb|mld_coarse_solve_| & \verb|character(len=*)|
-                         & 'BJAC' \ \ \ 'UMF' \ \ \ 'SLU' \ \ \ 'SLUDIST'
-                         & 'BJAC'
-                         & Only 'BJAC' and 'SLUDIST' can be used for distributed coarse matrix. 'BJAC' corresponds to some sweeps of a block-Jacobi solver, while 'SLUDIST' corresponds
-to the use of the external package SuperLU\_Dist~\cite{SUPERLUDIST}, version 2.0, for distributed sparse factorization and solve.  \\
-\verb|mld_coarse_subsolve_| & \verb|character(len=*)|
-                         & 'ILU' \ \ \ 'MILU' \ \ \ 'ILUT' \ \ \ 'UMF' \ \ \ 'SLU'
-                         & 'UMF'
-                         & available solver for diagonal local blocks of the coarse matrix, when 'BJAC' is used as coarse solver\\
-\verb|mld_coarse_sweeps_|& \verb|integer|                         
-                         & any number $> 0$
-                         & 4
-                         & number of Block-Jacobi sweeps when 'BJAC' is used as coarse solver    \\
-\verb|mld_coarse_fillin_| & \verb|integer|
-                         & any number $\ge 0$
-                         & 0
-                         & fill-in level in incomplete factorization of local diagonal blocks of the coarse matrix, when 'BJAC' is used as coarse solver and 'ILU' or 'MILU' is used as local solver  \textbf{MODIFICA NOME PARAM. NEL SW} \\
-\verb|mld_coarse_thresh_| & \verb|real|
-                         & any number $\ge 0.$
-                         & 0.
-                         & drop tolerance in incomplete factorization of local diagonal blocks of the coarse matrix, when 'BJAC' is used as coarse solver and 'ILUT' is used as local solver   \\ \hline
-\end{tabular}
-\end{center}
-\caption{Parameters defining the coarse-space correction at the coarsest
-level.\label{tab:p_coarse}} 
-\end{sidewaystable}
-
-
-\clearpage
-\subsection{Subroutine mld\_precbld\label{sec:precbld}}
-
-\begin{center}
-\verb|mld_precbld(a,desc_a,p,info)|\\
-\end{center}
-
-\noindent
-This routine builds the preconditioner according to the requirements made by
-the user through the routines \verb|mld_precinit| and \verb|mld_precset|.
-
-\subsubsection*{Arguments}
-
-\begin{tabular}{p{1.2cm}p{11.5cm}}
-\verb|a|      & \verb|type(psb_|\emph{x}\verb|spmat_type), intent(in)|. \\
-              & The sparse matrix structure containing the local part of the
-                matrix to be preconditioned. Note that \emph{x} must be chosen according
-                to the real/complex, single/double precision version of MLD2P4 under use.
-                See the PSBLAS User's Guide for details \cite{PSBLASGUIDE}.\\
-\verb|desc_a| & \verb|type(psb_desc_type), intent(in)|. \\
-              & The communication descriptor of a. See the PSBLAS User's Guide for
-                details \cite{PSBLASGUIDE}.\\
-\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(inout)|.\\
-              & The preconditioner data structure. Note that \emph{x} must be chosen according
-                to the real/complex, single/double precision version of MLD2P4 under use.\\
-\verb|info|   & \verb|integer, intent(out)|.\\
-              & Error code. See Section~\ref{sec:errors} for details.\\
-\end{tabular}
-
-\clearpage
-\subsection{Subroutine mld\_precaply\label{sec:precaply}}
-
-\begin{center}
-\verb|mld_precaply(p,x,y,desc_a,info)|\\
-\verb|mld_precaply(p,x,y,desc_a,info,trans,work)|\\
-\end{center}
-
-\noindent
-This routine computes $y = op(M^{-1})\, x$, where $M$ is a previously built
-preconditioner, stored in the \verb|p| data structure, and $op$
-denotes the preconditioner itself or its transpose, according to
-the value of \verb|trans|.
-Note that, when MLD2P4 is used with a Krylov solver from PSBLAS,
-\verb|mld_precaply| is called within the PSBLAS routine \verb|mld_krylov|
-and hence is completely transparent to the user.
-
-\subsubsection*{Arguments}
-
-\begin{tabular}{p{1.2cm}p{11.5cm}}
-\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(inout)|.\\
-              & The preconditioner data structure, containing the local part of $M$.
-                Note that \emph{x} must be chosen according
-                to the real/complex, single/double precision version of MLD2P4 under use.\\
-\verb|x|      & \emph{type}\verb|(|\emph{kind\_parameter}\verb|), dimension(:), intent(in)|.\\
-              & The local part of the vector $x$. Note that \emph{type} and   
-                \emph{kind\_parameter} must be chosen according
-                to the real/complex, single/double precision version of MLD2P4 under use.\\
-\verb|y|      & \emph{type}\verb|(|\emph{kind\_parameter}\verb|), dimension(:), intent(out)|.\\
-              & The local part of the vector $y$. Note that \emph{type} and
-                \emph{kind\_parameter} must be chosen according
-                to the real/complex, single/double precision version of MLD2P4 under use.\\
-\verb|desc_a| & \verb|type(psb_desc_type), intent(in)|. \\
-              & The communication descriptor associated to the matrix to be
-                preconditioned.\\
-\verb|info|   & \verb|integer, intent(out)|.\\
-              & Error code. See Section~\ref{sec:errors} for details.\\
-\verb|trans|  & \verb|character(len=1), optional, intent(in).|\\
-              & If \verb|trans| = \verb|'N','n'| then $op(M^{-1}) = M^{-1}$;
-                if \verb|trans| = \verb|'T','t'| then $op(M^{-1}) = M^{-T}$
-                (transpose of $M^{-1})$;  if \verb|trans| = \verb|'C','c'| then $op(M^{-1}) = M^{-C}$
-                (conjugate transpose of $M^{-1})$.\\
-\verb|work|  & \emph{type}\verb|(|\emph{kind\_parameter}\verb|), dimension(:), optional, target|.\\
-             & Workspace. Its size should be at
-               least \verb|4 * psb_cd_get_local_cols(desc_a)| (see the PSBLAS User's Guide).
-               Note that \emph{type} and \emph{kind\_parameter} must be chosen according
-               to the real/complex, single/double precision version of MLD2P4 under use.\\
-\end{tabular}
-
-
-\subsection{Subroutine mld\_precfree\label{sec:precfree}}
-
-\begin{center}
-\verb|mld_precfree(p,info)|\\
-\end{center}
-
-\noindent
-This routine deallocates the preconditioner data structure.
-
-\subsubsection*{Arguments}
-
-\begin{tabular}{p{1.2cm}p{11.5cm}}
-\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(inout)|.\\
-              & The preconditioner data structure. Note that \emph{x} must be chosen according
-                to the real/complex, single/double precision version of MLD2P4 under use.\\
-\verb|info|   & \verb|integer, intent(out)|.\\
-              & Error code. See Section~\ref{sec:errors} for details.\\
-\end{tabular}
-
-
-\subsection{Subroutine mld\_precdescr\label{sec:precdescr}}
-
-\begin{center}
-\verb|mld_precdescr(p,iout)|\\
-\end{center}
-
-\noindent
-This routine prints a description of the preconditioner
-to a file.
-
-\subsubsection*{Arguments}
-
-\begin{tabular}{p{1.2cm}p{10.6cm}}
-\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(in)|.\\
-              & The preconditioner data structure. Note that \emph{x} must be chosen according
-                to the real/complex, single/double precision version of MLD2P4 under use.\\
-\verb|iout|   & \verb|integer, intent(in), optional|.\\
-              & The id of the file where the preconditioner description
-                will be printed, default is  standard output.\\
-\end{tabular}
-
-%%% Local Variables: 
-%%% mode: latex
-%%% TeX-master: "userguide"
-%%% End: 
+\section{User Interface\label{sec:userinterface}}
+\markboth{\textsc{MLD2P4 User's and Reference Guide}}
+         {\textsc{\ref{sec:userinterface} User Interface}}
+
+The basic user interface of MLD2P4 consists of six routines. The four routines \verb|mld_precinit|,
+\verb|mld_precset|, \verb|mld_precbld| and \verb|mld_precaply| encapsulate all the functionalities
+for the setup and the application of any one-level and multi-level
+preconditioner implemented in the package.
+The routine \verb|mld_precfree| deallocates the preconditioner data structure, while
+\verb|mld_prec\-descr| prints a description of the preconditioner setup by the user.
+
+For each routine, the same user interface is overloaded with
+respect to the real/complex case and the single/double precision;
+arguments with appropriate data types must be passed to the routine,
+i.e.
+\begin{itemize}
+\item the sparse matrix data structure, containing the matrix to be
+  preconditioned, must be of type \verb|mld_|\emph{x}\verb|spmat_type|
+	with \emph{x} = \verb|s| for real single precision, \emph{x} = \verb|d|
+	for real double precision, \emph{x} = \verb|c| for complex single precision,
+	\emph{x} = \verb|z| for complex double precision;
+\item the preconditioner data structure must be of type
+  \verb|mld_|\emph{x}\verb|prec_type|, with \emph{x} =    
+  \verb|s|, \verb|d|, \verb|c|, \verb|z|, according to the sparse
+  matrix data structure;
+\item the arrays containing the vectors $v$ and $w$ involved in
+  the preconditioner application $w=M^{-1}v$ must be of type   
+  \emph{type}\verb|(|\emph{kind\_parameter}\verb|)|, with \emph{type} =
+  \verb|real|, \verb|complex| and \emph{kind\_parameter} = \verb|kind(1.e0)|,
+  \verb|kind(1.d0)|, according to the sparse matrix and preconditioner
+  data structure; note that the PSBLAS module \verb|psb_base_mod|
+  provides the constants \verb|psb_spk_|
+  = \verb|kind(1.e0)| and \verb|psb_dpk_| = \verb|kind(1.d0)|;
+\item real parameters defining the preconditioner must be declared
+  according to the precision of the previous data structures
+  (see Section~\ref{sec:precset}).
+\end{itemize}
+A description of each routine is given in the remainder of this section.
+
+
+\subsection{Subroutine mld\_precinit\label{sec:precinit}}
+
+\begin{center}
+\verb|mld_precinit(p,ptype,info)| \\
+\verb|mld_precinit(p,ptype,info,nlev)| \\
+\end{center}
+
+\noindent
+This routine allocates and initializes the preconditioner data structure,
+according to the preconditioner type chosen by the user.
+
+\subsubsection*{Arguments}
+
+\begin{tabular}{p{1.2cm}p{10.5cm}}
+\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(inout)|.\\
+              & The preconditioner data structure. Note that \emph{x}
+                must be chosen according to the real/complex, single/double
+                precision version of MLD2P4 under use.\\
+\verb|ptype|  & \verb|character(len=*), intent(in)|.\\
+              & The type of preconditioner. Its values are specified in Table~\ref{tab:precinit}.\\
+                Note that the uppercase and lowercase letters can be used indifferently.\\
+\verb|info|   & \verb|integer, intent(out)|.\\
+              & Error code. If no error, 0 is returned. See Section~\ref{sec:errors} for details.\\
+\verb|nlev|   & \verb|integer, optional, intent(in)|.\\
+              & The number of levels of the multilevel preconditioner.
+                If \verb|nlev| is not present and \verb|ptype|=\verb|'ML'|, \verb|'ml'|, 
+                then \verb|nlev|=2 is assumed. Otherwise, \verb|nlev| is ignored.\\
+\end{tabular}
+
+
+\subsection{Subroutine mld\_precset\label{sec:precset}}
+
+\begin{center}
+\verb|mld_precset(p,what,val,info)|\\
+\end{center}
+
+\noindent
+This routine sets the parameters defining the preconditioner. More
+precisely, the parameter identified by \verb|what| is assigned the value
+contained in \verb|val|.
+
+\subsubsection*{Arguments}
+
+\begin{tabular}{p{1.2cm}p{10.5cm}}
+\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(inout)|.\\
+              & The preconditioner data structure. Note that \emph{x} must
+                be chosen according to the real/complex, single/double precision
+                 version of MLD2P4 under use.\\
+\verb|what|   & \verb|integer, intent(in)|. \\
+              & The number identifying the parameter to be set.
+                A mnemonic constant has been associated to each of these
+                numbers, as reported in Tables~\ref{tab:p_type}-\ref{tab:p_coarse}.\\
+\verb|val |   & \verb|integer| \emph{or} \verb|character(len=*)| \emph{or}
+                \verb|real(psb_spk_)| \emph{or} \verb|real(psb_dpk_)|,
+                \verb|intent(in)|.\\
+              & The value of the parameter to be set. The list of allowed
+                values and the corresponding data types is given in
+                Tables~\ref{tab:p_type}-\ref{tab:p_coarse}.
+                When the value is of type \verb|character(len=*)|,
+                uppercase and lowercase letters can be used indifferently.\\
+\verb|info|   & \verb|integer, intent(out)|.\\
+              & Error code. If no error, 0 is returned. See Section~\ref{sec:errors}
+                for details.\\
+%
+%\verb|ilev|   & \verb|integer, optional, intent(in)|.\\
+%              & For the multilevel preconditioner, the level at which the
+%                preconditioner parameter has to be set.
+%                The levels are numbered in increasing
+%                order starting from the finest one, i.e.\ level 1 is the finest level.
+%                If \verb|ilev| is not present, the parameter identified by \verb|what|
+%                is set at all the appropriate levels (see Table~\ref{tab:params}).
+\end{tabular}
+
+\ \\
+A variety of (one-level and multi-level) preconditioners can be obtained
+by a suitable setting of the preconditioner parameters. These parameters
+can be logically divided into four groups, i.e.\ parameters defining
+\begin{enumerate}
+	\item the type of multi-level preconditioner;
+	\item the one-level preconditioner used as smoother;
+	\item the aggregation algorithm;
+	\item the coarse-space correction at the coarsest level.
+\end{enumerate}
+A list of the parameters that can be set, along with their allowed and
+default values, is given in Tables~\ref{tab:p_type}-\ref{tab:p_coarse}.
+For a better understanding of the meaning of the parameters, the user
+is referred to Section \ref{sec:background}.
+%
+%Note that the routine allows to set different features of the
+%preconditioner at each level through the use of \verb|ilev|.
+%This should be done by users with experience in the field of
+%multi-level preconditioners. Non-expert users are recommended
+%to call \verb| mld_precset| without specifying \verb|ilev|.
+
+
+\begin{sidewaystable}
+\begin{center}
+\begin{tabular}{|l|l|p{2cm}|l|p{7cm}|}
+\hline
+\verb|what|              & \textsc{data type}        &  \verb|val|      &  \textsc{default}  &
+\textsc{comments} \\ \hline
+%\multicolumn{5}{|c|}{\emph{type of the multi-level preconditioner}}\\ \hline
+\verb|mld_ml_type_|      & \verb|character(len=*)|
+                         & \texttt{'ADD'} \ \ \ \texttt{'MULT'}   
+                         & \texttt{'MULT'}
+                         & basic multi-level framework: additive or multiplicative
+                           among the levels (always additive inside a level)         \\  
+\verb|mld_smoother_type_|& \verb|character(len=*)|
+                         & \texttt{'DIAG'} \ \ \ \texttt{'BJAC'} \ \ \ \texttt{'AS'}
+                         & \texttt{'AS'}
+                         & basic one-level preconditioner (i.e.\ smoother) of the
+                           multi-level preconditioner: diagonal, block Jacobi,
+                           AS \\
+\verb|mld_smoother_pos_| & \verb|character(len=*)|
+                         & \texttt{'PRE'} \ \ \ \texttt{'POST'} \ \ \ \texttt{'TWOSIDE'}
+                         & \texttt{'POST'}
+                         & ``position'' of the smoother: pre-smoother, post-smoother, 
+                           pre- and post-smoother \\
+\hline
+\end{tabular}
+\end{center}
+\caption{Parameters defining the type of multi-level preconditioner.
+\label{tab:p_type}}                       
+\end{sidewaystable}
+                   
+\begin{sidewaystable}
+\begin{center}
+\begin{tabular}{|l|l|p{2.6cm}|l|p{7cm}|}
+\hline
+\verb|what|              & \textsc{data type}        &  \verb|val|      &  \textsc{default}  &
+\textsc{comments} \\ \hline
+%\multicolumn{5}{|c|}{\emph{basic one-level preconditioner (smoother)}} \\ \hline
+\verb|mld_sub_ovr_|       & \verb|integer|
+                         & any integer number $\ge 0$
+                         & 1
+                         & Number of overlap layers. \\
+\verb|mld_sub_restr_|    & \verb|character(len=*)|
+                         & \texttt{'HALO'} \ \ \ \ \ \texttt{'NONE'}
+                         & \texttt{'HALO'}
+                         & Type of restriction operator:
+                           \texttt{'HALO'} for taking into account the overlap, \texttt{'NONE'} 
+                           for neglecting it. \\
+\verb|mld_sub_prol_|     & \verb|character(len=*)|
+                         & \texttt{'SUM'} \ \ \ \ \ \texttt{'NONE'}
+                         & \texttt{'NONE'}
+                         & Type of prolongator operator:
+                           \texttt{'SUM'} for adding the contributions from the overlap, \texttt{'NONE'}
+                           for neglecting them.   \\
+\verb|mld_sub_solve_|    & \verb|character(len=*)|
+                         & \texttt{'ILU'} \ \ \ \ \ \texttt{'MILU'} \ \ \ \ \ \texttt{'ILUT'} \ \ \ \ \
+                           \texttt{'UMF'} \ \ \ \ \ \texttt{'SLU'}
+                         & \texttt{'UMF'}
+                         & Local solver: ILU($p$), MILU($p$), ILU($p,t$), LU from UMFPACK, LU from SuperLU,
+                           plus triangular solve. \\    
+\verb|mld_sub_fillin_|   & \verb|integer|
+                         & any integer number $\ge 0$
+                         & 0
+                         & Fill-in level $p$ of the incomplete LU factorizations. \\
+\verb|mld_sub_thresh_|   & \verb|real(|\emph{kind\_parameter}\verb|)|
+                         & any real number $\ge 0$
+                         & \texttt{0.e0} (or \texttt{0.d0})
+                         & Drop tolerance $t$ in the ILU($p,t$) factorization. \\
+\verb|mld_sub_ren_|      & \verb|character(len=*)|
+                         & \texttt{'RENUM\_NONE'}, \texttt{'RENUM\_GLOBAL'} %, \texttt{'RENUM_GPS'}
+                         & \texttt{'RENUM\_NONE'}
+                         & Row and column reordering of the local submatrices: no reordering,
+                           reordering according to the global numbering of the rows and columns of
+                           the whole matrix. \\
+\hline
+\end{tabular}
+\end{center}
+\caption{Parameters defining the one-level preconditioner used as smoother.
+\label{tab:p_smoother}}  
+\end{sidewaystable}                     
+                   
+\begin{sidewaystable}
+\begin{center}
+\begin{tabular}{|l|l|p{2.6cm}|l|p{7cm}|}
+\hline
+\verb|what|              & \textsc{data type}        &  \verb|val|      &  \textsc{default}  &
+\textsc{comments} \\ \hline
+%\multicolumn{5}{|c|}{\emph{aggregation algorithm}} \\ \hline
+\verb|mld_aggr_alg_|     & \verb|character(len=*)|
+                         & \texttt{'DEC'}
+                         & \texttt{'DEC'}
+                         & Aggregation algorithm. Currently, only the decoupled aggregation is available. \\
+\verb|mld_aggr_kind_|    & \verb|character(len=*)|
+                         & \texttt{'SMOOTH'} \ \ \ \ \ \texttt{'RAW'}
+                         & \texttt{'SMOOTH'}
+                         & Type of aggregation: smoothed or raw, i.e.\ using the tentative prolongator. \\
+\verb|mld_aggr_thresh_|  & \verb|real(|\emph{kind\_parameter}\verb|)|
+                         & any real number $\in [0, 1]$
+                         & \texttt{0.e0} (or \texttt{0.d0})
+                         & The threshold $\theta$ in the aggregation algorithm. \\
+\verb|mld_aggr_eig_|     & \verb|character(len=*)|
+                         & \texttt{'A\_NORMI'}
+                         & \texttt{'A\_NORMI'}
+                         & Estimate of the maximum eigenvalue of $D^{-1}A$
+                           for the smoothed aggregation. Currently, only the infinity norm of
+                           the matrix is available. \\
+\hline
+\end{tabular}
+\end{center}
+\caption{Parameters defining the aggregation algorithm.
+\label{tab:p_aggregation}} 
+\end{sidewaystable}
+                     
+\begin{sidewaystable}
+\begin{center}
+\begin{tabular}{|l|l|p{2.6cm}|l|p{7cm}|}
+\hline
+\verb|what|              & \textsc{data type}        &  \verb|val|      &  \textsc{default}  &
+\textsc{comments} \\ \hline
+%\multicolumn{5}{|c|}{\emph{coarse-space correction at the coarsest level}}\\ \hline
+\verb|mld_coarse_mat_|   & \verb|character(len=*)|
+                         & \texttt{'DISTR'} \ \ \ \ \ \texttt{'REPL'}
+                         & \texttt{'DISTR'}
+                         & Coarsest matrix: distributed among the processors or replicated on each of them.    \\
+\verb|mld_coarse_solve_| & \verb|character(len=*)|
+                         & \texttt{'BJAC'} \ \ \ \ \ \texttt{'UMF'} \ \ \ \ \ \ \ \ \texttt{'SLU'}
+                           \ \ \ \ \ \texttt{'SLUDIST'}
+                         & \texttt{'BJAC'}
+                         & Solver used at the coarsest level: block Jacobi, sequential LU from UMFPACK,
+                           sequential LU from SuperLU, or distributed LU from SuperLU\_Dist.
+                           If the coarsest matrix is distributed, only \texttt{'BJAC'} and \texttt{'SLUDIST'}
+                           can be chosen; if it is replicated, only \emph{'BJAC'} or \texttt{'SLUDIST'} can
+                           be selected. \\
+\verb|mld_coarse_subsolve_| & \verb|character(len=*)|
+                         & \texttt{'ILU'} \ \ \ \ \ \ \ \texttt{'MILU'} \ \ \ \ \ \ \ \ \texttt{'ILUT'}
+                           \ \ \ \ \ \ \ \texttt{'UMF'} \ \ \ \ \ \ \ \texttt{'SLU'}
+                         & \texttt{'UMF'}
+                         & Solver for the diagonal blocks of the coarse matrix, in case the block Jacobi solver
+                           is chosen as coarsest-level solver: ILU($p$), MILU($p$), ILU($p,t$), LU from UMFPACK,
+                           LU from SuperLU, plus triangular solve. \\
+\verb|mld_coarse_sweeps_|& \verb|integer|                         
+                         & any integer number $> 0$
+                         & \texttt{4}
+                         & Number of Block-Jacobi sweeps when 'BJAC' is used as coarsest-level solver. \\
+\verb|mld_coarse_fillin_| & \verb|integer|
+                         & any integer number $\ge 0$
+                         & \texttt{0}
+                         & Fill-in level $p$ of the incomplete LU factorizations. \\
+\verb|mld_coarse_thresh_| & \verb|real(|\emph{kind\_parameter}\verb|)|
+                         & any real number $\ge 0$
+                         & \texttt{0.d0} (or \texttt{0.e0})
+                         & Drop tolerance $t$ in the ILU($p,t$) factorization. \\
+\hline
+\end{tabular}
+\end{center}
+\caption{Parameters defining the coarse-space correction at the coarsest
+level.\label{tab:p_coarse}} 
+\end{sidewaystable}
+
+
+\clearpage
+\subsection{Subroutine mld\_precbld\label{sec:precbld}}
+
+\begin{center}
+\verb|mld_precbld(a,desc_a,p,info)|\\
+\end{center}
+
+\noindent
+This routine builds the preconditioner according to the requirements made by
+the user through the routines \verb|mld_precinit| and \verb|mld_precset|.
+
+\subsubsection*{Arguments}
+
+\begin{tabular}{p{1.2cm}p{10.5cm}}
+\verb|a|      & \verb|type(psb_|\emph{x}\verb|spmat_type), intent(in)|. \\
+              & The sparse matrix structure containing the local part of the
+                matrix to be preconditioned. Note that \emph{x} must be chosen according
+                to the real/complex, single/double precision version of MLD2P4 under use.
+                See the PSBLAS User's Guide for details \cite{PSBLASGUIDE}.\\
+\verb|desc_a| & \verb|type(psb_desc_type), intent(in)|. \\
+              & The communication descriptor of \verb|a|. See the PSBLAS User's Guide for
+                details \cite{PSBLASGUIDE}.\\
+\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(inout)|.\\
+              & The preconditioner data structure. Note that \emph{x} must be chosen according
+                to the real/complex, single/double precision version of MLD2P4 under use.\\
+\verb|info|   & \verb|integer, intent(out)|.\\
+              & Error code. If no error, 0 is returned. See Section~\ref{sec:errors} for details.\\
+\end{tabular}
+
+\clearpage
+\subsection{Subroutine mld\_precaply\label{sec:precaply}}
+
+\begin{center}
+\verb|mld_precaply(p,x,y,desc_a,info)|\\
+\verb|mld_precaply(p,x,y,desc_a,info,trans,work)|\\
+\end{center}
+
+\noindent
+This routine computes $y = op(M^{-1})\, x$, where $M$ is a previously built
+preconditioner, stored into \verb|p|, and $op$
+denotes the preconditioner itself or its transpose, according to
+the value of \verb|trans|.
+Note that, when MLD2P4 is used with a Krylov solver from PSBLAS,
+\verb|mld_precaply| is called within the PSBLAS routine \verb|mld_krylov|
+and hence it is completely transparent to the user.
+
+\subsubsection*{Arguments}
+
+\begin{tabular}{p{1.2cm}p{10.5cm}}
+\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(inout)|.\\
+              & The preconditioner data structure, containing the local part of $M$.
+                Note that \emph{x} must be chosen according
+                to the real/complex, single/double precision version of MLD2P4 under use.\\
+\verb|x|      & \emph{type}\verb|(|\emph{kind\_parameter}\verb|), dimension(:), intent(in)|.\\
+              & The local part of the vector $x$. Note that \emph{type} and   
+                \emph{kind\_parameter} must be chosen according
+                to the real/complex, single/double precision version of MLD2P4 under use.\\
+\verb|y|      & \emph{type}\verb|(|\emph{kind\_parameter}\verb|), dimension(:), intent(out)|.\\
+              & The local part of the vector $y$. Note that \emph{type} and
+                \emph{kind\_parameter} must be chosen according
+                to the real/complex, single/double precision version of MLD2P4 under use.\\
+\verb|desc_a| & \verb|type(psb_desc_type), intent(in)|. \\
+              & The communication descriptor associated to the matrix to be
+                preconditioned.\\
+\verb|info|   & \verb|integer, intent(out)|.\\
+              & Error code. If no error, 0 is returned. See Section~\ref{sec:errors} for details.\\
+\verb|trans|  & \verb|character(len=1), optional, intent(in).|\\
+              & If \verb|trans| = \verb|'N','n'| then $op(M^{-1}) = M^{-1}$;
+                if \verb|trans| = \verb|'T','t'| then $op(M^{-1}) = M^{-T}$
+                (transpose of $M^{-1})$;  if \verb|trans| = \verb|'C','c'| then $op(M^{-1}) = M^{-C}$
+                (conjugate transpose of $M^{-1})$.\\
+\verb|work|  & \emph{type}\verb|(|\emph{kind\_parameter}\verb|), dimension(:), optional, target|.\\
+             & Workspace. Its size should be at
+               least \verb|4 * psb_cd_get_local_cols(desc_a)| (see the PSBLAS User's Guide).
+               Note that \emph{type} and \emph{kind\_parameter} must be chosen according
+               to the real/complex, single/double precision version of MLD2P4 under use.\\
+\end{tabular}
+
+
+\subsection{Subroutine mld\_precfree\label{sec:precfree}}
+
+\begin{center}
+\verb|mld_precfree(p,info)|\\
+\end{center}
+
+\noindent
+This routine deallocates the preconditioner data structure.
+
+\subsubsection*{Arguments}
+
+\begin{tabular}{p{1.2cm}p{10.5cm}}
+\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(inout)|.\\
+              & The preconditioner data structure. Note that \emph{x} must be chosen according
+                to the real/complex, single/double precision version of MLD2P4 under use.\\
+\verb|info|   & \verb|integer, intent(out)|.\\
+              & Error code. If no error, 0 is returned. See Section~\ref{sec:errors} for details.\\
+\end{tabular}
+
+
+\subsection{Subroutine mld\_precdescr\label{sec:precdescr}}
+
+\begin{center}
+\verb|mld_precdescr(p,info,iout)|\\
+\end{center}
+
+\noindent
+This routine prints a description of the preconditioner
+to the standard output or to a file.
+
+\subsubsection*{Arguments}
+
+\begin{tabular}{p{1.2cm}p{10.6cm}}
+\verb|p|      & \verb|type(mld_|\emph{x}\verb|prec_type), intent(in)|.\\
+              & The preconditioner data structure. Note that \emph{x} must be chosen according
+                to the real/complex, single/double precision version of MLD2P4 under use.\\
+\verb|info|   & \verb|integer, intent(out)|.\\
+              & Error code. If no error, 0 is returned. See Section~\ref{sec:errors} for details.\\
+\verb|iout|   & \verb|integer, intent(in), optional|.\\
+              & The id of the file where the preconditioner description
+                will be printed; the default is the standard output.\\
+\end{tabular}
+
+%%% Local Variables: 
+%%% mode: latex
+%%% TeX-master: "userguide"
+%%% End: 
diff --git a/mlprec/mld_cprecinit.f90 b/mlprec/mld_cprecinit.f90
index ea1c9585..deb22ce7 100644
--- a/mlprec/mld_cprecinit.f90
+++ b/mlprec/mld_cprecinit.f90
@@ -117,12 +117,12 @@ subroutine mld_cprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:) = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_noprec_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_f_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_noprec_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_f_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
@@ -133,12 +133,12 @@ subroutine mld_cprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)              = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_diag_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_f_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0 
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_diag_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_f_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0 
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
@@ -149,14 +149,14 @@ subroutine mld_cprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)              = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
   case ('AS')
@@ -166,15 +166,15 @@ subroutine mld_cprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)               = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)      = mld_as_ 
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)      = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)      = psb_halo_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)       = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)        = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)          = 1
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)    = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_)  = 1
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_as_ 
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_halo_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 1
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
 
   case ('ML')
@@ -189,61 +189,62 @@ subroutine mld_cprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)                   = 0
-    p%baseprecv(ilev_)%rprcparm(:)                   = szero
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)      = mld_as_ 
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)      = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)      = psb_halo_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)       = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)        = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)          = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)    = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_)  = 1
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%rprcparm(:)                    = szero
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_as_ 
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_halo_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
     if (nlev_ == 1) return 
 
     do ilev_ = 2, nlev_ -1 
       if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
       if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
       if (info /= 0) return
-      p%baseprecv(ilev_)%iprcparm(:)                  = 0
-      p%baseprecv(ilev_)%rprcparm(:)                  = szero
-      p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+      p%baseprecv(ilev_)%iprcparm(:)                    = 0
+      p%baseprecv(ilev_)%rprcparm(:)                    = szero
+      p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
       p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-      p%baseprecv(ilev_)%iprcparm(mld_ml_type_)       = mld_mult_ml_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)      = mld_dec_aggr_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)     = mld_smooth_prol_
-      p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)    = mld_distr_mat_
+      p%baseprecv(ilev_)%iprcparm(mld_ml_type_)         = mld_mult_ml_
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)        = mld_dec_aggr_
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)       = mld_smooth_prol_
+      p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)      = mld_distr_mat_
       p%baseprecv(ilev_)%iprcparm(mld_smoother_pos_)    = mld_post_smooth_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)      = mld_max_norm_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)        = mld_max_norm_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
       p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
-      p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)     = 4.e0/3.e0         
+      p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)       = 4.e0/3.e0         
     end do
     ilev_ = nlev_
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)                  = 0
-    p%baseprecv(ilev_)%rprcparm(:)                  = szero
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%rprcparm(:)                    = szero
+    p%baseprecv(ilev_)%iprcparm(mld_coarse_solve_)    = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)      = mld_distr_mat_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-    p%baseprecv(ilev_)%iprcparm(mld_ml_type_)       = mld_mult_ml_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)      = mld_dec_aggr_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)     = mld_smooth_prol_
-    p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)    = mld_distr_mat_
+    p%baseprecv(ilev_)%iprcparm(mld_ml_type_)         = mld_mult_ml_
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)        = mld_dec_aggr_
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)       = mld_smooth_prol_
     p%baseprecv(ilev_)%iprcparm(mld_smoother_pos_)    = mld_post_smooth_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)      = mld_max_norm_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)        = mld_max_norm_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 4
-    p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)     = 4.e0/3.e0         
+    p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)       = 4.e0/3.e0         
 
   case default
     write(0,*) name,': Warning: Unknown preconditioner type request "',ptype,'"'
diff --git a/mlprec/mld_cprecset.f90 b/mlprec/mld_cprecset.f90
index c3b1f442..9ac8c29e 100644
--- a/mlprec/mld_cprecset.f90
+++ b/mlprec/mld_cprecset.f90
@@ -244,6 +244,8 @@ subroutine mld_cprecseti(p,what,val,info,ilev)
         select case (val) 
         case(mld_umf_, mld_slu_)
           p%baseprecv(nlev_)%iprcparm(mld_coarse_mat_)  = mld_repl_mat_
+        case(mld_bjac_)
+          p%baseprecv(nlev_)%iprcparm(mld_sub_solve_)   = mld_ilu_n_
         case default
           p%baseprecv(nlev_)%iprcparm(mld_coarse_mat_)  = mld_distr_mat_
         end select
diff --git a/mlprec/mld_dprecinit.f90 b/mlprec/mld_dprecinit.f90
index c4bdb2c1..a7dfdf2d 100644
--- a/mlprec/mld_dprecinit.f90
+++ b/mlprec/mld_dprecinit.f90
@@ -117,12 +117,12 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:) = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_noprec_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_f_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_noprec_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_f_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
@@ -133,12 +133,12 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)              = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_diag_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_f_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0 
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_diag_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_f_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0 
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
@@ -149,14 +149,14 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)              = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
   case ('AS')
@@ -166,15 +166,15 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)               = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)      = mld_as_ 
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)      = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)      = psb_halo_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)       = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)        = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)          = 1
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)    = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_)  = 1
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_as_ 
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_halo_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 1
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
 
   case ('ML')
@@ -189,61 +189,62 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)                   = 0
-    p%baseprecv(ilev_)%rprcparm(:)                   = dzero
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)      = mld_as_ 
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)      = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)      = psb_halo_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)       = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)        = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)          = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)    = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_)  = 1
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%rprcparm(:)                    = dzero
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_as_ 
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_halo_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
     if (nlev_ == 1) return 
 
     do ilev_ = 2, nlev_ -1 
       if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
       if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
       if (info /= 0) return
-      p%baseprecv(ilev_)%iprcparm(:)                  = 0
-      p%baseprecv(ilev_)%rprcparm(:)                  = dzero
-      p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+      p%baseprecv(ilev_)%iprcparm(:)                    = 0
+      p%baseprecv(ilev_)%rprcparm(:)                    = dzero
+      p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
       p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-      p%baseprecv(ilev_)%iprcparm(mld_ml_type_)       = mld_mult_ml_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)      = mld_dec_aggr_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)     = mld_smooth_prol_
-      p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)    = mld_distr_mat_
+      p%baseprecv(ilev_)%iprcparm(mld_ml_type_)         = mld_mult_ml_
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)        = mld_dec_aggr_
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)       = mld_smooth_prol_
+      p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)      = mld_distr_mat_
       p%baseprecv(ilev_)%iprcparm(mld_smoother_pos_)    = mld_post_smooth_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)      = mld_max_norm_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)        = mld_max_norm_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
       p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
-      p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)     = 4.d0/3.d0         
+      p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)       = 4.d0/3.d0         
     end do
     ilev_ = nlev_
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)                  = 0
-    p%baseprecv(ilev_)%rprcparm(:)                  = dzero
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%rprcparm(:)                    = dzero
+    p%baseprecv(ilev_)%iprcparm(mld_coarse_solve_)    = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)      = mld_distr_mat_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-    p%baseprecv(ilev_)%iprcparm(mld_ml_type_)       = mld_mult_ml_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)      = mld_dec_aggr_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)     = mld_smooth_prol_
-    p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)    = mld_distr_mat_
+    p%baseprecv(ilev_)%iprcparm(mld_ml_type_)         = mld_mult_ml_
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)        = mld_dec_aggr_
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)       = mld_smooth_prol_
     p%baseprecv(ilev_)%iprcparm(mld_smoother_pos_)    = mld_post_smooth_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)      = mld_max_norm_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)        = mld_max_norm_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 4
-    p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)     = 4.d0/3.d0         
+    p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)       = 4.d0/3.d0         
 
   case default
     write(0,*) name,': Warning: Unknown preconditioner type request "',ptype,'"'
diff --git a/mlprec/mld_dprecset.f90 b/mlprec/mld_dprecset.f90
index 62264096..a2107360 100644
--- a/mlprec/mld_dprecset.f90
+++ b/mlprec/mld_dprecset.f90
@@ -244,6 +244,8 @@ subroutine mld_dprecseti(p,what,val,info,ilev)
         select case (val) 
         case(mld_umf_, mld_slu_)
           p%baseprecv(nlev_)%iprcparm(mld_coarse_mat_)  = mld_repl_mat_
+        case(mld_bjac_)
+          p%baseprecv(nlev_)%iprcparm(mld_sub_solve_)   = mld_ilu_n_
         case default
           p%baseprecv(nlev_)%iprcparm(mld_coarse_mat_)  = mld_distr_mat_
         end select
diff --git a/mlprec/mld_sprecinit.f90 b/mlprec/mld_sprecinit.f90
index 38637309..7e5046ca 100644
--- a/mlprec/mld_sprecinit.f90
+++ b/mlprec/mld_sprecinit.f90
@@ -117,12 +117,12 @@ subroutine mld_sprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:) = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_noprec_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_f_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_noprec_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_f_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
@@ -133,12 +133,12 @@ subroutine mld_sprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)              = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_diag_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_f_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0 
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_diag_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_f_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0 
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
@@ -149,14 +149,14 @@ subroutine mld_sprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)              = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
   case ('AS')
@@ -166,15 +166,15 @@ subroutine mld_sprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)               = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)      = mld_as_ 
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)      = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)      = psb_halo_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)       = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)        = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)          = 1
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)    = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_)  = 1
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_as_ 
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_halo_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 1
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
 
   case ('ML')
@@ -189,61 +189,62 @@ subroutine mld_sprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)                   = 0
-    p%baseprecv(ilev_)%rprcparm(:)                   = szero
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)      = mld_as_ 
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)      = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)      = psb_halo_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)       = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)        = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)          = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)    = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_)  = 1
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%rprcparm(:)                    = szero
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_as_ 
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_halo_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
     if (nlev_ == 1) return 
 
     do ilev_ = 2, nlev_ -1 
       if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
       if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
       if (info /= 0) return
-      p%baseprecv(ilev_)%iprcparm(:)                  = 0
-      p%baseprecv(ilev_)%rprcparm(:)                  = szero
-      p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+      p%baseprecv(ilev_)%iprcparm(:)                    = 0
+      p%baseprecv(ilev_)%rprcparm(:)                    = szero
+      p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
       p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-      p%baseprecv(ilev_)%iprcparm(mld_ml_type_)       = mld_mult_ml_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)      = mld_dec_aggr_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)     = mld_smooth_prol_
-      p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)    = mld_distr_mat_
+      p%baseprecv(ilev_)%iprcparm(mld_ml_type_)         = mld_mult_ml_
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)        = mld_dec_aggr_
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)       = mld_smooth_prol_
+      p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)      = mld_distr_mat_
       p%baseprecv(ilev_)%iprcparm(mld_smoother_pos_)    = mld_post_smooth_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)      = mld_max_norm_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)        = mld_max_norm_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
       p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
-      p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)     = 4.e0/3.e0         
+      p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)       = 4.e0/3.e0         
     end do
     ilev_ = nlev_
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)                  = 0
-    p%baseprecv(ilev_)%rprcparm(:)                  = szero
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%rprcparm(:)                    = szero
+    p%baseprecv(ilev_)%iprcparm(mld_coarse_solve_)    = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)      = mld_distr_mat_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-    p%baseprecv(ilev_)%iprcparm(mld_ml_type_)       = mld_mult_ml_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)      = mld_dec_aggr_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)     = mld_smooth_prol_
-    p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)    = mld_distr_mat_
+    p%baseprecv(ilev_)%iprcparm(mld_ml_type_)         = mld_mult_ml_
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)        = mld_dec_aggr_
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)       = mld_smooth_prol_
     p%baseprecv(ilev_)%iprcparm(mld_smoother_pos_)    = mld_post_smooth_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)      = mld_max_norm_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)        = mld_max_norm_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 4
-    p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)     = 4.e0/3.e0         
+    p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)       = 4.e0/3.e0         
 
   case default
     write(0,*) name,': Warning: Unknown preconditioner type request "',ptype,'"'
diff --git a/mlprec/mld_sprecset.f90 b/mlprec/mld_sprecset.f90
index bf22fab4..a7dc898a 100644
--- a/mlprec/mld_sprecset.f90
+++ b/mlprec/mld_sprecset.f90
@@ -244,6 +244,8 @@ subroutine mld_sprecseti(p,what,val,info,ilev)
         select case (val) 
         case(mld_umf_, mld_slu_)
           p%baseprecv(nlev_)%iprcparm(mld_coarse_mat_)  = mld_repl_mat_
+        case(mld_bjac_)
+          p%baseprecv(nlev_)%iprcparm(mld_sub_solve_)   = mld_ilu_n_
         case default
           p%baseprecv(nlev_)%iprcparm(mld_coarse_mat_)  = mld_distr_mat_
         end select
diff --git a/mlprec/mld_zprecinit.f90 b/mlprec/mld_zprecinit.f90
index c20a6ae0..29241ce9 100644
--- a/mlprec/mld_zprecinit.f90
+++ b/mlprec/mld_zprecinit.f90
@@ -117,12 +117,12 @@ subroutine mld_zprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:) = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_noprec_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_f_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_noprec_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_f_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
@@ -133,12 +133,12 @@ subroutine mld_zprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)              = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_diag_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_f_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0 
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_diag_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_f_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0 
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
@@ -149,14 +149,14 @@ subroutine mld_zprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)              = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
   case ('AS')
@@ -166,15 +166,15 @@ subroutine mld_zprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)               = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)      = mld_as_ 
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)      = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)      = psb_halo_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)       = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)        = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)          = 1
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)    = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_)  = 1
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_as_ 
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_halo_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 1
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
 
 
   case ('ML')
@@ -189,61 +189,62 @@ subroutine mld_zprecinit(p,ptype,info,nlev)
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)                   = 0
-    p%baseprecv(ilev_)%rprcparm(:)                   = dzero
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)      = mld_as_ 
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)      = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)      = psb_halo_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)       = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)        = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)          = 0
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)    = 0
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_)  = 1
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%rprcparm(:)                    = dzero
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_as_ 
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_halo_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
     if (nlev_ == 1) return 
 
     do ilev_ = 2, nlev_ -1 
       if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
       if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
       if (info /= 0) return
-      p%baseprecv(ilev_)%iprcparm(:)                  = 0
-      p%baseprecv(ilev_)%rprcparm(:)                  = dzero
-      p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+      p%baseprecv(ilev_)%iprcparm(:)                    = 0
+      p%baseprecv(ilev_)%rprcparm(:)                    = dzero
+      p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
       p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-      p%baseprecv(ilev_)%iprcparm(mld_ml_type_)       = mld_mult_ml_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)      = mld_dec_aggr_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)     = mld_smooth_prol_
-      p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)    = mld_distr_mat_
+      p%baseprecv(ilev_)%iprcparm(mld_ml_type_)         = mld_mult_ml_
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)        = mld_dec_aggr_
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)       = mld_smooth_prol_
+      p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)      = mld_distr_mat_
       p%baseprecv(ilev_)%iprcparm(mld_smoother_pos_)    = mld_post_smooth_
-      p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)      = mld_max_norm_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-      p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+      p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)        = mld_max_norm_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+      p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
       p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1
-      p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)     = 4.d0/3.d0         
+      p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)       = 4.d0/3.d0         
     end do
     ilev_ = nlev_
     if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
     if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info)
     if (info /= 0) return
-    p%baseprecv(ilev_)%iprcparm(:)                  = 0
-    p%baseprecv(ilev_)%rprcparm(:)                  = dzero
-    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)     = mld_bjac_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)     = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)      = psb_none_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)       = 0
+    p%baseprecv(ilev_)%iprcparm(:)                    = 0
+    p%baseprecv(ilev_)%rprcparm(:)                    = dzero
+    p%baseprecv(ilev_)%iprcparm(mld_coarse_solve_)    = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)       = mld_ilu_n_
+    p%baseprecv(ilev_)%iprcparm(mld_smoother_type_)   = mld_bjac_
+    p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)      = mld_distr_mat_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_restr_)       = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_prol_)        = psb_none_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_ren_)         = 0
     p%baseprecv(ilev_)%iprcparm(mld_sub_ovr_)         = 0
-    p%baseprecv(ilev_)%iprcparm(mld_ml_type_)       = mld_mult_ml_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)      = mld_dec_aggr_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)     = mld_smooth_prol_
-    p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_)    = mld_distr_mat_
+    p%baseprecv(ilev_)%iprcparm(mld_ml_type_)         = mld_mult_ml_
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_)        = mld_dec_aggr_
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_)       = mld_smooth_prol_
     p%baseprecv(ilev_)%iprcparm(mld_smoother_pos_)    = mld_post_smooth_
-    p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)      = mld_max_norm_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_solve_)     = mld_ilu_n_
-    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)   = 0
+    p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_)        = mld_max_norm_
+    p%baseprecv(ilev_)%iprcparm(mld_sub_fillin_)      = 0
     p%baseprecv(ilev_)%iprcparm(mld_smoother_sweeps_) = 4
-    p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)     = 4.d0/3.d0         
+    p%baseprecv(ilev_)%rprcparm(mld_aggr_damp_)       = 4.d0/3.d0         
 
   case default
     write(0,*) name,': Warning: Unknown preconditioner type request "',ptype,'"'
diff --git a/mlprec/mld_zprecset.f90 b/mlprec/mld_zprecset.f90
index 01fc511b..2d58862d 100644
--- a/mlprec/mld_zprecset.f90
+++ b/mlprec/mld_zprecset.f90
@@ -244,6 +244,8 @@ subroutine mld_zprecseti(p,what,val,info,ilev)
         select case (val) 
         case(mld_umf_, mld_slu_)
           p%baseprecv(nlev_)%iprcparm(mld_coarse_mat_)  = mld_repl_mat_
+        case(mld_bjac_)
+          p%baseprecv(nlev_)%iprcparm(mld_sub_solve_)   = mld_ilu_n_
         case default
           p%baseprecv(nlev_)%iprcparm(mld_coarse_mat_)  = mld_distr_mat_
         end select