main.dox

/*!	\mainpage \authors <a href="https://www.cs.ubc.ca/~greif/">Chen Greif</a>, Shiwen He, <a href="https://cs.stanford.edu/people/paulliu/">Paul Liu</a>
*
*	\tableofcontents
*
*	\section intro_sec Overview
*
	\b sym-ildl is a C++ package for solving and producing fast incomplete factorizations of symmetric indefinite or skew-symmetric matrices. Given an \f$n\times n\f$ symmetric indefinite or skew-symmetric matrix \f$\mathbf{A}\f$, this package produces an incomplete \f$\mathbf{LDL^{T}}\f$ factorization. Prior to factorization, \b sym-ildl first scales the matrix to be equilibrated in the max-norm <a href="#refs">[2]</a>, and then preorders the matrix using either the Reverse Cuthill-McKee (RCM) algorithm or the Approximate Minimum Degree algorithm (AMD) <a href="#refs">[1]</a>. To maintain stability, the user can use Bunch-Kaufman or rook partial pivoting during the factorization process. The factorization produced is of the form 
	\f[
		\mathbf{P^{T}SASP \approx LDL^{T}}.
	\f]
	where \f$\mathbf{P}\f$ is a permutation matrix, \f$\mathbf{S}\f$ a scaling matrix, and \f$\mathbf{L}\f$ and \f$\mathbf{D}\f$ are the unit lower triangular and block diagonal factors respectively. The user can also optionally solve the given linear system, using \b sym-ildl's incomplete factorization to precondition the built-in solver or using the full factorization as a direct solver.
	
	This package is based on and extends an incomplete factorization approach proposed by Li and Saad <a href="#refs">[3]</a> (which itself builds on Li, Saad, and Chow <a href="#refs">[4]</a>).
*    
*   \section dlinks Download links
*       <i>[Latest release: 11/15/2015]</i>
*       \subsection matlab_dlinks MATLAB MEX Files
*           - <a href="https://github.com/inutard/matrix-factor/raw/master/release/win64/ildl-2015a.zip">For MATLAB 2015a and above on 64-bit Windows 8.1</a>
*           - <a href="https://github.com/inutard/matrix-factor/raw/master/release/linux/ildl-2014a.zip">For MATLAB 2014a and above on 64-bit OpenSUSE 13.1</a>
*       \subsection binary_dlinks Prebuilt binaries
*           - <a href="https://github.com/inutard/matrix-factor/raw/master/release/win64/ldl_driver.exe">For 64-bit Windows 8.1</a>
*           - <a href="https://github.com/inutard/matrix-factor/raw/master/release/linux/ldl_driver">For 64-bit OpenSUSE 13.2</a>
*       \subsection source_dlinks The source code
*           - <a href="https://github.com/inutard/matrix-factor/archive/master.zip">Latest stable release (gcc 4.7 and above)</a>
*
*
*	\section why_use Why use sym-ildl?
		- <b>It's easy to use</b>: We offer precompiled versions of sym-ildl on both Windows and Linux, as well as for use within MATLAB. If you want to use sym-ildl as a C++ library, there is no need to compile anything. You can just use the header files right away. sym-ildl is defined entirely in the header files.
		- <b>It's reliable</b>: sym-ildl includes several methods to improve stability of the preconditioner. These include different kinds of equilibration, reordering, as well partial pivoting. 
		- <b>It's not just a preconditioner</b>: sym-ildl implements a direct solver as well as several iterative solvers that are integrated with the preconditioner, allowing you to easily and quickly solve your linear system.
*	
*
*	\section quick_start Quick Start
*
	To begin using the package, first download the files above or compile the code hosted at <a href="https://github.com/inutard/matrix-factor">https://github.com/inutard/matrix-factor</a>. The GitHub repository contains the most up to date source code as well as an "experimental" branch that we release new features to. The package works under most Unix distributions as well as Cygwin under Windows. The default compiler used is \c g++, simply type \c make at the command line to compile the entire package. In addition to \subpage ldl_driver "usage as a standalone program", the package also has a \subpage matlab_mex "Matlab interface".
	
	For using sym-ildl as a software library, see the API <a href="http://www.cs.ubc.ca/~inutard/html/classsolver.html">here</a>. sym-ildl is a header only library, so one only needs to include \c source/solver.h for everything to work.

	\subsection ldl_driver Using the package as a standalone program
	The compiled program \c ldl_driver takes in through the command line one parameter (the matrix) as well as several additional optional parameters.
	
	The format of execution is: 
	\code 
		./ldl_driver -filename=[matrix-name.mtx] [-rhs_file=rhs-name.mtx ...]
	\endcode
	
	The parameters can be given in any order, and will use a default value when not specified.
	
	A description of each of these parameters can be accessed by typing 
	\code 
		./ldl_driver --help 
	\endcode
	
	Generally speaking, the code operates in two phases: generating the preconditioner (factorization) and solving the linear system (solver). The factorization parameters can be used to specify how the LDL' preconditioner is built. If a right hand side is specified, the built-in solver attempts to solve the given linear system. The solver takes the preconditioner and uses it in either preconditioned MINRES or a direct solve (full factorization with Bunch-Kaufman or rook pivoting).
    
    For convenience, the parameters are listed below. 
    
    \subsubsection fact_param Factorization parameters
    
	\param filename The filename of the matrix to be loaded. Several test matrices exist in the test_matrices folder. All matrices loaded are required to be in matrix market (.mtx) form.
	
	\param fill Controls memory usage. Each column of \f$\mathbf{L}\f$ is guaranteed to have fewer than \f$fill\cdot nnz(\mathbf{A})/n\f$ elements. Each column of \f$\mathbf{D}\f$ has at most 2 elements. When this argument is not given, the default value for \c fill is <c>3.0</c>.
	
	\param tol Controls agressiveness of dropping. In each column k, elements less than \f$tol \cdot \left|\left|\mathbf{L}_{k+1:n,k}\right|\right|_1\f$ are dropped. The default value for \c tol is <c>0.001</c>.
	
	\param pp_tol A parameter to aggressiveness of Bunch-Kaufman pivoting (BKP). When pp_tol is 0, there is no partial pivoting. Values between 0 and 1 vary the number of pivots of BKP makes. When pp_tol is equal to 1, standard BKP is used. The pp_tol parameter is ignored if the pivot parameter is set to 'rook'. See the \b pivot parameters for more details. The default is 'rook'.

	\param reordering Determines what sort of preordering will be used on the matrix. Choices are 'amd', 'rcm', and 'none'. The default is 'amd'.
	
    \param inplace Indicates whether the factorization should be performed in-place, leading to roughly a 33% saving in memory. This memory comes out of extra information used in the solver. If the solver is needed, then \c inplace should not be used. \c y indicates yes, \c n indicates no. The default flag is \c n.	
    
    \param pivot Indicates the pivoting algorithm used. Choices are 'rook' and 'bunch'. If \c rook is used, the \c pp_tol parameter is ignored. The default is 'rook'.
    
    \param save Indicates whether the output matrices should be saved. \c y indicates yes, \c n indicates no. The default flag is \c y. All matrices are saved in matrix market (.mtx) form. The matrices are saved into a folder named \c output_matrices. There are five saved files: <c>outA.mtx, outL.mtx, outD.mtx, outS.mtx</c>, and \c outP.mtx. \c outB.mtx is the matrix \f$\mathbf{B=P^{T}SASP}\f$. The rest of the outputs should be clear from the description above.
    
    \param display Indicates whether the output matrices should be displayed to the command line, used for debugging purposes. \c y indicates yes, \c n indicates no. The default flag is \c n.	For this parameter to appear, sym-ildl must be compiled with \c SYM_ILDL_DEBUG defined.
    
    \subsubsection solver_param Solver parameters
    \param solver This chooses the solver for \b sym-ildl if given a right hand side. The options are 'sqmr', 'minres', 'full', and 'none'. If the 'full' solver is chosen, the full factorization is produced and a straightforward direct solve is done. Setting this parameter to \c y overrides all other solver parameters as well as -fill and -tol (since we will no longer produce an incomplete factorization). The default solver is 'sqmr'.
    
    \param max_iters Number of iterations that the builtin iterative solvers (SQMR or MINRES) can use. The default is \c -1 (i.e. iterative solver is not applied). The output solution is written to <c>output_matrices\outsol.mtx</c>.
    
    \param solver_tol Relative tolerance for the builtin iterative solvers. When the iterate x satisfies ||Ax-b||/||b|| < \c solver_tol, the iterative solver is terminated. The default is \c 1e-6.
    
    \param rhs_file The filename of the right hand size we want to solve. All right hand sides loaded are required to be in matrix market (.mtx) form. If no right hand sides are given, only the preconditioner is generated.
    
    \par
	Typically, the \c pivot, \c equil, and \c reordering parameters are best left to the default options.
	
	\par Examples:
	Suppose we wish to factor the \c aug3dcqp matrix stored in <c>test_matrices/aug3dcqp.mtx</c>. Using the parameters described above, the execution of the program may go something like this:	
	\code
		./ldl_driver -filename=test_matrices/aug3dcqp.mtx -fill=3.0 tol=0.001 -save=y
		
		Load succeeded. File test_matrices/aug3dcqp.mtx was loaded.
		A is 35543 by 35543 with 128115 non-zeros.
		  Equilibration:                0.047 seconds.
		  AMD:                          0.047 seconds.
		  Permutation:                  0.047 seconds.
		  Factorization (BK pivoting):  0.109 seconds.
		Total time:     0.250 seconds.
		L is 35543 by 35543 with 162160 non-zeros.
		Saving matrices...
		Save complete.
		Factorization Complete. All output written to /output_matrices directory.
	\endcode	
	The code above factors the \c aug3dcqp.mtx matrix (<c>fill=3.0, tol=0.001</c>) from the \c test_matrices folder and saves the outputs. The time it took to pre-order and equilibrate the matrix (0.047s) as well as the actual factorization (0.109s) are also given.
	
	\par
	For convenience, we may use all optional arguments:
	\code
		./ldl_driver -filename=test_matrices/aug3dcqp.mtx
		
		Load succeeded. File test_matrices/aug3dcqp.mtx was loaded.
		A is 35543 by ...
	\endcode
	The code above would use the default arguments <c>-fill=3.0 -tol=0.001 -reordering=amd -save=y</c>. In general, we may give \c ldl_driver the arguments in any order, and omit any number of them (except the \c filename argument).
	
    \par Solving a linear system
	To actually solve the given linear system, simply supply a right hand side file (with the \c rhs_file argument) and a maximum number of solver iterations. When no right hand side is specified (but a solver iteration is), \b sym-ildl assumes a right hand side of all 1's for debugging purposes:
	\code
		./ldl_driver -filename=test_matrices/aug3dcqp.mtx -max_iters=300 -max_tol=1e-6 -solver=minres
        
        Load succeeded. File test_matrices/aug3dcqp.mtx was loaded.
        A is 35543 by 35543 with 77829 non-zeros.
        Right hand side has 35543 entries.
          Equilibration:                       0.000 seconds.
          AMD:                                  0.015 seconds.
          Permutation:                          0.047 seconds.
          Factorization (Rook pivoting):        0.047 seconds.
        Total time:     0.109 seconds.
        L is 35543 by 35543 with 162160 non-zeros.

        Solving matrix with MINRES...
        The estimated condition number of the matrix is 2.202256e+00.
        MINRES took 18 iterations and got down to relative residual 8.627871e-07.
        Solve time:             0.141 seconds.

        Solution saved to output_matrices/outsol.mtx.
        Saving matrices...
        Save complete.
        Factorization Complete. All output written to /output_matrices directory.
        
	\endcode
    
	\subsection matlab_mex Using sym-ildl within MATLAB
	If everything is compiled correctly, simply open MATLAB in the package directory. The \c startup.m script adds all necessary paths to MATLAB upon initiation. The program can now be called by its function handle, \c ildl.
	
	\c ildl takes in seven arguments, six of them being optional. A full description of the parameters can be displayed by typing 
	\code
		help ildl
	\endcode
	
	For convenience, the parameters are described below:
	\param A The matrix to be factored.
	
	\param fill Controls memory usage. Each column is guaranteed to have fewer than \f$fill\cdot nnz(\mathbf{A})/n\f$ elements. When this argument is not given, the default value for \c fill is <c>3.0</c>.
	
	\param tol Controls aggressiveness of dropping. In each column k, elements less than \f$tol \cdot \left|\left|\mathbf{L}_{k+1:n,k}\right|\right|_1\f$ are dropped. The default value for \c tol is <c>0.001</c>.

	\param reordering Determines what sort of pre-ordering will be used on the matrix. Choices are 'amd', 'rcm', and 'none'. The default is 'amd'.

	\param equil Determines if matrix is to be equilibriated (in the max norm) before anything. This parameter must be one of 'bunch' or 'none'. Default: 'bunch'
	
    \param pivot_type Chooses the pivoting scheme used during the factorization. This parameter must be one of 'rook' or 'bunch'. Tbe default is 'rook'.
    
	\param pp_tol Threshold parameter for Bunch-Kaufman pivoting (BKP). When pp_tol >= 1, full BKP is used. When pp_tol is 0, there is no partial pivoting. As pp_tol increases from 0 to 1, we smoothly switch from no pivoting to full BKP. Low values of pp_tol can be useful as an aggressive pivoting process may damage and permute any special structure present in the input matrix. The default value is 1.0. When rook pivoting is used, this parameter has no effect.
    
	As with the standalone executable, the function has five outputs: <c>L, D, p, S,</c> and \c B:
	\return \b L Unit lower triangular factor of \f$\mathbf{P^{T}SASP}\f$.
	\return \b D Block diagonal factor (consisting of 1x1 and 2x2 blocks) of \f$\mathbf{P^{T}SASP}\f$.
	\return \b p Permutation vector containing permutations done to \f$\mathbf{A}\f$.
	\return \b S Diagonal scaling matrix that equilibrates \f$\mathbf{A}\f$ in the max-norm.
	\return \b B Permuted and scaled matrix \f$\mathbf{B=P^{T}SASP}\f$ after factorization.
*
*	\par Examples:
	Before we begin, let's first generate some symmetric indefinite matrices:
	\code
		>> B = sparse(gallery('uniformdata',100,0));
		>> A = [speye(100) B; B' sparse(100, 100)];
	\endcode
	The \c A generated is a special type of matrix called a saddle-point matrix. These matrices are indefinite and arise often in optimzation problems. Note that A must be a MATLAB \b sparse matrix.
	
	\par
	To factor the matrix, we supply \c ildl with the parameters described above:
	\code
		>> [L, D, p, S, B] = ildl(A, 1.0, 0.001);
		Equilibration:	0.001 seconds.
		AMD:		    0.000 seconds.
		Permutation:	0.001 seconds.
		Factorization:	0.013 seconds.
		...
	\endcode
	As we can see above, \c ildl will supply some timing information to the console when used. The reordering time is the time taken to equilibrate and pre-order the matrix. The factorization time is the time it took to factor and pivot the matrix with partial pivoting.
	
	\par
	We may also take advantage of the optional parameters and simply feed \c ildl only one parameter:
	\code
		>> [L, D, p, S, B] = ildl(A);
		Equilibration:	0.001 seconds.
		AMD:		0.001 seconds.
		...
	\endcode	
	As specified above, the default values of <c>fill=3.0</c>, <c>tol=0.001</c>,  <c>pivot_type='rook'</c>, <c>pp_tol=1.0</c>, and <c>reordering=amd</c> are used.
	
	\subsubsection gmres_ildl Using sym-ildl with an iterative solver in MATLAB
	To use the factorization as a preconditioner, we must apply the permutation and scaling matrices returned by \c ildl. For convenience here is a complete example of how \c ildl can be used with GMRES:
	
	\par
	\code
		% Generate a test matrix
		B = sparse(gallery('uniformdata',100,0));
		A = [speye(100) B; B' sparse(100, 100)];
		
		% Run sym-ildl
		[L, D, p, S, B] = ildl(A, 1.0, 0.001);

		% Create a right hand side of all 1s
		rhs = ones(size(A,1), 1);
		
		% Solve A*x = e, or equivalently, P'SASP*y = P'S*e,
		% where y = P'S^(-1)*x. Since B = P'SASP, we can use give it to GMRES for convenience.
		% Create the new right hand side (scaled and permuted)
		new_rhs = S*rhs;
		new_rhs = new_rhs(p);
		
		% Create our preconditioner. LDLt(x) = (LDL')^(-1) * x.
		LDLt = @(x) L'\(D\(L\x));
		
		% Run GMRES(60), stopping when the relative residual is below 1e-6 or when 3*60 total iterations have passed.
		y = gmres(B, new_rhs, 60, 1e-8, 3, LDLt);
		
		% Get the solution x
		z(p) = y;
		x = S*z';
		
		% Sanity check: Is it the close to the direct solve?
		fprintf('This number should be small: %f\n', norm(x - A\rhs));
	\endcode
	
	\section skew_usage Usage for skew-symmetric matrices
	
    When the matrix is skew-symmetric, almost all documentation above still applies. The only difference is that the executable is \c skew_ldl_driver instead of \c ldl_driver. The skew functionality of \b sym-ildl can be found in the experimental branch of <a href="https://github.com/inutard/matrix-factor">https://github.com/inutard/matrix-factor</a>.
    
	\subsection skew_ldl_driver As a standalone program
	
	Let's factor the \c m3dskew50 matrix stored in <c>test_matrices/skew/m3dskew50.mtx</c>. As before, this is as simple as:
	\par
	\code
		./skew_ldl_driver -filename=test_matrices/skew/m2dskew50.mtx
		
		Load succeeded. File test_matrices/m2dskew50.mtx was loaded.
		A is 125000 by 125000 with 735000 non-zeros.
		Equilibration:  0.016 seconds.
		AMD:            0.203 seconds.
		...
	\endcode
	
	As in the symmetric case, we used the default values for the parameter we did not specify. A description of each of these parameters can be accessed by typing 
	\code 
		./skew_ldl_driver --help 
	\endcode
	
	\subsection skew_matlab_mex Within MATLAB
	Within MATLAB, using \b sym-ildl is even easier. As in the symmetric case, the command \c ildl can be used. Everything remains the same as the symmetric case, as \c ildl automatically detects whether the input is symmetric or skew-symmetric.
	
	Let's first generate a skew-symmetric matrix for testing:
	\par
	\code
		>> B = sparse(gallery('uniformdata',100,0));
		>> A = B-B';
	\endcode
	
	Since B is a matrix of random values between 0 and 1, A is almost certainly non-singular. Now we can call \c ildl exactly as before:
	\par
	\code
		>> [L, D, p, S, B] = ildl(A, 1.0, 0.001);
		Equilibration:	0.000 seconds.
        AMD:		    0.000 seconds.
        Permutation:	0.000 seconds.
        Factorization:	0.005 seconds.
        Total time:	    0.005 seconds.
        L is 100 by 100 with 4815 non-zeros.
	\endcode

	Finally, helpful information on the parameters for <c>ildl</c> can be found by typing:
	\code
		help ildl
	\endcode

*	\section contribute_sec Licensing and How to contribute
*
    \b sym-ildl is open source, and we're always looking for new contributions! The entire codebase is freely accessible at <a href="https://github.com/inutard/matrix-factor">https://github.com/inutard/matrix-factor</a>. We also use the MIT Licence, which essentially allow free use of this software in any way you want (see <a href="http://opensource.org/licenses/MIT">here</a> for more details). Simply send us a pull request on GitHub to contribute.
*
*
*   \section citations Citing this code
    If you have found our code useful, please consider citing us! A technical report (pdf), a BibTeX entry, and a link to the code are available here:
*    
    <em>
	SYM-ILDL: Incomplete LDL^{T} Factorization of Symmetric Indefinite and Skew-Symmetric Matrices. 
	Chen Greif, Shiwen He, Paul Liu. 
	UBC Computer Science Technical Report, 2016.
	</em> <a href="http://arxiv.org/pdf/1505.07589v1">[pdf]</a> <a href="http://dblp.uni-trier.de/rec/bibtex/journals/corr/GreifHL15">[bib]</a>
*    
*	\section refs References
    
	-#	J. A. George and J. W-H. Liu, <em>Computer Solution of Large Sparse Positive Definite Systems</em>, Prentice-Hall, 1981.
	-#	J. R. Bunch, <em>Equilibration of Symmetric Matrices in the Max-Norm</em>, JACM, 18 (1971), pp. 566-572.
	-#	N. Li and Y. Saad, <em>Crout versions of the ILU factorization with pivoting for sparse symmetric matrices</em>, ETNA, 20 (2006), pp. 75-85.
	-#	N. Li, Y. Saad, and E. Chow, <em>Crout versions of ILU for general sparse matrices</em>, SISC, 25 (2003), pp. 716-728.

*	
*/