
\documentclass{article}

\author{
Aart J.C. Bik
\thanks{ 
Computer Science Department, Indiana University, Lindley Hall 215,
Bloomington, Indiana 47405-4101, USA ({\tt ajcbik@cs.indiana.edu}).
}
\and Dennis B. Gannon
}

\title{A Note on Native Level~1 BLAS in Java}

\date{} 

\input{mystyle.tex}

\addtolength{\topmargin} {-2.0cm}
\addtolength{\textheight}{+4.0cm}
\addtolength{\textwidth}     {+2.0cm}
\addtolength{\oddsidemargin} {-1.0cm}
\addtolength{\evensidemargin}{-1.0cm}

\begin{document}

\maketitle 

\section{Introduction} 
Portability of the Java programming language~\cite{javalang}
is obtained by compiling Java programs into 
architectural neutral instructions (bytecode)
of the JVM (Java Virtual Machine)~\cite{javavm},
rather than into native machine code.  Bytecode runs on any platform
that supports an implementation of the JVM.
Although the interpretation of bytecode is substantially faster than
the interpretation of most high level languages,
still a performance penalty must be paid for this portability.
Clearly, for Java applications that are computational intensive,
it would be desirable to reduce this performance penalty
without sacrificing the portability of the language.
One approach, for example, is 
to optimize the Java bytecode~\cite{budimlic,cierniak},
either at compile-time (where a
machine independent bytecode to bytecode optimization 
is added as additional phase  to the Java compiler), or
at run-time (where optimizations that require
knowledge of the target machine are applied before execution).
Some implementations of the JVM further improve
performance  by means of JITC (just-in-time compilation),
where, at run-time, bytecode is compiled into native machine code.

In this note, we explore a way to speedup certain
operations in Java programs using ideas from the mathematical software community.
Here, it has been widely accepted
that adopting a set of basic routines for problems in linear
algebra can help in improving the clarity, portability,
modularity, maintenance, robustness, and even the efficiency of 
mathematical software.  The most well-known example of such a set of routines
is formed by the Basic Linear Algebra Subprograms~\cite[ch5]{dong}.
The original set of vector-vector operations is
now commonly referred to as Level~1 BLAS~\cite{blas1b,blas1a}.
The set has been extended to Level~2 BLAS~\cite{blas2a,blas2b} and Level~3 
BLAS~\cite{blas3a,blas3b} to provide more opportunities to exploit vector
processing facilities for matrix-vector operations
and memory hierarchies or parallelism
for matrix-matrix operations, respectively.
Once an efficient implementation of BLAS
is available, new mathematical software
can be easily build on top of the primitives.

Obviously,  a similar approach can be taken for Java by
extending the Java API (Application Programming Interface) with 
an appropriate set of mathematical primitives.
In first instance, a Java implementation can be provided 
for all these primitives to preserve the portability of all Java programs 
in which the mathematical primitives are used. 
On a particular machine, however,
the performance of all Java software that uses these
primitives is simply improved by providing native implementations 
of the mathematical primitives.
Although providing a broad range of highly optimized
mathematical primitives would offer the best potential to exploit all characteristic
of a particular target machine, this approach would
also require the most programming efforts to port the mathematical
primitives in the API to different machines. Therefore, in this
research note, we explore the potential of extending the
API with straightforward  native implementations of Level~1 BLAS only. We will
see that this extension alone already can improve
performance substantially, while combining these native Level~1 BLAS
with multi-threading in Java may even provide a simple and portable way 
to outperform compiled serial C code on multi-processors.

\section{Native BLAS}\label{native}
In Java, if the keyword {\tt native} appears as part of
a method definition (without an implementation), then
this implies that the method is implemented in another language. 
A set of primitives from BLAS, for instance,
can be defined using the following class 
that provides a shared library loader that will load
the actual implementation when required
as well as definitions of all primitives in the set:

\code
\begin{tabular}{l}
\verb"class Blas {" \\
\verb" " \\
\verb"  // Shared Library Loader" \\
\verb" " \\
\verb"  static {" \\
\verb!    System.loadLibrary("blas");!  \\
\verb"  }" \\
\verb" " \\
\verb"  // Level 1 BLAS" \\
\verb" " \\
\verb"  native static double ddot(int n, double[] x, int xoff, int incx," \\
\verb"                                                double[] y, int yoff, int incy);" \\
\verb"  native static void daxpy(int n, double alpha, double[] x, int xoff, int incx," \\
\verb"                                                double[] y, int yoff, int incy);" \\
\verb" " \\
\verb"  ..." \\
\verb"}" \\
\end{tabular}
\edoc

