\section{Manipulating and Solving Sparse Systems}
\label{sec:main}
This section is divided into two parts.  First, we describe how
to set up the matrix and preconditioner for parallel solution.
Second, we describe how to solve the linear systems after this setup
has taken place.

\subsection{Manipulation and Setup}
The data structure that \BS\ uses to represent
sparse matrices is {\tt typedef BSpar\_mat}.\findex{typedef BSpar_mat}
The routine that converts the user's data structure {\tt BSspmat}
to {\tt BSpar\_mat} is {\tt BSmain\_perm()}.
This routine is called as follows.
\begin{verbatim}
  BSpar_mat *BSmain_perm(BSprocinfo *procinfo,BSspmat *A);
\end{verbatim}
This routine colors and permutes the sparse matrix to create a new version
of the sparse matrix appropriate for parallel computation.  
A large number of options are available 
that are set through the context variable {\tt procinfo},
as described in \S\ref{sec:context}.
The user's sparse matrix is not changed permanently by this
routine, but may be manipulated and restored during execution.

\subsubsection{Symmetric and Nonsymmetric Matrices}
\BS\ has {\em two} versions of its internal data structure based
on whether
(1) the matrix nonzeros are symmetric \sindex{symmetric matrix}
and an incomplete Cholesky\sindex{incomplete Cholesky}
factor is to be computed, 
or (2) the matrix nonzeros are to be treated as if they
are nonsymmetric\sindex{nonsymmetric matrix}
(even if the matrix values are actually symmetric) and
an incomplete ILU\sindex{incomplete LU} factor is to be computed.

The routine 
{\tt BSset\_mat\_icc\_storage()}\findex{BSset_mat_icc_storage()}
must be called {\em before} {\tt BSmain\_perm()} to inform 
\BS\ whether internal incomplete 
Cholesky storage or the incomplete LU storage is to be used.
The calling sequence is
\findex{BSset_mat_icc_storage()}
\begin{verbatim}
   void BSset_mat_icc_storage(BSspmat *A,int storage);
\end{verbatim}
where the incomplete Cholesky storage will be used
if {\tt storage} is \true, and
the incomplete LU otherwise.
Note that if the matrix is
nonsymmetric and yet the incomplete Cholesky storage
is specified, incorrect results will be obtained.
\BS\ assumes that the matrix is symmetric if 
the incomplete Cholesky option is used to minimize
the required storage.

The routine {\tt BSset\_mat\_symmetric()} must be called to
inform \BS\ whether the matrix is symmetric (even if
you use the nonsymmetric ILU internal storage format).
The arguments to this routine are
\findex{BSset_mat_symmetric()}
\begin{verbatim}
   void  BSset_mat_symmetric(BSspmat *A,int sym);
\end{verbatim}
where {\tt sym} is \true\ if the matrix is symmetric
and \false\ if the matrix is nonsymmetric.

In Table~\ref{tab:sym}
we summarize the possible combination of options
for the two functions\linebreak {\tt BSset\_mat\_icc\_storage()}
and {\tt BSset\_mat\_symmetric()}.
The tricky choice is whether to use ILU for a 
symmetric matrix.
If the system is very indefinite, you may be forced to
this approach to be able to compute an incomplete factor
without breakdown.
However, remember that this approach requires essentially 
twice the storage of using the incomplete Cholesky option.
\begin{table} [htb]
\caption{\label{tab:sym}Consequences of choosing incomplete
Cholesky or incomplete LU for symmetric or nonsymmetric matrices}
\begin{center}
{\renewcommand{\arraystretch}{1.2}
\begin{tabular}{|c|c|c|} \hline
{\tt BSset\_mat\_symmetric()} & {\tt BSset\_mat\_icc\_storage()} & Consequences\\ \hline\hline
                            &                  & Positive definite \\
\true                       &  \true           & preconditioner  \\ \hline
                            &                  & ILU computed, twice the storage \\
\true                       &  \false           & of incomplete Cholesky  \\ \hline
                            &                  & Incorrect \\
\false                      &  \true           & results  \\ \hline
                            &                  & ILU computed, only\\
\false                      &  \false           & correct choice \\ \hline
\end{tabular}
}
%\begin{tabular}{c|c|c|} \cline{2-3}
%          & Incomplete & Incomplete  \\ 
%          & Cholesky   & LU          \\ \cline{2-3} \cline{2-3} \hline
%\multicolumn{1}{|c|}{Symmetric}    & Positive definite  & Twice the storage of \\
%\multicolumn{1}{|c|}{matrix}       & preconditioner     & incomplete Cholesky \\ \hline
%\multicolumn{1}{|c|}{Nonsymmetric} & Incorrect          & Only \\
%\multicolumn{1}{|c|}{matrix}       & results            & choice  \\\hline
%\end{tabular}
\end{center}
\end{table}

