paper.tex

这是Yousef Saad编写的矩阵运算的Fortran软件包（A basic tool-kit for sparse matrix computations (Version 2),包含常见的排序

\subsubsection{The Variable Block Row format (VBR)}
In many applications, matrices are blocked, but the blocks are not all
the same size.  These so-called variable block matrices arise from the
discretization of systems of partial differential equations where there
is a varying number of equations at each grid point.  Like in the Block
Sparse Row (BSR) format, all entries of nonzero blocks (blocks which
contain any nonzeros) are stored, even if their value is zero.  Also
like the BSR format, there is significant savings in integer pointer
overhead in the data structure.

Variable block generalizations can be made to many matrix storage
formats.  The Variable Block Row (VBR) format is a generalization of the
Compressed Sparse Row (CSR) format, and is similar to the variable
block format used at the University of Waterloo, and one currently
proposed in the Sparse BLAS toolkit.

In the VBR format, the $IA$ and $JA$ arrays of the CSR format store the
sparsity structure of the blocks.  The entries in each block are stored
in $A$ in column-major order so that each block may be passed as a small
dense matrix to a Fortran subprogram.  The block row and block column
partitionings are stored in the vectors {\em KVSTR} and {\em KVSTC},
by storing the first row or column number of each block row or column
respectively.  In most applications, the block row and column partitionings
will be conformal, and the same array may be used in the programs.
Finally, integer pointers to the beginning of each block in $A$ are stored
in the array $KA$.

$IA$ contains pointers to the beginning of each block row in $JA$ and $KA$.
Thus $IA$ has length equal to the number of block rows (plus one to mark
the end of the matrix), and $JA$ has length equal to the number of nonzero
blocks.  $KA$ has the same length as $JA$ plus one to mark the end of
the matrix.  {\em KVSTR} and {\em KVSTC} have length equal to the number
of block rows and columns respectively, and $A$ has length equal to the
number of nonzeros in the matrix.  The following
figure shows the VBR format applied to a small matrix.

This version of Sparskit has a number of routines to support the variable
block matrix format.  CSRVBR and VBRCSR convert between the VBR and CSR
formats; VBRINFO prints some elementary information about the block
structure of a matrix in VBR format; AMUXV performs a matrix-vector product
with a matrix in VBR format; CSRKVSTR and CSRKVSTC are used
to determine row and column block partitionings of a matrix in CSR format,
and KVSTMERGE is used to combine row and column partitionings to achieve
a conformal partitioning.

\begin{figure}[htb]
%\vspace{4.8in}
%\centerline{\psfig{figure=vbrpic.eps,width=6.2in}}
\includegraphics[width=6.2in]{vbrpic}
%\special{psfile=vbrpic.eps vscale = 100 hscale = 100}
%\special{psfile=vbrpic.eps vscale = 100 hscale = 100 voffset= -420}
\caption {A $6 \times 8$ sparse matrix and its storage vectors.}
\end{figure}

\subsection{The FORMATS conversion module}
It is important to note that there is no need to have a subroutine for
each pair of data structures, since all we need is to be able to
convert any format to the standard row-compressed format and then back
to any other format.  There are currently 32 different conversion
routines in this module all of which are devoted to converting from
one data structure into another.

The naming mechanism adopted is to use a 6-character name for each of
the subroutines, the first 3 for the input format and the last 3 for
the output format.  Thus COOCSR performs the conversion from the
coordinate format to the Compressed Sparse Row format.  However it was
necessary to break the naming rule in one exception.  We needed a
version of COOCSR that is in-place, i.e., which can take the input
matrix, and convert it directly into a CSR format by using very little
additional work space.  This routine is called COICSR.  Each of
the formats has a routine to translate it to the CSR format and a routine
to convert back to it from the CSR format.  The only exception is that
a CSCCSR routine is not necessary since
the conversion from Column Sparse format to Sparse Row format  can be
performed with the same routine CSRCSC.  This is essentially a transposition
operation.

