The R Project SVN R-packages

Rev

Blame | Last modification | View Log | Download | RSS feed 

\documentclass[11pt]{article}

\newcommand{\m}[1]{{\bf{#1}}}       % for matrices and vectors
\newcommand{\tr}{^{\sf T}}          % transpose

\topmargin 0in
\textheight 9in
\oddsidemargin 0pt
\evensidemargin 0pt
\textwidth 6.5in

%------------------------------------------------------------------------------
\begin{document}
%------------------------------------------------------------------------------

\title{AMD Version 1.2 User Guide}
\author{Patrick R. Amestoy\thanks{ENSEEIHT-IRIT,
2 rue Camichel 31017 Toulouse, France.
email: amestoy@enseeiht.fr.  http://www.enseeiht.fr/$\sim$amestoy.}
\and Timothy A. Davis\thanks{
Dept.~of Computer and Information Science and Engineering,
Univ.~of Florida, Gainesville, FL, USA.
email: davis@cise.ufl.edu.
http://www.cise.ufl.edu/$\sim$davis.
This work was supported by the National
Science Foundation, under grants ASC-9111263, DMS-9223088, and CCR-0203270.
Portions of the work were done while on sabbatical at Stanford University
and Lawrence Berkeley National Laboratory (with funding from Stanford
University and the SciDAC program).
}
\and Iain S. Duff\thanks{Rutherford Appleton Laboratory, Chilton, Didcot, 
Oxon OX11 0QX, England. email: i.s.duff@rl.ac.uk.  
http://www.numerical.rl.ac.uk/people/isd/isd.html.
This work was supported by the EPSRC under grant GR/R46441.
}}

\date{Aug. 30, 2005}
\maketitle

%------------------------------------------------------------------------------
\begin{abstract}
AMD is a set of routines that implements the approximate minimum degree ordering
algorithm to permute sparse matrices prior to
numerical factorization.
There are versions written in both C and Fortran 77.
A MATLAB interface is included.
\end{abstract}
%------------------------------------------------------------------------------

Technical report TR-04-002 (revised), CISE Department, University of Florida,
Gainesville, FL, 2005.

AMD Version 1.2, Copyright\copyright 2005 by Timothy A.
Davis, Patrick R. Amestoy, and Iain S. Duff.  All Rights Reserved.
AMD is available under alternate licences; contact T. Davis for details.

{\bf AMD License:}
    Your use or distribution of AMD or any modified version of
    AMD implies that you agree to this License.

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
    USA

    Permission is hereby granted to use or copy this program under the
    terms of the GNU LGPL, provided that the Copyright, this License,
    and the Availability of the original version is retained on all copies.
    User documentation of any code that uses this code or any modified
    version of this code must cite the Copyright, this License, the
    Availability note, and "Used by permission." Permission to modify
    the code and to distribute modified code is granted, provided the
    Copyright, this License, and the Availability note are retained,
    and a notice that the code was modified is included.

{\bf Availability:}
    http://www.cise.ufl.edu/research/sparse/amd

{\bf Acknowledgments:}

    This work was supported by the National Science Foundation, under
    grants ASC-9111263 and DMS-9223088 and CCR-0203270.
    The conversion to C, the addition of the elimination tree
    post-ordering, and the handling of dense rows and columns
    were done while Davis was on sabbatical at
    Stanford University and Lawrence Berkeley National Laboratory.

%------------------------------------------------------------------------------
\newpage
\section{Overview}
%------------------------------------------------------------------------------

AMD is a set of routines for preordering a sparse matrix prior to
numerical factorization.  It uses an approximate minimum degree ordering
algorithm \cite{AmestoyDavisDuff96,AmestoyDavisDuff04}
to find a permutation matrix $\m{P}$
so that the Cholesky factorization $\m{PAP}\tr=\m{LL}\tr$ has fewer
(often much fewer) nonzero entries than the Cholesky factorization of $\m{A}$.
The algorithm is typically much faster than other ordering methods
and  minimum degree ordering
algorithms that compute an exact degree \cite{GeorgeLiu89}.
Some methods, such as approximate deficiency
\cite{RothbergEisenstat98} and graph-partitioning based methods
\cite{Chaco,KarypisKumar98e,PellegriniRomanAmestoy00,schu:01}
can produce better orderings, depending on the matrix.

The algorithm starts with an undirected graph representation of a
symmetric sparse matrix $\m{A}$.  Node $i$ in the graph corresponds to row
and column $i$ of the matrix, and there is an edge $(i,j)$ in the graph if
$a_{ij}$ is nonzero.
The degree of a node is initialized to the number of off-diagonal nonzeros
in row $i$, which is the size of the set of nodes
adjacent to $i$ in the graph.

The selection of a pivot $a_{ii}$ from the diagonal of $\m{A}$ and the first
step of Gaussian elimination corresponds to one step of graph elimination.
Numerical fill-in causes new nonzero entries in the matrix
(fill-in refers to
nonzeros in $\m{L}$ that are not in $\m{A}$).
Node $i$ is eliminated and edges are added to its neighbors
so that they form a clique (or {\em element}).  To reduce fill-in,
node $i$ is selected as the node of least degree in the graph.
This process repeats until the graph is eliminated.

The clique is represented implicitly.  Rather than listing all the
new edges in the graph, a single list of nodes is kept which represents
the clique.  This list corresponds to the nonzero pattern of the first
column of $\m{L}$.  As the elimination proceeds, some of these cliques
become subsets of subsequent cliques, and are removed.   This graph
can be stored in place, that is
using the same amount of memory as the original graph.

The most costly part of the minimum degree algorithm is the recomputation
of the degrees of nodes adjacent to the current pivot element.
Rather than keep track of the exact degree, the approximate minimum degree
algorithm finds an upper bound on the degree that is easier to compute.
For nodes of least degree, this bound tends to be tight.  Using the
approximate degree instead of the exact degree leads to a substantial savings
in run time, particularly for very irregularly structured matrices.
It has no effect on the quality of the ordering.

In the C version of AMD, the elimination phase is followed by an
elimination tree post-ordering.  This has no effect on fill-in, but
reorganizes the ordering so that the subsequent numerical factorization is
more efficient.  It also includes a pre-processing phase in which nodes of
very high degree are removed (without causing fill-in), and placed last in the
permutation $\m{P}$.  This reduces the run time substantially if the matrix
has a few rows with many nonzero entries, and has little effect on the quality
of the ordering.
The C version operates on the
symmetric nonzero pattern of $\m{A}+\m{A}\tr$, so it can be given
an unsymmetric matrix, or either the lower or upper triangular part of
a symmetric matrix.

The two Fortran versions of AMD are essentially identical to two versions of
the AMD algorithm discussed in an earlier paper \cite{AmestoyDavisDuff96}
(approximate minimum external degree, both with and without aggressive
absorption).
For a discussion of the long history of the minimum degree algorithm,
see \cite{GeorgeLiu89}.

%------------------------------------------------------------------------------
\section{Availability}
%------------------------------------------------------------------------------

In addition to appearing as a Collected Algorithm of the ACM,
AMD Version 1.2 is available at http://www.cise.ufl.edu/research/sparse.
The Fortran version is available as the routine {\tt MC47} in HSL
(formerly the Harwell Subroutine Library) \cite{hsl:2002}.

%------------------------------------------------------------------------------
\section{Using AMD in MATLAB}
%------------------------------------------------------------------------------

To use AMD in MATLAB, you must first compile the AMD mexFunction.
Just type {\tt make} in the Unix system shell, while in the {\tt AMD}
directory.  You can also type {\tt amd\_make} in MATLAB, while in the
{\tt AMD/MATLAB} directory.  Place the {\tt AMD/MATLAB} directory in your
MATLAB path.  This can be done on any system with MATLAB, including Windows.
See Section~\ref{Install} for more details on how to install AMD.

The MATLAB statement {\tt p=amd(A)} finds a permutation vector {\tt p} such
that the Cholesky factorization {\tt chol(A(p,p))} is typically sparser than
{\tt chol(A)}.
If {\tt A} is unsymmetric, {\tt amd(A)} is identical to {\tt amd(A+A')}
(ignoring numerical cancellation).
If {\tt A} is not symmetric positive definite,
but has substantial diagonal entries and a mostly symmetric nonzero pattern,
then this ordering is also suitable for LU factorization.  A partial pivoting
threshold may be required to prevent pivots from being selected off the
diagonal, such as the statement {\tt [L,U,P] = lu (A (p,p), 0.1)}.
Type {\tt help lu} for more details.
The statement {\tt [L,U,P,Q] = lu (A (p,p))} in MATLAB 6.5 is
not suitable, however, because it uses UMFPACK Version 4.0 and thus
does not attempt to select pivots from the diagonal.  UMFPACK Version 4.1
uses several strategies, including a symmetric pivoting strategy, and
will give you better results if you want to factorize an unsymmetric matrix
of this type.  Refer to the UMFPACK User Guide for more details, at
http://www.cise.ufl.edu/research/sparse/umfpack.

The AMD mexFunction is much faster than the built-in MATLAB symmetric minimum
degree ordering methods, SYMAMD and SYMMMD.  Its ordering quality is
comparable to SYMAMD, and better than SYMMMD
\cite{DavisGilbertLarimoreNg04}.

An optional input argument can be used to modify the control parameters for
AMD (aggressive absorption, dense row/column handling, and printing of
statistics).  An optional output
argument provides statistics on the ordering, including an analysis of the
fill-in and the floating-point operation count for a subsequent factorization.
For more details (once AMD is installed),
type {\tt help amd} in the MATLAB command window.

%------------------------------------------------------------------------------
\section{Using AMD in a C program}
\label{Cversion}
%------------------------------------------------------------------------------

The C-callable AMD library consists of seven user-callable routines and one
include file.  There are two versions of each of the routines, with
{\tt int} and {\tt long} integers.
The routines with prefix
{\tt amd\_l\_} use {\tt long} integer arguments; the others use
{\tt int} integer arguments.  If you compile AMD in the standard
ILP32 mode (32-bit {\tt int}'s, {\tt long}'s, and pointers) then the versions
are essentially identical.  You will be able to solve problems using up to 2GB
of memory.  If you compile AMD in the standard LP64 mode, the size of an
{\tt int} remains 32-bits, but the size of a {\tt long} and a pointer both get
promoted to 64-bits.

The following routines are fully described in Section~\ref{Primary}:

\begin{itemize}
\item {\tt amd\_order}
({\tt long} version: {\tt amd\_l\_order})
    {\footnotesize
    \begin{verbatim}
    #include "amd.h"
    int n, Ap [n+1], Ai [nz], P [n] ;
    double Control [AMD_CONTROL], Info [AMD_INFO] ;
    int result = amd_order (n, Ap, Ai, P, Control, Info) ;
    \end{verbatim}
    }
    Computes the approximate minimum degree ordering of an $n$-by-$n$ matrix
    $\m{A}$.  Returns a permutation vector {\tt P} of size {\tt n}, where
    {\tt P[k] = i} if row and column {\tt i} are the {\tt k}th row and
    column in the permuted matrix.
    This routine allocates its own memory of size $1.2e+9n$ integers,
    where $e$ is the number of nonzeros in $\m{A}+\m{A}\tr$.
    It computes statistics about the matrix $\m{A}$, such as the symmetry of
    its nonzero pattern, the number of nonzeros in $\m{L}$,
    and the number of floating-point operations required for Cholesky and LU
    factorizations (which are returned in the {\tt Info} array).
    The user's input matrix is not modified.
    It returns {\tt AMD\_OK} if successful, {\tt AMD\_INVALID} if
    the matrix is invalid, or {\tt AMD\_OUT\_OF\_MEMORY} if out of memory.

\item {\tt amd\_defaults}
({\tt long} version: {\tt amd\_l\_defaults})
    {\footnotesize
    \begin{verbatim}
    #include "amd.h"
    double Control [AMD_CONTROL] ;
    amd_defaults (Control) ;
    \end{verbatim}
    }
    Sets the default control parameters in the {\tt Control} array.  These can
    then be modified as desired before passing the array to the other AMD
    routines.

\item {\tt amd\_control}
({\tt long} version: {\tt amd\_l\_control})
    {\footnotesize
    \begin{verbatim}
    #include "amd.h"
    double Control [AMD_CONTROL] ;
    amd_control (Control) ;
    \end{verbatim}
    }
    Prints a description of the control parameters, and their values.

\item {\tt amd\_info}
({\tt long} version: {\tt amd\_l\_info})
    {\footnotesize
    \begin{verbatim}
    #include "amd.h"
    double Info [AMD_INFO] ;
    amd_info (Info) ;
    \end{verbatim}
    }
    Prints a description of the statistics computed by AMD, and their values.

\item {\tt amd\_preprocess}
({\tt long} version: {\tt amd\_l\_preprocess})
    {\footnotesize
    \begin{verbatim}
    #include "amd.h"
    int n, Ap [n+1], Ai [nz], Rp [n+1], Ri [nz] ;
    int result = amd_preprocess (n, Ap, Ai, Rp, Ri) ;
    \end{verbatim}
    }
    Removes duplicate entries and sorts each column of its input $\m{A}$,
    and returns the nonzero pattern of the transpose, $\m{R}=\m{A}\tr$.
    It returns the transpose because this is the simplest way to sort
    a matrix and remove duplicate entries.  Either $\m{A}$ or $\m{A}\tr$
    can be passed to {\tt amd\_order} with little effect on the
    ordering (except for minor tie-breaking changes).

\item {\tt amd\_valid}
({\tt long} version: {\tt amd\_valid})
    {\footnotesize
    \begin{verbatim}
    #include "amd.h"
    int n, Ap [n+1], Ai [nz] ;
    int result = amd_valid (n, n, Ap, Ai) ;
    \end{verbatim}
    }
    Returns {\tt TRUE} (1) if the matrix is valid as input to {\tt amd\_order},
    {\tt FALSE} (0) otherwise.  For {\tt amd\_order}, the matrix must
    also be square.  The first two arguments are the number of rows and the
    number of columns of the matrix.  For its use in AMD, these must both
    equal {\tt n}.

\item {\tt amd\_2}
({\tt long} version: {\tt amd\_l2})
    AMD ordering kernel.  It is faster than {\tt amd\_order}, and
    can be called by the user, but it is difficult to use.
    It does not check its inputs for errors.
    It does not require the columns of its input matrix to be sorted,
    but it destroys the matrix on output.  Additional workspace must be passed.
    Refer to the source file {\tt AMD/Source/amd\_2.c} for a description.

\end{itemize}

The nonzero pattern of the matrix $\m{A}$ is represented in compressed column
form.
For an $n$-by-$n$ matrix $\m{A}$ with {\tt nz} nonzero entries, the
representation consists of two arrays: {\tt Ap} of size {\tt n+1} and {\tt Ai}
of size {\tt nz}.  The row indices of entries in column {\tt j} are stored in
    {\tt Ai[Ap[j]} $\ldots$ {\tt Ap[j+1]-1]}.
For {\tt amd\_order},
no duplicate row indices may be present, and the row indices in any given
column must be sorted in ascending order.
The matrix is 0-based, and thus
row indices must be in the range {\tt 0} to {\tt n-1}.
The first entry {\tt Ap[0]} must be zero.
The total number of entries in the matrix is thus {\tt nz = Ap[n]}.

The matrix must be square, but it does not need to be symmetric.
The {\tt amd\_order} routine constructs the nonzero pattern of
$\m{B} = \m{A}+\m{A}\tr$ (without forming $\m{A}\tr$ explicitly),
and then orders the matrix $\m{B}$.  Thus, either the
lower triangular part of $\m{A}$, the upper triangular part,
or any combination may be passed.  The transpose $\m{A}\tr$ may also be
passed to {\tt amd\_order}.
The diagonal entries may be present, but are ignored.

The input to {\tt amd\_order} must have sorted columns because it uses
an in-place algorithm to construct $\m{A}+\m{A}\tr$ without first constructing
$\m{A}\tr$.  This saves memory, but places an additional restriction on
the input.  If the input matrix has columns with unsorted and/or duplicate
row indices, it is not valid as input to {\tt amd\_order}.  To handle this
case, the {\tt amd\_preprocess} routine is provided.  It sorts, transposes,
and removes duplicate entries from its input matrix, returning its result
as another compressed-column matrix $\m{R}$ which can then be passed to
{\tt amd\_order}.

%------------------------------------------------------------------------------
\subsection{Control parameters}
\label{control_param}
%------------------------------------------------------------------------------

Control parameters are set an optional {\tt Control} array.
It is optional in the sense that if
a {\tt NULL} pointer is passed for the {\tt Control} input argument,
then default control parameters are used.
%
\begin{itemize}
\item {\tt Control[AMD\_DENSE]} (or {\tt Control(1)} in MATLAB):
controls the threshold for ``dense''
rows/columns.  A dense row/column in $\m{A}+\m{A}\tr$
can cause AMD to spend significant time
in ordering the matrix.  If {\tt Control[AMD\_DENSE]} $\ge 0$,
rows/columns with
more than {\tt Control[AMD\_DENSE]} $\sqrt{n}$ entries are ignored during
the ordering, and placed last in the output order.  The default
value of {\tt Control[AMD\_DENSE]} is 10.  If negative, no rows/columns
are treated as ``dense.''  Rows/columns with 16 or fewer off-diagonal
entries are never considered ``dense.''
%
\item {\tt Control[AMD\_AGGRESSIVE]} (or {\tt Control(2)} in MATLAB):
controls whether or not to use
aggressive absorption, in which a prior element is absorbed into the current
element if it is a subset of the current element, even if it is not
adjacent to the current pivot element (refer
to \cite{AmestoyDavisDuff96,AmestoyDavisDuff04}
for more details).  The default value is nonzero,
which means that aggressive absorption will be performed.  This nearly always
leads to a better ordering (because the approximate degrees are more
accurate) and a lower execution time.  There are cases where it can
lead to a slightly worse ordering, however.  To turn it off, set
{\tt Control[AMD\_AGGRESSIVE]} to 0.
%
\end{itemize}

Statistics are returned in the {\tt Info} array
(if {\tt Info} is {\tt NULL}, then no statistics are returned).
Refer to {\tt amd.h} file, for more details
(14 different statistics are returned, so the list is not included here).

New in version 1.2 is the ability to redefine the memory manager and
{\tt printf} routine.
Refer to Section~\ref{Primary} for more details.
Refer to the source file {\tt AMD/MATLAB/amd\_mex.c} for an example.

%------------------------------------------------------------------------------
\subsection{Sample C program}
%------------------------------------------------------------------------------

The following program, {\tt amd\_demo.c}, illustrates the basic use of AMD.
See Section~\ref{Synopsis} for a short description
of each calling sequence.

{\footnotesize
\begin{verbatim}
#include <stdio.h>
#include "amd.h"

int n = 5 ;
int Ap [ ] = { 0,   2,       6,       10,  12, 14} ;
int Ai [ ] = { 0,1, 0,1,2,4, 1,2,3,4, 2,3, 1,4   } ;
int P [5] ;

int main (void)
{
    int k ;
    (void) amd_order (n, Ap, Ai, P, (double *) NULL, (double *) NULL) ;
    for (k = 0 ; k < n ; k++) printf ("P [%d] = %d\n", k, P [k]) ;
    return (0) ;
}
\end{verbatim}
}

The {\tt Ap} and {\tt Ai} arrays represent the binary matrix
\[
\m{A} = \left[
\begin{array}{rrrrr}
 1 &  1 &  0 &  0 &  0 \\
 1 &  1 &  1 &  0 &  1 \\
 0 &  1 &  1 &  1 &  0 \\
 0 &  0 &  1 &  1 &  0 \\
 0 &  1 &  1 &  0 &  1 \\
\end{array}
\right].
\]
The diagonal entries are ignored.
%
AMD constructs the pattern of $\m{A}+\m{A}\tr$,
and returns a permutation vector of $(0, 3, 1, 4, 2)$.
%
Since the matrix is unsymmetric but with a mostly symmetric nonzero
pattern, this would be a suitable permutation for an LU factorization of a
matrix with this nonzero pattern and whose diagonal entries are not too small.
The program uses default control settings and does not return any statistics
about the ordering, factorization, or solution ({\tt Control} and {\tt Info}
are both {\tt (double *) NULL}).  It also ignores the status value returned by
{\tt amd\_order}.

More example programs are included with the AMD package.
The {\tt amd\_demo.c} program provides a more detailed demo of AMD.
Another example is the AMD mexFunction, {\tt amd\_mex.c}.

%------------------------------------------------------------------------------
\subsection{A note about zero-sized arrays}
%------------------------------------------------------------------------------

AMD uses several user-provided arrays of size {\tt n} or {\tt nz}.
Either {\tt n} or {\tt nz} can be zero.
If you attempt to {\tt malloc} an array of size zero,
however, {\tt malloc} will return a null pointer which AMD will report
as invalid.  If you {\tt malloc} an array of
size {\tt n} or {\tt nz} to pass to AMD, make sure that you handle the
{\tt n} = 0 and {\tt nz = 0} cases correctly.

%------------------------------------------------------------------------------
\section{Synopsis of C-callable routines}
\label{Synopsis}
%------------------------------------------------------------------------------

The matrix $\m{A}$ is {\tt n}-by-{\tt n} with {\tt nz} entries.

{\footnotesize
\begin{verbatim}
#include "amd.h"
int n, status, Ap [n+1], Ai [nz], P [n], Rp [n+1], Ri [nz] ;
double Control [AMD_CONTROL], Info [AMD_INFO] ;
amd_defaults (Control) ;
status = amd_order (n, Ap, Ai, P, Control, Info) ;
amd_control (Control) ;
amd_info (Info) ;
amd_preprocess (n, Ap, Ai, Rp, Ri) ;
status = amd_valid (n, n, Ap, Ai) ;
\end{verbatim}
}

The {\tt amd\_l\_*} routines are identical, except that all {\tt int}
arguments become {\tt long}:

{\footnotesize
\begin{verbatim}
#include "amd.h"
long n, status, Ap [n+1], Ai [nz], P [n], Rp [n+1], Ri [nz] ;
double Control [AMD_CONTROL], Info [AMD_INFO] ;
amd_l_defaults (Control) ;
status = amd_l_order (n, Ap, Ai, P, Control, Info) ;
amd_l_control (Control) ;
amd_l_info (Info) ;
amd_l_preprocess (n, Ap, Ai, Rp, Ri) ;
status = amd_l_valid (n, n, Ap, Ai) ;
\end{verbatim}
}

%------------------------------------------------------------------------------
\section{Using AMD in a Fortran program}
%------------------------------------------------------------------------------

Two Fortran versions of AMD are provided.  The {\tt AMD} routine computes the
approximate minimum degree ordering, using aggressive absorption.  The
{\tt AMDBAR} routine is identical, except that it does not perform aggressive
absorption.  The {\tt AMD} routine is essentially identical to the HSL
routine {\tt MC47B/BD}.
Note that earlier versions of the Fortran
{\tt AMD} and {\tt AMDBAR} routines included an {\tt IOVFLO} argument,
which is no longer present.

In contrast to the C version, the Fortran routines require a symmetric
nonzero pattern, with no diagonal entries present although the {\tt MC47A/AD}
wrapper in HSL allows duplicates, ignores out-of-range entries, and only
uses entries from the upper triangular part of the matrix.  Although we
have an experimental Fortran code for treating ``dense'' rows, the Fortran
codes in this release do not treat
``dense'' rows and columns of $\m{A}$ differently, and thus their run time
can be high if there are a few dense rows and columns in the matrix.
They do not perform a post-ordering of the elimination tree,
compute statistics on the ordering, or check the validity of their input
arguments. These facilities are provided by {\tt MC47A/AD} and other
subroutines from HSL.
Only one {\tt integer}
version of each Fortran routine is provided.  
Both Fortran routines overwrite the user's input
matrix, in contrast to the C version.  
%
The C version does not return the elimination or assembly tree.
The Fortran version returns an assembly tree;
refer to the User Guide for details.
The following is the syntax of the {\tt AMD} Fortran routine.
The {\tt AMDBAR} routine is identical except for the routine name.

{\footnotesize
\begin{verbatim}
        INTEGER N, IWLEN, PFREE, NCMPA, IW (IWLEN), PE (N), DEGREE (N), NV (N),
     $          NEXT (N), LAST (N), HEAD (N), ELEN (N), W (N), LEN (N)
        CALL AMD (N, PE, IW, LEN, IWLEN, PFREE, NV, NEXT,
     $          LAST, HEAD, ELEN, DEGREE, NCMPA, W)
        CALL AMDBAR (N, PE, IW, LEN, IWLEN, PFREE, NV, NEXT,
     $          LAST, HEAD, ELEN, DEGREE, NCMPA, W)
\end{verbatim}
}

The input matrix is provided to {\tt AMD} and {\tt AMDBAR}
in three arrays, {\tt PE}, of size {\tt N},
{\tt LEN}, of size {\tt N}, and {\tt IW}, of size {\tt IWLEN}.  The size of
{\tt IW} must be at least {\tt NZ+N}.  The recommended size is
{\tt 1.2*NZ + N}.
On input, the indices of nonzero entries in row {\tt I} are stored in {\tt IW}.
{\tt PE(I)} is the index in {\tt IW} of the start of row {\tt I}.
{\tt LEN(I)} is the number of entries in row {\tt I}.
The matrix is 1-based, with row and column indices in the range 1 to {\tt N}.
Row {\tt I} is contained in
{\tt IW (PE(I)} $\ldots \:$ {\tt PE(I) + LEN(I) - 1)}.
The diagonal entries must not be present.  The indices within each row must
not contain any duplicates, but they need not be sorted.  The rows
themselves need not be in any particular order, and there may be empty space
between the rows.  If {\tt LEN(I)} is zero, then there are no off-diagonal
entries in row {\tt I}, and {\tt PE(I)} is ignored.  The integer
{\tt PFREE} defines what part of {\tt IW} contains the user's input matrix,
which is held in {\tt IW(1}~$\ldots~\:${\tt PFREE-1)}.
The contents of {\tt IW} and {\tt LEN} are undefined on output,
and {\tt PE} is modified to contain information about the ordering.

As the algorithm proceeds, it modifies the {\tt IW} array, placing the
pattern of the partially eliminated matrix in
{\tt IW(PFREE} $\ldots \:${\tt IWLEN)}.
If this space is exhausted, the space is compressed.
The number of compressions performed on the {\tt IW} array is
returned in the scalar {\tt NCMPA}.  The value of {\tt PFREE} on output is the
length of {\tt IW} required for no compressions to be needed.

The output permutation is returned in the array {\tt LAST}, of size {\tt N}.
If {\tt I=LAST(K)}, then {\tt I} is the {\tt K}th row in the permuted
matrix.  The inverse permutation is returned in the array {\tt ELEN}, where
{\tt K=ELEN(I)} if {\tt I} is the {\tt K}th row in the permuted matrix.
On output, the {\tt PE} and {\tt NV} arrays hold the assembly tree,
a supernodal elimination tree that represents the relationship between
columns of the Cholesky factor $\m{L}$.
If {\tt NV(I)} $> 0$, then {\tt I} is a node in the assembly
tree, and the parent of {\tt I} is {\tt -PE(I)}.  If {\tt I} is a root of
the tree, then {\tt PE(I)} is zero.  The value of {\tt NV(I)} is the
number of entries in the corresponding column of $\m{L}$, including the
diagonal.
If {\tt NV(I)} is zero, then {\tt I} is a non-principal node that is
not in the assembly tree.  Node {\tt -PE(I)} is the parent of node {\tt I}
in a subtree, the root of which is a node in the assembly tree.  All nodes
in one subtree belong to the same supernode in the assembly tree.
The other size {\tt N} arrays
({\tt DEGREE}, {\tt HEAD}, {\tt NEXT}, and {\tt W}) are used as workspace,
and are not defined on input or output.

If you want to use a simpler user-interface and compute the elimination
tree post-ordering, you should be able to call the C routines {\tt amd\_order}
or {\tt amd\_l\_order} from a Fortran program.   Just be sure to take into
account the 0-based indexing in the {\tt P}, {\tt Ap}, and {\tt Ai} arguments
to {\tt amd\_order} and {\tt amd\_l\_order}.  A sample interface is provided
in the files {\tt AMD/Demo/amd\_f77cross.f} and
{\tt AMD/Demo/amd\_f77wrapper.c}.  To compile the {\tt amd\_f77cross} program,
type {\tt make cross} in the {\tt AMD/Demo} directory.  The
Fortran-to-C calling conventions are highly non-portable, so this example
is not guaranteed to work with your compiler C and Fortran compilers.
The output of {\tt amd\_f77cross} is in {\tt amd\_f77cross.out}.

%------------------------------------------------------------------------------
\section{Sample Fortran main program}
%------------------------------------------------------------------------------

The following program illustrates the basic usage of the Fortran version of AMD.
The {\tt AP} and {\tt AI} arrays represent the binary matrix
\[
\m{A} = \left[
\begin{array}{rrrrr}
 1 &  1 &  0 &  0 &  0 \\
 1 &  1 &  1 &  0 &  1 \\
 0 &  1 &  1 &  1 &  1 \\
 0 &  0 &  1 &  1 &  0 \\
 0 &  1 &  1 &  0 &  1 \\
\end{array}
\right]
\]
in a conventional 1-based column-oriented form,
except that the diagonal entries are not present.
The matrix has the same as nonzero pattern of $\m{A}+\m{A}\tr$ in the C
program, in Section~\ref{Cversion}.
The output permutation is $(4, 1, 3, 5, 2)$.
It differs from the permutation returned by the C routine {\tt amd\_order}
because a post-order of the elimination tree has not yet been performed.

{\footnotesize
\begin{verbatim}
        INTEGER N, NZ, J, K, P, IWLEN, PFREE, NCMPA
        PARAMETER (N = 5, NZ = 10, IWLEN = 17)
        INTEGER AP (N+1), AI (NZ), LAST (N), PE (N), LEN (N), ELEN (N),
     $      IW (IWLEN), DEGREE (N), NV (N), NEXT (N), HEAD (N), W (N)
        DATA AP / 1, 2,     5,     8,  9,  11/
        DATA AI / 2, 1,3,5, 2,4,5, 3,  2,3   /
C       load the matrix into the AMD workspace
        DO 10 J = 1,N
            PE (J) = AP (J)
            LEN (J) = AP (J+1) - AP (J)
10      CONTINUE
        DO 20 P = 1,NZ
            IW (P) = AI (P)
20      CONTINUE
        PFREE = NZ + 1
C       order the matrix (destroys the copy of A in IW, PE, and LEN)
        CALL AMD (N, PE, IW, LEN, IWLEN, PFREE, NV, NEXT, LAST, HEAD,
     $      ELEN, DEGREE, NCMPA, W)
        DO 60 K = 1, N
            PRINT 50, K, LAST (K)
50          FORMAT ('P (',I2,') = ', I2)
60      CONTINUE
        END
\end{verbatim}
}

The {\tt Demo} directory contains an example of how the C version
may be called from a Fortran program, but this is highly non-portable.
For this reason, it is placed in the {\tt Demo} directory, not in the
primary {\tt Source} directory.

%------------------------------------------------------------------------------
\section{Installation}
\label{Install}
%------------------------------------------------------------------------------

The following discussion assumes you have the {\tt make} program, either in
Unix, or in Windows with Cygwin.

System-dependent configurations are in the {\tt AMD/Make}
directory.  You can edit the {\tt Make.include}
file in that directory to customize the compilation.  The default
settings will work on most systems.
Sample configuration files are provided
for Linux, Sun Solaris, SGI IRIX, IBM AIX, and the DEC/Compaq Alpha.

To compile and install the C-callable AMD library,
go to the {\tt AMD} directory and type {\tt make}.
The library will be placed in {\tt AMD/Lib/libamd.a}.
Two demo programs of the AMD ordering routine will be compiled and tested in
the {\tt AMD/Demo} directory.
The outputs of these demo programs will then be compared with output
files in the distribution.  The AMD mexFunction for
use in MATLAB will also be compiled.  If you do not have MATLAB
type {\tt make lib} instead.

To compile and install the Fortran-callable AMD library,
go to the {\tt AMD} directory and type {\tt make fortran}.
The library will be placed in {\tt AMD/Lib/libamdf77.a}.
A demo program will be compiled and tested in the {\tt AMD/Demo} directory.
The output will be compared with an output file in the distribution.

Typing {\tt make clean} will remove all but the final compiled libraries
and demo programs.  Typing {\tt make purge} removes all files not in the
original distribution.
If you compile AMD and then later change the {\tt Make.include}
file or your system-specific configuration file such as {\tt Make.linux},
then you should type {\tt make purge} and then {\tt make} to recompile.

Here are the various parameters that you can control in your
{\tt Make.include} file:

\begin{itemize}
\item {\tt CC = } your C compiler, such as {\tt cc}.
\item {\tt RANLIB = } your system's {\tt ranlib} program, if needed.
\item {\tt CFLAGS = } optimization flags, such as {\tt -O}.
\item {\tt LIB = } your libraries, such as {\tt -lm} or {\tt -lblas}.
\item {\tt RM =} the command to delete a file.
\item {\tt MV =} the command to rename a file.
\item {\tt MEX =} the command to compile a MATLAB mexFunction.
\item {\tt F77 =} the command to compile a Fortran program (optional).
\item {\tt F77FLAGS =} the Fortran compiler flags (optional).
\item {\tt F77LIB =} the Fortran libraries (optional).
\end{itemize}

The {\tt Make.include} includes some definitions regarding the BLAS.
This is so that AMD and UMFPACK (which requires AMD) can share
the same configuration files.  If you wish to use AMD only, then
you can ignore any references to the BLAS (the -DNBLAS compile flag).

When you compile your program that uses the C-callable AMD library,
you need to add the {\tt AMD/Lib/libamd.a} library
and you need to tell your compiler to look in the
{\tt AMD/Include} directory for include
files.   To compile a Fortran program that calls the Fortran AMD library,
you need to add the {\tt AMD/Lib/libamdf77.a} library.
See {\tt AMD/Demo/Makefile} for an example.

If all you want to use is the AMD mexFunction in MATLAB, you can skip
the use of the {\tt make} command entirely.  Simply type
{\tt amd\_make} in MATLAB while in the {\tt AMD/MATLAB} directory.
This works on any system with MATLAB, including Windows.

If you are including AMD as a subset of a larger library and do not want
to link the C standard I/O library, or if you simply do not need to use
them, you can safely remove the {\tt amd\_control.c} and {\tt amd\_info.c}
files.  Similarly, if you use default parameters (or define your
own {\tt Control} array), then you can exclude the {\tt amd\_defaults.c}
file.  The {\tt amd\_preprocess.c} file is optional as well, if you
can ensure that the input matrix to {\tt amd\_order} is always sorted
and has no duplicate entries.
Each of these files contains the user-callable routines of the same
name.  None of these auxiliary routines are directly called by
{\tt amd\_order}.
The {\tt amd\_dump.c} file contains debugging routines
that are neither used nor compiled unless debugging is enabled.
The {\tt amd\_internal.h} file must be edited to enable debugging;
refer to the instructions in that file.  Thus, it too can be excluded
if compiled into a larger production program or library.
The bare minimum files required to use just {\tt amd\_order} are
{\tt amd.h} in the {\tt Include} directory,
and
{\tt amd\_1.c},
{\tt amd\_2.c},
{\tt amd\_aat.c},
{\tt and\_order.c},
{\tt amd\_postorder.c},
{\tt amd\_post\_tree.c},
{\tt amd\_valid.c},
and
{\tt amd\_internal.h},
in the {\tt Source} directory.

%------------------------------------------------------------------------------
\newpage
\section{The AMD routines}
\label{Primary}
%------------------------------------------------------------------------------

The file {\tt AMD/Include/amd.h} listed below
describes each user-callable routine in the C version of AMD,
and gives details on their use.

{\footnotesize
\begin{verbatim}

/* ========================================================================= */
/* === AMD:  approximate minimum degree ordering =========================== */
/* ========================================================================= */

/* ------------------------------------------------------------------------- */
/* AMD Version 1.2, Copyright (c) 2005 by Timothy A. Davis,                  */
/* Patrick R. Amestoy, and Iain S. Duff.  See ../README.txt for License.     */
/* email: davis at cise.ufl.edu    CISE Department, Univ. of Florida.        */
/* web: http://www.cise.ufl.edu/research/sparse/amd                          */
/* ------------------------------------------------------------------------- */

/* AMD finds a symmetric ordering P of a matrix A so that the Cholesky
 * factorization of P*A*P' has fewer nonzeros and takes less work than the
 * Cholesky factorization of A.  If A is not symmetric, then it performs its
 * ordering on the matrix A+A'.  Two sets of user-callable routines are
 * provided, one for "int" integers and the other for "long" integers.
 *
 * The method is based on the approximate minimum degree algorithm, discussed
 * in Amestoy, Davis, and Duff, "An approximate degree ordering algorithm",
 * SIAM Journal of Matrix Analysis and Applications, vol. 17, no. 4, pp.
 * 886-905, 1996.  This package can perform both the AMD ordering (with
 * aggressive absorption), and the AMDBAR ordering (without aggressive
 * absorption) discussed in the above paper.  This package differs from the
 * Fortran codes discussed in the paper:
 *
 *      (1) it can ignore "dense" rows and columns, leading to faster run times
 *      (2) it computes the ordering of A+A' if A is not symmetric
 *      (3) it is followed by a depth-first post-ordering of the assembly tree
 *          (or supernodal elimination tree)
 *
 * For historical reasons, the Fortran versions, amd.f and amdbar.f, have
 * been left (nearly) unchanged.  They compute the identical ordering as
 * described in the above paper.
 */

#ifndef AMD_H
#define AMD_H

/* make it easy for C++ programs to include AMD */
#ifdef __cplusplus
extern "C" {
#endif

/* get the definition of size_t: */
#include <stddef.h>

int amd_order (             /* returns 0 if OK, negative value if error */
    int n,                  /* A is n-by-n.  n must be >= 0. */
    const int Ap [ ],       /* column pointers for A, of size n+1 */
    const int Ai [ ],       /* row indices of A, of size nz = Ap [n] */
    int P [ ],              /* output permutation, of size n */
    double Control [ ],     /* input Control settings, of size AMD_CONTROL */
    double Info [ ]         /* output Info statistics, of size AMD_INFO */
) ;

long amd_l_order (          /* see above for description of arguments */
    long n,
    const long Ap [ ],
    const long Ai [ ],
    long P [ ],
    double Control [ ],
    double Info [ ]
) ;

/* Input arguments (not modified):
 *
 *      n: the matrix A is n-by-n.
 *      Ap: an int/long array of size n+1, containing the column pointers of A.
 *      Ai: an int/long array of size nz, containing the row indices of A,
 *          where nz = Ap [n].
 *      Control:  a double array of size AMD_CONTROL, containing control
 *          parameters.  Defaults are used if Control is NULL.
 *
 * Output arguments (not defined on input):
 *
 *      P: an int/long array of size n, containing the output permutation. If
 *          row i is the kth pivot row, then P [k] = i.  In MATLAB notation,
 *          the reordered matrix is A (P,P).
 *      Info: a double array of size AMD_INFO, containing statistical
 *          information.  Ignored if Info is NULL.
 *
 * On input, the matrix A is stored in column-oriented form.  The row indices
 * of nonzero entries in column j are stored in Ai [Ap [j] ... Ap [j+1]-1].
 * The row indices must appear in ascending order in each column, and there
 * must not be any duplicate entries.  Row indices must be in the range 0 to
 * n-1.  Ap [0] must be zero, and thus nz = Ap [n] is the number of nonzeros
 * in A.  The array Ap is of size n+1, and the array Ai is of size nz = Ap [n].
 * The matrix does not need to be symmetric, and the diagonal does not need to
 * be present (if diagonal entries are present, they are ignored except for
 * the output statistic Info [AMD_NZDIAG]).  The arrays Ai and Ap are not
 * modified.  This form of the Ap and Ai arrays to represent the nonzero
 * pattern of the matrix A is the same as that used internally by MATLAB.
 * If you wish to use a more flexible input structure, please see the
 * umfpack_*_triplet_to_col routines in the UMFPACK package, at
 * http://www.cise.ufl.edu/research/sparse/umfpack, or use the amd_preprocess
 * routine discussed below.
 *
 * Restrictions:  n >= 0.  Ap [0] = 0.  Ap [j] <= Ap [j+1] for all j in the
 *      range 0 to n-1.  nz = Ap [n] >= 0.  For all j in the range 0 to n-1,
 *      and for all p in the range Ap [j] to Ap [j+1]-2, Ai [p] < Ai [p+1] must
 *      hold.  Ai [0..nz-1] must be in the range 0 to n-1.  To avoid integer
 *      overflow, (2.4*nz + 8*n) < INT_MAX / sizeof (int) for must hold for the
 *      "int" version. (2.4*nz + 8*n) < LONG_MAX / sizeof (long) must hold
 *      for the "long" version.  Finally, Ai, Ap, and P must not be NULL.  If
 *      any of these restrictions are not met, AMD returns AMD_INVALID.
 *
 * AMD returns:
 *
 *      AMD_OK if the matrix is valid and sufficient memory can be allocated to
 *          perform the ordering.
 *
 *      AMD_OUT_OF_MEMORY if not enough memory can be allocated.
 *
 *      AMD_INVALID if the input arguments n, Ap, Ai are invalid, or if P is
 *          NULL.
 *
 * The AMD routine first forms the pattern of the matrix A+A', and then
 * computes a fill-reducing ordering, P.  If P [k] = i, then row/column i of
 * the original is the kth pivotal row.  In MATLAB notation, the permuted
 * matrix is A (P,P), except that 0-based indexing is used instead of the
 * 1-based indexing in MATLAB.
 *
 * The Control array is used to set various parameters for AMD.  If a NULL
 * pointer is passed, default values are used.  The Control array is not
 * modified.
 *
 *      Control [AMD_DENSE]:  controls the threshold for "dense" rows/columns.
 *          A dense row/column in A+A' can cause AMD to spend a lot of time in
 *          ordering the matrix.  If Control [AMD_DENSE] >= 0, rows/columns
 *          with more than Control [AMD_DENSE] * sqrt (n) entries are ignored
 *          during the ordering, and placed last in the output order.  The
 *          default value of Control [AMD_DENSE] is 10.  If negative, no
 *          rows/columns are treated as "dense".  Rows/columns with 16 or
 *          fewer off-diagonal entries are never considered "dense".
 *
 *      Control [AMD_AGGRESSIVE]: controls whether or not to use aggressive
 *          absorption, in which a prior element is absorbed into the current
 *          element if is a subset of the current element, even if it is not
 *          adjacent to the current pivot element (refer to Amestoy, Davis,
 *          & Duff, 1996, for more details).  The default value is nonzero,
 *          which means to perform aggressive absorption.  This nearly always
 *          leads to a better ordering (because the approximate degrees are
 *          more accurate) and a lower execution time.  There are cases where
 *          it can lead to a slightly worse ordering, however.  To turn it off,
 *          set Control [AMD_AGGRESSIVE] to 0.
 *
 *      Control [2..4] are not used in the current version, but may be used in
 *          future versions.
 *
 * The Info array provides statistics about the ordering on output.  If it is
 * not present, the statistics are not returned.  This is not an error
 * condition.
 * 
 *      Info [AMD_STATUS]:  the return value of AMD, either AMD_OK,
 *          AMD_OUT_OF_MEMORY, or AMD_INVALID.
 *
 *      Info [AMD_N]: n, the size of the input matrix
 *
 *      Info [AMD_NZ]: the number of nonzeros in A, nz = Ap [n]
 *
 *      Info [AMD_SYMMETRY]:  the symmetry of the matrix A.  It is the number
 *          of "matched" off-diagonal entries divided by the total number of
 *          off-diagonal entries.  An entry A(i,j) is matched if A(j,i) is also
 *          an entry, for any pair (i,j) for which i != j.  In MATLAB notation,
 *              S = spones (A) ;
 *              B = tril (S, -1) + triu (S, 1) ;
 *              symmetry = nnz (B & B') / nnz (B) ;
 *
 *      Info [AMD_NZDIAG]: the number of entries on the diagonal of A.
 *
 *      Info [AMD_NZ_A_PLUS_AT]:  the number of nonzeros in A+A', excluding the
 *          diagonal.  If A is perfectly symmetric (Info [AMD_SYMMETRY] = 1)
 *          with a fully nonzero diagonal, then Info [AMD_NZ_A_PLUS_AT] = nz-n
 *          (the smallest possible value).  If A is perfectly unsymmetric
 *          (Info [AMD_SYMMETRY] = 0, for an upper triangular matrix, for
 *          example) with no diagonal, then Info [AMD_NZ_A_PLUS_AT] = 2*nz
 *          (the largest possible value).
 *
 *      Info [AMD_NDENSE]: the number of "dense" rows/columns of A+A' that were
 *          removed from A prior to ordering.  These are placed last in the
 *          output order P.
 *
 *      Info [AMD_MEMORY]: the amount of memory used by AMD, in bytes.  In the
 *          current version, this is 1.2 * Info  [AMD_NZ_A_PLUS_AT] + 9*n
 *          times the size of an integer.  This is at most 2.4nz + 9n.  This
 *          excludes the size of the input arguments Ai, Ap, and P, which have
 *          a total size of nz + 2*n + 1 integers.
 *
 *      Info [AMD_NCMPA]: the number of garbage collections performed.
 *
 *      Info [AMD_LNZ]: the number of nonzeros in L (excluding the diagonal).
 *          This is a slight upper bound because mass elimination is combined
 *          with the approximate degree update.  It is a rough upper bound if
 *          there are many "dense" rows/columns.  The rest of the statistics,
 *          below, are also slight or rough upper bounds, for the same reasons.
 *          The post-ordering of the assembly tree might also not exactly
 *          correspond to a true elimination tree postordering.
 *
 *      Info [AMD_NDIV]: the number of divide operations for a subsequent LDL'
 *          or LU factorization of the permuted matrix A (P,P).
 *
 *      Info [AMD_NMULTSUBS_LDL]:  the number of multiply-subtract pairs for a
 *          subsequent LDL' factorization of A (P,P).
 *
 *      Info [AMD_NMULTSUBS_LU]:  the number of multiply-subtract pairs for a
 *          subsequent LU factorization of A (P,P), assuming that no numerical
 *          pivoting is required.
 *
 *      Info [AMD_DMAX]:  the maximum number of nonzeros in any column of L,
 *          including the diagonal.
 *
 *      Info [14..19] are not used in the current version, but may be used in
 *          future versions.
 */    

/* ------------------------------------------------------------------------- */
/* AMD preprocess */
/* ------------------------------------------------------------------------- */

/* amd_preprocess: sorts, removes duplicate entries, and transposes the
 * nonzero pattern of a column-form matrix A, to obtain the matrix R.
 *
 * Alternatively, you can consider this routine as constructing a row-form
 * matrix from a column-form matrix.  Duplicate entries are allowed in A (and
 * removed in R). The columns of R are sorted.  Checks its input A for errors.
 *
 * On input, A can have unsorted columns, and can have duplicate entries.
 * Ap [0] must still be zero, and Ap must be monotonically nondecreasing.
 * Row indices must be in the range 0 to n-1.
 *
 * On output, if this routine returns AMD_OK, then the matrix R is a valid
 * input matrix for AMD_order.  It has sorted columns, with no duplicate
 * entries in each column.  Since AMD_order operates on the matrix A+A', it
 * can just as easily use A or A', so the transpose has no significant effect
 * (except for minor tie-breaking, which can lead to a minor effect in the
 * quality of the ordering).  As an example, compare the output of amd_demo.c
 * and amd_demo2.c.
 *
 * This routine transposes A to get R because that's the simplest way to
 * sort and remove duplicate entries from a matrix.
 *
 * Allocates 2*n integer work arrays, and free's them when done.
 *
 * If you wish to call amd_order, but do not know if your matrix has unsorted
 * columns or duplicate entries, then you can use the following code, which is
 * fairly efficient.  amd_order will not allocate any internal matrix until
 * it checks that the input matrix is valid, so the method below is memory-
 * efficient as well.  This code snippet assumes that Rp and Ri are already
 * allocated, and are the same size as Ap and Ai respectively.

    result = amd_order (n, p, Ap, Ai, Control, Info) ;
    if (result == AMD_INVALID)
    {
        if (amd_preprocess (n, Ap, Ai, Rp, Ri) == AMD_OK)
        {
            result = amd_order (n, p, Rp, Ri, Control, Info) ;
        }
    }

 * amd_preprocess will still return AMD_INVALID if any row index in Ai is out
 * of range or if the Ap array is invalid.  These errors are not corrected by
 * amd_preprocess since they represent a more serious error that should be
 * flagged with the AMD_INVALID error code.
 *
 * You may also call amd_valid directly (see below).
 */ 

int amd_preprocess
(
    int n,
    const int Ap [ ],
    const int Ai [ ],
    int Rp [ ],
    int Ri [ ]
) ;

long amd_l_preprocess
(
    long n,
    const long Ap [ ],
    const long Ai [ ],
    long Rp [ ],
    long Ri [ ]
) ;

/* Input arguments (not modified):
 *
 *      n: the matrix A is n-by-n.
 *      Ap: an int/long array of size n+1, containing the column pointers of A.
 *      Ai: an int/long array of size nz, containing the row indices of A,
 *          where nz = Ap [n].
 *      The nonzero pattern of column j of A is in Ai [Ap [j] ... Ap [j+1]-1].
 *      Ap [0] must be zero, and Ap [j] <= Ap [j+1] must hold for all j in the
 *      range 0 to n-1.  Row indices in Ai must be in the range 0 to n-1.
 *      The row indices in any one column need not be sorted, and duplicates
 *      may exist.
 *
 * Output arguments (not defined on input):
 *
 *      Rp: an int/long array of size n+1, containing the column pointers of R.
 *      Ri: an int/long array of size rnz, containing the row indices of R,
 *          where rnz = Rp [n].  Note that Rp [n] will be less than Ap [n] if
 *          duplicates appear in A.  In general, Rp [n] <= Ap [n].
 *      The data structure for R is the same as A, except that each column of
 *      R contains sorted row indices, and no duplicates appear in any column.
 *
 * amd_preprocess returns:
 *
 *      AMD_OK if the matrix A is valid and sufficient memory can be allocated
 *          to perform the preprocessing.
 *
 *      AMD_OUT_OF_MEMORY if not enough memory can be allocated.
 *
 *      AMD_INVALID if the input arguments n, Ap, Ai are invalid, or if Rp or
 *          Ri are NULL.
 */

/* ------------------------------------------------------------------------- */
/* direct interface to AMD */
/* ------------------------------------------------------------------------- */

/* This is the primary AMD ordering routine.  It is not meant to be
 * user-callable because of its restrictive inputs and because it destroys
 * the user's input matrix.  It does not check its inputs for errors, either.
 * However, if you can work with these restrictions it can be faster than
 * amd_order and use less memory (assuming that you can create your own copy
 * of the matrix for AMD to destroy).  Refer to AMD/Source/amd_2.c for a
 * description. */

void amd_2 (
    int n,
    int Pe [ ],
    int Iw [ ],
    int Len [ ],
    int iwlen,
    int pfree,
    int Nv [ ],
    int Next [ ], 
    int Last [ ],
    int Head [ ],
    int Elen [ ],
    int Degree [ ],
    int W [ ],
    double Control [ ],
    double Info [ ]
) ;

void amd_l2 (
    long n,
    long Pe [ ],
    long Iw [ ],
    long Len [ ],
    long iwlen,
    long pfree,
    long Nv [ ],
    long Next [ ], 
    long Last [ ],
    long Head [ ],
    long Elen [ ],
    long Degree [ ],
    long W [ ],
    double Control [ ],
    double Info [ ]
) ;

/* ------------------------------------------------------------------------- */
/* amd_valid */
/* ------------------------------------------------------------------------- */

/* Returns TRUE (1) if the matrix is valid as input to amd_order, FALSE (0)
 * otherwise.  For amd_order, the matrix must also be square. */

int amd_valid
(
    int n_row,              /* # of rows */
    int n_col,              /* # of columns */
    const int Ap [ ],       /* column pointers, of size n_col+1 */
    const int Ai [ ]        /* row indices, of size Ap [n_col] */
) ;

long amd_l_valid
(
    long n_row,
    long n_col,
    const long Ap [ ],
    const long Ai [ ]
) ;

/* ------------------------------------------------------------------------- */
/* AMD memory manager and printf routines */
/* ------------------------------------------------------------------------- */

/* The user can redefine these to change the malloc, free, and printf routines
 * that AMD uses. */

extern void *(*amd_malloc) (size_t) ;               /* pointer to malloc */
extern void (*amd_free) (void *) ;                  /* pointer to free */
extern void *(*amd_realloc) (void *, size_t) ;      /* pointer to realloc */
extern void *(*amd_calloc) (size_t, size_t) ;       /* pointer to calloc */
extern int (*amd_printf) (const char *, ...) ;      /* pointer to printf */

/* ------------------------------------------------------------------------- */
/* AMD Control and Info arrays */
/* ------------------------------------------------------------------------- */

/* amd_defaults:  sets the default control settings */
void amd_defaults   (double Control [ ]) ;
void amd_l_defaults (double Control [ ]) ;

/* amd_control: prints the control settings */
void amd_control    (double Control [ ]) ;
void amd_l_control  (double Control [ ]) ;

/* amd_info: prints the statistics */
void amd_info       (double Info [ ]) ;
void amd_l_info     (double Info [ ]) ;

#define AMD_CONTROL 5       /* size of Control array */
#define AMD_INFO 20         /* size of Info array */

/* contents of Control */
#define AMD_DENSE 0         /* "dense" if degree > Control [0] * sqrt (n) */
#define AMD_AGGRESSIVE 1    /* do aggressive absorption if Control [1] != 0 */

/* default Control settings */
#define AMD_DEFAULT_DENSE 10.0      /* default "dense" degree 10*sqrt(n) */
#define AMD_DEFAULT_AGGRESSIVE 1    /* do aggressive absorption by default */

/* contents of Info */
#define AMD_STATUS 0        /* return value of amd_order and amd_l_order */
#define AMD_N 1             /* A is n-by-n */
#define AMD_NZ 2            /* number of nonzeros in A */ 
#define AMD_SYMMETRY 3      /* symmetry of pattern (1 is sym., 0 is unsym.) */
#define AMD_NZDIAG 4        /* # of entries on diagonal */
#define AMD_NZ_A_PLUS_AT 5  /* nz in A+A' */
#define AMD_NDENSE 6        /* number of "dense" rows/columns in A */
#define AMD_MEMORY 7        /* amount of memory used by AMD */
#define AMD_NCMPA 8         /* number of garbage collections in AMD */
#define AMD_LNZ 9           /* approx. nz in L, excluding the diagonal */
#define AMD_NDIV 10         /* number of fl. point divides for LU and LDL' */
#define AMD_NMULTSUBS_LDL 11 /* number of fl. point (*,-) pairs for LDL' */
#define AMD_NMULTSUBS_LU 12  /* number of fl. point (*,-) pairs for LU */
#define AMD_DMAX 13          /* max nz. in any column of L, incl. diagonal */

/* ------------------------------------------------------------------------- */
/* return values of AMD */
/* ------------------------------------------------------------------------- */

#define AMD_OK 0                /* success */
#define AMD_OUT_OF_MEMORY -1    /* malloc failed */
#define AMD_INVALID -2          /* input arguments are not valid */

/* ========================================================================== */
/* === AMD version ========================================================== */
/* ========================================================================== */

/* AMD Version 1.2 and later will include the following definitions.
 * As an example, to test if the version you are using is 1.2 or later:
 *
 * #ifdef AMD_VERSION
 *      if (AMD_VERSION >= AMD_VERSION_CODE (1,2)) ...
 * #endif
 *
 * This also works during compile-time:
 *
 *      #if defined(AMD_VERSION) && (AMD_VERSION >= AMD_VERSION_CODE (1,2))
 *          printf ("This is version 1.2 or later\n") ;
 *      #else
 *          printf ("This is an early version\n") ;
 *      #endif
 *
 * Versions 1.1 and earlier of AMD do not include a #define'd version number.
 */

#define AMD_DATE "Aug. 30, 2005"
#define AMD_VERSION_CODE(main,sub) ((main) * 1000 + (sub))
#define AMD_MAIN_VERSION 1
#define AMD_SUB_VERSION 2
#define AMD_VERSION AMD_VERSION_CODE(AMD_MAIN_VERSION,AMD_SUB_VERSION)

#ifdef __cplusplus
}
#endif

#endif
\end{verbatim}
}

%------------------------------------------------------------------------------
\newpage
% References
%------------------------------------------------------------------------------

\bibliographystyle{plain}
\bibliography{AMD_UserGuide}

\end{document}