\subsubsection{Fast Reorderings for Different Matrix Nonzero Values}
Often you may wish to solve more than one linear system
with the same nonzero structure but different nonzero values.
In this case the context should be set to ``retain'' important
intermediate data structures by using the function
{\tt BSctx\_set\_rt()}\findex{BSctx_set_rt()} described 
in \S\ref{sec:context}.
If {\tt BSmain\_perm()} has already been called with the ``retain''
parameter set to true, you can call {\tt BSmain\_reperm()} to
permute a matrix with the same structure. 
The routine {\tt BSmain\_reperm()}\findex{BSmain_reperm()}
is called with the following arguments.
\begin{verbatim}
  void BSmain_reperm(BSprocinfo *procinfo,BSspmat *A,BSpar_mat *pA);
\end{verbatim}

\subsubsection{Diagonal Scaling of the Linear System}
After calling {\tt BSmain\_perm()}, the matrix then can be scaled 
diagonally by calling {\tt BSscale\_diag()}.\findex{BSscale_diag()}
The function arguments are the following.
\begin{verbatim}
  void BSscale_diag(BSpar_mat *A,FLOAT *sc_diag,BSprocinfo *procinfo);
\end{verbatim}
Once the matrix has been scaled, \BS\ automatically
solves the scaled system
\begin{equation} \nonumber
  (D^{-1}AD^{-1}) (Dx) = (D^{-1}b)~,
\end{equation}
where $D$ is the diagonal matrix whose values are
the square root of the absolute value of the 
vector {\tt sc\_diag} (usually chosen as the diagonal
of $A$).
Thus, if the unscaled matrix has a negative diagonal entry,
it will be $-1$ after scaling.

\subsubsection{The Communication Data Structure {\tt BScomm}}
Prior to either factoring or solving the matrix, the communication 
patterns used by \BS\ must be created.
For the factorization phase, this pattern is compiled by calling
the routine
{\tt BSsetup\_factor()}.\findex{BSsetup_factor()}
\begin{verbatim}
  BScomm *BSsetup_factor(BSpar_mat *A,BSprocinfo *procinfo);
\end{verbatim}
For matrix solution, it is compiled by calling
{\tt BSsetup\_forward()}.\findex{BSsetup_forward()}
\begin{verbatim}
  BScomm *BSsetup_forward(BSpar_mat *A,BSprocinfo *procinfo);
\end{verbatim}
Both routines return a communication pattern of
{\tt typedef BScomm}.\findex{typedef BScomm}
The communication patterns may be freed by
calling {\tt BSfree\_comm()}.\findex{BSfree_comm()}
\begin{verbatim}
  void BSfree_comm(BScomm *comm_ptr);
\end{verbatim}

\subsubsection{Computing an Incomplete Factorization}
If an incomplete factor is to be computed, a copy of 
the matrix must be made by using the routine 
{\tt BScopy\_par\_mat()}.\findex{BScopy_par_mat()}
\begin{verbatim}
  BSpar_mat *BScopy_par_mat(BSpar_mat *A);
\end{verbatim}
The copy is necessary because the iterative solver will
require both the initial matrix and the factorization.
As an aside, we note that (for the incomplete Cholesky storage
scheme) the copy of the sparse matrix shares the clique storage 
space with the matrix from which it is copied.
For more information see the {\tt man} page on
{\tt BScopy\_par\_mat()}.\findex{BScopy_par_mat()}

The incomplete factorization is computed by calling
the routine {\tt BSfactor()}.\findex{BSfactor()}
\begin{verbatim}
  int BSfactor(BSpar_mat *A,BScomm *comm,BSprocinfo *procinfo);