Considerable effort has been put at attempting to make the conversion
routines in-place, i.e., in allowing some or all of 
the output arrays to be the same as the input arrays. 
The purpose is to save storage whenever possible without
sacrificing performance. The added flexibility can be 
very convenient in some situations. 
When the additional coding complexity to  permit
the routine to be in-place was not too high this was always done.
If the subroutine is in-place this is clearly indicated in the
documentation. As mentioned above, we found it necessary in one instance
to provide both the in-place version as well as the regular version:
 COICSR is an in-place version of the COOCSR routine.
We would also like to add that other routines that avoid the
CSR format for some of the more important data structures may
eventually be included. For now, there is only one such routine
\footnote{Contributed by E. Rothman from Cornell University.}
namely, COOELL. 

\subsection{Internal format used in SPARSKIT}
Most of the routines in SPARSKIT  use internally the Compressed
Sparse Row format. The selection of 
the CSR mode has been motivated by several factors. Simplicity,
generality, and widespread use are certainly the most 
important ones. However, it has often been argued 
that the column scheme may
have been a better choice. One argument in this favor is that
vector machines usually give a better performance for such 
operations as matrix vector by multiplications for matrices
stored in CSC format. In fact for parallel machines which
have a low overhead in loop synchronization (e.g., the Alliants),
the situation is reversed, see \cite{Saad-Boeing} for details.
For almost any argument in favor of one scheme there seems
to be an argument in favor of the other. Fortunately,
the difference provided in functionality is rather minor.
For example the subroutine APLB to add two matrices in CSR format,
described in Section 5.1, can actually be also used to add 
two matrices in CSC format, since the data structures
are identical. Several such subroutines can be used for both
schemes, by pretending that the input matrices
are stored in CSR mode whereas in fact they are 
stored in CSC mode. 

\section{Manipulation routines} The module UNARY   
of SPARSKIT consists of a  number of  utilities to  manipulate and perform
basic operations with sparse matrices. The  following sections 
give an overview of this part of the package.