Here, method definitions of DDOT
($d \leftarrow \vec{x}\,^T \vec{y}$ in double precision)
and DAXPY ($\vec{y} \leftarrow \alpha \vec{x} + \vec{y}$
in double precision) are shown.\footnote{
Because Java does not support passing of 
arbitrary sub-arrays, offsets into arrays are added
as additional parameters (cf. the Java implementation
of BLAS in~\cite{java_linpack}).}
Alternatively, as stated earlier, Java implementations of the primitives
could be given in first instance.

\plaat{h}
\centerline{\resizebox{7cm}{!}{\includegraphics{nat.ps}}}
\taalp{Integration of Native Methods}{nat}

In figure~\ref{nat},  the integration
of native methods in Java is illustrated~\cite{native}.
Using JDK1.0.2, first the source file {\tt Blas.java} is compiled
into {\tt Blas.class} using the Java compiler {\tt javac}.
Subsequently, the tool {\tt javah}
is used to generate a stubs file {\tt Blas.c} and
a header file {\tt Blas.h}. The latter file contains C prototypes
for all native methods in the class.  The actual implementation
of these methods is given in a file {\tt BlasImp.c}.
Method {\tt ddot()}, for example,
can be implemented in a straightforward
manner in C as follows, where
the macro {\tt unhand()} is used to dereference an object handle:

\code
\begin{tabular}{l}
\verb!#include <StubPreamble.h>! \\
\verb!#include "Blas.h"! \\
\verb" " \\
\verb"double Blas_ddot(struct HBlas *this, long n," \\
\verb"                 HArrayOfDouble *x_, long xoff, long incx," \\
\verb"                 HArrayOfDouble *y_, long yoff, long incy) {" \\
\verb"  double *x = unhand(x_)->body;" \\
\verb"  double *y = unhand(y_)->body;" \\
\verb"  double  d = 0;" \\
\verb" " \\
\verb"  if (n > 0) {" \\
\verb"    if ((incx == 1) && (incy == 1)) { /* Unit strides */" \\
\verb"      int i;" \\
\verb"      for (i=0;i < n; i++)" \\
\verb"         d += x[xoff+i]*y[yoff+i];" \\
\verb"    }" \\
\verb"    else {                        /* Non-Unit strides */" \\
\verb"     ..." \\
\verb"    }" \\
\verb"  }" \\
\verb"  return d;" \\
\verb"}" \\
\end{tabular}
\edoc

Eventually, the files {\tt Blas.c} and {\tt BlasImp.c} are
compiled into a shared library. The way in which this is done
and the naming convention of shared libraries depends
on the target architecture. Under Solaris, for instance,
the shared library is created as follows:

\code
\verb"cc -G -O  -I $JAVA_HOME/include -I $JAVA_HOME/include/solaris Blas.c BlasImp.c -o libblas.so"
\edoc

As illustrated in the second picture in figure~\ref{nat},
integrating native methods in JDK1.1 is slightly different.
Only  a header file {\tt Blas.h} which contains the appropriate
prototypes is generated using  {\tt javah -jni}. The actual implementation
of methods is, again, given in a file {\tt BlasImp.c}.
Method {\tt ddot()}, for example,
can be implemented in C as follows,
where functions {\tt GetDoubleArrayElements()} 
and {\tt ReleaseDoubleArrayElements()} 
are used to obtain and release
a double array that cannot be moved by the garbage collector.
Eventually, the whole file is
compiled into a shared library: 

\code
\begin{tabular}{l}
\verb!#include "Blas.h"! \\
\verb"JNIEXPORT jdouble JNICALL Java_Blas_ddot" \\
\verb"    (JNIEnv *env, jclass class, jint n, jdoubleArray x_, jint xoff, jint incx," \\
\verb"                                        jdoubleArray y_, jint yoff, jint incy) {" \\
\verb"  jdouble *x = (*env)->GetDoubleArrayElements(env, x_, 0);" \\
\verb"  jdouble *y = (*env)->GetDoubleArrayElements(env, y_, 0);" \\
\verb"  jdouble  d = 0;" \\
\verb" " \\
\verb"  if (n > 0) {" \\
\verb"    if ((incx == 1) && (incy == 1)) { /* Unit strides */" \\
\verb"      int i;" \\
\verb"      for (i=0;i < n; i++)" \\
\verb"        d += x[xoff+i]*y[yoff+i];" \\
\verb"    }" \\
\verb"    else {                         /* Non-Unit strides */" \\
\verb"      ..." \\
\verb"    } " \\
\verb"  }" \\
\verb"  (*env)->ReleaseDoubleArrayElements(env, x_, x, 0);" \\
\verb"  (*env)->ReleaseDoubleArrayElements(env, y_, y, 0);" \\
\verb"  return d;" \\
\verb"}" \\
\end{tabular}
\edoc