\end{verbatim}
The matrix {\tt A} will be overwritten with the incomplete
factorization.
If the factorization is successful the routine returns
{\tt 0}; otherwise a negative integer is returned whose
absolute value is the row number of the color (less one)
where the failure occurred.
Note that the kind of factorization computed is determined
by the internal representation of the matrix produced by
{\tt BSmain\_perm()}.
For example, consider the following code segment.

{\small
\begin{verbatim}
   BSset_mat_symmetric(A,TRUE);
   BSset_mat_icc_storage(A,FALSE);
   pA = BSmain_perm(procinfo,A); CHKERR(0);
   /* diagonally scale the matrix */
   BSscale_diag(pA,pA->diag,procinfo); CHKERR(0);
   /* set up the communication structure for triangular matrix solution */
   Acomm = BSsetup_forward(pA,procinfo); CHKERR(0);
   /* get a copy of the sparse matrix */
   f_pA = BScopy_par_mat(pA); CHKERR(0);
   /* set up a communication structure for factorization */
   f_comm = BSsetup_factor(f_pA,procinfo); CHKERR(0);
   /* compute the incomplete LU factorization */
   ierr = BSfactor(f_pA,f_comm,procinfo);
\end{verbatim}
}
In this example the call to 
{\tt BSset\_mat\_symmetric()}\findex{BSset_mat_symmetric()}
specifies the matrix as being symmetric.
The internal storage format for the incomplete LU
factorization is chosen prior to the call to 
{\tt BSmain\_perm()} by calling the routine
{\tt BSset\_mat\_icc\_storage()}.\findex{BSset_mat_icc_storage()}
Thus, the call to {\tt BSfactor()} will compute
the incomplete LU factorization rather than
the incomplete Cholesky factorization even though
the matrix is specified as symmetric.
Also, note that the matrix has been diagonally
scaled before the factorization with the call
\findex{BSscale_diag()}
to {\tt BSscale\_diag()}.

\subsubsection{What to Do If the Factorization Fails}
An attempt to compute the incomplete Cholesky factorization 
of a positive definite matrix can fail if a zero
or negative diagonal is encountered during the factorization.
We note that the incomplete factorization has been shown to exist 
only if the matrix is diagonally dominant or in several other
special cases; in general, there is no guarantee that it will exist.
In the case of failure the matrix must be recopied and the 
factorization retried.
We recommend diagonally scaling the matrix and using 
the following loop to accomplish this task.

{\small
\begin{verbatim}
    my_alpha = 1.0;
    /* get a copy of the sparse matrix */
    f_pA = BScopy_par_mat(pA);
    /* factor the matrix until successful */
    while (BSfactor(f_pA,f_comm,procinfo) != 0) {
        /* recopy just the nonzero values */
        BScopy_nz(pA,f_pA);
        /* increment the diagonal shift */
        my_alpha += 0.1;
        BSset_diag(f_pA,my_alpha,procinfo);
    }
\end{verbatim}
}
In this code segment, we are shifting the diagonal
of the matrix by 0.1 every time the factorization fails.
Other strategies are certainly possible and could easily be
implemented by the user.
The routine {\tt BSset\_diag()}\findex{BSset_diag()}
is used to change the entire diagonal to {\tt my\_alpha}.
\begin{verbatim}
  void BSset_diag(BSpar_mat *A,FLOAT my_alpha,BSprocinfo *procinfo);
\end{verbatim}

\subsection{Solving the Linear System}
Once the parallel matrix and the communication structures
have been created, we are ready to solve the sparse linear system.
One of two routines can be called to do this: 
(1) {\tt BSpar\_solve()},\findex{BSpar_solve()} for
either symmetric or nonsymmetric linear systems; and
(2) {\tt BSpar\_isolve()},\findex{BSpar_isolve()} for
symmetric indefinite systems, especially those involved
in ``shift and invert'' strategies.

{\tt BSpar\_solve()} can be used repeatedly to solve systems 
of linear equations with one or with multiple right-hand sides.
The calling arguments are as follows.\findex{BSpar_solve()}
\begin{verbatim}
  int BSpar_solve(BSpar_mat *A, BSpar_mat *fact_A, BScomm *comm_A,
      FLOAT *rhs, FLOAT *x, FLOAT *residual, BSprocinfo *procinfo);
