\documentclass[12pt,oneside]{book} \usepackage{amsmath,amssymb,amsthm,latexsym,paralist,comment} \usepackage{graphicx} \usepackage{psfrag} \usepackage{soul} \usepackage{amsmath} \usepackage{multirow} \usepackage{algorithm} \usepackage{algpseudocode} \usepackage{graphicx} \usepackage{subcaption} \usepackage{listings} \usepackage[dvipsnames]{xcolor} \usepackage[hidelinks]{hyperref} \usepackage{framed} \usepackage{mdframed} \hypersetup{ colorlinks = true, citecolor = blue, linkcolor = blue, urlcolor = Maroon } \usepackage[letterpaper,margin=1in]{geometry} \theoremstyle{definition} \newcommand{\N}{\mathbf{N}} \newcommand{\R}{\mathbf{R}} \newcommand{\Z}{\mathbf{Z}} \newcommand{\import}{\textcolor{red}{\textbf{**IMPORTANT**}}} % Whether to include SPEX update or not % If update is set to true, then all of the text on updates is shown % if update is set to false, then all of the text on updates is hidden \usepackage{ifthen} \newboolean{update} \setboolean{update}{false} \newcommand{\update}[1]{\ifthenelse{\boolean{update}}{#1}{}} \usepackage[Glenn]{fncychap} %Fancy Chapter Headers \usepackage{cprotect} % ------------------------------------------------------ % Centered Tabular Environment, with default small size % ------------------------------------------------------ % Example 1 (small): \begin{SizedCenteredTabular}|l|l|l|} % Example 2 (scriptsize): \begin{SizedCenteredTabular}[\scriptsize]{|l|l|l|} \newenvironment{SizedCenteredTabular}[2][\small] { #1\begin{center}\begin{tabular}{#2} }{ \end{tabular}\end{center} } \begin{document} \thispagestyle{empty} \begin{center}\begin{large} \phantom{.}\\[1in] \textbf{User Guide for the SPEX Software Package} \\ \vspace{5mm} Version 3.1, March, 2024 % VERSION \vspace{20mm} Jinhao Chen, Timothy A. Davis, Christopher Lourenco, Lorena Mejia-Domenzain, Erick Moreno-Centeno \\ Texas A\&M University and US Naval Academy \vspace{20mm} Contact Information: Contact Chris Lourenco, \href{mailto:chrisjlourenco@gmail.com}{chrisjlourenco@gmail.com}, \href{mailto:lourenco@usna.edu}{lourenco@usna.edu}, or Tim Davis, \href{mailto:timdavis@aldenmath.com}{timdavis@aldenmath.com}, \href{mailto:davis@tamu.edu}{davis@tamu.edu}, \href{DrTimothyAldenDavis@gmail.com}{DrTimothyAldenDavis@gmail.com} \end{large}\end{center} %------------------------------------------------------------------------------- \newpage \tableofcontents %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- \chapter{SPEX Overview}\vspace{-0.75in} %------------------------------------------------------------------------------- SPEX is a software package comprising several state-of-the-art SParse EXact linear algebra routines. It currently consists of the following: \begin{description} \item[SPEX Utilities] Utility and auxiliary functions for all SPEX routines: interface to the GMP/MPFR library, memory management functions, the \verb|SPEX_matrix|, \linebreak \verb|SPEX_factorization|, and \verb|SPEX_symbolic_analysis| data structures, and various functions that are auxiliary to the factorization and solve functions. Please refer to Chapter \ref{ch:Util} for further details. \item[SPEX LU] Sparse exact left-looking LU factorization to solve the linear system $A \mathbf{x} = \mathbf{b}$. The solution time is proportional to the arithmetic work in the bit-complexity model; which is asymptotically efficient. Please refer to Chapter \ref{ch:LeftLU} for further details. \item[SPEX Cholesky] Sparse exact left-looking and up-looking Cholesky factorizations to solve the symmetric positive definite (SPD) linear system $A \mathbf{x} = \mathbf{b}$. The solution time is proportional to the arithmetic work in the bit-complexity model; this is an asymptotically efficient complexity bound. Please refer to Chapter \ref{ch:Chol} for further details. \item[SPEX Backslash] Routines to exactly solve the system $A \mathbf{x} = \mathbf{b}$ using either LU or Cholesky factorization. This is the simplest way to access the SPEX software package. Please refer to Chapter \ref{ch:Backslash} for further details. \update{ \item[SPEX Update] Sparse exact factorization update. Currently consists of LU and Cholesky column replacement updates and Cholesky rank 1 update. Please refer to Chapter \ref{ch:Update} for further details.} \end{description} \noindent \textbf{Location:} \url{https://github.com/clouren/SPEX} and \url{www.suitesparse.com}\\ \noindent \textbf{Required Packages:} SPEX depends on the following packages: \begin{itemize} \item GNU GMP \cite{granlund2015gnu} and MPFR \cite{fousse2007mpfr} libraries. Distributed under the LGPL3 and GPL2 and can be acquired and installed from \url{https://gmplib.org/} and \url{http://www.mpfr.org/}, respectively. \item CMake \cite{hoffman2003cmake}, available under a BSD 3-clause license. May be independently obtained at \url{https://cmake.org}. \item AMD \cite{amestoy1996approximate,amestoy2004algorithmamd}, available under a BSD 3-clause license and distributed along with SPEX. May be independently obtained at \url{www.suitesparse.com} \item COLAMD \cite{davis2004column,davis2004algorithmcolamd}, available under a BSD 3-clause license and distributed along with SPEX. May be independently obtained at \url{www.suitesparse.com} %COMMENTED OUT because (1) no license, (2) distributed in SPEX, and (3) consistency with paper %\item \verb|SuiteSparse_config| \cite{davis2020suitesparse}, no license restrictions and distributed along with SPEX. May be independently obtained at \url{www.suitesparse.com} \end{itemize} %------------------------------------------------------------------------------- \chapter{Setting up SPEX}\vspace{-0.75in} %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- \section{Licensing} \label{s:util:licensing} %------------------------------------------------------------------------------- \textbf{Copyright:} The copyright of this software is held by Christopher Lourenco, Jinhao Chen, Lorena Mejia-Domenzain, Erick Moreno-Centeno, and Timothy A. Davis.\\ \noindent \textbf{Contact Info:} Chris Lourenco, \href{mailto:chrisjlourenco@gmail.com}{chrisjlourenco@gmail.com} \href{mailto:lourenco@usna.edu}{lourenco@usna.edu}, or Tim Davis, \href{mailto:timdavis@aldenmath.com}{timdavis@aldenmath.com}, \href{DrTimothyAldenDavis@gmail.com}{DrTimothyAldenDavis@gmail.com}, or \href{mailto:davis@tamu.edu}{davis@tamu.edu}\\ \noindent \textbf{License:} This software package is dual licensed under the GNU General Public License version 2 or the GNU Lesser General Public License version 3. Details of this license are in \verb|SPEX/License/license.txt|. For alternative licenses, please contact the authors. %------------------------------------------------------------------------------- \section{Installation} \label{s:util:install} %------------------------------------------------------------------------------- Installation of SPEX requires the \verb|cmake| utility in Linux, MacOS, and Windows. With the appropriate compiler and version of cmake, typing \verb|make| under the main directory will compile AMD, COLAMD, and SPEX to their respective \verb|build| folder. All shared library files can be found in the top level \verb|build| folder. To further install the libraries onto your computer, simply type \verb|make install|. Thereafter, to use the code inside of your program, precede your code with \newline \verb|#include "SPEX.h"|. SPEX is also distributed with MATLAB and Python interfaces. Note that these interfaces have been thoroughly tested in Linux and are tuned to work ``out of the box'' on these types of machines. However, if the end user wishes to utilize the MATLAB or Python interfaces within a MacOS or Windows system, they may require additional library linkage in order to function properly. For example, on the MacOS, MATLAB R2022 does not currently support binaries compiled on ARM architecture, thus the code would have to be compiled with x-86. To install the MATLAB interface, navigate to the \verb|SPEX/MATLAB| folder from the MATLAB command window and type \verb|spex_mex_install| which will install the MATLAB interfaces to all SPEX packages. These packages can then be used outside of the \verb|SPEX/MATLAB| folder by using the MATLAB addpath tool. The Python interface does not need any additional installation, but does require the \texttt{Numpy}, \texttt{SciPy}, and \texttt{ctypes} libraries. Note that \chapter{General SPEX Data Structures and Macros}\vspace{-0.75in} The following macros/data structures are defined in \verb|SPEX.h| and are used in all SPEX functions. %------------------------------------------------------------------------------- %\cprotect\section{\verb|SPEX_VERSION|: the software package version} \section{\texttt{SPEX\_VERSION}: the software package version} %------------------------------------------------------------------------------- SPEX defines the following strings with \verb|#define|. Refer to the \verb|SPEX.h| file for details. \begin{SizedCenteredTabular}{ll} \hline Macro & Purpose \\ \hline \verb|SPEX_VERSION| & Current version of the code (as a string)\\ \verb|SPEX_VERSION_MAJOR| & Major version of the code \\ \verb|SPEX_VERSION_MINOR| & Minor version of the code \\ \verb|SPEX_VERSION_SUB| & Sub version of the code \\ \hline \end{SizedCenteredTabular} %------------------------------------------------------------------------------- %\cprotect\section{\verb|SPEX_info|: status codes returned by SPEX} \section{\texttt{SPEX\_info}: status codes returned by SPEX} \label{ss:SPEX_info} %------------------------------------------------------------------------------- Most SPEX functions return their status to the caller as their return value, an enumerated type called \verb|SPEX_info|. All current possible values for \verb|SPEX_info| are listed as follows: \begin{SizedCenteredTabular}{rll} \hline 0 & \verb|SPEX_OK|& The function was successfully executed.\\ \hline -1 & \verb|SPEX_OUT_OF_MEMORY|& Out of memory\\ \hline -2 & \verb|SPEX_SINGULAR|& The input matrix $A$ is exactly singular.\\ \hline -3 & \verb|SPEX_INCORRECT_INPUT|& One or more input arguments are incorrect.\\ \hline -4 & \verb|SPEX_NOTSPD| & The input matrix is not SPD (thus can't use Cholesky) \\ \hline -5 & \verb|SPEX_INCORRECT_ALGORITHM| & The algorithm is not compatible with the factorization \\ \hline -6 & \verb|SPEX_PANIC| & SPEX environment error \\ \hline \end{SizedCenteredTabular} %------------------------------------------------------------------------------- %\cprotect\section{\verb|SPEX_pivot|: enum for pivoting schemes} \section{\texttt{SPEX\_pivot}: enum for pivoting schemes}\label{ss:SPEX_pivot} %------------------------------------------------------------------------------- There are six available pivoting schemes provided in SPEX that can be selected with the \verb|SPEX_options| structure. If the matrix is non-singular (in an exact sense), then the pivot is always nonzero, and is chosen as the {\em smallest} nonzero entry, with the smallest magnitude. This may seem counter-intuitive, but selecting a small nonzero pivot leads to smaller growth in the number of digits in the entries of \verb|L| and \verb|U|. This choice does not lead to any kind of numerical inaccuracy, since SPEX is guaranteed to find an exact roundoff-error free factorization of a non-singular matrix (unless it runs out of memory), for any nonzero pivot choice. The pivot tolerance for two of the pivoting schemes is specified by the \verb|tol| component in \verb|SPEX_options|. The pivoting schemes are as follows: \begin{SizedCenteredTabular}{llp{4in}}\hline 0 & \verb|SPEX_SMALLEST| & The $k$-th pivot is selected as the smallest entry in the $k$-th column.\\ \hline 1 & \verb|SPEX_DIAGONAL| & The $k$-th pivot is selected as the diagonal entry. If the diagonal entry is zero, this method instead selects the smallest pivot in the column.\\ \hline 2 & \verb|SPEX_FIRST_NONZERO| & The $k$-th pivot is selected as the first eligible nonzero in the column. \\ \hline 3 & \verb|SPEX_TOL_SMALLEST| & The $k$-th pivot is selected as the diagonal entry if the diagonal is within a specified tolerance of the smallest entry in the column. Otherwise, the smallest entry in the $k$-th column is selected. This is the default pivot selection strategy. \\ \hline 4 & \verb|SPEX_TOL_LARGEST| & The $k$-th pivot is selected as the diagonal entry if the diagonal is within a specified tolerance of the largest entry in the column. Otherwise, the largest entry in the $k$-th column is selected. \\ \hline 5 & \verb|SPEX_LARGEST| & The $k$-th pivot is selected as the largest entry in the $k$-th column. \\ \hline \end{SizedCenteredTabular} %------------------------------------------------------------------------------- %\cprotect\section{\verb|SPEX_preorder|} \label{ss:SPEX_preorder} \section{\texttt{SPEX\_preorder}} \label{ss:SPEX_preorder} %------------------------------------------------------------------------------- The SPEX Library provides three ordering schemes: no ordering, COLAMD, and AMD. In LU factorization, the ordering is applied only to the columns, that is this ordering gives the matrix Q. In Cholesky factorizations, the ordering is applied to both the rows and columns, that is the ordering gives the matrices P and Q. \begin{SizedCenteredTabular}{llp{4in}} \hline 1 & \verb|SPEX_NO_ORDERING| & No pre-ordering is performed on the matrix $A$, that is $Q = I$. \\ \hline 2 & \verb|SPEX_COLAMD| & The rows and/or columns of $A$ are permuted prior to factorization using the COLAMD \cite{davis2004algorithmcolamd} ordering. This is recommended for LU factorization. \\ \hline 3 & \verb|SPEX_AMD| & The rows and/or columns of $A$ are permuted prior to the factorization using the the AMD \cite{amestoy2004algorithmamd}. This is recommended for Cholesky factorization. \\ \hline \end{SizedCenteredTabular} %------------------------------------------------------------------------------- %\cprotect\section{\verb|SPEX_factorization_algorithm|} \label{ss:SPEX_factorization_algorithm} \section{\texttt{SPEX\_factorization\_algorithm}} \label{ss:SPEX_factorization_algorithm} %------------------------------------------------------------------------------- This code tells SPEX which factorization is being used. Importantly, this is only used within a given solver. That is, this code is only used within LU/Cholesky factorization codes themselves. This is \textbf{NOT} used in the SPEX Backslash routines as that code selects the type of factorization using its own logic. \begin{SizedCenteredTabular}{llp{4in}} \hline 1 & \verb|SPEX_LU_LEFT| & Left-looking LU factorization \\ \hline 2 & \verb|SPEX_CHOL_LEFT| & Left-looking Chokesy factorization\\ \hline 3 & \verb|SPEX_CHOL_UP| & Up-looking Cholesky factorization \\ \hline \end{SizedCenteredTabular} %------------------------------------------------------------------------------- %\cprotect\section{ \verb|SPEX_options| structure} \section{\texttt{SPEX\_options} structure} \label{ss:SPEX_options_struct} %------------------------------------------------------------------------------- The \verb|SPEX_options| struct stores key command parameters for various functions used in the SPEX package. The \verb|SPEX_options* option| struct contains the following components: \begin{itemize} \item \verb|option->pivot|: An enum \verb|SPEX_pivot| type which controls the type of pivoting used. Default value: \verb|SPEX_SMALLEST| (3). \item \verb|option->order|: An enum \verb|SPEX_preorder| type which controls what column ordering is used. Default value: \verb|SPEX_COLAMD| for LU and \verb|SPEX_AMD| for Cholesky. \item \verb|option->tol|: A \verb|double| tolerance for the tolerance-based pivoting scheme, i.e., \newline \verb|SPEX_TOL_SMALLEST| or \verb|SPEX_TOL_LARGEST|. \verb|option->tol| must be in the range of $(0,1]$. Default value: 1 meaning that the diagonal entry will be selected if it has the same magnitude as the smallest entry in the $k$ the column. \item \verb|option->print_level|: An \verb|int| which controls the amount of output: 0: print nothing, 1: just errors, 2: terse, with basic stats from COLAMD/AMD and SPEX, 3: all, with matrices and results. Default value: 0. \item \verb|option->prec|: An \verb|int32_t| which specifies the precision used for multiple precision floating point numbers, (i.e., MPFR). This can be any integer larger than \verb|MPFR_PREC_MIN| (value of 1 in MPFR 4.0.2 and 2 in some legacy versions) and smaller than \verb|MPFR_PREC_MAX| (usually the largest possible integer available in your system). Default value: 128 (quad precision). \item \verb|option->round|: A \verb|mpfr_rnd_t| which determines the type of MPFR rounding to be used by SPEX. This is a parameter of the MPFR library. The options for this parameter are: \begin{itemize} \item \verb|MPFR_RNDN|: Round to nearest (roundTiesToEven in IEEE 754-2008) \item \verb|MPFR_RNDZ|: Round toward zero (roundTowardZero in IEEE 754-2008) \item \verb|MPFR_RNDU|: Round toward plus infinity (roundTowardPositive in IEEE 754-2008) \item \verb|MPFR_RNDD|: Round toward minus infinity (roundTowardNegative in IEEE 754-2008) \item \verb|MPFR_RNDA|: Round away from zero \item \verb|MPFR_RNDF|: Faithful rounding. This is not stable. \end{itemize} \noindent Refer to the MPFR User Guide available at \url{https://www.mpfr.org/mpfr-current/mpfr.pdf} for details on the MPFR rounding style and any other utilized MPFR convention. Default value: \verb|MPFR_RNDN|. \item \verb|option->algo|: A \verb|SPEX_factorization_algorithm| which indicates which type of factorization is being used. \end{itemize} All SPEX routines except basic memory management routines in Sections \ref{ss:SPEX_finalize}-\ref{ss:SPEX_calloc} and \verb|SPEX_options| allocation routine in \ref{ss:create_default_options} require \verb|option| as an input argument. The construction of the \verb|option| struct can be avoided by passing \verb|NULL| for the default settings. Otherwise, the following functions create and destroy a \verb|SPEX_options| structure: \begin{SizedCenteredTabular}{lp{2.5in}l} \hline Function/Macro Name & Description & Section \\ \hline \verb|SPEX_create_default_options| & create and return \verb|SPEX_options| pointer with default parameters upon successful allocation & \ref{ss:create_default_options} \\ \hline \verb|SPEX_FREE| & destroy \verb|SPEX_options| structure & \ref{ss:SPEX_free} \\ \hline \end{SizedCenteredTabular} %------------------------------------------------------------------------------- %\cprotect\section{\verb|SPEX_vector|} \label{ss:SPEX_vector} \section{\texttt{SPEX\_vector}} \label{ss:SPEX_vector} %------------------------------------------------------------------------------- SPEX vector is a compressed sparse vector data structure which will be used for SPEX dynamic CSC matrices. This struct is not used in SPEX version 3.0 and its funcionality will be fully developed in a future release of SPEX; however the struct is provided here so that future versions of SPEX have backward compatibility. \update{ This is only used publicly when calling the functions in SPEX Update to construct the vector to modify original matrix A, (either \verb|w| for $A'=A+\sigma ww^T$ in rank-1 update/downdate or \verb|vk| to be swapped with \verb|A->v[k]| in the update for column replacement). } This is \textbf{NOT} intended to be used for building any n-by-1 vector (e.g., the right-hand-side vector b in Ax=b), which should be considered as a n-by-1 \verb|SPEX_matrix|. This struct contains the following components: \begin{itemize} \item \verb|vector->nz|: The number of explicit entries in the vector. Data Type: \verb|int64_t|. \item \verb|vector->nzmax|: The size of the \verb|i| and \verb|x| arrays. Note that \verb|nz| $\le$ \verb|nzmax|. Data Type: \verb|int64_t|. \item \verb|vector->i|: An array of size \verb|nzmax| containing the row indices of all explicit entries in the vector. The last \verb|(nzmax-nz)| entries are undefined. Data Type: \verb|int64_t*|. \item \verb|vector->x|: An array of size \verb|nzmax| containing the numeric values of all explicit entries in the vector. The last \verb|(nzmax-nz)| entries are undefined. Data Type: \verb|mpz_t*|. \item \verb|vector->scale|: Scaling parameter. The actual value of the $k$-th nonzero should be computed as \verb|x[k]*scale|. Both \verb|x[k]*scale| and \verb|x[k]/mpq_denref(scale)| must be integer for all entries, where \verb|mpq_denref(scale)| is a GMP macro that gives the denominator of \verb|scale|. This is used to skip explicit update(s) for a column/row of the factorization matrix, when all entries are to be multiplied with the same scaling factor(s). Data Type: \verb|mpq_t|. \end{itemize} In the current release, the \verb|SPEX_vector| is only used as a part of the \verb|SPEX_matrix| struct and is always a NULL pointer. \update{ The SPEX package has a set of functions to allocate, destroy and reallocate a SPEX vector, \verb|SPEX_vector|, as shown in the following table. \begin{SizedCenteredTabular}{lll} \hline Function Name & Description & Section \\ \hline \verb|SPEX_vector_allocate| & allocate \verb|SPEX_vector| with \verb|nzmax| entries & \ref{ss:spex_vector_allocate} \\ \hline \verb|SPEX_vector_realloc| & reallocate \verb|SPEX_vector| with \verb|new_size| entries & \ref{ss:spex_vector_realloc} \\ \hline \verb|SPEX_vector_free| & destroy a \verb|SPEX_vector| and free its allocated & \ref{ss:spex_vector_free} \\& memory&\\ \hline \end{SizedCenteredTabular} }%End Update %------------------------------------------------------------------------------- %\cprotect\section{The \verb|SPEX_matrix| structure} \label{ss:SPEX_matrix} \section{The \texttt{SPEX\_matrix} structure} \label{ss:SPEX_matrix} %------------------------------------------------------------------------------- SPEX operates on matrices stored in any of the 16 different matrix formats: 15 of which are combinations of matrix formats and entry data-types: $\{$Static Compressed Sparse Column (CSC), triplet, dense$\} \times \{$ \verb|mpz_t|, \verb|mpq_t|, \verb|mpfr_t|, \verb|int64_t|, or \verb|double|$\}$, and the 16th of which is the dynamic CSC matrix with \verb|mpz_t| entries. Using the SPEX matrix copy function, a matrix of any given form and data-type can be copied and converted into a matrix of any one of the 16 matrix-form and data-type combinations. Most routines require the matrix to be in CSC form with \verb|mpz_t| (i.e., arbitrary-sized integer) data type. This data structure stores the matrix $A$ as a sequence of three arrays: \begin{itemize} \item \verb|A->p|: Column pointers; an array of size \verb|n+1|. The row indices of column $j$ are located in positions \verb|A->p[j]| to \verb|A->p[j+1]-1| of the array \verb|A->i|. Data type: \verb|int64_t|. \item \verb|A->i|: Row indices; an array of size equal to the number of entries in the matrix. The entry \verb|A->i[k]| is the row index of the $k$th nonzero in the matrix. Data type: \verb|int64_t|. \item \verb|A->x|: Numeric entries. The entry \verb|A->x[k]| is the numeric value of the $k$th nonzero in the matrix. The array \verb|A->x| has a union type and must be accessed via a suffix according to the type of \verb|A|. For details, please refer to Section~\ref{ss:SPEX_matrix}. \end{itemize} An example matrix $A$ with \verb|mpz_t| type is stored as follows (note that indexing is zero based as per the C convention). \[A = \begin{bmatrix} 1 & 0 & 0 & 1 \\ 2 & 0 & 4 & 12 \\ 7 & 1 & 1 & 1 \\ 0 & 2 & 3 & 0 \\ \end{bmatrix}\] \begin{verbatim} A->p = [0, 3, 5, 8, 11] A->i = [0, 1, 2, 2, 3, 1, 2, 3, 0, 1, 2] A->x.mpz = [1, 2, 7, 1, 2, 4, 1, 3, 1, 12, 1] \end{verbatim} For example, the last column appears in positions 8 to 10 of \verb|A->i| and \verb|A->x.mpz|, with row indices 0, 1, and 2, and values $a_{03}=1$, $a_{13}=12$, and $a_{23}=1$. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_kind|: enum for matrix formats} \subsection{\texttt{SPEX\_kind}: enum for matrix formats} \label{ss:SPEX_kind} %------------------------------------------------------------------------------- The SPEX library provides four available matrix formats: sparse CSC (compressed sparse column), sparse triplet, dense and sparse dynamic CSC. \begin{SizedCenteredTabular}{llp{4in}} \hline 0 & \verb|SPEX_CSC| & Matrix is in compressed sparse column format. \\ \hline 1 & \verb|SPEX_TRIPLET| & Matrix is in sparse triplet format. \\ \hline 2 & \verb|SPEX_DENSE| & Matrix is in dense format.\\ \hline 3 & \verb|SPEX_DYNAMIC_CSC| & Matrix is in dynamic CSC format. \\ \hline \end{SizedCenteredTabular} %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_type|: enum for data types of matrix entries} \subsection{\texttt{SPEX\_type}: enum for data types of matrix entries} \label{ss:SPEX_type} %------------------------------------------------------------------------------- The SPEX library provides five data types for matrix entries: \verb|mpz_t|, \verb|mpq_t|, \verb|mpfr_t|, \verb|int64_t| and \verb|double|. \begin{SizedCenteredTabular}{llp{4in}} \hline 0 & \verb|SPEX_MPZ| & Matrix entries are in \verb|mpz_t| type: an integer of arbitrary size. \\ \hline 1 & \verb|SPEX_MPQ| & Matrix entries are in \verb|mpq_t| type: a rational number with arbitrary-sized integer numerator and denominator. \\ \hline 2 & \verb|SPEX_MPFR| & Matrix entries are in \verb|mpfr_t| type: a floating-point number of arbitrary precision. \\ \hline 3 & \verb|SPEX_INT64| & Matrix entries are in \verb|int64_t| type. \\ \hline 4 & \verb|SPEX_FP64| & Matrix entries are in \verb|double| type. \\ \hline \end{SizedCenteredTabular} %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_matrix| structure} \subsection{\texttt{SPEX\_matrix} structure} %------------------------------------------------------------------------------- A matrix \verb|SPEX_matrix *A| has the following components: \begin{itemize} \item \verb|A->kind|: Indicating the kind of matrix A: CSC, triplet, dense or dynamic CSC. Data Type: \verb|SPEX_kind|. \item \verb|A->type|: Indicating the type of entries in matrix A: \verb|mpz_t|, \verb|mpq_t|, \verb|mpfr_t|, \verb|int64_t| or \verb|double|. Data Type: \verb|SPEX_type|. \item \verb|A->m|: Number of rows in the matrix. Data Type: \verb|int64_t|. \item \verb|A->n|: Number of columns in the matrix. Data Type: \verb|int64_t|. \item \verb|A->scale|: A scaling parameter for matrix of \verb|mpz_t| type. For all matrices whose entries are stored in data type other than \verb|mpz_t|, SPEX assumes and maintains \verb|A->scale = 1|. This is used to ensure that entry can be represented as an integer in an \verb|mpz_t| matrix if these entries are converted from non-integer type data (such as double, variable precision floating point, or rational). Data Type: \verb|mpq_t|. \item \verb|A->nzmax|: The allocated size of the vectors \verb|A->i|, \verb|A->j| and \verb|A->x|. Note that \verb|A->nzmax| $\geq$ \verb|nnz(A)|, where \verb|nnz(A)| is the return value of \verb|SPEX_matrix_nnz(A,option)|. Data Type: \verb|int64_t|. \item \verb|A->nz|: The number of nonzeros in the matrix $A$, if $A$ is a triplet matrix (ignored for matrices in CSC, dense or dynamic CSC formats). Data Type: \verb|int64_t|. \item \verb|A->p|: An array of size \verb|A->n|$+1$ which contains column pointers of $A$, if $A$ is a CSC matrix (\verb|NULL| for matrices in triplet or dense formats). Data Type: \verb|int64_t*|. \item \verb|A->p_shallow|: A boolean indicating whether \verb|A->p| is shallow. A {\em shallow} pointer is one that refers to a component of another matrix or data structure. If \verb|A->p| is shallow, then it should not be modified as part of the \verb|A| matrix, and it is not freed if \verb|A| is freed. Data Type: \verb|bool|. \item \verb|A->i|: An array of size \verb|A->nzmax| which contains the row indices of the nonzeros in $A$, if $A$ is a CSC or triplet matrix (\verb|NULL| for dense matrices). The matrix is zero-based, so row indices are in the range of $[0,$ \verb|A->m|$-1]$. Data Type: \verb|int64_t*|. \item \verb|A->i_shallow|: A boolean indicating whether \verb|A->i| is shallow. Data Type: \verb|bool|. \item \verb|A->j|: An array of size \verb|A->nzmax| which contains the column indices of the nonzeros in $A$, if $A$ is a triplet matrix (\verb|NULL| for matrices in CSC or dense formats). The matrix is zero-based, so column indices are in the range of $[0,$ \verb|A->n|$-1]$. Data Type: \verb|int64_t*|. \item \verb|A->j_shallow|: A boolean indicating whether \verb|A->j| is shallow. Data Type: \verb|bool|. \item \verb|A->x|: An array of size \verb|A->nzmax| which contains the numeric values of the matrix. This array is a union, and must be accessed via one of: \verb|A->x.mpz|, \verb|A->x.mpq|, \verb|A->x.mpfr|, \verb|A->x.int64|, or \verb|A->x.fp64|, depending on the \verb|A->type| parameter. Data Type: \verb|union|. \item \verb|A->x_shallow|: A boolean indicating whether \verb|A->x| is shallow. Data Type: \verb|bool|. \item \verb|A->v|: If the matrix is a \verb|SPEX_DYNAMIC_CSC| this is an array of size \verb|A->n|, each of which is a dynamic column vector. Data Type: \verb|SPEX_vector**|. Always NULL in SPEX 3.0 \end{itemize} Specifically, for different kinds of \verb|A| of size \verb|A->m| $\times$ \verb|A->n| with \verb|nz| nonzero entries, its components are defined as: \begin{itemize} \item (0) \verb|SPEX_CSC|: A sparse matrix in CSC (compressed sparse column) format. \verb|A->p| is an \verb|int64_t| array of size \verb|A->n|+1, \verb|A->i| is an \verb|int64_t| array of size \verb|A->nzmax| (with $nz$ $\le$ \verb|A->nzmax|), and \verb|A->x.TYPE| is an array of size \verb|A->nzmax| of matrix entries (\verb'TYPE' is one of \verb|mpz|, \verb|mpq|, \verb|mpfr|, \verb|int64|, or \verb|fp64|). The row indices of column $j$ appear in \verb|A->i [A->p [j] ... A->p [j+1]-1]|, and the values appear in the same locations in \verb|A->x.TYPE|. The \verb|A->j| array is \verb|NULL|. \verb|A->nz| is ignored; the number of entries in \verb|A| is given by \verb|A->p [A->n]|. Row indices need not be sorted in each column, but duplicates cannot appear. \item (1) \verb|SPEX_TRIPLET|: A sparse matrix in triplet format. \verb|A->i| and \verb|A->j| are both \verb|int64_t| arrays of size \verb|A->nzmax|, and \verb|A->x.TYPE| is an array of values of the same size. The $k$th tuple has row index \verb|A->i [k]|, column index \verb|A->j [k]|, and value \verb|A->x.TYPE [k]|, with 0 $\le$ $k <$ \verb|A->nz|. The \verb|A->p| array is \verb|NULL|. Triplets can be unsorted, but duplicates cannot appear. \item (2) \verb|SPEX_DENSE|: A dense matrix. The integer arrays \verb|A->p|, \verb|A->i|, and \verb|A->j| are all \verb|NULL|. \verb|A->x.TYPE| is a pointer to an array of size \verb|A->m|*\verb|A->n|, stored in column-oriented format. The value of $A(i,j)$ is \verb|A->x.TYPE [p]| with \verb|p| = $i + j*$\verb|A->m|. \verb|A->nz| is ignored; the number of entries in \verb|A| is \verb|A->m| $\times$ \verb|A->n|. \item (3) \verb|SPEX_DYNAMIC_CSC|: Currently unused \update{ A sparse matrix in dynamic CSC format with the number of nonzeros in each column changing independently and dynamically, which is only used in the SPEX update functions. The matrix is held as an array of \verb|A->n| SPEX vectors, one per column. Each column is held as a \verb|SPEX_vector|, containing \verb|mpz_t| values and its own scale factor. For this kind, \verb|A->nzmax|, \verb|A->nz|, \verb|A->p|, \verb|A->i|, \verb|A->x| and \verb|A->*_shallow| are ignored and pointers \verb|p|, \verb|i| and \verb|x| are remained as NULL pointers. To access entries in column $j$, \verb|A->v[j]->i[0 ... A->v[j]->nz-1]| give the row indices of all nonzeros, and the \verb|mpz_t| values of these entries appear in the same locations in \verb|A->v[j]->x|. \verb|A->v[j]->nzmax| is the max number of nonzeros allocated for the $j$th column, while \verb|A->v[j]->nz| is the number of existing nonzeros in the $j$th column. Therefore, the total number of existing nonzeros is computed as $\sum_{j=0}^{n-1}$(\verb|A->v[j]->nz|).} \end{itemize} \verb|A| may contain shallow components, \verb|A->p|, \verb|A->i|, \verb|A->j|, and \verb|A->x|. For example, if \verb|A->p_shallow| is true, then a non-\verb|NULL| \verb|A->p| is a pointer to a read-only array, and the \verb|A->p| array is not freed by \verb|SPEX_matrix_free|. If \verb|A->p| is \verb|NULL| (for a triplet or dense matrix), then \verb|A->p_shallow| has no effect. %removed %To simplify the access the entries in \verb|A|, SPEX package provides the %following macros (Note that the \verb|TYPE| parameter in the macros is one of: %\verb|mpz|, \verb|mpq|, \verb|mpfr|, \verb|int64| or \verb|fp64|): %\begin{itemize} %\item %\verb|SPEX_1D(A,k,TYPE)|: used to access the $k$th entry in % \verb|SPEX_matrix *A| using 1D linear addressing for % any matrix kind (CSC, triplet or dense), in any type % with \verb|TYPE| specified corresponding % %\item %\verb|SPEX_2D(A,i,j,TYPE)|: used to access the $(i,j)$th entry in a dense % \verb|SPEX_matrix *A|. % %\end{itemize} The SPEX package has a set of functions to allocate, copy(convert), query and destroy a SPEX matrix, \verb|SPEX_matrix|, as shown in the following table. \begin{SizedCenteredTabular}{lp{2.5in}l} \hline Function Name & Description & Section \\ \hline \verb|SPEX_matrix_allocate| & allocate a $m$-by-$n$ \verb|SPEX_matrix| & \ref{s:user:matrix_allocate} \\ \hline \verb|SPEX_matrix_free| & destroy a \verb|SPEX_matrix| and free its allocated memory & \ref{s:user:matrix_free} \\ \hline \verb|SPEX_matrix_copy| & make a copy of a matrix, into another kind and/or type & \ref{s:user:matrix_copy} \\ \hline \verb|SPEX_matrix_nnz| & get the number of entries in a matrix & \ref{s:user:matrix_nnz} \\ \hline \verb|SPEX_matrix_check| & check the validity of a matrix and print it & \ref{s:user:matrix_check} \\ \hline \end{SizedCenteredTabular} %------------------------------------------------------------------------------- %\cprotect\section{The \verb|SPEX_symbolic_analysis| struct} \label{s:spex_symbolic_analysis} \section{The \texttt{SPEX\_symbolic\_analysis} struct} \label{s:spex_symbolic_analysis} %------------------------------------------------------------------------------- The symbolic analysis structure handles all preorderings and graphical struture information for each factorization within SPEX. First, section \ref{ss:spex_factorization_kind} discusses an enum for the type of factorization and next section \ref{ss:SPEX_symbolic_struct} discusses the components of this data structure. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_factorization_kind|: enum for kind of factorization} \subsection{\texttt{SPEX\_factorization\_kind}: enum for kind of factorization} \label{ss:spex_factorization_kind} %------------------------------------------------------------------------------- The SPEX library currently provides two types of factorizations: LU and Cholesky. The value \verb|SPEX_QR_FACTORIZATION| is reserved for future development. \begin{SizedCenteredTabular}{lll} \hline 0 & \verb|SPEX_LU_FACTORIZATION| & LU factorization is being used \\ \hline 1 & \verb|SPEX_CHOLESKY_FACTORIZATION| & Cholesky factorization is being used\\ \hline 2 & \verb|SPEX_QR_FACTORIZATION| & QR factorization is being used \\ & & (reserved for future use)\\ \hline \end{SizedCenteredTabular} %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_symbolic_analysis| Data Structure} \subsection{\texttt{SPEX\_symbolic\_analysis} Data Structure} \label{ss:SPEX_symbolic_struct} %------------------------------------------------------------------------------- A symbolic analysis \verb|SPEX_symbolic_analysis *S| has the following components: \begin{itemize} \item \verb|S->kind|: Indicating the kind of factorization either LU or Cholesky. Data type: \verb|SPEX_factorization_kind| \item \verb|S->P_perm|: Row permutation for Cholesky and LU factorization. Data type: \verb|int64_t*| \item \verb|S->Pinv_perm|: Inverse row permutation for Cholesky and LU factorization. Data type: \verb|int64_t*| \item \verb|S->Q_perm|: Column permutation for LU factorization. This is always \verb|NULL| and ignored for Cholesky factorization since its row and column permutations are the same. Data type: \verb|int64_t*| \item \verb|S->Qinv_perm|: Inverse column permutation for LU factorization. This is always \verb|NULL| and ignored for Cholesky factorization since its inverse row and column permutations are the same. Data type: \verb|int64_t*| \item \verb|S->lnz|: Approximate number of nonzeros in $L$. In LU factorization, this is a crude estimate based on either AMD or COLAMD. In Cholesky factorization, if AMD is used, this is the exact number of nonzeros in $L$. Data type: \verb|int64_t| \item \verb|S->unz|: Approximate number of nonzeros in $U$. In LU factorization, this is a crude estimate based on either AMD or COLAMD. In Cholesky factorization this is not used. Data type: \verb|int64_t| \item \verb|S->parent|: This is the elimination tree of the input matrix for Cholesky factorization. This is always \verb|NULL| for LU factorization. Data type: \verb|int64_t*| \item \verb|S->cp|: Column pointers of L for Cholesky factorization. This is always \verb|NULL| for LU factorization. Data type: \verb|int64_t*| \end{itemize} This data type is constructed when analysis is called in the appropriate factorizations. See sections \ref{ss:spex_lu_analyze} and \ref{ss:spex_chol_analyze} for further details. To free this data structure, the function \verb|SPEX_symbolic_analysis_free| is used and discussed further in section \ref{s:spex_symbolic_analysis_helper}. %------------------------------------------------------------------------------- %\cprotect\section{The \verb|SPEX_factorization| data structure} \section{The \texttt{SPEX\_factorization} data structure} %------------------------------------------------------------------------------- The \verb|SPEX_factorization| object holds an LU or Cholesky numerical factorization. \update{in either non-updatable or updatable form} The introduction of this structure is one of the largest API update for SPEX 2.0, as the components of all factorizations are now held in this structure instead of being carried around by the user. The components of the factorization structure are accessible to the user application. However, they should only be modified by calling SPEX methods. Changing them directly can lead to undefined behavior. The components of a \verb|SPEX_factorization* F| are as follows: \begin{itemize} \item \verb|F->kind|: Indicating the kind of factorization either LU or Cholesky. Data type: \verb|SPEX_factorization_kind| \item \verb|F->updatable|: a flag that indicates whether the factorization is in an updatable format. Reserved for future development. Data type: \verb|bool| \item \verb|F->scale_for_A|: Scaling factor of the input matrix $A$. As discussed in section \ref{ss:SPEX_matrix}, all matrices in SPEX are integral, thus, if A must be scaled the scaling factor applied is stored here. Data type: \verb|mpq_t| \item \verb|F->L|: The lower triangular matrix for either LU or Cholesky factorization. Data type: \verb|SPEX_matrix*| \item \verb|F->U|: The upper triangular matrix for LU factorization. This is always \verb|NULL| for Cholesky factorization. Data type: \verb|SPEX_matrix*| \item \verb|F->Q|: The matrix for (future) QR factorization. Provided here so that future versions of SPEX have backward compatibility. Data type: \verb|SPEX_matrix*| \item \verb|F->R|: The right triangular matrix for (future) QR factorization. Provided here so that future versions of SPEX have backward compatibility. Data type: \verb|SPEX_matrix*| \item \verb|F->rhos|: An $n \times 1$ dense matrix containing the pivot values used for LU or Cholesky factorization. Data type: \verb|SPEX_matrix*| \item \verb|F->P_perm|: Row permutation of the LU or Cholesky factors. Data type: \verb|int64_t*| \item \verb|F->Pinv_perm|: Inverse row permutation of the LU or Cholesky factors. Data type: \verb|int64_t*| \item \verb|F->Q_perm|: Column permutation of the LU factors. This is \verb|NULL| and ignored for Cholesky factorization. Data type: \verb|int64_t*| \item \verb|F->Qinv_perm|: Inverse column permutation of the LU factors. This is \verb|NULL| and ignored for Cholesky factorization. Data type: \verb|int64_t*| \end{itemize} \update{ As mentioned in the above description, one or more of these components could be \verb|NULL| for certain factorization. For example, \verb|F->U|, \verb|F->Q_perm| and \verb|F->Qinv_perm| are \verb|NULL| for Cholesky factorization. Moreover, \verb|F->Qinv_perm| can be \verb|NULL| for a non-updatable LU factorization with \verb|F->updatable=false|, but it will be generated when the LU factorization converted to the updatable form. Aside from that \verb|Qinv_perm| will be generated for the updatable LU factorization, the only difference between non-updatable and updatable factorizations is the way that $L$ (and $U$ if exists) is stored. Specifically, a updatable factorization must meet the following conditions. \begin{itemize} \item Both \verb|F->L| and \verb|F->U| are \verb|SPEX_DYNAMIC_CSC| matrices. Notably, \verb|F->U| in the updatable factorization is actually the transpose of $U$ (or equivalently, $U$ is stored as compressed sparse row (CSR) form instead), since $U$ will be updated one row at a time. \item $A = LDU$, which means $L$ and $U$ are properly permuted. (Recall that $PAQ = LDU$ or $PAP^T = LDL^T$ holds for static factorization). That is, for a updatable factorization, the rows of $L$ are in the same order as the rows of $A$, while the $j$-th column of $L$ (\verb|F->L->v[j]|) contains the $j$-th pivot, which would be \verb|F->L->v[j]->x[0]|, (i.e., \verb|F->L->v[j]->i[0] == F->P_perm[j]|); the columns of $U$ (i.e., the rows of $U^T$) are in the same order as the columns of $A$, while the $j$-th row of $U$ (i.e., the $j$-th column of $U^T$) (\verb|F->U->v[j]|) contains the $j$-th pivot, which would be \verb|F->U->v[j]->x[0]|, (i.e., \verb|F->U->v[j]->i[0] == F->Q_perm[j]|). \end{itemize} Due to these non-trivial conditions, users cannot simply perform \verb|SPEX_matrix_copy| to obtain \verb|F->L| and/or \verb|F->U| in \verb|SPEX_DYNAMIC_CSC SPEX_MPZ| format and claim that it is updatable, and vice versa. To correctly convert the factorization, user should call \verb|SPEX_factorization_convert| to perform in-place conversion for a given \verb|F| to either updatable or non-updatable as specified. } A \verb|SPEX_factorization| is constructed by the appropriate factorizations (see sections \ref{ss:spex_left_lu_factorize} and \ref{ss:spex_chol_factorize} for further details). To free this data structure, the function \verb|SPEX_factorization_free| is used and discussed further in section \ref{ss:spex_factorization_free}.\update{, and updated by appropriate \verb|SPEX_Update_*| functions (see sections \ref{ss:SPEX_Update_Chol_Rank1} and \ref{ss:SPEX_Update_LU_ColRep} for further details). In addition, the SPEX package has a set of functions to convert, query and destroy a SPEX factorization, \verb|SPEX_factorization|, as shown in the following table.} \update{ \begin{SizedCenteredTabular}{lp{2.5in}l} \hline Function Name & Description & Section \\ \hline \verb|SPEX_factorization_check| & check the validity of a factorization and print it & \ref{ss:spex_factorization_check} \\ \hline \verb|SPEX_factorization_convert| & convert a factorization to non-/updatable as specified & \ref{ss:spex_factorization_convert} \\ \hline \verb|SPEX_factorization_free| & destroy a \verb|SPEX_factorization| and free its allocated memory & \ref{ss:spex_factorization_free} \\ \hline \end{SizedCenteredTabular} } \update{It should be noted that all solvers (except \verb|SPEX_Update_(t)solve|) and functions that create factorization return non-updatable factorization with \verb|F->L| (and \verb|F->U| if exists) in \verb|SPEX_CSC SPEX_MPZ| form. On the other hand, all \verb|SPEX_Update_*| functions require and output updatable factorization with \verb|F->L| (and \verb|F->U| if exists) in \verb|SPEX_DYNAMIC_CSC SPEX_MPZ| form. These \verb|SPEX_Update_*| functions check the input factorization, convert it (if not updatable) automatically, and output updatable factorization.} %------------------------------------------------------------------------------- \chapter{SPEX Utilities}\vspace{-0.75in} \label{ch:Util} %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- \section{Overview} \label{s:util:overview} %------------------------------------------------------------------------------- SPEX Util contains utility and auxiliary functions for the SPEX factorizations. Additionally, SPEX Util provides a wrapper class for the GNU Multiple Precision Arithmetic (GMP) \cite{granlund2015gnu} and GNU Multiple Precision Floating Point Reliable (MPFR) \cite{fousse2007mpfr} libraries that prevent memory leaks and improve the overall stability of these external libraries. SPEX Util is written in ANSI C. %------------------------------------------------------------------------------- \section{Managing the SPEX environment} \label{s:user:setup} %------------------------------------------------------------------------------- Either \verb|SPEX_initialize| or \verb|SPEX_initialize_expert| (but not both) must be called prior to using any other SPEX functions. Otherwise, all SPEX user-callable functions would return \verb|SPEX_PANIC|. \verb|SPEX_finalize| must be called as the last SPEX function. Note that if a user is working in a multi threaded environment then only one user thread should call the \verb|SPEX_initialize| and \verb|SPEX_finalize| functions. Subsequent SPEX sessions can be restarted after a call to \verb|SPEX_finalize|, by calling either \verb|SPEX_initialize| or \verb|SPEX_initialize_expert| (but not both), followed by a final call to \verb|SPEX_finalize| when finished. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_initialize|: initialize the working environment} \subsection{\texttt{SPEX\_initialize}: initialize the working environment} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_initialize ( void ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_initialize| initializes the working environment for SPEX functions. SPEX utilizes a specialized memory management scheme in order to prevent potential memory failures caused by GMP and MPFR libraries. Either this function or \verb|SPEX_initialize_expert| must be called prior to using any other function in the library. Returns \verb|SPEX_PANIC| if SPEX has already been initialized, or \verb|SPEX_OK| if successful. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_initialize_expert|: initialize environment (expert version)}\label{ss:SPEX_initialize_expert} \subsection{\texttt{SPEX\_initialize\_expert}: initialize environment (expert version)} \label{ss:SPEX_initialize_expert} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_initialize_expert ( void* (*MyMalloc) (size_t), // user-defined malloc void* (*MyCalloc) (size_t, size_t), // user-defined calloc void* (*MyRealloc) (void *, size_t), // user-defined realloc void (*MyFree) (void *) // user-defined free ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_initialize_expert| is the same as \verb|SPEX_initialize| except that it allows for a redefinition of custom memory functions that are used for SPEX and GMP/ MPFR. The four inputs to this function are pointers to four functions with the same signatures as the ANSI C \verb'malloc', \verb'calloc', \verb'realloc', and \verb'free' functions. That is: \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} #include void *malloc (size_t size) ; void *calloc (size_t nmemb, size_t size) ; void *realloc (void *ptr, size_t size) ; void free (void *ptr) ; \end{verbatim} } \end{mdframed} Returns \verb|SPEX_PANIC| if SPEX has already been initialized, or \verb|SPEX_OK| if successful. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_finalize|: free the working environment} \subsection{\texttt{SPEX\_finalize}: free the working environment} \label{ss:SPEX_finalize} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_finalize ( void ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_finalize| finalizes the working environment for SPEX library, and frees any internal workspace created by SPEX. It must be called as the last \verb|SPEX_*| function called, except that a subsequent call to \verb|SPEX_initialize*| may be used to start another SPEX session. Returns \verb|SPEX_PANIC| if SPEX has not been initialized, or \verb|SPEX_OK| if successful. \subsection{\texttt{SPEX\_thread\_initialize}: initialize working environment for a single thread} \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_thread_initialize ( void ); \end{verbatim} } \end{mdframed} \verb|SPEX_thread_initialize| initializes the working environment of SPEX for a single user thread. If the user is working in a multithreaded environment, they must call this function at the beginning of each user thread. Returns \verb|SPEX_OK| if successful or \verb|SPEX_PANIC| if SPEX was already initialized. This function is only required for a multithreaded user application that calls SPEX functions from threads other than the primary thread that called \verb'SPEX_initialize'. When the primary thread of the user application starts, it must call \verb'SPEX_initialize'. When the user application enters a parallel region (say with OpenMP) or creates its own threads with a threading library, each user thread must call \verb'SPEX_thread_initialize' when it starts, and \verb'SPEX_thread_finalize' when it finishes. An example usage can be found in the \verb'SPEX/Demo' folder in the \verb'spex_demo_threaded.c' main program. \subsection{\texttt{SPEX\_thread\_finalize}: finalize the working environment for a single thread} \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_thread_finalize ( void ); \end{verbatim} } \end{mdframed} \verb|SPEX_thread_finalize| finalizes the working environment and frees any internal workspace created by SPEX for a single user thread. If the user is working in a multithreaded environment, they must call this function at the end of each user thread. Returns \verb|SPEX_OK| if successful or \verb|SPEX_PANIC| if SPEX was not initialized. %------------------------------------------------------------------------------- \section{Memory Management} \label{s:user:memmanag} %------------------------------------------------------------------------------- The routines in this section are used to allocate and free memory for the data structures used in SPEX. By default, SPEX relies on the SuiteSparse memory management functions, \verb|SuiteSparse_malloc|, \verb|SuiteSparse_calloc|, \verb|SuiteSparse_realloc|, and \verb|SuiteSparse_free|. By default, those functions rely on the ANSI C \verb|malloc|, \verb|calloc|, \verb|realloc|, and \verb|free|, but this may be changed by initializing the SPEX environment with \verb|SPEX_initialize_expert|. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_calloc|: allocate initialized memory} \subsection{\texttt{SPEX\_calloc}: allocate initialized memory} \label{ss:SPEX_calloc} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} void *SPEX_calloc ( size_t nitems, // number of items to allocate size_t size // size of each item ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_calloc| allocates a block of memory for an array of \verb|nitems| elements, each of them \verb|size| bytes long, and initializes all its bits to zero. If any input is less than 1, it is treated as if equal to 1. If the function failed to allocate the requested block of memory, then a \verb|NULL| pointer is returned. Returns \verb|NULL| if SPEX has not been initialized. %------------------------------------------------------------------------------- \newpage %\cprotect\subsection{\verb|SPEX_malloc|: allocate uninitialized memory} \subsection{\texttt{SPEX\_malloc}: allocate uninitialized memory} \label{ss:SPEX_malloc} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} void *SPEX_malloc ( size_t size // size of memory space to allocate ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_malloc| allocates a block of \verb|size| bytes of memory, returning a pointer to the beginning of the block. The content of the newly allocated block of memory is not initialized, remaining with indeterminate values. If \verb|size| is less than 1, it is treated as if equal to 1. If the function fails to allocate the requested block of memory, then a \verb|NULL| pointer is returned. Returns \verb|NULL| if SPEX has not been initialized. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_realloc|: resize allocated memory} \subsection{\texttt{SPEX\_realloc}: resize allocated memory} \label{ss:SPEX_realloc} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} void *SPEX_realloc // pointer to reallocated block, or original block // if the realloc failed ( int64_t nitems_new, // new number of items int64_t nitems_old, // previous/old number of items size_t size_of_item, // size of each item void *p, // pointer to reallocate bool *ok // true if success, false on failure ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_realloc| is a wrapper for realloc. If \verb|p| is non-\verb|NULL| on input, it points to a previously allocated array of size \verb|nitems_old| $\times$ \verb|size_of_item|. The array is reallocated to be of size \verb|nitems_new| $\times$ \verb|size_of_item|. If \verb|p| is \verb|NULL| on input, then a new array of that size is allocated. On success, a pointer to the new array is returned. Returns \verb|ok| as \verb|false| if SPEX has not been initialized. If the reallocation fails, \verb|p| is not modified, and \verb|ok| is returned as \verb|false| to indicate that the reallocation failed. If the size decreases or remains the same, then the method always succeeds (\verb|ok| is returned as \verb|true|), unless SPEX has not been initialized. Typical usage: the following code fragment allocates an array of 10 \verb|int|'s, and then increases the size of the array to 20 \verb|int|'s. If the \verb|SPEX_malloc| succeeds but the \verb|SPEX_realloc| fails, then the array remains unmodified, of size 10. \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} int *p ; p = SPEX_malloc (10 * sizeof (int)) ; if (p == NULL) { error here ... } printf ("p points to an array of size 10 * sizeof (int)\n") ; bool ok ; p = SPEX_realloc (20, 10, sizeof (int), p, &ok) ; if (ok) printf ("p has size 20 * sizeof (int)\n") ; else printf ("realloc failed; p still has size 10 * sizeof (int)\n") ; SPEX_free (p) ; \end{verbatim} } \end{mdframed} %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_free|: free allocated memory} \newpage \subsection{\texttt{SPEX\_free}: free allocated memory} \label{ss:SPEX_free} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} void SPEX_free ( void *p // Pointer to memory space to free ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_free| frees the memory previously allocated by a call to \verb|SPEX_calloc|, \verb|SPEX_malloc|, or \verb|SPEX_realloc|. If \verb|p| is \verb|NULL| on input, then no action is taken (this is not an error condition). To guard against freeing the same memory space twice, the following macro \verb|SPEX_FREE| is provided, which calls \verb|SPEX_free| and then sets the freed pointer to \verb|NULL|. \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} #define SPEX_FREE(p) \ { \ SPEX_free (p) ; \ (p) = NULL ; \ } \end{verbatim} } \end{mdframed} No action is taken if SPEX has not been initialized. %------------------------------------------------------------------------------- %\cprotect\section{\verb|SPEX_options| helper functions} \label{ss:SPEX_options} \section{\texttt{SPEX\_options} helper function} \label{ss:SPEX_options} %------------------------------------------------------------------------------- The \verb|SPEX_options| structure contains numerous parameters that may be modified to change the behavior of the SPEX functions. Default values of these parameters will lead to good performance in most cases. The following helper functions are provided. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_create_default_options|: create default \verb|SPEX_options| structure} \subsection{\texttt{SPEX\_create\_default\_options}: create default \texttt{SPEX\_options} structure} \label{ss:create_default_options} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_options* SPEX_create_default_options ( void ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_create_default_options| creates and returns a pointer to a \verb|SPEX_options| struct with default parameters upon successful allocation, which are discussed in Section \ref{ss:SPEX_options_struct}. To safely free the \verb|SPEX_options* option| structure, simply use \newline \verb|SPEX_FREE(option)|. All functions that require \verb|SPEX_options *option| as an input argument can have a \verb'NULL' pointer passed instead. In this case, the default value of the corresponding command option is used. \update{ %-------------------------------------------------------------------------------%\cprotect\section{\verb|SPEX_vector| helper functions} \label{s:spex_vector_helper} \section{\texttt{SPEX\_vector} helper functions} \label{s:spex_vector_helper} %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_vector_allocate|: Allocate a spex vector} \label{ss:spex_vector_allocate} \subsection{\texttt{SPEX\_vector\_allocate}: Allocate a SPEX vector} \label{ss:spex_vector_allocate} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_vector_allocate ( SPEX_vector **v_handle, // vector to be allocated const int64_t nzmax, // number of nnz entries in v const SPEX_options *option ); \end{verbatim} } \end{mdframed} \verb|SPEX_vector_allocate| creates and initializes a \verb|SPEX_vector| of size \verb|nzmax|. On input, the \verb|SPEX_vector| \verb|*v| that \verb|v_handle| points to is \verb|NULL|. On output, \verb|SPEX_vector| \verb|*v| is allocated with \verb|v->nz|=0, \verb|v->nzmax|=\verb|nzmax|, \verb|v->scale|=1, \verb|v->i| being an \verb|int64_t| array of size \verb|nzmax| and \verb|v->x| being an array containing \verb|nzmax| initialized \verb|mpz_t| entries. %-------------------------------------------------------------------------------%\cprotect\subsection{\verb|SPEX_vector_realloc|: Reallocate a spex vector} \label{ss:spex_vector_realloc} \subsection{\texttt{SPEX\_vector\_realloc}: Reallocate a SPEX vector} \label{ss:spex_vector_realloc} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_vector_realloc ( SPEX_vector* v, // the vector to be expanded const int64_t new_size, // desired new size for v const SPEX_options *option ); \end{verbatim} } \end{mdframed} \verb|SPEX_vector_realloc| either expands or shrinks the given \verb|SPEX_vector *v| of size \verb|v->nzmax| so that its new size is \verb|new_size|. If \verb|new_size > v->nzmax|, all newly allocated \verb|mpz_t| entries are initialized. %-------------------------------------------------------------------------------%\cprotect\subsection{\verb|SPEX_vector_free|: Free a spex vector} \label{ss:spex_vector_free} \subsection{\texttt{SPEX\_vector\_free}: Free a SPEX vector} \label{ss:spex_vector_free} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_vector_free ( SPEX_vector **v_handle, // vector to be deleted const SPEX_options *option ); \end{verbatim} } \end{mdframed} \verb|SPEX_vector_free| frees the \verb|SPEX_vector *v| that \verb|v_handle| points to. On output, the memory associated with \verb|v| is freed and \verb|v| is set to \verb|NULL|. } %-------------------------------------------------------------------------------%\cprotect\section{\verb|SPEX_matrix| helper functions} \label{s:spex_matrix_functions} \newpage \section{\texttt{SPEX\_matrix} helper functions} \label{s:spex_matrix_functions} %------------------------------------------------------------------------------- These functions provide several utilities for a \verb|SPEX_matrix|. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_matrix_allocate|: allocate an $m$-by-$n$ \verb|SPEX_matrix|} \subsection{\texttt{SPEX\_matrix\_allocate}: allocate an m-by-n \texttt{SPEX\_matrix}} \label{s:user:matrix_allocate} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_matrix_allocate ( SPEX_matrix **A_handle, // matrix to allocate SPEX_kind kind, // CSC, triplet, dense or SPEX_DYNAMIC_CSC SPEX_type type, // mpz, mpq, mpfr, int64, or double int64_t m, // # of rows int64_t n, // # of columns int64_t nzmax, // max # of entries bool shallow, // if true, matrix is shallow. A->p, A->i, A->j, // A->x are all returned as NULL and must be set // by the caller. All A->*_shallow are returned // as true. Ignored for SPEX_DYNAMIC_CSC // kind matrix. bool init, // If true, and the data types are mpz, mpq, or // mpfr, the entries of A->x are initialized // (using the proper SPEX_mp*_init function). // If false, the mpz, mpq, and mpfr arrays are // allocated but not initialized. Meaningless // for data types FP64 or INT64. Ignored if kind // is SPEX_DYNAMIC_CSC or shallow is true. const SPEX_options *option ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_matrix_allocate| allocates memory space for a $m$-by-$n$ \verb|SPEX_matrix| whose kind (CSC, triplet, dense, or dynamic CSC) and data type (\verb|mpz|, \verb|mpq|, \verb|mpfr|, \verb|int64| or \verb|fp64|) is specified. On input, the SPEX matrix that \verb|A_handle| points to is \verb|NULL|. On output, \verb|A_handle| points to a SPEX matrix of specified type, kind and size. %Returns SPEX_PANIC if SPEX has not been initialized. For a CSC, triplet or dense matrix, if \verb|shallow| is true, all components (\verb|A->p|, \verb|A->i|, \verb|A->j|, \verb|A->x|) are returned as \verb|NULL|, and their shallow flags are all true. The pointers \verb|A->p|, \verb|A->i|, \verb|A->j|, and/or \verb|A->x| can then be assigned from arrays in the calling application. If \verb|shallow| is false, the appropriate individual arrays are allocated (via \verb|SPEX_calloc|). The second boolean parameter \verb|init| is used if the entries are \verb|mpz_t|, \verb|mpq_t|, or \verb|mpfr_t|. Specifically, if \verb|init| is true, the individual entries within \verb|A->x.TYPE| are initialized using the appropriate \verb|SPEX_mp*_init| function. Otherwise, if \verb|init| is false, the \verb|A->x.TYPE| array is allocated (via \verb|SPEX_calloc|) and left that way. They are not otherwise initialized, and attempting to access the values of these uninitialized entries will lead to undefined behavior. For a \verb|SPEX_DYNAMIC_CSC| matrix, \verb|type|, \verb|shallow| and \verb|init| are ignored (since it only allows \verb|mpz_t| entries). Moreover, each column of the returned \verb|SPEX_DYNAMIC_CSC| matrix will be allocated as \verb|SPEX_vector| with zero available entry. Additional reallocation for each column will be needed. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_matrix_free|: free a \verb|SPEX_matrix|} \newpage \subsection{\texttt{SPEX\_matrix\_free}: free a \texttt{SPEX\_matrix}} \label{s:user:matrix_free} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_matrix_free ( SPEX_matrix **A_handle, // matrix to free const SPEX_options *option ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_matrix_free| frees the \verb|SPEX_matrix *A|. Note that the input of the function is the pointer to the pointer of a \verb|SPEX_matrix| structure. This is because this function internally sets the pointer of a \verb|SPEX_matrix| to be \verb|NULL| to prevent potential segmentation fault that could be caused by double \verb|free|. %If default settings are desired, %\verb|option| can be input as \verb|NULL|. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_matrix_copy|: make a copy of a \verb|SPEX_matrix| with a potentially different matrix-format and data-type} \subsection{\texttt{SPEX\_matrix\_copy}: make a copy of a \texttt{SPEX\_matrix} with a potentially different matrix-format and data-type} \label{s:user:matrix_copy} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_matrix_copy ( SPEX_matrix **C_handle, // matrix to create (never shallow) // inputs, not modified: SPEX_kind C_kind, // C->kind: CSC, triplet, dense, or dynamic SPEX_type C_type, // C->type: mpz_t, mpq_t, mpfr_t, int64_t, or double const SPEX_matrix *A, // matrix to make a copy of (may be shallow) const SPEX_options *option ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_matrix_copy| makes a deep copy of a \verb|SPEX_matrix *A| as a new \verb|SPEX_matrix *C|, which can be any of the 16 matrix formats discussed in Section \ref{ss:SPEX_matrix}. That is, the new matrix \verb|C| can be exactly the same as \verb|A| or any other type or kind different than \verb|A|. On input, the SPEX matrix that \verb|C_handle| points to must be \verb|NULL| and will be ignored, and \verb|A| is a valid matrix that can be potentially shallow. On output, \verb|C_handle| points to the matrix \verb|C|, which is a copy of \verb|A| of kind \verb|kind| and type \verb|type|. Results are undefined for an invalid input matrix \verb|A|. Though all matrices generated from any SPEX user-callable functions are valid, they could become invalid when user directly modifies their component(s). To check the validity of the input matrix, call \verb|SPEX_matrix_check| (Section \ref{s:user:matrix_check}). %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_matrix_nnz|: get the number of entries in a \verb|SPEX_matrix|} \subsection{\texttt{SPEX\_matrix\_nnz}: get the number of entries in a \texttt{SPEX\_matrix}} \label{s:user:matrix_nnz} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_matrix_nnz // return # of entries in A, or -1 on error ( int64_t *nnz, const SPEX_matrix *A, // matrix to query const SPEX_options *option ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_matrix_nnz| returns an integer, \verb|nnz|, which is equal to the number of entries in a \verb|SPEX_matrix *A|. For details regarding how the number of entries is obtained for different kinds of matrices, refer to Section \ref{ss:SPEX_matrix}. For any matrix with invalid dimension(s), nnz is returned as -1. %If default settings are desired, \verb|option| can be input as \verb|NULL|. %Returns \verb|SPEX_PANIC| if the SPEX working environment has not been initialized (e.g. via %\verb|SPEX_initialize|). %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_matrix_check|: check and optionally print a \verb|SPEX_matrix|} \subsection{\texttt{SPEX\_matrix\_check}: check and optionally print a \texttt{SPEX\_matrix}} \label{s:user:matrix_check} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_matrix_check // returns a SPEX status code ( const SPEX_matrix *A, // matrix to check const SPEX_options* option // defines the print level ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_matrix_check| checks the validity of a \verb|SPEX_matrix *A| in any of the 16 matrix formats discussed in Section \ref{ss:SPEX_matrix}. %the 15 different matrix types (CSC, triplet, dense) $\times$ (\verb|mpz|, %\verb|mpq|, \verb|mpfr|, \verb|int64|, \verb|fp64|). In addition, it prints the matrix and any error found with proper print level specified by \verb|option->print_level|. Specifically, \verb|SPEX_matrix_check| prints nothing for \verb|print_level|=0 (default); or just errors for \verb|print_level|=1; or errors and terse output of the matrix for \verb|print_level|=2; or errors and detailed output of the matrix for \verb|print_level|=3. %(refer to Section \ref{ss:SPEX_options} for more details). As mentioned, if default settings are desired, \verb|option| can be input as \verb|NULL|. % Returns \verb|SPEX_PANIC| if SPEX has not been initialized. %-------------------------------------------------------------------------------%\cprotect\section{\verb|SPEX_symbolic_analysis| helper functions} \label{s:spex_symbolic_analysis_helper} \section{\texttt{SPEX\_symbolic\_analysis} helper function} \label{s:spex_symbolic_analysis_helper} %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_symbolic_analysis_free|: free a symbolic analysis struct} \subsection{\texttt{SPEX\_symbolic\_analysis\_free}: free a symbolic analysis struct} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_symbolic_analysis_free ( SPEX_symbolic_analysis **S_handle, // Structure to be deleted const SPEX_options *option ); \end{verbatim} } \end{mdframed} \verb|SPEX_symbolic_analysis_free| frees the memory of the \verb|SPEX_symbolic_analysis| \verb|*S| that \verb|S_handle| points to. On output, the symbolic analysis \verb|S| is set to \verb|NULL|. %-------------------------------------------------------------------------------%\cprotect\section{\verb|SPEX_factorization| helper functions} \label{s:spex_factorization_helper} \section{\texttt{SPEX\_factorization} helper functions} \label{s:spex_factorization_helper} %------------------------------------------------------------------------------- These functions provide several utilities for a \verb|SPEX_factorization| \update{ %-------------------------------------------------------------------------------%\cprotect\subsection{\verb|SPEX_factorization_check|: check correctness of a factorizations struct} \subsection{\texttt{SPEX\_factorization\_check}: check correctness of a factorizations struct} \label{ss:spex_factorization_check} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_factorization_check ( SPEX_factorization *F, // The factorization to check const SPEX_options* option ); \end{verbatim} } \end{mdframed} \verb|SPEX_factorization_check| checks if a given factorization \verb|F| is correctly formatted. Specifically, it checks the following 5 conditions: \begin{enumerate} \item All required components of \verb|F| are present (e.g., an array was not erroneously freed) \item The sizes of all matrices match (e.g., \verb|F->L| and \verb|F->U| are square matrices of the same dimension) \item \verb|F->L| (and \verb|F->U| if exists) is correctly formatted via \verb|SPEX_matrix_check| \item \verb|F->L|, \verb|F->U|, and \verb|F->rhos| have the correct and same pivot values. Additionally, when \verb|F->updatable=true|, \verb|F->L| (and \verb|F->U| if exists) is of \verb|SPEX_DYNAMIC_CSC|, and the $i$-th pivot is the first entry in the nonzero list of $i$-th vector of \verb|L| (and \verb|U| if exists), i.e., \verb|F->L->v[i]->i[0] == F->P_perm[i]| and \verb|F->U->v[i]->i[0] ==| \verb|F->Q_perm[i]| \item Each permutation has no repeated indices and maps from $[0..n-1]$ to $[0..n-1]$; \verb|P_perm| and \verb|Pinv_perm| are mutually inverse vectors, same applied to \verb|(Q_perm,| \verb|Qinv_perm)| if exists \end{enumerate} Similar to \verb|SPEX_matrix_check|, \verb|SPEX_factorization_check| also prints values of the factorization, together with any error found, with proper print level specified by \verb|option->print_level|. Specifically, \verb|SPEX_factorization_check| prints nothing for \verb|print_level|=0 (default); or just errors for \verb|print_level|=1; or errors and terse output of the factorization for \verb|print_level|=2; or errors and detailed output of the factorization for \verb|print_level|=3. As mentioned, if default settings are desired, \verb|option| can be input as \verb|NULL|. } %-------------------------------------------------------------------------------%\cprotect\subsection{\verb|SPEX_factorization_convert|: Convert between updatable and non updatable} \update{\subsection{\texttt{SPEX\_factorization\_convert}: Convert between updatable and non-updatable} \label{ss:spex_factorization_convert} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_factorization_convert ( SPEX_factorization *F, // The factorization to be converted bool updatable, // if true, make F updatable. false: make non-updatable const SPEX_options* option // Command options ); \end{verbatim} } \end{mdframed} \verb|SPEX_factorization_convert| performs in-place conversion between updatable and non-updatable factorization as specified by the boolean input argument \verb'updatable'. If \verb|F->updatable == updatable| holds upon input, this function does nothing. Otherwise, it performs the corresponding in-place conversion. In case of any error, the returned factorization should be considered as undefined. Upon input, \verb|F->L| (and \verb|F->U| if exists) must be non-shallow \verb|SPEX_CSC SPEX_MPZ| matrix for non-updatable (static) factorization (i.e., \verb|F->updatalbe == false|), otherwise, the input format is considered as incorrect and \verb|SPEX_INCORRECT_INPUT| is returned. Likewise, \verb|F->L| (and \verb|F->U| if exists) must be \verb|SPEX_DYNAMIC_CSC SPEX_MPZ| matrix for updatable factorization. All SPEX functions output factorization in either of these two formats and non-shallow. Therefore, these input requirements can be met easily if users do not try to modify any individual component of \verb|F|. } %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_factorization_free|: Free a SPEX factorization} \label{ss:spex_factorization_free} \subsection{\texttt{SPEX\_factorization\_free}: Free a SPEX factorization} \label{ss:spex_factorization_free} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_factorization_free ( SPEX_factorization **F_handle, // Structure to be deleted const SPEX_options *option ); \end{verbatim} } \end{mdframed} \verb|SPEX_factorization_free| frees the memory of the \verb|SPEX_factorization *F| that \verb|F_handle| points to, and sets \verb|F| to \verb|NULL|. %------------------------------------------------------------------------------- \newpage \section{Misc Utilty Functions} %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- \subsection{\texttt{SPEX\_version}: Return version of the code} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_version ( int version [3], // SPEX major, minor, and sub version char date [128] // date of this version ) \end{verbatim} } \end{mdframed} \verb|SPEX_version| returns the library version and date. The \verb'version' array contains the three version numbers that are available at compile-time \verb'#define''d values: \verb'SPEX_VERSION_MAJOR', \verb'SPEX_VERSION_MINOR', and \verb'SPEX_VERSION_SUB', in that order. The \verb'SPEX_version' function allows the user application to check which version of SPEX it has been linked with. The three \verb'#define''d values allow the user application to know which version of SPEX was used at compile-time, which might not be the same version that was linked later on. The \verb'date' is the string \verb'SPEX_DATE', in the form \verb'"Mar 31, 2023"' for example. The string is null-terminated. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_determine_symmetry|: Determine if a matrix is symmetric} \subsection{\texttt{SPEX\_determine\_symmetry}: Determine if a matrix is symmetric} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_determine_symmetry ( bool *is_symmetric, // true if symmetric SPEX_matrix* A, // Input matrix to be checked for symmetry const SPEX_options* option // Command options ); \end{verbatim} } \end{mdframed} \verb|SPEX_determine_symmetry| checks if \verb|A| is pattern and numerically symmetric. It first checks for pattern symmetry. If it is pattern symmetric, it is checked for numerical symmetry. If $A$ is a symmetric matrix, \verb|is_symmetric| is returned as \verb|true|. % removed %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_scale|: Scale a matrix by a constant} %\subsection{\texttt{SPEX\_scale}: Scale a matrix by a constant} %------------------------------------------------------------------------------- %\begin{mdframed}[userdefinedwidth=\textwidth] %{\footnotesize %\begin{verbatim} % SPEX_info SPEX_scale % ( % // Output % SPEX_matrix* x, % // Input % const mpq_t scaling_num, //numerator % const mpq_t scaling_den, //denominator % const SPEX_options* option // command options % ); %\end{verbatim} %} \end{mdframed} %This function scales the matrix \verb|x| by the term \verb|scaling_num|/\verb|scaling_den|. This is %primarily used during forward and backward solve routines. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_transpose|: Transpose a csc mpz matrix} \subsection{\texttt{SPEX\_transpose}: Transpose a CSC mpz matrix} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_transpose ( SPEX_matrix **C_handle, // C = A' SPEX_matrix *A, // Matrix to be transposed const SPEX_options *option ); \end{verbatim} } \end{mdframed} \verb|SPEX_transpose| sets $C = A^T$. Currently, it is only supported if $A$ is CSC and \verb|mpz_t|. Returns \verb|SPEX_OK| if successful otherwise returns the appropriate error code. %------------------------------------------------------------------------------- \newpage \section{\texttt{SPEX\_gmp}: SPEX wrapper functions for GMP/MPFR} %------------------------------------------------------------------------------- SPEX provides a wrapper class for all GMP and MPFR functions used by SPEX. The wrapper class provides error-handling for out-of-memory conditions that are not handled by the GMP and MPFR libraries. These wrapper functions are used inside all SPEX functions, wherever any GMP or MPFR functions are used. These functions may also be called by the end-user application. Each wrapped function has the same name as its corresponding GMP/MPFR function with the added prefix \verb|SPEX_|. For example, the default GMP function \verb|mpz_mul| is changed to \verb|SPEX_mpz_mul|. Each SPEX GMP/MPFR function returns \verb|SPEX_OK| if successful or the correct error code if not. The following table gives a brief list of each currently covered SPEX GMP/MPFR function. Each function is declared in \verb|SPEX.h| and defined in \verb|SPEX/SPEX_Util/Source/SPEX_gmp.c|. % \newpage \thispagestyle{empty} \begin{SizedCenteredTabular}[\scriptsize]{|l|l|l|} \hline {\bf MPFR Function} & \verb|SPEX_MPFR| {\bf Function} & {\bf Description} \\ \hline\hline \verb|n = mpfr_asprintf(&buff, fmt, ...)| & \verb|n = SPEX_mpfr_asprintf(&buff, fmt, ...)| & Print format to allocated string \\ \hline \verb|mpfr_free_str(buff)| & \verb|SPEX_mpfr_free_str(buff)| & Free string allocated by MPFR \\ \hline \verb|mpfr_init2(x, s)| & \verb|SPEX_mpfr_init2(x, s)| & Initialize x with \verb's' bits \\ \hline \verb|mpfr_set_prec(x, s)| & \verb|SPEX_mpfr_set_prec(x, s)| & Set x to contain \verb's' bits \\ \hline \verb|mpfr_clear(x)| & \verb|SPEX_mpfr_clear(x)| & Safely free \verb|mpfr_t| value \\ \hline \verb|mpfr_set_null(x)| & \verb|SPEX_mpfr_set_null(x)| & Initialize the (pointer) contents of a \verb|mpfr_t| value \\ \hline \verb|mpfr_set(x, y, rnd)| & \verb|SPEX_mpfr_set(x, y, rnd)| & $x = y$ \\ \hline \verb|mpfr_set_d(x, y, rnd)| & \verb|SPEX_mpfr_set_d(x, y, rnd)| & $x = y$ (double) \\ \hline \verb|mpfr_set_si(x, y, rnd)| & \verb|SPEX_mpfr_set_si(x, y, rnd)| & $x = y$ (\verb|int64_t|) \\ \hline \verb|mpfr_set_q(x, y, rnd)| & \verb|SPEX_mpfr_set_q(x, y, rnd)| & $x = y$ (\verb|mpq_t|) \\ \hline \verb|mpfr_set_z(x, y, rnd)| & \verb|SPEX_mpfr_set_z(x, y, rnd)| & $x = y$ (\verb|mpz_t|) \\ \hline \verb|r = mpfr_get_z(x, y, rnd)| & \verb|SPEX_mpfr_get_z(x, y, rnd)| & (\verb|mpz_t|) $x = y$\\ \hline \verb|mpfr_get_q(x, y)| & \verb|SPEX_mpfr_get_q(x, y, rnd)| & (\verb|mpq_t|) $x = y$\\ \hline \verb|x = mpfr_get_d(y, rnd)| & \verb|SPEX_mpfr_get_d(x, y, rnd)| & (double) $x = y$\\ \hline \verb|x = mpfr_get_si(y, rnd)| & \verb|SPEX_mpfr_get_si(x, y, rnd)| & (\verb|int64_t|) $x = y$\\ \hline \verb|mpfr_mul(x, y, z, rnd)| & \verb|SPEX_mpfr_mul(x, y, z, rnd)| & $x = y*z$ (\verb|mpfr_t|) \\ \hline \verb|mpfr_mul_d(x, y, z, rnd)| & \verb|SPEX_mpfr_mul_d(x, y, z, rnd)| & $x = y*z$ (double) \\ \hline \verb|mpfr_div_d(x, y, z, rnd)| & \verb|SPEX_mpfr_div_d(x, y, z, rnd)| & $x = y/z$ (double) \\ \hline \verb|mpfr_ui_pow_ui(x, y, z, rnd)| & \verb|SPEX_mpfr_ui_pow_ui(x, y, z, rnd)| & $x = y^z$ (\verb|uint64_t|) \\ \hline \verb|sgn = mpfr_sgn(x)| & \verb|SPEX_mpfr_sgn(sgn, x)| & $sgn =\text{sgn}(x)$ \\ \hline %\verb|mpfr_log2(x, y, rnd)| % & \verb|SPEX_mpfr_log2(x, y, rnd )| % & $x = \log_2 (y)$ \\ \hline \verb|mpfr_free_cache()| & \verb|SPEX_mpfr_free_cache()| & Free all caches and pools used by \\&&MPFR internally \\ \hline \end{SizedCenteredTabular} \begin{SizedCenteredTabular}[\scriptsize]{|l|l|l|} \hline {\bf GMP Function} & \verb|SPEX_GMP| {\bf Function} & {\bf Description} \\ \hline\hline \verb|n = gmp_fscanf(fp, fmt, ...)| & \verb|n = SPEX_gmp_fscanf(fp, fmt, ...)| & Read from file fp \\ \hline \verb|mpz_init(x)| & \verb|SPEX_mpz_init(x)| & Initialize x \\ \hline \verb|mpz_init2(x, size)| & \verb|SPEX_mpz_init2(x, size)| & Initialize x to size bits \\ \hline \verb|mpz_clear(x)| & \verb|SPEX_mpz_clear(x)| & Safely free \verb|mpz_t| value \\ \hline \verb|mpz_set(x, y)| & \verb|SPEX_mpz_set(x, y)| & $x = y$ (\verb|mpz_t|) \\ \hline \verb|mpz_set_null(x)| & \verb|SPEX_mpz_set_null(x)| & Initialize the (pointer) contents of a \verb|mpz_t| value \\ \hline \verb|mpz_set_ui(x, y)| & \verb|SPEX_mpz_set_ui(x, y)| & $x = y$ (\verb|uint64_t|) \\ \hline \verb|mpz_set_si(x, y)| & \verb|SPEX_mpz_set_si(x, y)| & $x = y$ (\verb|int64_t|) \\ \hline %\verb|mpz_swap(x, y)| % & \verb|SPEX_mpz_swap(x, y)| % & Swap the values of $x$ and $y$ \\ \hline %\verb|mpz_set_d(x, y)| % & \verb|SPEX_mpz_set_d(x, y)| % & $x = y$ (double)\\ \hline %\verb|mpz_set_q(x, y)| % & \verb|SPEX_mpz_set_q(x, y)| % & $x = y$ (\verb|mpz_t|) \\ \hline \verb|x = mpz_get_d(y)| & \verb|SPEX_mpz_get_d(x, y)| & (double) $x = y$\\ \hline \verb|x = mpz_get_si(y)| & \verb|SPEX_mpz_get_si(x, y)| & (\verb|int64_t|) $x = y$ \\ \hline \verb|mpz_mul(x, y, z)| & \verb|SPEX_mpz_mul(x, y, z)| & $x = y*z$ \\ \hline \verb|mpz_mul_si(x, y, z)| & \verb|SPEX_mpz_mul(x, y, z)| & $x = y*z (\verb|int64_t|)$ \\ \hline %\verb|mpz_add(x, y, z)| % & \verb|SPEX_mpz_add(x, y, z)| % & $x = y+z$ \\ \hline %\verb|mpz_addmul(x, y, z)| % & \verb|SPEX_mpz_addmul(x, y, z)| % & $x = x+y*z$ \\ \hline \verb|mpz_sub(x, y, z)| & \verb|SPEX_mpz_sub(x, y, z)| & $x = y-z$ \\ \hline \verb|mpz_submul(x, y, z)| & \verb|SPEX_mpz_submul(x, y, z)| & $x = x-y*z$ \\ \hline %\verb|mpz_fdiv_q(q, x, y)| % & \verb|SPEX_mpz_fdiv_q(q, x, y)| % & $q = \text{floor}(x/y)$ \\ \hline %\verb|mpz_cdiv_q(q, x, y)| % & \verb|SPEX_mpz_cdiv_q(q, x, y)| % & $q = \text{ceil}(x/y)$ \\ \hline \verb|mpz_cdiv_qr(q, r, x, y)| & \verb|SPEX_mpz_cdiv_qr(q, r, x, y)| & $q = \text{ceil}(x/y), r = x-q*y$ \\ \hline \verb|mpz_divexact(x, y, z)| & \verb|SPEX_mpz_divexact(x, y, z)| & $x = y/z$ \\ \hline \verb|gcd = mpz_gcd(x, y)| & \verb|SPEX_mpz_gcd(gcd, x, y)| & $gcd = \text{gcd}(x,y)$\\ \hline \verb|lcm = mpz_lcm(x, y)| & \verb|SPEX_mpz_lcm(lcm, x, y)| & $lcm = \text{lcm}(x,y)$ \\ \hline \verb|mpz_neg(x, y)| & \verb|SPEX_mpz_neg(x, y)| & $x = -y$ \\ \hline \verb|mpz_abs(x, y)| & \verb|SPEX_mpz_abs(x, y)| & $x = |y|$ \\ \hline \verb|r = mpz_cmp(x, y)| & \verb|SPEX_mpz_cmp(r, x, y)| & $r = \text{sgn}(x-y)$ \\ \hline %\verb|r = mpz_cmpabs(x, y)| % & \verb|SPEX_mpz_cmpabs(r, x, y)| % & $r = \text{sgn}(|x|-|y|)$\\ \hline \verb|r = mpz_cmp_ui(x, y)| & \verb|SPEX_mpz_cmp_ui(r, x, y)| & $r = \text{sgn}(x-y)$ (\verb|uint64_t|) \\ \hline \verb|r = mpz_cmpabs_ui(x, y)| & \verb|SPEX_mpz_cmpabs_ui(r, x, y)| & $r = \text{sgn}(|x|-|y|)$ (\verb|uint64_t|) \\ \hline \verb|sgn = mpz_sgn(x)| & \verb|SPEX_mpz_sgn(sgn, x)| & $sgn = \text{sgn}(x)$ \\ \hline \verb|size = mpz_sizeinbase(x, base)| & \verb|SPEX_mpz_sizeinbase(size, x, base)| & size of x in base \\ \hline \verb|mpq_init(x)| & \verb|SPEX_mpq_init(x)| & Initialize x \\ \hline \verb|mpq_set_null(x)| & \verb|SPEX_mpq_set_null(x)| & Initialize the (pointer) contents of a \verb|mpq_t| value \\ \hline \verb|mpq_clear(x)| & \verb|SPEX_mpq_clear(x)| & Safely free \verb|mpq_t| value \\ \hline \verb|mpq_set(x, y)| & \verb|SPEX_mpq_set(x, y)| & $x = y$ \\ \hline \verb|mpq_set_z(x, y)| & \verb|SPEX_mpq_set_z(x, y)| & $x = y$ (\verb|mpz|) \\ \hline \verb|mpq_set_d(x, y)| & \verb|SPEX_mpq_set_d(x, y)| & $x=y$ (double) \\ \hline \verb|mpq_set_ui(x, y, z)| & \verb|SPEX_mpq_set_ui(x, y, z)| & $x = y/z$ (\verb|uint64_t|/\verb|uint64_t|) \\ \hline \verb|mpq_set_si(x, y, z)| & \verb|SPEX_mpq_set_si(x, y, z)| & $x = y/z$ (\verb|int64_t|/\verb|uint64_t|) \\ \hline \verb|mpq_set_num(x, y)| & \verb|SPEX_mpq_set_num(x, y)| & $num(x) = y$ \\ \hline \verb|mpq_set_den(x, y)| & \verb|SPEX_mpq_set_den(x, y)| & $den(x) = y$ \\ \hline %\verb|mpq_canonicalize(x)| % & \verb|SPEX_mpq_canonicalize(x)| % & Canonicalize x \\ \hline %\verb|mpq_get_den(x, y)| % & \verb|SPEX_mpq_get_den(x, y)| % & $x = den(y)$ \\ \hline \verb|x = mpq_get_d(y)| & \verb|SPEX_mpq_get_d(x, y)| & (double) $x = y$ \\ \hline \verb|mpq_neg(x, y)| & \verb|SPEX_mpq_neg(x, y)| & $x = -y$ \\ \hline \verb|mpq_abs(x, y)| & \verb|SPEX_mpq_abs(x, y)| & $x = |y|$ \\ \hline \verb|mpq_add(x, y, z)| & \verb|SPEX_mpq_add(x, y, z)| & $x = y+z$ \\ \hline \verb|mpq_mul(x, y, z)| & \verb|SPEX_mpq_mul(x, y, z)| & $x = y*z$ \\ \hline \verb|mpq_div(x, y, z)| & \verb|SPEX_mpq_div(x, y, z)| & $x = y/z$ \\ \hline \verb|r = mpq_cmp(x, y)| & \verb|SPEX_mpq_cmp(r, x, y)| & $r = \text{sgn}(x-y)$ \\ \hline \verb|r = mpq_cmp_ui(x, n, d)| & \verb|SPEX_mpq_cmp_ui(r, x, n, d)| & $r = \text{sgn}(x-n/d)$ (\verb|uint64_t|/\verb|uint64_t|) \\ \hline \verb|sgn = mpq_sgn(x)| & \verb|SPEX_mpq_sgn(sgn, x)| & $sgn = \text{sgn}(x)$ \\ \hline \verb|r = mpq_equal(x, y)| & \verb|SPEX_mpq_equal(r, x, y)| & $r \neq 0$ if $x=y$, $r= 0$ if $x\neq y$ \\ \hline \end{SizedCenteredTabular} If additional GMP and MPFR functions are needed in the end-user application, this wrapper mechanism can be extended to those functions, which \ul{requires user to edit the source files} of the SPEX library, (i.e., both \verb|SPEX.h| and \verb|SPEX_gmp.c|). Below are instructions on how to do this. Given a GMP function \verb|void gmpfunc(TYPEa a, TYPEb b, ...)|, where \verb|TYPEa| and \verb|TYPEb| can be GMP type data (\verb|mpz_t|, \verb|mpq_t| and \verb|mpfr_t|, for example) or non-GMP type data (\verb|int|, \verb|double|, for example), and they need not to be the same. \pagebreak A wrapper for a new GMP or MPFR function can be created by following this outline: %\newpage \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_gmpfunc ( TYPEa a, TYPEb b, ... ) { // Start the GMP Wrappter // uncomment one of the following: // If this function is not modifying any GMP/MPFR type variable, use //SPEX_GMP_WRAPPER_START; // If this function is modifying mpz_t type (say TYPEa = mpz_t), use //SPEX_GMPZ_WRAPPER_START(a) ; // If this function is modifying two variables of mpz_t type (say // TYPEa = mpz_t, TYPEb = mpz_t), use //SPEX_GMPZ_WRAPPER_START2(a, b) ; // If this function is modifying mpq_t type (say TYPEa = mpq_t), use //SPEX_GMPQ_WRAPPER_START(a) ; // If this function is modifying mpfr_t type (say TYPEa = mpfr_t), use //SPEX_GMPFR_WRAPPER_START(a) ; // Call the GMP function gmpfunc(a,b,...) ; //Finish the wrapper and return ok if successful. SPEX_GMP_WRAPPER_FINISH; return SPEX_OK; } \end{verbatim} } \end{mdframed} \newpage Note that, other than \verb|SPEX_mpfr_fprintf|, \verb|SPEX_gmp_fprintf|, \verb|SPEX_gmp_printf| and \verb|SPEX_gmp_fscanf|, all of the wrapped GMP/MPFR functions always return \verb|SPEX_info| to the caller. Therefore, for some GMP/MPFR functions that have their own return value. For example, for \verb|int mpq_cmp(const mpq_t a, const mpq_t b)|, the return value becomes a parameter of the wrapped function. In general, a GMP/MPFR function in the form of \verb|TYPEr gmpfunc(TYPEa a, TYPEb b, ...)|, the wrapped function can be constructed as follows: \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_gmpfunc ( TYPEr *r, // return value of the GMP/MPFR function TYPEa a, TYPEb b, ... ) { // Start the GMP Wrappter //SPEX_GMP_WRAPPER_START; // Call the GMP function *r = gmpfunc(a,b,...) ; //Finish the wrapper and return ok if successful. SPEX_GMP_WRAPPER_FINISH; return SPEX_OK; } \end{verbatim} } \end{mdframed} \section{SPEX Helper Macros} \label{ss:SPEX_helper_macros} In addition to the functionality described in this section; SPEX offers several helper macros to increase ease for the end user application. The first two macros are a simple try/catch mechanism which can be used to wrap functions for error handling. The next two give an easy interface to access entries $(i,j)$ in a matrix. \subsection{\texttt{SPEX\_TRY} and \texttt{SPEX\_CATCH}} In a robust application, the return values from SPEX should be checked and properly handled in the case an error occurs. SPEX is written in C and thus it cannot rely on the try/catch mechanism of C++. Thus, \verb|SPEX_TRY| and \verb|SPEX_CHECK| aim to achieve this goal. We provide \verb|SPEX_TRY| and leave \verb|SPEX_CATCH| to the user to define. \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} #define SPEX_TRY(method) \ { \ SPEX_info info = (method) ; \ if (info != SPEX_OK) \ { \ SPEX_CATCH (info) ; \ } \ } \end{verbatim} } \end{mdframed} An example definition of a \verb|SPEX_CATCH| is below. This example assumes that the user needs to free a matrix and return an error code. \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} #define SPEX_CATCH(info) \ { \ SPEX_matrix_free (&A, NULL) ; \ fprintf (stderr, "SPEX failed: info %d, \ line %d, file %s\n", \ info, __LINE__, __FILE__) ; \ return (info) ; \ } \end{verbatim} } \end{mdframed} With this mechanism, the user can safely wrap any SPEX function which returns \verb|SPEX_info| with \verb|SPEX_TRY|. For example, one can wrap. \subsection{\texttt{SPEX\_1D}: Access matrix entries with 1D linear indexing.} \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} #define SPEX_1D(A,k,type) ((A)->x.type [k]) \end{verbatim} } \end{mdframed} This allows the $k$th entry of a matrix stored in any kind (CSC, triplet, dense) of any type (mpq, mpz, int64, double, int) to be returned. For example, to return the $n$th entry of a CSC matrix with \verb|mpz_t| data types, one would use \verb|SPEX_1D(A, n, mpz)|. \subsection{\texttt{SPEX\_2D}: Access dense matrix with 2D indexing.} \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} #define SPEX_2D(A,i,j,type) SPEX_1D (A, (i)+(j)*((A)->m), type) \end{verbatim} } \end{mdframed} This allows the $(i,j)$ entry of a dense matrix of any type (mpq, mpz, int64, double, int). For example to return the $(m,n)$ entry of a dense matrix with \verb|mpq_t| data types, one would use \verb|SPEX_2D(A, m, n, mpq)|. %------------------------------------------------------------------------------- \chapter{SPEX LU}\vspace{-0.75in} \label{ch:LeftLU} %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- \section{Overview} \label{s:LeftLU:intro} %------------------------------------------------------------------------------- SPEX LU is a software package designed to exactly solve unsymmetric sparse linear systems, $ A x = b$, where $A \in \mathbb{Q}^{n \times n}$, $b \in \mathbb{Q}^{n \times r}$, and $x \in \mathbb{Q}^{n \times r}$. This package performs a left-looking, roundoff-error-free (REF) LU factorization $P A Q = L D U$, where $L$ and $U$ are integer, $D$ is diagonal, and $P$ and $Q$ are row and column permutations, respectively. Note that, in order to solve a linear system, the matrix $D$ is never explicitly computed nor needed; thus this package uses only the matrices $L$ and $U$. The theory associated with this code is the Sparse Left-looking Integer-Preserving (SLIP) LU factorization \cite{lourenco2019exact}. Aside from solving sparse linear systems exactly, one of the key goals of this package is to provide a framework for other solvers to benchmark the reliability and stability of their linear solvers, as our final solution vector $x$ is \ul{guaranteed} to be exact. SPEX LU is written in ANSI C and is accompanied by a MATLAB interface. Version 1.1.2 of SPEX Left LU was published in ACM TOMS as: Lourenco, C., Chen, J., Moreno-Centeno, E., \& Davis, T. A. (2022). Algorithm 1021: SPEX Left LU, Exactly Solving Sparse Linear Systems via a Sparse Left-looking Integer-preserving LU Factorization. ACM Transactions on Mathematical Software (TOMS), 48(2), 1-23. %------------------------------------------------------------------------------- \section{Licensing} \label{s:LeftLU:licensing} %------------------------------------------------------------------------------- \textbf{Copyright:} The copyright of this software is held by Christopher Lourenco, Jinhao Chen, Erick Moreno-Centeno, and Timothy A. Davis.\\ \noindent \textbf{Contact Info:} Contact Chris Lourenco, \href{mailto:chrisjlourenco@gmail.com}{chrisjlourenco@gmail.com}, or Tim Davis, \href{mailto:timdavis@aldenmath.com}{timdavis@aldenmath.com}, \href{mailto:davis@tamu.edu}{davis@tamu.edu}, or \href{DrTimothyAldenDavis@gmail.com}{DrTimothyAldenDavis@gmail.com}\\ \noindent \textbf{License:} This software package is dual licensed under the GNU General Public License version 2 or the GNU Lesser General Public License version 3. Details of this license are in \verb|SPEX/License/license.txt|. For alternative licenses, please contact the authors. \section{Factorization and Solve Routines} To factorize and solve a linear system $A \mathbf{x} = \mathbf{b}$ via the SPEX Left LU factorization, a user must call analyze, factorize, and solve. The functions are explained below: %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_LU_analyze|: symbolic analysis for LU factorization} \label{ss:spex_lu_analyze} \subsection{\texttt{SPEX\_lu\_analyze}: symbolic analysis for LU factorization} \label{ss:spex_lu_analyze} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_lu_analyze ( SPEX_symbolic_analysis** S_handle, // symbolic analysis including // column perm. and nnz of L and U const SPEX_matrix *A, // Input matrix const SPEX_options *option // Control parameters, if NULL, use default ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_lu_analyze| performs symbolic analysis for the REF LU factorization. On input, the \verb|SPEX_symbolic_analysis *S| that \verb|S_handle| points to is undefined; \verb|A| must be a square matrix of \verb|SPEX_CSC| kind; and \verb|option| contains any command parameters (default settings are used if \verb|option| is input as \verb|NULL|). On output, \verb|S| contains the column preordering of \verb|A| and estimates on the number of nonzeros in $L$ and $U$. The type of ordering can be chosen with \verb|option->order|. It is suggested that COLAMD is used. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Left_LU_factorize|: Compute the LU factorization of A} \label{ss:spex_left_lu_factorize} \subsection{\texttt{SPEX\_lu\_factorize}: Compute the LU factorization of A} \label{ss:spex_left_lu_factorize} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_lu_factorize ( // output: SPEX_factorization **F_handle, // LU factorization // input: const SPEX_matrix *A, // matrix to be factored const SPEX_symbolic_analysis *S, // symbolic analysis const SPEX_options* option // command options ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_lu_factorize| performs the left-looking LU factorization. On input, the \linebreak \verb|SPEX_factorization *F| that \verb|F_handle| points to is undefined; \verb|A| must be a square matrix of \verb|SPEX_CSC SPEX_MPZ| format; \verb|S| is obtained from \verb|SPEX_lu_analyze| that contains the column ordering of \verb|A|; and \verb|option| contains any command parameters (default settings are used if \verb|option| is input as \verb|NULL|). On output, \verb|A|, \verb|S|, and \verb|option| are unmodified and \verb|F| contains the REF LU factorization of \verb|A|. %Returns \verb|SPEX_PANIC| if SPEX has not been initialized. Otherwise, if another If any error occurs, \verb|F| is returned as \verb|NULL|, and an appropriate error code is returned. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Left_LU_solve|: solve the linear system $Ax=b$} \newpage \subsection{\texttt{SPEX\_lu\_solve}: solve the linear system} \label{ss:SPEX_Left_LU_solve} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_lu_solve // solves the linear system LD^(-1)U x = b ( // Output SPEX_matrix **x_handle, // rational solution to the system // input/output: SPEX_factorization *F, // The LU factorization. // input: const SPEX_matrix *b, // right hand side vector const SPEX_options* option // Command options ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_lu_solve| obtains the solution of \verb|mpq_t| type to the linear system $Ax=b$ upon a successful factorization. This function may be called after a successful return from \linebreak \verb|SPEX_lu_factorize|. On input, \verb|SPEX_matrix *x| that \verb|x_handle| points to is undefined; \verb|F| must be a valid LU factorization; and \verb|b| must be a dense \verb|mpz_t| matrix with same number of rows as \verb|F->L|; Default settings are used if \verb|option| is input as \verb|NULL|. Upon successful completion, the function returns \verb|SPEX_OK|, and \verb|x| contains the solution of \verb|mpq_t| type with dense format to the linear system $Ax=b$. \update{ \verb|F| is mathematically unchanged on output. However, if \verb|F| is updatable on input, it is converted to non-updatable. If \verb|F| is already non-updatable, it is not modified.} In case of failure, \verb|x| is returned as \verb|NULL| and the appropriate error code is returned. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Left_LU_backslash|: solve $Ax=b$} \subsection{\texttt{SPEX\_lu\_backslash}: solve a linear system} \label{ss:SPEX_Left_LU_backslash} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_lu_backslash ( // Output SPEX_matrix **X_handle, // Final solution vector // Input SPEX_type type, // Type of output desired. Must be // SPEX_FP64, SPEX_MPFR, or SPEX_MPQ const SPEX_matrix* A, // Input matrix of SPEX_CSC SPEX_MPZ const SPEX_matrix* b, // Right hand side vector(s). Must be // SPEX_DENSE SPEX_MPZ const SPEX_options* option // Command options (Default if NULL) ) ; \end{verbatim} } \end{mdframed} \verb|SPEX_lu_backslash| solves the linear system $Ax=b$ and returns the solution as a dense matrix of \verb|mpq_t|, \verb|mpfr_t| or \verb|double| entries. This function performs symbolic analysis, factorization, and solving all in one line. It can be thought of as an exact version of MATLAB sparse backslash. On input, \verb|SPEX_matrix *x| that \verb|X_handle| points to is undefined. \verb|type| must be one of: \verb|SPEX_MPQ|, \verb|SPEX_MPFR| or \verb|SPEX_FP64| to specify the data type of the solution entries. \verb|A| should be a square CSC \verb|mpz_t| matrix while \verb|b| should be a dense \verb|mpz_t| matrix. In addition, \verb|A->m| should be equal to \verb|b->m|. Default settings are used if \verb|option| is input as \verb|NULL|. Upon successful completion, the function returns \verb|SPEX_OK|, and \verb|x| contains the solution of data type specified by \verb|type| to the linear system $Ax=b$. In case of failure, \verb|x| is returned as \verb|NULL| and the appropriate error code is returned. %Returns \verb|SPEX_PANIC| if SPEX has not been initialized. \begin{comment} %------------------------------------------------------------------------------- \cprotect\section{Using SPEX Left LU in C} \label{s:Using} %------------------------------------------------------------------------------- Using SPEX Left LU in C has three steps: \begin{enumerate} \item initialize and populate data structures, \item perform symbolic analysis, factorize the matrix $A$ and solve the linear system for each $b$ vector, and \item free all used memory and finalize. \end{enumerate} Step 1 is discussed in Section \ref{s:Using:init}. For Step 2, performing symbolic analysis and factorizing $A$ and solving the linear $A x =b$ can be done in one of two ways. If only the solution vector $x$ is required, SPEX Left LU provides a simple interface for this purpose which is discussed in Section \ref{s:Using:simple}. Alternatively, if the $L$ and $U$ factors are required, refer to Section \ref{s:Using:expert}. Finally, step 3 is discussed in Section \ref{s:Using:free}. For the remainder of this section, \verb|n| will indicate the dimension of $A$ (that is, $A \in \mathbb{Z}^{n \times n}$) and \verb|numRHS| will indicate the number of right hand side vectors being solved (that is, if \verb|numRHS|$= r$, then $b \in \mathbb{Z}^{n \times r}$). %------------------------------------------------------------------------------- \cprotect\subsection{SPEX Left LU initialization and population of data structures} \label{s:Using:init} %------------------------------------------------------------------------------- This section discusses how to initialize and populate the global data structures required for SPEX Left LU. %------------------------------------------------------------------------------- \subsubsection{Initializing the environment} %------------------------------------------------------------------------------- SPEX is built upon the GNU Multiple Precision Arithmetic (GMP) \cite{granlund2015gnu} and GNU Multiple Precision Floating Point Reliable (MPFR) \cite{fousse2007mpfr} libraries and provides wrappers to all GMP/MPFR functions it uses. This allows SPEX to properly handle memory management failures, which GMP/MPFR does not handle. To enable this mechanism, SPEX requires initialization. The following must be done before using any other SPEX function: \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_initialize ( ) ; // or SPEX_initialize_expert (...); if custom memory functions are desired \end{verbatim} } \end{mdframed} %------------------------------------------------------------------------------- \subsubsection{Initializing data structures} \label{ss:init} %------------------------------------------------------------------------------- SPEX assumes three specific input options for all functions. These are: \begin{itemize} \item \verb|SPEX_matrix* A| and \verb|SPEX_matrix *b|: \verb|A| contains the input coefficient matrix, while \verb|b| contains the right hand side vector(s) of the linear system $Ax=b$. \item \verb|SPEX_LU_analysis* S|: \verb|S| contains the column permutation used for $A$ as well as estimates of the number of nonzeros in $L$ and $U$. \item \verb|SPEX_options* option|: \verb|option| contains various control options for the factorization including column ordering used, pivot selection scheme, and others. For a full list of the contents of the \verb|SPEX_options| structure, refer to Section \ref{ss:SPEX_options}. If default settings are desired, \verb|option| can be set to \verb|NULL|. \end{itemize} %------------------------------------------------------------------------------- \subsubsection{Populating data structures} \label{ss:populate_Ab} %------------------------------------------------------------------------------- Of the three data structures discussed in Section~\ref{ss:init}, \verb|S| is constructed during symbolic analysis (Section \ref{s:SPEX_LU_analyze}), and \verb|option| is an optional parameter for selecting non-default parameters. Refer to Section \ref{ss:SPEX_options} for the contents of \verb|option|. SPEX allows the input numerical data for \verb|A| and \verb|b| to come in one of 5 types: \verb|int64_t|, \verb|double|, \verb|mpfr_t|, \verb|mpq_t|, and \verb|mpz_t|. Moreover, both \verb|A| and \verb|b| can be stored in CSC form, sparse triplet form or dense form. CSC form is discussed in Section \ref{s:util:overview}. The triplet form stores the contents of the matrix $A$ in three arrays \verb|i|, \verb|j|, and \verb|x| where the $k$th nonzero entry is stored as $A ( i[k], j[k]) = x[k]$. SPEX stores its dense matrices in in column-oriented format, that is, the $(i,j)$th entry in \verb|A| is \verb|A->x.TYPE[p]| with $p = i+j$*\verb|A->m|. If the data for matrices are in file format to be read, refer to \newline \verb|SPEX/SPEX/SPEX_Left_LU/Demo| \verb|/example2.c| on how to read in data and construct \verb|A| and \verb|b|. If the data for matrices are already stored in vectors corresponding to CSC form, sparse triplet form or dense form, allocate a shallow \verb|SPEX_matrix| and assign vectors accordingly, then use \verb|SPEX_matrix_copy| to get a \verb|SPEX_matrix| in the desired kind and type. For more details, refer to \verb|SPEX/SPEX/SPEX_Left_LU/Demo/example.c|. In a case when \verb|A| is available in format other than CSC \verb|mpz|, and/or \verb|b| is available in format other than dense \verb|mpz|, the following code snippet shows how to get \verb|A| and \verb|b| in a required format. {\small \begin{verbatim} /* Get the matrix A. Assume that A1 is stored in CSC form with mpfr_t entries, while b1 is stored in triplet form with mpq_t entries. (for A1 and b1 in any other form, the exact same code will work) */ SPEX_matrix *A, *b; // A is a copy of the A1. A is a CSC matrix with mpz_t entries SPEX_matrix_copy(&A, SPEX_CSC, SPEX_MPZ, A1, option); // b is a copy of the b1. b is a dense matrix with mpz_t entries. SPEX_matrix_copy(&b, SPEX_DENSE, SPEX_MPZ, b1, option); \end{verbatim} } %------------------------------------------------------------------------------- \cprotect\subsection{Simple SPEX Left LU routines for solving linear systems} \label{s:Using:simple} %------------------------------------------------------------------------------- After initializing the necessary data structures, SPEX obtains the solution to $Ax=b$ using the simple interface of SPEX Left LU, \verb|SPEX_Left_LU_backslash|. The \newline \verb|SPEX_Left_LU_backslash| function can return \verb|x| as \verb|double|, \verb|mpq_t|, or \verb|mpfr_t| with an associated precision. See Section \ref{ss:SPEX_Left_LU_backslash} for more details. The following code snippet shows how to get solution as a dense \verb|mpq_t| matrix. {\small \begin{verbatim} SPEX_matrix *x; SPEX_type my_type = SPEX_MPQ; // SPEX_MPQ, SPEX_MPFR, SPEX_FP64 SPEX_Left_LU_backslash(&x, my_type, A, b, option) ; \end{verbatim} } On successful return, this function returns \verb|SPEX_OK| (see Section \ref{ss:SPEX_info}). %------------------------------------------------------------------------------- \cprotect\subsection{Expert SPEX Left LU routines} \label{s:Using:expert} %------------------------------------------------------------------------------- If the $L$ and $U$ factors from the SPEX Left LU factorization of the matrix $A$ are required, the steps performed by \verb|SPEX_Left_LU_backslash| can be done with a sequence of calls to SPEX functions: \begin{enumerate} \item declare \verb|L|, \verb|U|, the solution matrix \verb|x|, and others, \item perform symbolic analysis, \item compute the factorization $PAQ = L D U$, \item solve the linear system $Ax =b$, and \item convert the final solution into the final desired form. \end{enumerate} \noindent These steps are discussed below, along with examples. %------------------------------------------------------------------------------- \subsubsection{Declare workspace} %------------------------------------------------------------------------------- Using SPEX in this form requires the intermediate variables be declared, such as \verb|L|, \verb|U|, etc. The following code snippet shows the detailed list. {\small \begin{verbatim} // A and b are in required type and ready to use SPEX_matrix *L = NULL; SPEX_matrix *U = NULL; SPEX_matrix *x = NULL; SPEX_matrix *rhos = NULL; int64_t* pinv = NULL; SPEX_LU_analysis* S = NULL; // option needs no declaration if default setting is desired // only declare option for further modification on default setting SPEX_options *option = SPEX_create_default_options(); \end{verbatim} } %------------------------------------------------------------------------------- \subsubsection{SPEX Left LU symbolic analysis} %------------------------------------------------------------------------------- The symbolic analysis phase of an LU factorization entails computing the symbolic column ordering and estimating the number of nonzeros in $L$ and $U$. This is performed by calling the following function: {\small \begin{verbatim} SPEX_LU_analyze (&S, A, option) ; \end{verbatim} } %------------------------------------------------------------------------------- \subsubsection{Computing the factorization} %------------------------------------------------------------------------------- The matrices \verb|L| and \verb|U|, the pivot sequence \verb|rhos|, and the row permutation \verb|pinv| are computed via the \verb|SPEX_Left_LU_factorize| function (Section \ref{s:LeftLU:SPEX_Left_LU_factorize}). Upon successful completion, this function returns \verb|SPEX_OK|. %------------------------------------------------------------------------------- \subsubsection{Solving the linear system} %------------------------------------------------------------------------------- After factorization, the next step is to solve the linear system and store the solution as a dense matrix \verb|x| with entries of rational number \verb|mpq_t|. This solution is done via the \verb|SPEX_Left_LU_solve| function (Section \ref{ss:SPEX_Left_LU_solve}). Upon successful completion, this function returns \verb|SPEX_OK|. In this step, \verb|option->check| can be set to \verb|true| to enable the solution check process as discussed in Section \ref{ss:SPEX_Left_LU_solve}. The process can verify that the solution vector x satisfies $Ax=b$ in perfect precision intended for debugging. This step is not needed, since the solution returned is guaranteed to be exact. It appears here simply as debugging tool, and as a verification that SPEX is computing its expected result. This test can fail only if it runs out of memory, or if there is a bug in the code (in which case, please notify the authors). Also, note that this process can be quite time consuming; thus it is not recommended to be used in general. %------------------------------------------------------------------------------- \subsubsection{Converting the solution vector to the final desired form} %------------------------------------------------------------------------------- Upon completion of the above routines, the solution to the linear system is in a dense \verb|mpq_t| matrix. SPEX allows this to be converted into any form of matrix in the set of (CSC, sparse triplet, dense) $\times$ (\verb|mpfr_t|, \verb|mpq_t|, \verb|double|) using \verb|SPEX_matrix_copy|. The following code snippet shows how to get solution as a dense \verb|double| matrix; since this involves a floating-point representation, the solution \verb|my_x| will no longer be exact, even though \verb|x| is the exact solution. {\small \begin{verbatim} SPEX_kind my_kind = SPEX_DENSE; // SPEX_CSC, SPEX_TRIPLET or SPEX_DENSE SPEX_type my_type = SPEX_FP64; // SPEX_MPQ, SPEX_MPFR, or SPEX_FP64 SPEX_matrix* my_x = NULL; // New output // Create copy which is stored as my_kind and my_type: SPEX_matrix_copy( &my_x, my_kind, my_type, x, option);\end{verbatim} } %------------------------------------------------------------------------------- \cprotect\subsection{Freeing memory} \label{s:Using:free} %------------------------------------------------------------------------------- As described in Section \ref{s:user:memmanag}, SPEX provides a number of functions/macros to free SPEX structures: \begin{itemize} \item \verb|SPEX_matrix*|: A \verb|SPEX_matrix* A| data structure can be freed with a call to \verb|SPEX_matrix_free(&A, NULL) ;| \item \verb|SPEX_LU_analysis*|: A \verb|SPEX_LU_analysis* S| data structure can be freed with a call to \verb|SPEX_LU_analysis_free(&S, NULL) ;| \item All others including \verb|SPEX_options*|: These data structures can be freed with a call to the macro \verb|SPEX_FREE()|, for example, \verb|SPEX_FREE(option)| for \newline \verb|SPEX_options* option|. \end{itemize} After all usage of the SPEX routines is finished, \verb|SPEX_finalize()| must be called (Section \ref{ss:SPEX_finalize}) to finalize usage of the library. %------------------------------------------------------------------------------- \cprotect\subsection{Examples} \label{s:Using:Examples} %------------------------------------------------------------------------------- The \verb|SPEX/SPEX/SPEX_Left_LU/Demo| folder contains three sample C codes which utilize SPEX. These files demonstrate the usage of SPEX as follows: \begin{itemize} \item \verb|example.c|: This example generates a random dense $50 \times 50$ matrix and a random dense $50 \times 1$ right hand side vector $b$ and solves the linear system. In this function, the \verb|SPEX_Left_LU_backslash| function is used; and the output is given as a double matrix. \item \verb|example2.c|: This example reads in a matrix stored in triplet format from the \verb|ExampleMats| folder. Additionally, it reads in a right hand side vector from this folder and solves the associated linear system via the \verb|SPEX_Left_LU_backslash| function, and, the solution is given as a matrix of rational numbers. \item \verb|spexlu_demo.c|: This example reads in a matrix and right hand side vector from a file and solves the linear system $A x = b$ using the techniques discussed in Section \ref{s:Using:expert}. This file also allows command line arguments (discussed in \verb|README.md|) and can be used to replicate the results from \cite{lourenco2019exact}. \end{itemize} \end{comment} %------------------------------------------------------------------------------- \chapter{SPEX Cholesky}\vspace{-0.75in} \label{ch:Chol} %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- \section{Overview} \label{s:Chol:intro} %------------------------------------------------------------------------------- SPEX Cholesky is a software package designed to exactly solve symmetric positive definite linear systems, $A x = b$ where $A \in \mathbb{Q}^{n \times n}$, $b \in \mathbb{Q}^{n \times r}$, and $x \in \mathbb{Q}^{n \times r}$. This package performs either a left-looking or up-looking sparse roundoff-error-free Cholesky factorization $P A P^T = L D L^T$ where $L$ is integer, and $P$ is the symmetric permutation. Note that, in order to solve a linear system, the matrix $D$ is never explicitly computed nor needed; thus this package uses only the matrix $L$. The theory associated with this code can be found at \cite{lourenco2022exactly}. SPEX Cholesky is written in ANSI C and is accompanied by MATLAB and Python interfaces. %------------------------------------------------------------------------------- \section{Licensing} \label{s:Chol:licensing} %------------------------------------------------------------------------------- \textbf{Copyright:} The copyright of this software is held by Christopher Lourenco, Lorena Mejia Domenzain, Jinhao Chen, Erick Moreno-Centeno, and Timothy A. Davis.\\ \noindent \textbf{Contact Info:} Contact Chris Lourenco, \href{mailto:chrisjlourenco@gmail.com}{chrisjlourenco@gmail.com}, or Tim Davis, \href{mailto:timdavis@aldenmath.com}{timdavis@aldenmath.com}, \href{mailto:davis@tamu.edu}{davis@tamu.edu}, or \href{DrTimothyAldenDavis@gmail.com}{DrTimothyAldenDavis@gmail.com}\\ \noindent \textbf{License:} This software package is dual licensed under the GNU General Public License version 2 or the GNU Lesser General Public License version 3. Details of this license are in \verb|SPEX/License/license.txt|. For alternative licenses, please contact the authors. \section{Factorization and Solve Routines} To factorize and solve a linear system $A \mathbf{x} = \mathbf{b}$ via the SPEX Cholesky factorization, a user must call analyze, factorize, and solve. The functions are explained below: %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Chol_analyze|: symbolic analysis for Cholesky factorization} \label{ss:spex_chol_analyze} \newpage \subsection{\texttt{SPEX\_cholesky\_analyze}: symbolic analysis for Cholesky factorization} \label{ss:spex_chol_analyze} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_cholesky_analyze ( // Output SPEX_symbolic_analysis** S_handle, // Symbolic analysis data structure // Input const SPEX_matrix* A, // Input matrix of SPEX_CSC const SPEX_options* option // Command options (Default if NULL) ); \end{verbatim} } \end{mdframed} \verb|SPEX_cholesky_analyze| performs symbolic analysis for the REF Cholesky factorization. On input, the \verb|SPEX_symbolic_analysis *S| that \verb|S_handle| points to is undefined; \verb|A| must be an SPD matrix of \verb|SPEX_CSC| kind; and \verb|option| contains any command parameters (default settings are used if \verb|option| is input as \verb|NULL|). On output, \verb|S| contains the row and column ordering of \verb|A|, the exact number of nonzeros in $L$, the elimination tree of \verb|A|, and the column pointers of $L$. The type of ordering can be chosen with \verb|option->order|. It is suggested that AMD is used. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Chol_factorize|: Compute the Cholesky factorization of A} \label{ss:spex_chol_factorize} \subsection{\texttt{SPEX\_cholesky\_factorize}: Compute the Cholesky factorization of A} \label{ss:spex_chol_factorize} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_cholesky_factorize ( // Output SPEX_factorization **F_handle, // Cholesky factorization struct //Input const SPEX_matrix* A, // CSC MPZ Matrix to be factored const SPEX_symbolic_analysis* S,// Symbolic analysis struct from // SPEX_Chol_analyze. const SPEX_options* option // command options, option->chol_type can be // either CHOL_UP (default) or CHOL_LEFT. ); \end{verbatim} } \end{mdframed} \verb|SPEX_cholesky_factorize| performs the REF Cholesky factorization via either the up-looking (default) or left-looking manner (specified by \verb|option->chol_type|). On input, the \verb|SPEX_factorization *F| that \verb|F_handle| points to is undefined; \verb|A| must be an SPD matrix of \verb|SPEX_CSC SPEX_MPZ| format; \verb|S| is obtained from \verb|SPEX_Chol_analyze| that contains the column/row ordering of \verb|A|; and \verb|option| contains any command parameters (default settings are used if \verb|option| is input as \verb|NULL|). On output, \verb|A|, \verb|S|, and \verb|option| are unmodified and \verb|F| contains the REF Cholesky factorization of \verb|A|. %Returns \verb|SPEX_PANIC| if SPEX has not been initialized. Otherwise, if another If error occurs, \verb|F| is returned as \verb|NULL|, and an appropriate error code is returned. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Chol_solve|: solve the linear system $Ax=b$} \newpage \subsection{\texttt{SPEX\_cholesky\_solve}: solve the linear system} \label{ss:SPEX_Chol_solve} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_cholesky_solve // solves the linear system LD^(-1)L^T x = b ( // Output SPEX_matrix** x_handle, // rational solution to the system. // input/output: SPEX_factorization *F, // The non-updatable Cholesky factorization. // input: const SPEX_matrix* b, // Right hand side vector const SPEX_options* option // command options ); \end{verbatim} } \end{mdframed} \verb|SPEX_cholesky_solve| obtains the solution of \verb|mpq_t| type to the linear system $Ax=b$ upon a successful factorization. This function may be called after a successful return from \verb|SPEX_Chol_factorize|. On input, \verb|SPEX_matrix *x| that \verb|x_handle| points to is undefined; \verb|F| must be a valid Cholesky factorization and \verb|b| must be dense \verb|mpz_t| with same number of rows as \verb|F->L|. Default settings are used if \verb|option| is input as \verb|NULL|. Upon successful completion, the function returns \verb|SPEX_OK|, and \verb|x| contains the solution of \verb|mpq_t| type with dense format to the linear system $Ax=b$. \update{ \verb|F| is mathematically unchanged on output. However, if \verb|F| is updatable on input, it is converted to non-updatable. If \verb|F| is already non-updatable, it is not modified.} In case of failure, \verb|x| is returned as \verb|NULL| and the appropriate error code is returned. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Chol_backslash|: solve $Ax=b$} \subsection{\texttt{SPEX\_cholesky\_backslash}: solve a linear system} \label{ss:SPEX_Chol_backslash} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_cholesky_backslash ( // Output SPEX_matrix** x_handle, // Final solution vector(s) // Input SPEX_type type, // Type of output desired. Must be // SPEX_FP64, SPEX_MPFR, or SPEX_MPQ const SPEX_matrix* A, // Input matrix of SPEX_CSC SPEX_MPZ const SPEX_matrix* b, // Right hand side vector(s). Must be // SPEX_DENSE SPEX_MPZ const SPEX_options* option // Command options (Default if NULL) ); \end{verbatim} } \end{mdframed} \verb|SPEX_cholesky_backslash| solves the linear system $Ax=b$ and returns the solution as a dense matrix of \verb|mpq_t|, \verb|mpfr_t| or \verb|double| entries. This function performs symbolic analysis, factorization, and solving all in one line. It can be thought of as an exact version of MATLAB sparse backslash for SPD matrices. If $A$ is not SPD, this function will not work and LU factorization must be used. On input, \verb|SPEX_matrix *x| that \verb|x_handle| points to is undefined. \verb|type| must be one of: \verb|SPEX_MPQ|, \verb|SPEX_MPFR| or \verb|SPEX_FP64| to specify the data type of the solution entries. \verb|A| should be a square CSC \verb|mpz_t| matrix while \verb|b| should be a dense \verb|mpz_t| matrix. In addition, \verb|A->m| should be equal to \verb|b->m|. Default settings are used if \verb|option| is input as \verb|NULL|. Upon successful completion, the function returns \verb|SPEX_OK|, and \verb|x| contains the solution of data type specified by \verb|type| to the linear system $Ax=b$. In case of failure, \verb|x| is returned as \verb|NULL| and the appropriate error code is returned. %Returns \verb|SPEX_PANIC| if SPEX has not been initialized. %------------------------------------------------------------------------------- \chapter{SPEX Backslash}\vspace{-0.75in} \label{ch:Backslash} %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- \section{Overview} \label{s:Backslash:intro} %------------------------------------------------------------------------------- SPEX Backslash is a software package designed to exactly solve sparse linear systems, $A x = b$ where $A \in \mathbb{Q}^{n \times n}$, $b \in \mathbb{Q}^{n \times r}$, and $x \in \mathbb{Q}^{n \times r}$. This package determines the appropraite factorization to apply based on the structure of the input matrix. SPEX Backslash is written in ANSI C and is accompanied by MATLAB and Python interfaces. %------------------------------------------------------------------------------- \section{Licensing} \label{s:Backslash:licensing} %------------------------------------------------------------------------------- \textbf{Copyright:} The copyright of this software is held by Christopher Lourenco, Lorena Mejia Domenzain, Jinhao Chen, Erick Moreno-Centeno, and Timothy A. Davis.\\ \noindent \textbf{Contact Info:} Contact Chris Lourenco, \href{mailto:chrisjlourenco@gmail.com}{chrisjlourenco@gmail.com}, or Tim Davis, \href{mailto:timdavis@aldenmath.com}{timdavis@aldenmath.com}, \href{mailto:davis@tamu.edu}{davis@tamu.edu}, or \href{DrTimothyAldenDavis@gmail.com}{DrTimothyAldenDavis@gmail.com}\\ \noindent \textbf{License:} This software package is dual licensed under the GNU General Public License version 2 or the GNU Lesser General Public License version 3. Details of this license are in \verb|SPEX/License/license.txt|. For alternative licenses, please contact the authors. %------------------------------------------------------------------------------- %\cprotect\section{\verb|SPEX_Backslash|: Exactly solve sparse linear systems} \newpage \section{\texttt{SPEX\_backslash}: Exactly solve sparse linear systems} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_backslash ( // Output SPEX_matrix **X_handle, // On output: Final solution vector // On input: undefined // Input SPEX_type type, // Type of output desired. Must be // SPEX_FP64, SPEX_MPFR, or SPEX_MPQ const SPEX_matrix* A, // Input matrix of SPEX_CSC SPEX_MPZ const SPEX_matrix* b, // Right hand side vector(s). Must be // SPEX_DENSE SPEX_MPZ const SPEX_options* option // Command options (Default if NULL) ); \end{verbatim} } \end{mdframed} \verb|SPEX_backslash| exactly solves the linear system $A \mathbf{x} = \mathbf{b}$ using the appropriate factorization. On input, \verb|SPEX_matrix *x| that \verb|X_handle| points to is undefined. \verb|type| must be one of: \verb|SPEX_MPQ|, \verb|SPEX_MPFR| or \verb|SPEX_FP64| to specify the data type of the solution entries. \verb|A| should be a square CSC \verb|mpz_t| matrix while \verb|b| should be a dense \verb|mpz_t| matrix. In addition, \verb|A->m| should be equal to \verb|b->m|. Default settings are used if \verb|option| is input as \verb|NULL|. This function first checks the symmetry of \verb|A|. If \verb|A| is numerically and pattern symmetric, SPEX Cholesky factorization is attempted. If the Cholesky factorization is successful, it is used to solve $A x = b$. Otherwise, LU factorization is used. Upon successful completion, the function returns \verb|SPEX_OK|, and \verb|x| contains the solution of data type specified by \verb|type| to the linear system $Ax=b$. If an error occurs, \verb|x| is freed and the appropriate error code is returned. \update{ \chapter{SPEX Update}\vspace{-0.75in} \label{ch:Update} %------------------------------------------------------------------------------- \section{Overview} \label{s:Update:intro} %------------------------------------------------------------------------------- SPEX Update is a software package designed to exactly update an available REF LU/Cholesky factorization for the sparse matrix $A$ when a low-rank modification is performed on $A$. The low-rank modification could be either column/row replacement or rank-1 update/downdate. Currently, SPEX provides functions for column-replacement LU update and rank-1 Cholesky update/downdate. The theoretical basis for this code is to be submitted as \begin{itemize} \item J. Chen, T. A. Davis, “Sparse Exact LU Update for Column Replacement,” SIAM Journal on Matrix Analysis and Applications. \item J. Chen, T. A. Davis, “Sparse Exact Rank-1 Cholesky Update,” SIAM Journal on Matrix Analysis and Applications. \end{itemize} SPEX Update is written in ANSI C. %------------------------------------------------------------------------------- \section{Licensing} \label{s:Update:licensing} %------------------------------------------------------------------------------- \textbf{Copyright:} The copyright of this software is held by Jinhao Chen, Christopher Lourenco, Lorena Mejia Domenzain, Erick Moreno-Centeno, and Timothy A. Davis.\\ \noindent \textbf{Contact Info:} Contact Jinhao Chen \href{mailto:cjh10644@hotmail.com}{cjh10644@hotmail.com}, or Chris Lourenco, \href{mailto:chrisjlourenco@gmail.com}{chrisjlourenco@gmail.com}, or Tim Davis, \href{mailto:timdavis@aldenmath.com}{timdavis@aldenmath.com}, \href{mailto:davis@tamu.edu}{davis@tamu.edu}, or \href{DrTimothyAldenDavis@gmail.com}{DrTimothyAldenDavis@gmail.com}\\ \noindent \textbf{License:} This software package is dual licensed under the GNU General Public License version 2 or the GNU Lesser General Public License version 3. Details of this license are in \verb|SPEX/License/license.txt|. For alternative licenses, please contact the authors. %------------------------------------------------------------------------------- \section{SPEX Factorization Update Functions} %------------------------------------------------------------------------------- SPEX update provides the following functions to perform factorization updates. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Update_Chol_Rank1|: rank 1 update to REF Cholesky} \subsection{\texttt{SPEX\_Update\_Chol\_Rank1}: rank-1 REF Cholesky update/downdate}\label{ss:SPEX_Update_Chol_Rank1} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_Update_Chol_Rank1 ( // input/output SPEX_factorization *F, // The SPEX Cholesky factorization of A. SPEX_matrix *w, // a n-by-1 SPEX_DYNAMIC_CSC matrix that // contains the vector to modify the original // matrix A, the resulting A is A+sigma*w*w^T. // input const int64_t sigma, // a nonzero scalar that determines whether // this is an update (sigma > 0) or downdate // (sigma < 0). const SPEX_options* option // Command options ); \end{verbatim} } \end{mdframed} \verb|SPEX_Update_Chol_Rank1| performs the rank-1 REF Cholesky update/downdate for a $n$-by-$n$ sparse matrix $A$. On input, \verb|F| contains the valid Cholesky factorization to be updated that consists of \verb|L|, \verb|rhos|, \verb|P_perm| and \verb|Pinv_perm|; \verb|w| must be a $n$-by-$1$ matrix of \verb|SPEX_DYNAMIC_CSC| kind; \verb|sigma| is a nonzero scalar, where \verb|sigma > 0| indicates update and \verb|sigma < 0| indicates downdate; and \verb|option| contains command options (default settings are used if \verb|option| is input as \verb|NULL|). This function requires that the rows of \verb|w| are in the same order as that of $A$ (and \verb|F->L| if \verb|F| is updatable). In addition, \verb|w| and \verb|SPEX_matrix *A| that contains $A$ are scaled with same factor, i.e., \verb|w->scale = A->scale|, while \verb|w->v[0]->scale = 1|. (Users can refer to section \ref{ss:SPEX_matrix} for the difference between \verb|w->scale| and \verb|w->v[0]->scale|.) As for \verb|F|, if it is input in a non-updatable format, it is first converted to updatable format. On success \verb|SPEX_OK| is returned, and \verb|F| contains the updatable REF Cholesky factorization of $A + \sigma w w^T$, and \verb|w| contains the solution $w_{out}$ to $LD^{-1}w_{out} = w_{in}$, where $A=LD^{-1}L^T$ and $w_{in}$ is the input matrix \verb|w|. It should be noted that $w_{out}$ is also the solution to $\bar{L}\bar{D}^{-1}w_{out} = w_{in}$, where $A+ \sigma w w^T=\bar{L}\bar{D}^{-1}\bar{L}^T$. Otherwise, if this function fails for any reason, the appropriate error code is returned and \verb|F| is undefined on output (since the modification for the factorization is done in place during the update process). This function does not require/use/modify the original matrix $A$. Therefore, if the updated matrix $\bar{A}$ is needed, users should implement it by their own {\bf BEFORE} using this function (since \verb|w| will be modified by this function). %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Update_LU_ColRep|: Column replacement update for REF LU} \subsection{\texttt{SPEX\_Update\_LU\_ColRep}: Column-replacement REF LU update}\label{ss:SPEX_Update_LU_ColRep} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_Update_LU_ColRep ( // input/output SPEX_factorization* F, // The SPEX LU factorization of A. // input SPEX_matrix *vk, // a n-by-1 SPEX_DYNAMIC_CSC matrix that // contains the column to be inserted. int64_t k, // The column index that vk will be inserted const SPEX_options *option// Command parameters ); \end{verbatim} } \end{mdframed} \verb|SPEX_Update_LU_ColRep| exactly updates the REF LU factorization of a $n$-by-$n$ sparse matrix $A$ when one column of $A$ is replaced. On input \verb|F| contains the valid LU factorization to be updated that consists of \verb|L, U, rhos, P_perm, Pinv_perm| and \verb|Q_perm| (and \verb|Qinv_perm| for an updatable \verb|F|); \verb|vk| is a $n$-by-$1$ \verb|SPEX_DYNAMIC_CSC| matrix that contains the column to be inserted; \verb|k| is the index of the column of $A$ to be replaced and thus must be in the range of $[0,n-1]$; and option contains command parameters (default settings are used if \verb|option| is input as \verb|NULL|). This function requires that the rows of \verb|vk| are in the same order as that of $A$ (and \verb|F->L| if \verb|F| is updatable). In addition, \verb|vk| and \verb|SPEX_matrix *A| that contains $A$ are scaled with same factor, i.e., \verb|vk->scale = A->scale|, while \verb|vk->v[0]->scale = 1|. (Users can refer to section \ref{ss:SPEX_matrix} for the difference between \verb|vk->scale| and \verb|vk->v[0]->scale|.) As for \verb|F|, if it is input in a non-updatable format, it is first converted to updatable format. On success \verb|SPEX_OK| is returned, and \verb|F| contains the updatable REF LU factorization of $\bar{A}$ which is obtained by replacing column $k$ of $A$ with \verb|vk|, while \verb|vk| is unmodified. Otherwise, if this function fails for any reason, the appropriate error code is returned and \verb|F| is undefined on output (since the modification for the factorization is done in place during the update process). This function does not require/use/modify the original matrix $A$. Therefore, if the updated matrix $\bar{A}$ is needed, users should call \verb|SPEX_Update_matrix_colrep| {\bf AFTER} using this function (since \verb|vk| will be modified by \verb|SPEX_Update_matrix_colrep|). %------------------------------------------------------------------------------- \section{SPEX Update Solve Routines} %------------------------------------------------------------------------------- When factorization is updated iteratively, it should not be converted back and forth between updatable and non-updatable formats. However, this happens when users use the default LU/Cholesky solves (i.e., \verb|SPEX_Left_LU_solve| or \verb|SPEX_Chol_solve|) after using either \verb|SPEX_Update_Chol_Rank1| or \verb|SPEX_Update_LU_ColRep|, since the default LU/Cholesky solves always convert factorization to non-updatable format while the SPEX Update functions require and output updatable factorization. To avoid this, SPEX Update provides functions to solve the linear equation $Ax=b$ based on the updatable LU/Cholesky factorization of $A$. In particular, both a normal solve and a transpose solve using the updatable factorizations are provided. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Update_solve|: Solve $Ax=b$ with the updatable factorization} \subsection{\texttt{SPEX\_Update\_solve}: Solve linear system with the updatable factorization}\label{ss:SPEX_Update_solve} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_Update_solve // solves Ax = b ( // Output SPEX_matrix **x_handle, // rational solution to the system // input/output SPEX_factorization *F, // The updatable LU/Cholesky factorization of A // input: const SPEX_matrix *b, // Right hand side vector const SPEX_options* option // Command options ); \end{verbatim} } \end{mdframed} \verb|SPEX_Update_solve| obtains the solution of \verb|mpq_t| type to the linear system $Ax=b$ upon a valid LU/Cholesky factorization of $A$. On input, \verb|SPEX_matrix *x| that \verb|x_handle| points to is undefined; \verb|F| must be a valid LU/Cholesky factorization; \verb|b| must be dense \verb|mpz_t| with same number of rows as \verb|F->L|; and default settings are used if \verb|option| is input as \verb|NULL|. Upon successful completion, the function returns \verb|SPEX_OK|, and the dense \verb|mpq_t| matrix \verb|x| contains the solution to the linear system $Ax=b$. \verb|F| is mathematically unchanged on output. However, if \verb|F| is non-updatable on input, it is converted to updatable. If \verb|F| is already updatable, it is not modified. In case of failure, \verb|x| is returned as \verb|NULL| and the appropriate error code is returned. %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Update_tsolve|: Solve $A^Tx=b$ with the updatable factorization} \subsection{\texttt{SPEX\_Update\_tsolve}: Solve transposed linear system with the updatable factorization} \label{ss:SPEX_Update_tsolve} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_Update_tsolve // solves A^T x = b ( // Output SPEX_matrix **x_handle, // rational solution to the system // input/output SPEX_factorization *F, // The updatable LU/Cholesky factorization of A // input: const SPEX_matrix *b, // Right hand side vector const SPEX_options* option // Command options ); \end{verbatim} } \end{mdframed} \verb|SPEX_Update_tsolve| obtains the solution of \verb|mpq_t| type to the linear system $A^Tx=b$ upon a valid LU/Cholesky factorization of $A$. On input, \verb|SPEX_matrix *x| that \verb|x_handle| points to is undefined; \verb|F| must be a valid LU/Cholesky factorization; \verb|b| must be dense \verb|mpz_t| with same number of rows as \verb|F->L|; and default settings are used if \verb|option| is input as \verb|NULL|. Upon successful completion, the function returns \verb|SPEX_OK|, and the dense \verb|mpq_t| matrix \verb|x| contains the solution to the linear system $A^Tx=b$. \verb|F| is mathematically unchanged on output. However, if \verb|F| is non-updatable on input, it is converted to updatable. If \verb|F| is already updatable, it is not modified. In case of failure, \verb|x| is returned as \verb|NULL| and the appropriate error code is returned. %------------------------------------------------------------------------------- \section {SPEX Update Helper Functions} %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- \subsection{\texttt{SPEX\_Update\_matrix\_colrep}: Swap the columns of a matrix} \label{ss:SPEX_Update_matrix_colrep} %------------------------------------------------------------------------------- \begin{mdframed}[userdefinedwidth=\textwidth] {\footnotesize \begin{verbatim} SPEX_info SPEX_Update_matrix_colrep// performs column replacement ( // input/output SPEX_matrix *A, // m-by-n target matrix of SPEX_DYNAMIC_CSC SPEX_matrix *vk, // m-by-1 SPEX_DYNAMIC_CSC matrix that contains // the column vector to replace the k-th column // of A // input int64_t k, // The column index that vk will be inserted const SPEX_options *option// Command parameters ); \end{verbatim} } \end{mdframed} \verb|SPEX_Update_matrix_colrep| simply serves as a helper function in the column-replacement scenario where users wish to obtain the matrix after one of its columns is replaced. It should be noted that this does {\bf NOT} update the factorization of \verb|A|, and thus users should call \verb|SPEX_Update_LU_ColRep| to update the REF LU factorization for column replacement. To be more specific, this function swaps column \verb|k| of the matrix \verb|A| with column of \verb|vk| (i.e., swaps \verb|A->v[k]| and \verb|vk->v[0]|). Therefore, both \verb|A| and \verb|vk| are modified by this function, which requires users to call this function {\bf BEFORE} using \verb|SPEX_Update_LU_ColRep|. On input, \verb|A| and \verb|vk| must be matrices of \verb|SPEX_DYNAMIC_CSC| kind with the same number of rows $m$; \verb|k| is the index of the column of $A$ to be replaced and thus must be in the range of $[0,m-1]$; and option contains command parameters (default settings are used if \verb|option| is input as \verb|NULL|). This function requires that the rows of \verb|vk| are in the same order as that of \verb|A|. In addition, \verb|vk| and \verb|A| are scaled with the same factor, i.e., \verb|vk->scale = A->scale|, while \verb|vk->v[0]->scale = 1|. (Users can refer to section \ref{ss:SPEX_matrix} for the difference between \verb|vk->scale| and \verb|vk->v[0]->scale|.) Different from \verb|SPEX_Update_LU_ColRep|, this function does not require \verb|A| to be a square matrix, as long as it has the same number of rows as \verb|vk|. On output, \verb|A->v[k]| and \verb|vk->v[0]| are simply swapped. } %------------------------------------------------------------------------------- \chapter{Using SPEX in MATLAB}\vspace{-0.75in} \label{s:Use:MATLAB} %------------------------------------------------------------------------------- The MATLAB interface of SPEX can be installed by navigating to the \verb|MATLAB| folder and typing \verb|spex_mex_install|. Doing so installs SPEX and allows the use of 3 mex functions \verb|spex_lu_backslash.m|, \verb|spex_cholesky_backslash.m|, and \verb|spex_backslash.m|. First, this section describes the \verb|option| struct in Section \ref{s:Use:MATLAB:setup}. The use of the factorization is discussed in Section \ref{s:Use:MATLAB:factor}. The \verb|SPEX/SPEX/MATLAB| folder must be in your MATLAB path. %------------------------------------------------------------------------------- \section{Optional parameter settings} \label{s:Use:MATLAB:setup} %------------------------------------------------------------------------------- The SPEX MATLAB interface includes an \verb|option| struct as in optional input parameter that modifies behavior. If this parameter is not provided, default parameter settings are used. The elements of the \verb'option' struct are listed below. Any fields not present in the struct are treated as their default values. \begin{itemize} \item \verb|option.pivot|: This parameter is a string that controls the pivoting scheme used. When selecting a pivot entry in a given column, the factorization method uses one of the following pivoting strategies. Note that importantly this is only valid for LU factorization: \begin{itemize} \item \verb|'smallest'|: (default) smallest pivot, \item \verb|'diagonal'|: diagonal pivot if possible, otherwise smallest pivot, \item \verb|'first'|: first nonzero pivot in each column, \item \verb|'tol smallest'|: diagonal pivot with a tolerance (\verb|option.tol|) for the smallest pivot, \item \verb|'tol largest'|: diagonal pivot with a tolerance (\verb|option.tol|) for the largest pivot, \item \verb|'largest'|: largest pivot. \end{itemize} \item \verb|option.order|: This parameter is a string controls the fill-reducing column preordering used. This is valid for either LU or Cholesky as Backslash will choose its own ordering. \begin{itemize} \item \verb|'none'|: no column ordering; factorize \verb'A' as-is. \item \verb|'colamd'|: COLAMD ordering (default for LU) \item \verb|'amd'|: AMD ordering (default for Cholesky) \end{itemize} \item \verb|option.tol|: This parameter determines the tolerance used if one of the threshold pivoting schemes is chosen. The default value is 1 and this parameter can take any value in the range $(0,1]$. This is only valid for LU factorization. \item \verb|option.solution|: a string determining how \verb|x| is to be returned: \begin{itemize} \item \verb|'double'|: \verb|x| is converted to a 64-bit floating-point approximate solution. This is the default. \item \verb|'vpa'|: \verb|x| is returned as a \verb|vpa| array with \verb|option.digits| digits (default is given by the MATLAB \verb|digits| function). The result may be inexact, if an entry in \verb|x| cannot be represented in the specified number of digits. To convert this \verb|x| to double, use \verb|x=double(x)|. \item \verb|'char'|: \verb|x| is returned as a cell array of strings, where \verb|x {i} =| \newline \verb|'numerator/denominator'| and both \verb|numerator| and \verb|denominator| are arbitrary-length strings of decimal digits. The result is always exact, although \verb|x| cannot be directly used in MATLAB for numerical calculations. It can be inspected or analyzed using MATLAB string manipulation. To convert \verb|x| to \verb|vpa|, use \verb|x=vpa(x)|. To convert \verb|x| to double, use \verb|x=double(vpa(x))|. \end{itemize} \item \verb|option.digits|: the number of decimal digits to use for \verb|x|, if \verb|option.solution| is \verb|'vpa'|. Must be in range 2 to $2^{29}$. \item \verb|option.print|: display the inputs and outputs (0: nothing (default), 1: just errors, 2: terse, 3: all). \end{itemize} %------------------------------------------------------------------------------- \section{SPEX m files for use} \label{s:Use:MATLAB:factor} %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- %\cprotect\subsection{\verb|SPEX_Left_LU_backslash.m|} \subsection{\texttt{spex\_lu\_backslash.m}} %------------------------------------------------------------------------------- The \verb|spex_lu_backslash.m| function solves the linear system $A x = b$ where $A \in \mathtt{R}^{n \times n}$, $x \in \mathtt{R}^{n \times m}$ and $b \in \mathtt{R}^{n \times m}$. The final solution vector(s) obtained via this function are exact prior to their conversion to double precision. This function expects as input a sparse matrix $A$ and dense set of right hand side vectors $b$. Optionally, \verb|option| struct can be passed in. Currently, there are 2 ways to use this function outlined below: \begin{itemize} \item \verb|x = spex_lu_backslash(A,b)| returns the solution to $A x = b$ using default settings. The solution vectors are more accurate than the solution obtained via \verb|x = A \ b|. The solution \verb|x| is returned as a MATLAB double matrix. \item \verb|x = spex_lu_backslash(A,b,option)| returns the solution to $A x = b$ using non-default settings from the \verb|option| struct. \end{itemize} If the result \verb|x| is held as a MATLAB double matrix, in conventional floating-point representation (\verb|double|), it is guaranteed to be exact only if the exact solution can be held in \verb|double| without modification. The solution \verb|x| may also be returned as a MATLAB \verb|vpa| array, or as a cell array of strings; See Section \ref{s:Use:MATLAB:setup} for details. %------------------------------------------------------------------------------- \subsection{\texttt{spex\_cholesky\_backslash.m}} %------------------------------------------------------------------------------- The \verb|spex_cholesky_backslash.m| function solves the linear system $A x = b$ where $A \in \mathtt{R}^{n \times n}$, $x \in \mathtt{R}^{n \times m}$ and $b \in \mathtt{R}^{n \times m}$. The final solution vector(s) obtained via this function are exact prior to their conversion to double precision. Note that A must be SPD otherwise this function returns an error. This function expects as input a sparse matrix $A$ and dense set of right hand side vectors $b$. Optionally, \verb|option| struct can be passed in. Currently, there are 2 ways to use this function outlined below: \begin{itemize} \item \verb|x = spex_cholesky_backslash(A,b)| returns the solution to $A x = b$ using default settings. The solution vectors are more accurate than the solution obtained via \verb|x = A \ b|. The solution \verb|x| is returned as a MATLAB double matrix. \item \verb|x = spex_cholesky_backslash(A,b,option)| returns the solution to $A x = b$ using non-default settings from the \verb|option| struct. \end{itemize} If the result \verb|x| is held as a MATLAB double matrix, in conventional floating-point representation (\verb|double|), it is guaranteed to be exact only if the exact solution can be held in \verb|double| without modification. The solution \verb|x| may also be returned as a MATLAB \verb|vpa| array, or as a cell array of strings; See Section \ref{s:Use:MATLAB:setup} for details. %------------------------------------------------------------------------------- \subsection{\texttt{spex\_backslash.m}} %------------------------------------------------------------------------------- The \verb|spex_backslash.m| function solves the linear system $A x = b$ where $A \in \mathtt{R}^{n \times n}$, $x \in \mathtt{R}^{n \times m}$ and $b \in \mathtt{R}^{n \times m}$. The final solution vector(s) obtained via this function are exact prior to their conversion to double precision. This function expects as input a sparse matrix $A$ and dense set of right hand side vectors $b$. Optionally, \verb|option| struct can be passed in. If A is numerically symmetric, it attempts a Cholesky factorization. If the Cholesky fails or if the matrix is not numerically symmetric it performs an LU factorization. Currently, there are 2 ways to use this function outlined below: \begin{itemize} \item \verb|x = spex_backslash(A,b)| returns the solution to $A x = b$ using default settings. The solution vectors are more accurate than the solution obtained via \verb|x = A \ b|. The solution \verb|x| is returned as a MATLAB double matrix. \item \verb|x = spex_backslash(A,b,option)| returns the solution to $A x = b$ using non-default settings from the \verb|option| struct. \end{itemize} If the result \verb|x| is held as a MATLAB double matrix, in conventional floating-point representation (\verb|double|), it is guaranteed to be exact only if the exact solution can be held in \verb|double| without modification. The solution \verb|x| may also be returned as a MATLAB \verb|vpa| array, or as a cell array of strings; See Section \ref{s:Use:MATLAB:setup} for details. %------------------------------------------------------------------------------- \subsection{\texttt{spex\_mex\_demo.m}} %------------------------------------------------------------------------------- This function provides a demo of the SPEX library. It shows the usage for an exact solution as well as error checking and tuning the parameters. The typical output of this function may be seen in the provided \verb|MATLAB/html| folder. %------------------------------------------------------------------------------- \chapter{Using SPEX in Python}\vspace{-0.75in} %------------------------------------------------------------------------------- The Python interface of SPEX can be installed by navigating to the Python folder and typing make. Doing so allows the use of the Python SPEX library. First, this section describes the \verb|Option| object in Section \ref{s:Python:option}. The use of SPEX to solve $Ax=b$ is discussed in Section \ref{s:Python:Funcs}. %------------------------------------------------------------------------------- \section{Optional parameter settings}\label{s:Python:option} %------------------------------------------------------------------------------- The SPEX Python interface includes an object as an optional input parameter that modifies behaviour. If this is not provided, default parameter settings are used. \begin{itemize} \item output: This parameter is a string that determines how the solution is to be returned \begin{itemize} \item \verb|'double'|: \verb|x| is converted to a 64-bit floating-point approximate solution. This is the default. \item \verb|'string'|: \verb|x| is returned as an array of strings. \end{itemize} \item ordering: This parameter is a string that controls the fill-reducing column preordering used. By default it is initialized as None, if this option is chosen, the solve functions use the appropriate default ordering (AMD for Cholesky and COLAMD for Left LU). \begin{itemize} \item \verb|'none'|: no column ordering; factorize $A$ as-is. \item \verb|'colamd'|: COLAMD ordering \item \verb|'amd'|: AMD ordering \end{itemize} \end{itemize} %------------------------------------------------------------------------------- \section{Functions in Python SPEX}\label{s:Python:Funcs} %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- \subsection{\texttt{lu\_backslash}} %------------------------------------------------------------------------------- The \verb|lu_backslash| function solves the linear system $Ax=b$ where$A \in \mathtt{R}^{n \times n}$, $x \in \mathtt{R}^{n \times 1}$ and $b \in \mathtt{R}^{n \times 1}$. The final solution vector(s) obtained via this function are exact prior to their conversion to double precision. The LU function expects as input a \verb|scipy| sparse matrix $A$ and a right hand side vector $b$. Optionally, \verb|option| object can be passed in. Currently, there are 2 ways to use this function outlined below: \begin{itemize} \item \verb|x=SPEX.lu_backslash(A,b)| returns the solution to $A x = b$ using default settings. The solution \verb|x| is returned as a \verb|numpy| double array. \item \verb|x=SPEX.lu_backslash(A,b,options)| returns the solution to $A x = b$ using non-default settings from the \verb|option| object. \end{itemize} If the result \verb|x| is held as a \verb|numpu| double array, in conventional floating-point representation (\verb|double|), it is guaranteed to be exact only if the exact solution can be held in \verb|double| without modification. The solution \verb|x| may also be returned as a list of strings; See Section \ref{s:Python:option} for details. %------------------------------------------------------------------------------- \subsection{\texttt{cholesky\_backslash}} %------------------------------------------------------------------------------- The \verb|cholesky_backslash| function solves the linear system $Ax=b$ where$A \in \mathtt{R}^{n \times n}$, $x \in \mathtt{R}^{n \times 1}$ and $b \in \mathtt{R}^{n \times 1}$. The final solution vector(s) obtained via this function are exact prior to their conversion to double precision. Note that $A$ must be symmetric positive definite. The Cholesky function expects as input a \verb|scipy| sparse matrix $A$ and a right hand side vector $b$. Optionally, \verb|option| object can be passed in. Currently, there are 2 ways to use this function outlined below: \begin{itemize} \item \verb|x=SPEX.cholesky_backslash(A,b)| returns the solution to $A x = b$ using default settings. The solution \verb|x| is returned as a \verb|numpy| double array. \item \verb|x=SPEX.cholesky_backslash(A,b,options)| returns the solution to $A x = b$ using non-default settings from the \verb|option| object. \end{itemize} If the result \verb|x| is held as a \verb|numpu| double array, in conventional floating-point representation (\verb|double|), it is guaranteed to be exact only if the exact solution can be held in \verb|double| without modification. The solution \verb|x| may also be returned as a list of strings; See Section \ref{s:Python:option} for details. %------------------------------------------------------------------------------- \subsection{\texttt{backslash}} %------------------------------------------------------------------------------- The \verb|backslash| function solves the linear system $Ax=b$ where$A \in \mathtt{R}^{n \times n}$, $x \in \mathtt{R}^{n \times 1}$ and $b \in \mathtt{R}^{n \times 1}$. The final solution vector(s) obtained via this function are exact prior to their conversion to double precision. Note that $A$ must be symmetric positive definite. The Backslash function expects as input a \verb|scipy| sparse matrix $A$ and a right hand side vector $b$. Optionally, \verb|option| object can be passed in. If A is numerically symmetric, it attempts a Cholesky factorization. If the Cholesky fails or if the matrix is not numerically symmetric it performs an LU factorization. Currently, there are 2 ways to use this function outlined below: \begin{itemize} \item \verb|x=SPEX.backslash(A,b)| returns the solution to $A x = b$ using default settings. The solution \verb|x| is returned as a \verb|numpy| double array. \item \verb|x=SPEX.backslash(A,b,options)| returns the solution to $A x = b$ using non-default settings from the \verb|option| object. \end{itemize} If the result \verb|x| is held as a \verb|numpu| double array, in conventional floating-point representation (\verb|double|), it is guaranteed to be exact only if the exact solution can be held in \verb|double| without modification. The solution \verb|x| may also be returned as a list of strings; See Section \ref{s:Python:option} for details. %------------------------------------------------------------------------------- \section{Demo} %------------------------------------------------------------------------------- There is a file that provides a demo of the SPEX library in Python \verb|demo.py|. It shows the usage for an exact solution as well as tuning the parameters. %------------------------------------------------------------------------------- % References %------------------------------------------------------------------------------- \bibliographystyle{siam} \bibliography{SPEX_UserGuide.bib} \end{document}