If native methods are used, care must be taken
with respect to the precision of the computed result,
because \emph{Java data types are mapped to the
nearest matching data types on the native language
side}.  

\section{Experiments}\label{experiments}
In this section, we present the results of a series
of experiments that have been conducted on
an IBM RS/6000 G30 with four PowerPC 604 processors 
using the AIX4.2 JDK1.0.2B (with JITC)
and the  AIX4.2 JDK1.1beta (without JITC),
a Sun with two 175MHz. ultra SPARC processors
using the Solaris 2.5 JDK1.0.2dp (with JITC),   and
an SGI Indy 4600 with one
133 MHZ IP22 Processor
using the IRIX6.2 JDK1.0.2 (without JITC).

All Java and C programs are compiled using the flag `{\tt -O}'.
Bytecode is interpreted using the flag `-noasyncgc' and, if available,
with JITC enabled.
We compare the performance of a pure Java implementation, a
Java implementation that uses native Level~1 BLAS
(and possibly parallelism by means of Java multi-threading),
and a compiled C implementation of
an $\Theta(N)$ DAXPY-operation, an
$\Theta(N^2)$ matrix times vector operation,
and an $\Theta(N^3)$ matrix times matrix operation.

Below we show the pure Java implementation 
of the DAXPY operation ($\vec{y} \leftarrow \alpha \vec{x} + \vec{y}$),
and the Java implementation that uses
a primitive from Level~1 BLAS:

\code
\begin{tabular}{l}
pure Java: \\
\verb" " \\
\verb"for (int i = 0; i < N; i++)" \\
\verb"  y[i] += alpha * x[i];" \\
\end{tabular}
\begin{tabular}{l}
Java + Level~1 BLAS: \\
\verb" " \\
\verb"Blas.daxpy(N, alpha, x, 0, 1, y, 0, 1);" \\
\verb" " \\
\end{tabular}
\edoc

The inner product Java implementation of
matrix times vector and a Java implementation that uses 
a call to Level~1 BLAS DDOT ($d \leftarrow \vec{x}\,^T \vec{y}$)
are shown below:

\code
\begin{tabular}{l}
pure Java: \\
\verb" " \\
\verb"for (int i = 0; i < N; i++)" \\
\verb"  for (int j = 0; j < N; j++)" \\
\verb"    b[i] += a[i][j] * x[j];" \\
\end{tabular}
\begin{tabular}{l}
Java + Level~1 BLAS: \\
\verb" " \\
\verb"for (int i = 0; i < N; i++)" \\
\verb"  b[i] = Blas.ddot(N, a[i], 0, 1, x, 0, 1);" \\
\verb" " \\
\end{tabular}
\edoc

Finally, the product of two matrices can be computed
using either the pure Java fragment shown below, or
a similar Java fragment that uses a call to DAXPY:

\code
\begin{tabular}{l}
pure Java: \\
\verb" " \\
\verb"for (int i = 0; i < N; i++)" \\
\verb"  for (int k = 0; k < N; k++)" \\
\verb"    for (int j = 0; j < N; j++)" \\
\verb"      c[i][j] += a[i][k] * b[k][j];" \\
\end{tabular}
\begin{tabular}{l}
Java + Level~1 BLAS: \\
\verb" " \\
\verb"for (int i = 0; i < N; i++)" \\
\verb"  for (int k = 0; k < N; k++)" \\
\verb"    Blas.daxpy(N, a[i][k], b[k], 0, 1, c[i], 0, 1);" \\
\verb" " \\
\end{tabular}
\edoc