\end{verbatim}
The routine returns the number of iterations taken by
the solver.
If the solver fails (for example, a negative curvature direction
is found in a conjugate gradient iteration), the solver returns
the negative of the iteration number.

The following code segment presents an example of
how {\tt BSpar\_solve()} could be called.

{\small
\begin{verbatim}
    BSctx_set_method(procinfo,GMRES);
    BSctx_set_pre(procinfo,PRE_ILU);
    BSctx_set_max_it(procinfo,50);
    BSctx_set_restart(procinfo,25);
    BSctx_set_guess(procinfo,TRUE);
    BSctx_set_tol(procinfo,1.0e-7);
    num_iter = BSpar_solve(pA,f_pA,Acomm,rhs,x,residual,procinfo); CHKERR();
\end{verbatim}
}
In this example the GMRES iterative method has been specified,
preconditioned by an incomplete LU factorization.
The preconditioner specified must agree with the
preconditioner that has been computed by 
{\tt BSfactor()}.\findex{BSfactor()}
The maximum number of iterations allowed the solver has been
specified as 50, and the number of storage vector allowed GMRES
before restart as 25.
In this example {\tt BSctx\_set\_guess()}\findex{BSctx_set_guess()}
has specified that the solver use an initial guess for {\tt x}
as the zero vector.
Finally, the relative residual tolerance, $\| Ax - b\| / \| b \|$,
has been set to {\tt 1.0e-7}.
Additional details on the arguments used can be found in 
the {\tt man} page.

A second solver,
{\tt BSpar\_isolve()},\findex{BSpar_isolve()} is set up to solve
the shifted system $(A-\sigma B) x = b$, where $A$ and $B$ are
symmetric matrices, $\sigma$ is a real constant, $x$ is the 
solution vector, and $b$ is the right-hand side.
\BS\ is set up to take advantage of $B$ being NULL or $\sigma$ being zero.
The calling sequence is the following.

{\small
\begin{verbatim}
  int BSpar_isolve(BSpar_mat *A, BSpar_mat *fact_A, BScomm *comm_A,
      BSpar_mat *B, BScomm *comm_B, FLOAT *in_rhs, FLOAT *out_x,
      FLOAT shift, FLOAT *residual, BSprocinfo *procinfo);
\end{verbatim}
}
{\tt BSpar\_isolve()} uses the SYMMLQ algorithm,
which requires that the preconditioner, if any, be positive definite.
Symmetric diagonal scaling is not possible for an indefinite matrix,
so one of the other preconditioners must be used.  The restriction that
the preconditioner be positive definite is too restrictive for many
problems, but we know of no general-purpose alternative to SYMMLQ that 
takes advantage of symmetry while allowing an indefinite preconditioner.
Another option would be to explicitly form $(A-\sigma B)$
and solve the resulting system with GMRES preconditioned with ILU
(at the expense of twice the memory).
The program examples {\tt grid7} and {\tt grid8}
compare these two options.
See \S\ref{sec:examples} for details.

If you wish simultaneously to solve for more than one right-hand
side, you must call the routine {\tt BSsetup\_block()}
\findex{BSsetup_block()}
to modify the communication structure to 
accommodate the multiple right-hand sides.
The arguments for this function are the following.
\begin{verbatim}
   void BSsetup_block(BSpar_mat *A,BScomm *comm,int block_size,
                      BSprocinfo *procinfo);
\end{verbatim}

\subsection{Freeing Matrices}
To free the parallel matrix created by {\tt BSmain\_perm()},
call the routine {\tt BSfree\_par\_mat()}\findex{BSfree_par_mat()}.
\begin{verbatim}
  void BSfree_par_mat(BSpar_mat *A);
\end{verbatim}
To free a copy of a parallel matrix created by
{\tt BScopy\_par\_mat()}, call the routine\linebreak
{\tt BSfree\_copy\_par\_mat()}\findex{BSfree_copy_par_mat()}.
\begin{verbatim}
  void BSfree_copy_par_mat(BSpar_mat *A);
\end{verbatim}

\newpage