\subsection{Miscellaneous operations with sparse matrices}
There are a large number of non-algebraic operations that
are commonly used when working with sparse matrices. 
A typical example is to transform $A$ into $B = P A Q $
where $P$ and $Q$ are two permutation matrices. 
Another example is to extract the lower triangular
part of $A$ or a given diagonal from $A$. Several other 
such `extraction' operations are supplied in SPARSKIT. 
Also provided is the transposition function. This may seem
as an unnecessary addition since the routine 
CSRCSC already does perform this function economically. However,
the new transposition provided is in-place, in that it may
transpose the matrix and overwrite the result on the original matrix,
thus saving memory usage. Since many of these manipulation routines
involve one matrix (as opposed to two in the basic linear algebra routines)
we created a module called  UNARY to include these subroutines.

Another set of subroutines that are sometimes useful 
are those involving  a `mask'. A mask
defines a given nonzero pattern and for all practical
purposes a mask matrix is a sparse matrix whose nonzero
entries are all ones (therefore there is no need to store 
its real values).  Sometimes it is useful
to extract from a given matrix $A$ the `masked' matrix according
to a mask $M$, i.e., to compute the matrix
$A \odot M$ , where $\odot $ denotes the element-wise matrix
product, and $M$ is some mask matrix.

\subsection{The module UNARY}
This module of SPARSKIT consists of a number of routines
to  perform some basic non-algebraic operations on a matrix.
The following is a list of the routines currently supported with
a brief explanation. 
%%There are many other routines which are not 
%% listed their inclusion still being debated.

\vskip .5in

\marg{SUBMAT}\disp{Extracts a square or
rectangular submatrix from a sparse matrix. Both the 
input and output  matrices are in CSR format.
The routine is in-place.} 

\marg{FILTER}\disp{Filters out elements from a
 matrix according to their magnitude. Both the input and
the output matrices are in CSR format.
The output matrix, is obtained from the
input matrix by removing all the elements that are smaller than
a  certain threshold. The threshold is computed for each row
according to one of three provided options.
The algorithm is in-place.}

\marg{FILTERM}\disp{Same as above, but for the MSR format.}

\marg{CSORT}\disp{Sorts the elements of a matrix stored in
CSR format in increasing order of the column numbers. }

\marg{ TRANSP  }\disp{ This is an in-place transposition routine,
i.e., it can be viewed as an in-place version of the CSRCSC 
routine in FORMATS.
One notable disadvantage of  TRANSP is that unlike CSRCSC it 
does not sort the
nonzero elements in increasing number of the column positions.} 

\marg{ COPMAT  }\disp{Copy of a matrix into another 
matrix (both stored CSR).}

\marg{MSRCOP}\disp{Copies a matrix in MSR format into a matrix in MSR format.}
\marg{ GETELM }\disp{Function returning 
the value of $a_{ij}$ for any pair $(i,j)$. Also returns address
of the element in arrays $A, JA$. }

\marg{ GETDIA  }\disp{ Extracts a specified diagonal from a matrix.
An option is provided to transform the input matrix so that the 
extracted diagonal is zeroed out in input matrix. Otherwise the diagonal 
is extracted and the input matrix remains untouched.}

\marg{ GETL    }\disp{This subroutine extracts the
lower triangular part of a matrix, including the main diagonal.
The algorithm is in-place.}

\marg{ GETU    }\disp{Extracts the upper triangular part of
a matrix. Similar to GETL.}

\marg{ LEVELS  }\disp{ 
Computes the level scheduling data structure for lower
triangular matrices, see \cite{Anderson-Saad}.} 

\marg{ AMASK   }\disp{ Extracts    $ C = A \odot M  $,
i.e., performs the mask operation. 
This routine computes a sparse matrix from an input matrix $A$ by 
extracting only the elements in $A$, where the corresponding
elements of $M$ are nonzero. The mask matrix $M$, is a
sparse matrix in CSR format without the real values, i.e.,
only the integer arrays of the CSR format are passed. }

\marg{ CPERM   }\disp{ Permutes the columns of a matrix, i.e.,
computes the matrix $B = A Q$ where $Q$ is a permutation matrix. }

\marg{RPERM}\disp{ Permutes the rows of a matrix, i.e.,
computes the matrix $B = P A$ where $P$ is a permutation matrix. }

\marg{DPERM}\disp{Permutes the rows and columns of
a matrix, i.e., computes $B = P A Q$
given two permutation matrices  $ P$ and $ Q$. This routine gives a
special treatment to the common case where $Q=P^T$.} 

\marg{DPERM2}\disp{General submatrix permutation/extraction routine.}

\marg{DMPERM}\disp{Symmetric permutation of row and column (B=PAP') in MSR format}

\marg{DVPERM}\disp{ Performs an in-place permutation of a real vector, i.e.,
performs $x := P x $, where $P$ is a permutation matrix. }

\marg{IVPERM}\disp{ Performs an in-place permutation of an
integer vector.} 

\marg{RETMX}\disp{ Returns the maximum absolute value in each row  
of an input matrix.  } 

\marg{DIAPOS}\disp{ Returns the positions in the arrays $A$ and 
$ JA$ of the  diagonal elements, for a matrix stored in CSR format. } 

\marg{ EXTBDG  }\disp{ Extracts the main diagonal blocks of a matrix. 
The output is a rectangular matrix of dimension $N \times NBLK$,
containing the $N/NBLK$ blocks, 
in which $NBLK$ is the block-size (input).}

\marg{ GETBWD  }\disp{ Returns bandwidth information on a matrix. This
subroutine returns the bandwidth of the lower part and the upper part
of a given matrix. May be used to determine these two parameters
for converting a matrix into the BND format.}

\marg{ BLKFND  }\disp{ Attempts to find the block-size of a matrix
stored in CSR format. One restriction is that the zero elements in 
each block if there are any are assumed to be represented as nonzero 
elements in the data structure for the $A$ matrix, with zero values. }

\marg{ BLKCHK  }\disp{ Checks whether a given integer is the block 
size of A.  This routine is called by BLKFND. Same restriction as
above.}