The performance of the fragments for matrix times vector and matrix times matrix
with Level~1 BLAS may be further improved by parallelization of
the outermost {\tt i}-loop. In earlier work~\cite{bikgana,javar}, we have
shown how parallel loops can be expressed in Java  by means
of multi-threading.  In this manner the transformed
Java program remains portable (i.e. the parallelized version
still runs on uni-processors with only a slight overhead),
while the programming efforts of the parallelization are reduced substantially
with respect to exploiting parallelism in a native language.
The versions with a parallel outermost loop are labeled as `parallel' 
in subsequent figures. 
The C versions of the three mathematical operations
are similar to the pure Java implementations.

In figures~\ref{ibm_daxpy_0}--\ref{sgi_mm_0},
we present the results of all experiments.
On the IBM, we see that with JITC, providing Level~1 BLAS primitives only
suffices to obtain performance that is close to the performance
of native C code. Moreover, because IBM's implementation
of the JVM supports the actual parallel execution 
of threads, the parallel versions with native Level~1 BLAS
even outperform serial C code. Without JITC, however, the
$\Theta(N^3)$ operation still suffers from much overhead,
and here it would probably be desirable to provide
primitives from Level~2 BLAS as well. On the Sun,
the performance of Java with native methods is frequently even slightly better than
the performance of straightforward
native C code. Obviously, Sun's implementation of the JVM does not
support the actual parallel execution of threads.
Finally, on the SGI without JITC,  the Java performance of
$\Theta(N^3)$ matrix times matrix operations 
is substantially less than the performance of compiled C code.
Obviously, on this uni-processor,
no speedup can be expected from loop parallelization.

\section{Conclusions}\label{conclusions}
In this research note, we have explored the potential of extending the Java
Application Programming Interface (API) with a 
native  implementation of Level~1 BLAS to improve the performance of 
certain mathematical operations in Java programs.
Because, in first instance, a Java implementation of the
primitives can be provided,
the portability of Java is maintained.

We have shown that for implementations of the 
JVM that support JITC,
providing straightforward native Level~1 BLAS already suffices to
obtain performance that is close to the performance
of straightforward implementations in native C code for $\Theta(N)$,
$\Theta(N^2)$, and $\Theta(N^3)$ mathematical operations.
In case JITC is not supported, the performance difference
of the latter operations is more profound, and here it would
be desirable to  provide higher level primitives as well.
In addition, we have shown that, on multi-processors,
combining native Level~1 BLAS 
with the multi-threading mechanism
of Java may provide a simple and portable way to obtain a
Java program that runs faster than compiled serial C code.

In conclusion, although it is obvious that more speedup
can be expected by  providing highly optimized
native implementations of a broad range of mathematical primitives in which
characteristics of the target machine
(such as multiple processors, memory hierarchies or 
vector processing facilities) are fully exploited at native
level, providing straightforward native implementations of Level~1 BLAS
alone already can help in improving the performance of
mathematical operations in Java. In all cases where native
methods are used, however, care must be taken with respect
to the precision of the computed result, because Java data types are mapped to the
nearest matching data types on the native language side.  

{\bf Acknowledgments}
This project is supported by DARPA under contract ARPA F19628-94-C-0057
through a subcontract from Syracuse University. 

\footnotesize

\nocite{carpenter,fox1}
\bibliographystyle{plain}  
\bibliography{biblio}

\clearpage

\twocolumn

\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{ibm_daxpy_0.ps}}}
\taalp{DAXPY on the IBM}{ibm_daxpy_0}
\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{ibm_mv_0.ps}}}
\taalp{Matrix x Vector on the IBM}{ibm_mv_0}
\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{ibm_mm_0.ps}}}
\taalp{Matrix x Matrix on the IBM}{ibm_mm_0}

\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{ibm_daxpy_1.ps}}}
\taalp{DAXPY on the IBM}{ibm_daxpy_1}
\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{ibm_mv_1.ps}}}
\taalp{Matrix x Vector on the IBM}{ibm_mv_1}
\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{ibm_mm_1.ps}}}
\taalp{Matrix x Matrix on the IBM}{ibm_mm_1}

\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{sol_daxpy_0.ps}}}
\taalp{DAXPY on the Sun}{sol_daxpy_0}
\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{sol_mv_0.ps}}}
\taalp{Matrix x Vector on the Sun}{sol_mv_0}
\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{sol_mm_0.ps}}}
\taalp{Matrix x Matrix on the Sun}{sol_mm_0}

\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{sgi_daxpy_0.ps}}}
\taalp{DAXPY on the SGI}{sgi_daxpy_0}
\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{sgi_mv_0.ps}}}
\taalp{Matrix x Vector on the SGI}{sgi_mv_0}
\plaat{g}
\centerline{\resizebox{7cm}{!}{\includegraphics{sgi_mm_0.ps}}}
\taalp{Matrix x Matrix on the SGI}{sgi_mm_0}

\end{document}
