增加GKLib

third_party/METIS/.gitignore

@@ -57,6 +57,5 @@ graphs/*.part.*
 graphs/*.iperm
 graphs/*.epart.*
 graphs/*.npart.*
-GKlib
 .svn/
 

third_party/METIS/GKlib/CMakeLists.txt

+cmake_minimum_required(VERSION 2.8)
+project(GKlib C)
+
+get_filename_component(abs "." ABSOLUTE)
+set(GKLIB_PATH ${abs})
+unset(abs)
+include(GKlibSystem.cmake)
+
+include_directories(".")
+add_library(GKlib STATIC ${GKlib_sources})
+if(UNIX)
+  target_link_libraries(GKlib m)
+endif(UNIX)
+
+include_directories("test")
+add_subdirectory("test")
+
+install(TARGETS GKlib
+  ARCHIVE DESTINATION lib/${LINSTALL_PATH}
+  LIBRARY DESTINATION lib/${LINSTALL_PATH})
+install(FILES ${GKlib_includes} DESTINATION include/${HINSTALL_PATH})

third_party/METIS/GKlib/GKlib.h

+/*
+ * GKlib.h
+ * 
+ * George's library of most frequently used routines
+ *
+ * $Id: GKlib.h 14866 2013-08-03 16:40:04Z karypis $
+ *
+ */
+
+#ifndef _GKLIB_H_
+#define _GKLIB_H_ 1
+
+#define GKMSPACE
+
+#if defined(_MSC_VER)
+#define __MSC__
+#endif
+#if defined(__ICC)
+#define __ICC__
+#endif
+
+
+#include "gk_arch.h" /*!< This should be here, prior to the includes */
+
+
+/*************************************************************************
+* Header file inclusion section
+**************************************************************************/
+#include <stddef.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <stdio.h>
+#include <memory.h>
+#include <errno.h>
+#include <ctype.h>
+#include <math.h>
+#include <float.h>
+#include <time.h>
+#include <string.h>
+#include <unistd.h>
+#include <limits.h>
+#include <signal.h>
+#include <setjmp.h>
+#include <assert.h>
+#include <sys/stat.h>
+
+#if defined(__WITHPCRE__)
+  #include <pcreposix.h>
+#else
+  #if defined(USE_GKREGEX)
+    #include "gkregex.h"
+  #else
+    #include <regex.h>
+  #endif /* defined(USE_GKREGEX) */
+#endif /* defined(__WITHPCRE__) */
+
+
+
+#if defined(__OPENMP__) 
+#include <omp.h>
+#endif
+
+
+
+
+#include <gk_types.h>
+#include <gk_struct.h>
+#include <gk_externs.h>
+#include <gk_defs.h>
+#include <gk_macros.h>
+#include <gk_getopt.h>
+
+#include <gk_mksort.h>
+#include <gk_mkblas.h>
+#include <gk_mkmemory.h>
+#include <gk_mkpqueue.h>
+#include <gk_mkpqueue2.h>
+#include <gk_mkrandom.h>
+#include <gk_mkutils.h>
+
+#include <gk_proto.h>
+
+
+#endif  /* GKlib.h */
+
+

third_party/METIS/GKlib/GKlibSystem.cmake

+# Helper modules.
+include(CheckFunctionExists)
+include(CheckIncludeFile)
+
+# Setup options.
+option(GDB "enable use of GDB" OFF)
+option(ASSERT "turn asserts on" OFF)
+option(ASSERT2 "additional assertions" OFF)
+option(DEBUG "add debugging support" OFF)
+option(GPROF "add gprof support" OFF)
+option(OPENMP "enable OpenMP support" OFF)
+option(PCRE "enable PCRE support" OFF)
+option(GKREGEX "enable GKREGEX support" OFF)
+option(GKRAND "enable GKRAND support" OFF)
+
+
+# Add compiler flags.
+if(MSVC)
+  set(GKlib_COPTS "/Ox")
+  set(GKlib_COPTIONS "-DWIN32 -DMSC -D_CRT_SECURE_NO_DEPRECATE -DUSE_GKREGEX")
+elseif(MINGW)
+  set(GKlib_COPTS "-DUSE_GKREGEX")
+else()
+  set(GKlib_COPTIONS "-DLINUX -D_FILE_OFFSET_BITS=64")
+endif(MSVC)
+if(CYGWIN)
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -DCYGWIN")
+endif(CYGWIN)
+if(CMAKE_COMPILER_IS_GNUCC)
+# GCC opts.
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -std=c99 -fno-strict-aliasing")
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -march=native")
+  if(NOT MINGW)
+      set(GKlib_COPTIONS "${GKlib_COPTIONS} -fPIC")
+  endif(NOT MINGW)
+# GCC warnings.
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -Werror -Wall -pedantic -Wno-unused-function -Wno-unused-but-set-variable -Wno-unused-variable -Wno-unknown-pragmas")
+elseif(${CMAKE_C_COMPILER_ID} MATCHES "Sun")
+# Sun insists on -xc99.
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -xc99")
+endif(CMAKE_COMPILER_IS_GNUCC)
+
+# Intel compiler
+if(${CMAKE_C_COMPILER_ID} MATCHES "Intel")
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -xHost -std=c99")
+endif()
+
+# Find OpenMP if it is requested.
+if(OPENMP)
+  include(FindOpenMP)
+  if(OPENMP_FOUND)
+    set(GKlib_COPTIONS "${GKlib_COPTIONS} -D__OPENMP__ ${OpenMP_C_FLAGS}")
+  else()
+    message(WARNING "OpenMP was requested but support was not found")
+  endif(OPENMP_FOUND)
+endif(OPENMP)
+
+
+# Add various definitions.
+if(GDB)
+  set(GKlib_COPTS "${GKlib_COPTS} -g")
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -Werror")
+else()
+  set(GKlib_COPTS "-O3")
+endif(GDB)
+
+
+if(DEBUG)
+  set(GKlib_COPTS "-g")
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -DDEBUG")
+endif(DEBUG)
+
+if(GPROF)
+  set(GKlib_COPTS "-pg")
+endif(GPROF)
+
+if(NOT ASSERT)
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -DNDEBUG")
+endif(NOT ASSERT)
+
+if(NOT ASSERT2)
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -DNDEBUG2")
+endif(NOT ASSERT2)
+
+
+# Add various options
+if(PCRE)
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -D__WITHPCRE__")
+endif(PCRE)
+
+if(GKREGEX)
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -DUSE_GKREGEX")
+endif(GKREGEX)
+
+if(GKRAND)
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -DUSE_GKRAND")
+endif(GKRAND)
+
+
+# Check for features.
+check_include_file(execinfo.h HAVE_EXECINFO_H)
+if(HAVE_EXECINFO_H)
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -DHAVE_EXECINFO_H")
+endif(HAVE_EXECINFO_H)
+
+check_function_exists(getline HAVE_GETLINE)
+if(HAVE_GETLINE)
+  set(GKlib_COPTIONS "${GKlib_COPTIONS} -DHAVE_GETLINE")
+endif(HAVE_GETLINE)
+
+
+# Custom check for TLS.
+if(MSVC)
+   set(GKlib_COPTIONS "${GKlib_COPTIONS} -D__thread=__declspec(thread)")
+
+  # This if checks if that value is cached or not.
+  if("${HAVE_THREADLOCALSTORAGE}" MATCHES "^${HAVE_THREADLOCALSTORAGE}$")
+    try_compile(HAVE_THREADLOCALSTORAGE
+      ${CMAKE_BINARY_DIR}
+      ${GKLIB_PATH}/conf/check_thread_storage.c)
+    if(HAVE_THREADLOCALSTORAGE)
+      message(STATUS "checking for thread-local storage - found")
+    else()
+      message(STATUS "checking for thread-local storage - not found")
+    endif()
+  endif()
+  if(NOT HAVE_THREADLOCALSTORAGE)
+    set(GKlib_COPTIONS "${GKlib_COPTIONS} -D__thread=")
+  endif()
+endif()
+
+# Finally set the official C flags.
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${GKlib_COPTIONS} ${GKlib_COPTS}")
+
+# Find GKlib sources.
+file(GLOB GKlib_sources ${GKLIB_PATH}/*.c)
+file(GLOB GKlib_includes ${GKLIB_PATH}/*.h)

third_party/METIS/GKlib/LICENSE.txt

+
+Copyright & License Notice
+---------------------------
+
+Copyright 1995-2018, Regents of the University of Minnesota
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or 
+implied. See the License for the specific language governing 
+permissions and limitations under the License.
+

third_party/METIS/GKlib/Makefile

+# Configuration options.
+cc       = gcc
+prefix   = ~/local
+openmp   = not-set
+gdb      = not-set
+assert   = not-set
+assert2  = not-set
+debug    = not-set
+gprof    = not-set
+pcre     = not-set
+gkregex  = not-set
+gkrand   = not-set
+
+
+# Basically proxies everything to the builddir cmake.
+cputype = $(shell uname -m | sed "s/\\ /_/g")
+systype = $(shell uname -s)
+
+BUILDDIR = build/$(systype)-$(cputype)
+
+# Process configuration options.
+CONFIG_FLAGS = -DCMAKE_VERBOSE_MAKEFILE=1
+ifneq ($(gdb), not-set)
+    CONFIG_FLAGS += -DGDB=$(gdb)
+endif
+ifneq ($(assert), not-set)
+    CONFIG_FLAGS += -DASSERT=$(assert)
+endif
+ifneq ($(assert2), not-set)
+    CONFIG_FLAGS += -DASSERT2=$(assert2)
+endif
+ifneq ($(debug), not-set)
+    CONFIG_FLAGS += -DDEBUG=$(debug)
+endif
+ifneq ($(gprof), not-set)
+    CONFIG_FLAGS += -DGPROF=$(gprof)
+endif
+ifneq ($(openmp), not-set)
+    CONFIG_FLAGS += -DOPENMP=$(openmp)
+endif
+ifneq ($(pcre), not-set)
+    CONFIG_FLAGS += -DPCRE=$(pcre)
+endif
+ifneq ($(gkregex), not-set)
+    CONFIG_FLAGS += -DGKREGEX=$(pcre)
+endif
+ifneq ($(gkrand), not-set)
+    CONFIG_FLAGS += -DGKRAND=$(pcre)
+endif
+ifneq ($(prefix), not-set)
+    CONFIG_FLAGS += -DCMAKE_INSTALL_PREFIX=$(prefix)
+endif
+ifneq ($(cc), not-set)
+    CONFIG_FLAGS += -DCMAKE_C_COMPILER=$(cc)
+endif
+
+define run-config
+mkdir -p $(BUILDDIR)
+cd $(BUILDDIR) && cmake $(CURDIR) $(CONFIG_FLAGS)
+endef
+
+all clean install: $(BUILDDIR)
+	make -C $(BUILDDIR) $@
+
+uninstall:
+	 xargs rm < $(BUILDDIR)/install_manifest.txt
+
+$(BUILDDIR):
+	$(run-config)
+
+config: distclean
+	$(run-config)
+
+distclean:
+	rm -rf $(BUILDDIR)
+
+remake:
+	find . -name CMakeLists.txt -exec touch {} ';'
+
+.PHONY: config distclean all clean install uninstall remake

third_party/METIS/GKlib/README.md

+# GKlib
+A library of various helper routines and frameworks used by many of the lab's software
+
+## Build requirements
+ - CMake 2.8, found at http://www.cmake.org/, as well as GNU make. 
+
+Assuming that the above are available, two commands should suffice to 
+build the software:
+```
+make config 
+make
+```
+
+## Configuring the build
+It is primarily configured by passing options to make config. For example:
+```
+make config cc=icc
+```
+
+would configure it to be built using icc.
+
+Configuration options are:
+```
+cc=[compiler]     - The C compiler to use [default: gcc]
+prefix=[PATH]     - Set the installation prefix [default: ~/local]
+openmp=set        - To build a version with OpenMP support
+```
+
+
+## Building and installing
+To build and install, run the following
+```
+make
+make install
+```
+
+By default, the library file, header file, and binaries will be installed in
+```
+~/local/lib
+~/local/include
+~/local/bin
+```
+
+## Other make commands
+    make uninstall 
+         Removes all files installed by 'make install'.
+   
+    make clean 
+         Removes all object files but retains the configuration options.
+   
+    make distclean 
+         Performs clean and completely removes the build directory.
+
+

third_party/METIS/GKlib/b64.c

+/*! 
+\file  b64.c
+\brief This file contains some simple 8bit-to-6bit encoding/deconding routines
+
+Most of these routines are outdated and should be converted using glibc's equivalent
+routines.
+
+\date   Started 2/22/05
+\author George
+\version\verbatim $Id: b64.c 10711 2011-08-31 22:23:04Z karypis $ \endverbatim
+
+\verbatim 
+$Copyright$ 
+$License$
+\endverbatim
+
+*/
+
+
+#include "GKlib.h"
+
+#define B64OFFSET       48      /* This is the '0' number */
+
+
+/******************************************************************************
+* Encode 3 '8-bit' binary bytes as 4 '6-bit' characters
+*******************************************************************************/
+void encodeblock(unsigned char *in, unsigned char *out)
+{
+  out[0] = (in[0] >> 2);
+  out[1] = (((in[0] & 0x03) << 4) | (in[1] >> 4));
+  out[2] = (((in[1] & 0x0f) << 2) | (in[2] >> 6));
+  out[3] = (in[2] & 0x3f);
+
+  out[0] += B64OFFSET;
+  out[1] += B64OFFSET;
+  out[2] += B64OFFSET;
+  out[3] += B64OFFSET;
+
+//  printf("%c %c %c %c %2x %2x %2x %2x %2x %2x %2x\n", out[0], out[1], out[2], out[3], out[0], out[1], out[2], out[3], in[0], in[1], in[2]);
+}
+
+/******************************************************************************
+* Decode 4 '6-bit' characters into 3 '8-bit' binary bytes
+*******************************************************************************/
+void decodeblock(unsigned char *in, unsigned char *out)
+{   
+  in[0] -= B64OFFSET;
+  in[1] -= B64OFFSET;
+  in[2] -= B64OFFSET;
+  in[3] -= B64OFFSET;
+
+  out[0] = (in[0] << 2 | in[1] >> 4);
+  out[1] = (in[1] << 4 | in[2] >> 2);
+  out[2] = (in[2] << 6 | in[3]);
+}
+
+
+/******************************************************************************
+* This function encodes an input array of bytes into a base64 encoding. Memory
+* for the output array is assumed to have been allocated by the calling program
+* and be sufficiently large. The output string is NULL terminated.
+*******************************************************************************/
+void GKEncodeBase64(int nbytes, unsigned char *inbuffer, unsigned char *outbuffer)
+{
+  int i, j;
+
+  if (nbytes%3 != 0)
+    gk_errexit(SIGERR, "GKEncodeBase64: Input buffer size should be a multiple of 3! (%d)\n", nbytes);
+
+  for (j=0, i=0; i<nbytes; i+=3, j+=4) 
+    encodeblock(inbuffer+i, outbuffer+j);
+
+//printf("%d %d\n", nbytes, j);
+  outbuffer[j] = '\0';
+}
+
+
+
+/******************************************************************************
+* This function decodes an input array of base64 characters into their actual
+* 8-bit codes. Memory * for the output array is assumed to have been allocated 
+* by the calling program and be sufficiently large. The padding is discarded.
+*******************************************************************************/
+void GKDecodeBase64(int nbytes, unsigned char *inbuffer, unsigned char *outbuffer)
+{
+  int i, j;
+
+  if (nbytes%4 != 0)
+    gk_errexit(SIGERR, "GKDecodeBase64: Input buffer size should be a multiple of 4! (%d)\n", nbytes);
+
+  for (j=0, i=0; i<nbytes; i+=4, j+=3) 
+    decodeblock(inbuffer+i, outbuffer+j);
+}
+

third_party/METIS/GKlib/blas.c

+/*!
+\file blas.c
+\brief This file contains GKlib's implementation of BLAS-like routines
+
+The BLAS routines that are currently implemented are mostly level-one.
+They follow a naming convention of the type gk_[type][name], where
+[type] is one of c, i, f, and d, based on C's four standard scalar
+datatypes of characters, integers, floats, and doubles.
+
+These routines are implemented using a generic macro template,
+which is used for code generation.
+
+\date   Started 9/28/95
+\author George
+\version\verbatim $Id: blas.c 14330 2013-05-18 12:15:15Z karypis $ \endverbatim
+*/
+
+#include <GKlib.h>
+
+
+
+/*************************************************************************/
+/*! Use the templates to generate BLAS routines for the scalar data types */
+/*************************************************************************/
+GK_MKBLAS(gk_c,   char,     int)
+GK_MKBLAS(gk_i,   int,      int)
+GK_MKBLAS(gk_i32, int32_t,  int32_t)
+GK_MKBLAS(gk_i64, int64_t,  int64_t)
+GK_MKBLAS(gk_z,   ssize_t,  ssize_t)
+GK_MKBLAS(gk_zu,  size_t,   size_t)
+GK_MKBLAS(gk_f,   float,    float)
+GK_MKBLAS(gk_d,   double,   double)
+GK_MKBLAS(gk_idx, gk_idx_t, gk_idx_t)
+
+
+
+

third_party/METIS/GKlib/conf/check_thread_storage.c

+extern __thread int x;
+
+int main(int argc, char **argv) {
+  return 0;
+}

third_party/METIS/GKlib/csr.c

+/*!
+ * \file 
+ *
+ * \brief Various routines with dealing with CSR matrices
+ *
+ * \author George Karypis
+ * \version\verbatim $Id: csr.c 21044 2017-05-24 22:50:32Z karypis $ \endverbatim
+ */
+
+#include <GKlib.h>
+
+#define OMPMINOPS       50000
+
+/*************************************************************************/
+/*! Allocate memory for a CSR matrix and initializes it 
+    \returns the allocated matrix. The various fields are set to NULL.
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_Create()
+{
+  gk_csr_t *mat=NULL;
+
+  if ((mat = (gk_csr_t *)gk_malloc(sizeof(gk_csr_t), "gk_csr_Create: mat")))
+    gk_csr_Init(mat);
+
+  return mat;
+}
+
+
+/*************************************************************************/
+/*! Initializes the matrix 
+    \param mat is the matrix to be initialized.
+*/
+/*************************************************************************/
+void gk_csr_Init(gk_csr_t *mat)
+{
+  memset(mat, 0, sizeof(gk_csr_t));
+  mat->nrows = mat->ncols = 0;
+}
+
+
+/*************************************************************************/
+/*! Frees all the memory allocated for matrix.
+    \param mat is the matrix to be freed.
+*/
+/*************************************************************************/
+void gk_csr_Free(gk_csr_t **mat)
+{
+  if (*mat == NULL)
+    return;
+  gk_csr_FreeContents(*mat);
+  gk_free((void **)mat, LTERM);
+}
+
+
+/*************************************************************************/
+/*! Frees only the memory allocated for the matrix's different fields and
+    sets them to NULL.
+    \param mat is the matrix whose contents will be freed.
+*/    
+/*************************************************************************/
+void gk_csr_FreeContents(gk_csr_t *mat)
+{
+  gk_free((void *)&mat->rowptr, &mat->rowind, &mat->rowval, 
+      &mat->rowids, &mat->rlabels, &mat->rmap,
+      &mat->colptr, &mat->colind, &mat->colval, 
+      &mat->colids, &mat->clabels, &mat->cmap,
+      &mat->rnorms, &mat->cnorms, &mat->rsums, &mat->csums, 
+      &mat->rsizes, &mat->csizes, &mat->rvols, &mat->cvols, 
+      &mat->rwgts, &mat->cwgts, 
+          LTERM);
+}
+
+
+/*************************************************************************/
+/*! Returns a copy of a matrix.
+    \param mat is the matrix to be duplicated.
+    \returns the newly created copy of the matrix.
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_Dup(gk_csr_t *mat)
+{
+  gk_csr_t *nmat;
+
+  nmat = gk_csr_Create();
+
+  nmat->nrows  = mat->nrows;
+  nmat->ncols  = mat->ncols;
+
+  /* copy the row structure */
+  if (mat->rowptr)
+    nmat->rowptr = gk_zcopy(mat->nrows+1, mat->rowptr, 
+                            gk_zmalloc(mat->nrows+1, "gk_csr_Dup: rowptr"));
+  if (mat->rowids)
+    nmat->rowids = gk_icopy(mat->nrows, mat->rowids, 
+                            gk_imalloc(mat->nrows, "gk_csr_Dup: rowids"));
+  if (mat->rlabels)
+    nmat->rlabels = gk_icopy(mat->nrows, mat->rlabels, 
+                            gk_imalloc(mat->nrows, "gk_csr_Dup: rlabels"));
+  if (mat->rnorms)
+    nmat->rnorms = gk_fcopy(mat->nrows, mat->rnorms, 
+                            gk_fmalloc(mat->nrows, "gk_csr_Dup: rnorms"));
+  if (mat->rsums)
+    nmat->rsums = gk_fcopy(mat->nrows, mat->rsums, 
+                            gk_fmalloc(mat->nrows, "gk_csr_Dup: rsums"));
+  if (mat->rsizes)
+    nmat->rsizes = gk_fcopy(mat->nrows, mat->rsizes, 
+                            gk_fmalloc(mat->nrows, "gk_csr_Dup: rsizes"));
+  if (mat->rvols)
+    nmat->rvols = gk_fcopy(mat->nrows, mat->rvols, 
+                            gk_fmalloc(mat->nrows, "gk_csr_Dup: rvols"));
+  if (mat->rwgts)
+    nmat->rwgts = gk_fcopy(mat->nrows, mat->rwgts, 
+                            gk_fmalloc(mat->nrows, "gk_csr_Dup: rwgts"));
+  if (mat->rowind)
+    nmat->rowind = gk_icopy(mat->rowptr[mat->nrows], mat->rowind, 
+                            gk_imalloc(mat->rowptr[mat->nrows], "gk_csr_Dup: rowind"));
+  if (mat->rowval)
+    nmat->rowval = gk_fcopy(mat->rowptr[mat->nrows], mat->rowval, 
+                            gk_fmalloc(mat->rowptr[mat->nrows], "gk_csr_Dup: rowval"));
+
+  /* copy the col structure */
+  if (mat->colptr)
+    nmat->colptr = gk_zcopy(mat->ncols+1, mat->colptr, 
+                            gk_zmalloc(mat->ncols+1, "gk_csr_Dup: colptr"));
+  if (mat->colids)
+    nmat->colids = gk_icopy(mat->ncols, mat->colids, 
+                            gk_imalloc(mat->ncols, "gk_csr_Dup: colids"));
+  if (mat->clabels)
+    nmat->clabels = gk_icopy(mat->ncols, mat->clabels, 
+                            gk_imalloc(mat->ncols, "gk_csr_Dup: clabels"));
+  if (mat->cnorms)
+    nmat->cnorms = gk_fcopy(mat->ncols, mat->cnorms, 
+                            gk_fmalloc(mat->ncols, "gk_csr_Dup: cnorms"));
+  if (mat->csums)
+    nmat->csums = gk_fcopy(mat->ncols, mat->csums, 
+                            gk_fmalloc(mat->ncols, "gk_csr_Dup: csums"));
+  if (mat->csizes)
+    nmat->csizes = gk_fcopy(mat->ncols, mat->csizes, 
+                            gk_fmalloc(mat->ncols, "gk_csr_Dup: csizes"));
+  if (mat->cvols)
+    nmat->cvols = gk_fcopy(mat->ncols, mat->cvols, 
+                            gk_fmalloc(mat->ncols, "gk_csr_Dup: cvols"));
+  if (mat->cwgts)
+    nmat->cwgts = gk_fcopy(mat->ncols, mat->cwgts, 
+                            gk_fmalloc(mat->ncols, "gk_csr_Dup: cwgts"));
+  if (mat->colind)
+    nmat->colind = gk_icopy(mat->colptr[mat->ncols], mat->colind, 
+                            gk_imalloc(mat->colptr[mat->ncols], "gk_csr_Dup: colind"));
+  if (mat->colval)
+    nmat->colval = gk_fcopy(mat->colptr[mat->ncols], mat->colval, 
+                            gk_fmalloc(mat->colptr[mat->ncols], "gk_csr_Dup: colval"));
+
+  return nmat;
+}
+
+
+/*************************************************************************/
+/*! Returns a submatrix containint a set of consecutive rows.
+    \param mat is the original matrix.
+    \param rstart is the starting row.
+    \param nrows is the number of rows from rstart to extract.
+    \returns the row structure of the newly created submatrix.
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_ExtractSubmatrix(gk_csr_t *mat, int rstart, int nrows)
+{
+  ssize_t i;
+  gk_csr_t *nmat;
+
+  if (rstart+nrows > mat->nrows)
+    return NULL;
+
+  nmat = gk_csr_Create();
+
+  nmat->nrows  = nrows;
+  nmat->ncols  = mat->ncols;
+
+  /* copy the row structure */
+  if (mat->rowptr)
+    nmat->rowptr = gk_zcopy(nrows+1, mat->rowptr+rstart, 
+                              gk_zmalloc(nrows+1, "gk_csr_ExtractSubmatrix: rowptr"));
+  for (i=nrows; i>=0; i--)
+    nmat->rowptr[i] -= nmat->rowptr[0];
+  ASSERT(nmat->rowptr[0] == 0);
+
+  if (mat->rowids)
+    nmat->rowids = gk_icopy(nrows, mat->rowids+rstart, 
+                            gk_imalloc(nrows, "gk_csr_ExtractSubmatrix: rowids"));
+  if (mat->rnorms)
+    nmat->rnorms = gk_fcopy(nrows, mat->rnorms+rstart, 
+                            gk_fmalloc(nrows, "gk_csr_ExtractSubmatrix: rnorms"));
+
+  if (mat->rsums)
+    nmat->rsums = gk_fcopy(nrows, mat->rsums+rstart, 
+                            gk_fmalloc(nrows, "gk_csr_ExtractSubmatrix: rsums"));
+
+  ASSERT(nmat->rowptr[nrows] == mat->rowptr[rstart+nrows]-mat->rowptr[rstart]);
+  if (mat->rowind)
+    nmat->rowind = gk_icopy(mat->rowptr[rstart+nrows]-mat->rowptr[rstart], 
+                            mat->rowind+mat->rowptr[rstart], 
+                            gk_imalloc(mat->rowptr[rstart+nrows]-mat->rowptr[rstart],
+                                       "gk_csr_ExtractSubmatrix: rowind"));
+  if (mat->rowval)
+    nmat->rowval = gk_fcopy(mat->rowptr[rstart+nrows]-mat->rowptr[rstart], 
+                            mat->rowval+mat->rowptr[rstart], 
+                            gk_fmalloc(mat->rowptr[rstart+nrows]-mat->rowptr[rstart],
+                                       "gk_csr_ExtractSubmatrix: rowval"));
+
+  return nmat;
+}
+
+
+/*************************************************************************/
+/*! Returns a submatrix containing a certain set of rows.
+    \param mat is the original matrix.
+    \param nrows is the number of rows to extract.
+    \param rind is the set of row numbers to extract.
+    \returns the row structure of the newly created submatrix.
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_ExtractRows(gk_csr_t *mat, int nrows, int *rind)
+{
+  ssize_t i, ii, j, nnz;
+  gk_csr_t *nmat;
+
+  nmat = gk_csr_Create();
+
+  nmat->nrows = nrows;
+  nmat->ncols = mat->ncols;
+
+  for (nnz=0, i=0; i<nrows; i++)  
+    nnz += mat->rowptr[rind[i]+1]-mat->rowptr[rind[i]];
+
+  nmat->rowptr = gk_zmalloc(nmat->nrows+1, "gk_csr_ExtractPartition: rowptr");
+  nmat->rowind = gk_imalloc(nnz, "gk_csr_ExtractPartition: rowind");
+  nmat->rowval = gk_fmalloc(nnz, "gk_csr_ExtractPartition: rowval");
+
+  nmat->rowptr[0] = 0;
+  for (nnz=0, j=0, ii=0; ii<nrows; ii++) {
+    i = rind[ii];
+    gk_icopy(mat->rowptr[i+1]-mat->rowptr[i], mat->rowind+mat->rowptr[i], nmat->rowind+nnz);
+    gk_fcopy(mat->rowptr[i+1]-mat->rowptr[i], mat->rowval+mat->rowptr[i], nmat->rowval+nnz);
+    nnz += mat->rowptr[i+1]-mat->rowptr[i];
+    nmat->rowptr[++j] = nnz;
+  }
+  ASSERT(j == nmat->nrows);
+
+  return nmat;
+}
+
+
+/*************************************************************************/
+/*! Returns a submatrix corresponding to a specified partitioning of rows.
+    \param mat is the original matrix.
+    \param part is the partitioning vector of the rows.
+    \param pid is the partition ID that will be extracted.
+    \returns the row structure of the newly created submatrix.
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_ExtractPartition(gk_csr_t *mat, int *part, int pid)
+{
+  ssize_t i, j, nnz;
+  gk_csr_t *nmat;
+
+  nmat = gk_csr_Create();
+
+  nmat->nrows = 0;
+  nmat->ncols = mat->ncols;
+
+  for (nnz=0, i=0; i<mat->nrows; i++) {
+    if (part[i] == pid) {
+      nmat->nrows++;
+      nnz += mat->rowptr[i+1]-mat->rowptr[i];
+    }
+  }
+
+  nmat->rowptr = gk_zmalloc(nmat->nrows+1, "gk_csr_ExtractPartition: rowptr");
+  nmat->rowind = gk_imalloc(nnz, "gk_csr_ExtractPartition: rowind");
+  nmat->rowval = gk_fmalloc(nnz, "gk_csr_ExtractPartition: rowval");
+
+  nmat->rowptr[0] = 0;
+  for (nnz=0, j=0, i=0; i<mat->nrows; i++) {
+    if (part[i] == pid) {
+      gk_icopy(mat->rowptr[i+1]-mat->rowptr[i], mat->rowind+mat->rowptr[i], nmat->rowind+nnz);
+      gk_fcopy(mat->rowptr[i+1]-mat->rowptr[i], mat->rowval+mat->rowptr[i], nmat->rowval+nnz);
+      nnz += mat->rowptr[i+1]-mat->rowptr[i];
+      nmat->rowptr[++j] = nnz;
+    }
+  }
+  ASSERT(j == nmat->nrows);
+
+  return nmat;
+}
+
+
+/*************************************************************************/
+/*! Splits the matrix into multiple sub-matrices based on the provided
+    color array.
+    \param mat is the original matrix.
+    \param color is an array of size equal to the number of non-zeros
+           in the matrix (row-wise structure). The matrix is split into
+           as many parts as the number of colors. For meaningfull results,
+           the colors should be numbered consecutively starting from 0.
+    \returns an array of matrices for each supplied color number.
+*/
+/**************************************************************************/
+gk_csr_t **gk_csr_Split(gk_csr_t *mat, int *color)
+{
+  ssize_t i, j;
+  int nrows, ncolors;
+  ssize_t *rowptr;
+  int *rowind;
+  float *rowval;
+  gk_csr_t **smats;
+
+  nrows  = mat->nrows;
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+  rowval = mat->rowval;
+
+  ncolors = gk_imax(rowptr[nrows], color, 1)+1;
+
+  smats = (gk_csr_t **)gk_malloc(sizeof(gk_csr_t *)*ncolors, "gk_csr_Split: smats");
+  for (i=0; i<ncolors; i++) {
+    smats[i] = gk_csr_Create();
+    smats[i]->nrows  = mat->nrows;
+    smats[i]->ncols  = mat->ncols;
+    smats[i]->rowptr = gk_zsmalloc(nrows+1, 0, "gk_csr_Split: smats[i]->rowptr"); 
+  }
+
+  for (i=0; i<nrows; i++) {
+    for (j=rowptr[i]; j<rowptr[i+1]; j++) 
+      smats[color[j]]->rowptr[i]++;
+  }
+  for (i=0; i<ncolors; i++) 
+    MAKECSR(j, nrows, smats[i]->rowptr);
+
+  for (i=0; i<ncolors; i++) {
+    smats[i]->rowind = gk_imalloc(smats[i]->rowptr[nrows], "gk_csr_Split: smats[i]->rowind"); 
+    smats[i]->rowval = gk_fmalloc(smats[i]->rowptr[nrows], "gk_csr_Split: smats[i]->rowval"); 
+  }
+
+  for (i=0; i<nrows; i++) {
+    for (j=rowptr[i]; j<rowptr[i+1]; j++) {
+      smats[color[j]]->rowind[smats[color[j]]->rowptr[i]] = rowind[j];
+      smats[color[j]]->rowval[smats[color[j]]->rowptr[i]] = rowval[j];
+      smats[color[j]]->rowptr[i]++;
+    }
+  }
+
+  for (i=0; i<ncolors; i++) 
+    SHIFTCSR(j, nrows, smats[i]->rowptr);
+
+  return smats;
+}
+
+
+/**************************************************************************/
+/*! Determines the format of the CSR matrix based on the extension.
+    \param filename is the name of the file.
+    \param the user-supplied format.
+    \returns the type. The extension of the file directly maps to the
+           name of the format.
+*/
+/**************************************************************************/
+int gk_csr_DetermineFormat(char *filename, int format)
+{
+  if (format != GK_CSR_FMT_AUTO)
+    return format;
+
+  format = GK_CSR_FMT_CSR;
+  char *extension = gk_getextname(filename);
+
+  if (!strcmp(extension, "csr"))
+    format = GK_CSR_FMT_CSR;
+  else if (!strcmp(extension, "ijv"))
+    format = GK_CSR_FMT_IJV;
+  else if (!strcmp(extension, "cluto"))
+    format = GK_CSR_FMT_CLUTO;
+  else if (!strcmp(extension, "metis"))
+    format = GK_CSR_FMT_METIS;
+  else if (!strcmp(extension, "binrow"))
+    format = GK_CSR_FMT_BINROW;
+  else if (!strcmp(extension, "bincol"))
+    format = GK_CSR_FMT_BINCOL;
+  else if (!strcmp(extension, "bijv"))
+    format = GK_CSR_FMT_BIJV;
+
+  gk_free((void **)&extension, LTERM);
+
+  return format;
+}
+
+
+/**************************************************************************/
+/*! Reads a CSR matrix from the supplied file and stores it the matrix's 
+    forward structure.
+    \param filename is the file that stores the data.
+    \param format is either GK_CSR_FMT_METIS, GK_CSR_FMT_CLUTO, 
+           GK_CSR_FMT_CSR, GK_CSR_FMT_BINROW, GK_CSR_FMT_BINCOL 
+           specifying the type of the input format. 
+           The GK_CSR_FMT_CSR does not contain a header
+           line, whereas the GK_CSR_FMT_BINROW is a binary format written 
+           by gk_csr_Write() using the same format specifier.
+    \param readvals is either 1 or 0, indicating if the CSR file contains
+           values or it does not. It only applies when GK_CSR_FMT_CSR is
+           used.
+    \param numbering is either 1 or 0, indicating if the numbering of the 
+           indices start from 1 or 0, respectively. If they start from 1, 
+           they are automatically decreamented during input so that they
+           will start from 0. It only applies when GK_CSR_FMT_CSR is
+           used.
+    \returns the matrix that was read.
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_Read(char *filename, int format, int readvals, int numbering)
+{
+  ssize_t i, k, l;
+  size_t nfields, nrows, ncols, nnz, fmt, ncon;
+  size_t lnlen;
+  ssize_t *rowptr;
+  int *rowind, *iinds, *jinds, ival;
+  float *rowval=NULL, *vals, fval;
+  int readsizes, readwgts;
+  char *line=NULL, *head, *tail, fmtstr[256];
+  FILE *fpin;
+  gk_csr_t *mat=NULL;
+
+  format = gk_csr_DetermineFormat(filename, format);
+
+  if (!gk_fexists(filename)) 
+    gk_errexit(SIGERR, "File %s does not exist!\n", filename);
+
+  switch (format) {
+    case GK_CSR_FMT_BINROW:
+      mat = gk_csr_Create();
+
+      fpin = gk_fopen(filename, "rb", "gk_csr_Read: fpin");
+      if (fread(&(mat->nrows), sizeof(int32_t), 1, fpin) != 1)
+        gk_errexit(SIGERR, "Failed to read the nrows from file %s!\n", filename);
+      if (fread(&(mat->ncols), sizeof(int32_t), 1, fpin) != 1)
+        gk_errexit(SIGERR, "Failed to read the ncols from file %s!\n", filename);
+      mat->rowptr = gk_zmalloc(mat->nrows+1, "gk_csr_Read: rowptr");
+      if (fread(mat->rowptr, sizeof(ssize_t), mat->nrows+1, fpin) != mat->nrows+1)
+        gk_errexit(SIGERR, "Failed to read the rowptr from file %s!\n", filename);
+      mat->rowind = gk_imalloc(mat->rowptr[mat->nrows], "gk_csr_Read: rowind");
+      if (fread(mat->rowind, sizeof(int32_t), mat->rowptr[mat->nrows], fpin) != mat->rowptr[mat->nrows])
+        gk_errexit(SIGERR, "Failed to read the rowind from file %s!\n", filename);
+      if (readvals == 1) {
+        mat->rowval = gk_fmalloc(mat->rowptr[mat->nrows], "gk_csr_Read: rowval");
+        if (fread(mat->rowval, sizeof(float), mat->rowptr[mat->nrows], fpin) != mat->rowptr[mat->nrows])
+          gk_errexit(SIGERR, "Failed to read the rowval from file %s!\n", filename);
+      }
+
+      gk_fclose(fpin);
+      return mat;
+
+      break;
+
+    case GK_CSR_FMT_BINCOL:
+      mat = gk_csr_Create();
+
+      fpin = gk_fopen(filename, "rb", "gk_csr_Read: fpin");
+      if (fread(&(mat->nrows), sizeof(int32_t), 1, fpin) != 1)
+        gk_errexit(SIGERR, "Failed to read the nrows from file %s!\n", filename);
+      if (fread(&(mat->ncols), sizeof(int32_t), 1, fpin) != 1)
+        gk_errexit(SIGERR, "Failed to read the ncols from file %s!\n", filename);
+      mat->colptr = gk_zmalloc(mat->ncols+1, "gk_csr_Read: colptr");
+      if (fread(mat->colptr, sizeof(ssize_t), mat->ncols+1, fpin) != mat->ncols+1)
+        gk_errexit(SIGERR, "Failed to read the colptr from file %s!\n", filename);
+      mat->colind = gk_imalloc(mat->colptr[mat->ncols], "gk_csr_Read: colind");
+      if (fread(mat->colind, sizeof(int32_t), mat->colptr[mat->ncols], fpin) != mat->colptr[mat->ncols])
+        gk_errexit(SIGERR, "Failed to read the colind from file %s!\n", filename);
+      if (readvals) {
+        mat->colval = gk_fmalloc(mat->colptr[mat->ncols], "gk_csr_Read: colval");
+        if (fread(mat->colval, sizeof(float), mat->colptr[mat->ncols], fpin) != mat->colptr[mat->ncols])
+          gk_errexit(SIGERR, "Failed to read the colval from file %s!\n", filename);
+      }
+
+      gk_fclose(fpin);
+      return mat;
+
+      break;
+
+
+    case GK_CSR_FMT_IJV:
+      gk_getfilestats(filename, &nrows, &nnz, NULL, NULL);
+
+      if (readvals == 1 && 3*nrows != nnz)
+        gk_errexit(SIGERR, "Error: The number of numbers (%zd %d) in the input file is not a multiple of 3.\n", nnz, readvals);
+      if (readvals == 0 && 2*nrows != nnz)
+        gk_errexit(SIGERR, "Error: The number of numbers (%zd %d) in the input file is not a multiple of 2.\n", nnz, readvals);
+
+      nnz = nrows;
+      numbering = (numbering ? - 1 : 0);
+
+      /* read the data into three arrays */
+      iinds = gk_i32malloc(nnz, "iinds");
+      jinds = gk_i32malloc(nnz, "jinds");
+      vals  = (readvals ? gk_fmalloc(nnz, "vals") : NULL);
+
+      fpin = gk_fopen(filename, "r", "gk_csr_Read: fpin");
+      for (nrows=0, ncols=0, i=0; i<nnz; i++) {
+        if (readvals) {
+          if (fscanf(fpin, "%d %d %f", &iinds[i], &jinds[i], &vals[i]) != 3)
+            gk_errexit(SIGERR, "Error: Failed to read (i, j, val) for nnz: %zd.\n", i);
+        }
+        else {
+          if (fscanf(fpin, "%d %d", &iinds[i], &jinds[i]) != 2)
+            gk_errexit(SIGERR, "Error: Failed to read (i, j) value for nnz: %zd.\n", i);
+        }
+        iinds[i] += numbering;
+        jinds[i] += numbering;
+
+        if (nrows < iinds[i])
+          nrows = iinds[i];
+        if (ncols < jinds[i])
+          ncols = jinds[i];
+      }
+      nrows++;
+      ncols++;
+      gk_fclose(fpin);
+
+      /* convert (i, j, v) into a CSR matrix */
+      mat = gk_csr_Create();
+      mat->nrows = nrows;
+      mat->ncols = ncols;
+      rowptr = mat->rowptr = gk_zsmalloc(nrows+1, 0, "rowptr");
+      rowind = mat->rowind = gk_i32malloc(nnz, "rowind");
+      if (readvals)
+        rowval = mat->rowval = gk_fmalloc(nnz, "rowval");
+
+      for (i=0; i<nnz; i++)
+        rowptr[iinds[i]]++;
+      MAKECSR(i, nrows, rowptr);
+
+      for (i=0; i<nnz; i++) {
+        rowind[rowptr[iinds[i]]] = jinds[i];
+        if (readvals)
+          rowval[rowptr[iinds[i]]] = vals[i];
+        rowptr[iinds[i]]++;
+      }
+      SHIFTCSR(i, nrows, rowptr);
+
+      gk_free((void **)&iinds, &jinds, &vals, LTERM);
+
+      return mat;
+
+      break;
+
+    case GK_CSR_FMT_BIJV:
+      mat = gk_csr_Create();
+
+      fpin = gk_fopen(filename, "rb", "gk_csr_Read: fpin");
+
+      if (fread(&(mat->nrows), sizeof(int32_t), 1, fpin) != 1)
+        gk_errexit(SIGERR, "Failed to read the nrows from file %s!\n", filename);
+      if (fread(&(mat->ncols), sizeof(int32_t), 1, fpin) != 1)
+        gk_errexit(SIGERR, "Failed to read the ncols from file %s!\n", filename);
+      if (fread(&nnz, sizeof(size_t), 1, fpin) != 1)
+        gk_errexit(SIGERR, "Failed to read the nnz from file %s!\n", filename);
+      if (fread(&readvals, sizeof(int32_t), 1, fpin) != 1)
+        gk_errexit(SIGERR, "Failed to read the readvals from file %s!\n", filename);
+
+      /* read the data into three arrays */
+      iinds = gk_i32malloc(nnz, "iinds");
+      jinds = gk_i32malloc(nnz, "jinds");
+      vals  = (readvals ? gk_fmalloc(nnz, "vals") : NULL);
+
+      for (i=0; i<nnz; i++) {
+        if (fread(&(iinds[i]), sizeof(int32_t), 1, fpin) != 1)
+          gk_errexit(SIGERR, "Failed to read iinds[i] from file %s!\n", filename);
+        if (fread(&(jinds[i]), sizeof(int32_t), 1, fpin) != 1)
+          gk_errexit(SIGERR, "Failed to read jinds[i] from file %s!\n", filename);
+        if (readvals) {
+          if (fread(&(vals[i]), sizeof(float), 1, fpin) != 1)
+            gk_errexit(SIGERR, "Failed to read vals[i] from file %s!\n", filename);
+        }
+        //printf("%d %d\n", iinds[i], jinds[i]);
+      }
+      gk_fclose(fpin);
+
+      /* convert (i, j, v) into a CSR matrix */
+      rowptr = mat->rowptr = gk_zsmalloc(mat->nrows+1, 0, "rowptr");
+      rowind = mat->rowind = gk_i32malloc(nnz, "rowind");
+      if (readvals)
+        rowval = mat->rowval = gk_fmalloc(nnz, "rowval");
+
+      for (i=0; i<nnz; i++)
+        rowptr[iinds[i]]++;
+      MAKECSR(i, mat->nrows, rowptr);
+
+      for (i=0; i<nnz; i++) {
+        rowind[rowptr[iinds[i]]] = jinds[i];
+        if (readvals)
+          rowval[rowptr[iinds[i]]] = vals[i];
+        rowptr[iinds[i]]++;
+      }
+      SHIFTCSR(i, mat->nrows, rowptr);
+
+      gk_free((void **)&iinds, &jinds, &vals, LTERM);
+
+      return mat;
+
+      break;
+
+
+    /* the following are handled by a common input code, that comes after the switch */
+
+    case GK_CSR_FMT_CLUTO:
+      fpin = gk_fopen(filename, "r", "gk_csr_Read: fpin");
+      do {
+        if (gk_getline(&line, &lnlen, fpin) <= 0)
+          gk_errexit(SIGERR, "Premature end of input file: file:%s\n", filename);
+      } while (line[0] == '%');
+
+      if (sscanf(line, "%zu %zu %zu", &nrows, &ncols, &nnz) != 3)
+        gk_errexit(SIGERR, "Header line must contain 3 integers.\n");
+
+      readsizes = 0;
+      readwgts  = 0;
+      readvals  = 1;
+      numbering = 1;
+
+      break;
+
+    case GK_CSR_FMT_METIS:
+      fpin = gk_fopen(filename, "r", "gk_csr_Read: fpin");
+      do {
+        if (gk_getline(&line, &lnlen, fpin) <= 0)
+          gk_errexit(SIGERR, "Premature end of input file: file:%s\n", filename);
+      } while (line[0] == '%');
+
+      fmt = ncon = 0;
+      nfields = sscanf(line, "%zu %zu %zu %zu", &nrows, &nnz, &fmt, &ncon);
+      if (nfields < 2)
+        gk_errexit(SIGERR, "Header line must contain at least 2 integers (#vtxs and #edges).\n");
+
+      ncols = nrows;
+      nnz *= 2;
+
+      if (fmt > 111)
+        gk_errexit(SIGERR, "Cannot read this type of file format [fmt=%zu]!\n", fmt);
+
+      sprintf(fmtstr, "%03zu", fmt%1000);
+      readsizes = (fmtstr[0] == '1');
+      readwgts  = (fmtstr[1] == '1');
+      readvals  = (fmtstr[2] == '1');
+      numbering = 1;
+      ncon      = (ncon == 0 ? 1 : ncon);
+
+      break;
+
+    case GK_CSR_FMT_CSR:
+      readsizes = 0;
+      readwgts  = 0;
+
+      gk_getfilestats(filename, &nrows, &nnz, NULL, NULL);
+
+      if (readvals == 1 && nnz%2 == 1)
+        gk_errexit(SIGERR, "Error: The number of numbers (%zd %d) in the input file is not even.\n", nnz, readvals);
+      if (readvals == 1)
+        nnz = nnz/2;
+      fpin = gk_fopen(filename, "r", "gk_csr_Read: fpin");
+
+      break;
+
+    default:
+      gk_errexit(SIGERR, "Unknown csr format.\n");
+      return NULL;
+  }
+
+  mat = gk_csr_Create();
+
+  mat->nrows = nrows;
+
+  rowptr = mat->rowptr = gk_zmalloc(nrows+1, "gk_csr_Read: rowptr");
+  rowind = mat->rowind = gk_imalloc(nnz, "gk_csr_Read: rowind");
+  if (readvals != 2)
+    rowval = mat->rowval = gk_fsmalloc(nnz, 1.0, "gk_csr_Read: rowval");
+
+  if (readsizes)
+    mat->rsizes = gk_fsmalloc(nrows, 0.0, "gk_csr_Read: rsizes");
+
+  if (readwgts)
+    mat->rwgts = gk_fsmalloc(nrows*ncon, 0.0, "gk_csr_Read: rwgts");
+
+  /*----------------------------------------------------------------------
+   * Read the sparse matrix file
+   *---------------------------------------------------------------------*/
+  numbering = (numbering ? -1 : 0);
+  for (ncols=0, rowptr[0]=0, k=0, i=0; i<nrows; i++) {
+    do {
+      if (gk_getline(&line, &lnlen, fpin) == -1)
+        gk_errexit(SIGERR, "Premature end of input file: file while reading row %d\n", i);
+    } while (line[0] == '%');
+
+    head = line;
+    tail = NULL;
+
+    /* Read vertex sizes */
+    if (readsizes) {
+#ifdef __MSC__
+      mat->rsizes[i] = (float)strtod(head, &tail);
+#else
+      mat->rsizes[i] = strtof(head, &tail);
+#endif
+      if (tail == head)
+        gk_errexit(SIGERR, "The line for vertex %zd does not have size information\n", i+1);
+      if (mat->rsizes[i] < 0)
+        errexit("The size for vertex %zd must be >= 0\n", i+1);
+      head = tail;
+    }
+
+    /* Read vertex weights */
+    if (readwgts) {
+      for (l=0; l<ncon; l++) {
+#ifdef __MSC__
+        mat->rwgts[i*ncon+l] = (float)strtod(head, &tail);
+#else
+        mat->rwgts[i*ncon+l] = strtof(head, &tail);
+#endif
+        if (tail == head)
+          errexit("The line for vertex %zd does not have enough weights "
+                  "for the %d constraints.\n", i+1, ncon);
+        if (mat->rwgts[i*ncon+l] < 0)
+          errexit("The weight vertex %zd and constraint %zd must be >= 0\n", i+1, l);
+        head = tail;
+      }
+    }
+
+   
+    /* Read the rest of the row */
+    while (1) {
+      ival = (int)strtol(head, &tail, 0);
+      if (tail == head) 
+        break;
+      head = tail;
+      
+      if ((rowind[k] = ival + numbering) < 0)
+        gk_errexit(SIGERR, "Error: Invalid column number %d at row %zd.\n", ival, i);
+
+      ncols = gk_max(rowind[k], ncols);
+
+      if (readvals == 1) {
+#ifdef __MSC__
+        fval = (float)strtod(head, &tail);
+#else
+	fval = strtof(head, &tail);
+#endif
+        if (tail == head)
+          gk_errexit(SIGERR, "Value could not be found for column! Row:%zd, NNZ:%zd\n", i, k);
+        head = tail;
+
+        rowval[k] = fval;
+      }
+      k++;
+    }
+    rowptr[i+1] = k;
+  }
+
+  if (format == GK_CSR_FMT_METIS) {
+    ASSERT(ncols+1 == mat->nrows);
+    mat->ncols = mat->nrows;
+  }
+  else {
+    mat->ncols = ncols+1;
+  }
+
+  if (k != nnz)
+    gk_errexit(SIGERR, "gk_csr_Read: Something wrong with the number of nonzeros in "
+                       "the input file. NNZ=%zd, ActualNNZ=%zd.\n", nnz, k);
+
+  gk_fclose(fpin);
+
+  gk_free((void **)&line, LTERM);
+
+  return mat;
+}
+
+
+/**************************************************************************/
+/*! Writes the row-based structure of a matrix into a file.
+    \param mat is the matrix to be written,
+    \param filename is the name of the output file.
+    \param format is one of: GK_CSR_FMT_CLUTO, GK_CSR_FMT_CSR, 
+           GK_CSR_FMT_BINROW, GK_CSR_FMT_BINCOL, GK_CSR_FMT_BIJV.
+    \param writevals is either 1 or 0 indicating if the values will be 
+           written or not. This is only applicable when GK_CSR_FMT_CSR
+           is used.
+    \param numbering is either 1 or 0 indicating if the internal 0-based 
+           numbering will be shifted by one or not during output. This 
+           is only applicable when GK_CSR_FMT_CSR is used.
+*/
+/**************************************************************************/
+void gk_csr_Write(gk_csr_t *mat, char *filename, int format, int writevals, int numbering)
+{
+  ssize_t i, j;
+  int32_t edge[2];
+  FILE *fpout;
+
+  format = gk_csr_DetermineFormat(filename, format);
+
+  switch (format) {
+    case GK_CSR_FMT_METIS:
+      if (mat->nrows != mat->ncols || mat->rowptr[mat->nrows]%2 == 1)
+        gk_errexit(SIGERR, "METIS output format requires a square symmetric matrix.\n");
+
+      if (filename)
+        fpout = gk_fopen(filename, "w", "gk_csr_Write: fpout");
+      else
+        fpout = stdout; 
+
+      fprintf(fpout, "%d %zd\n", mat->nrows, mat->rowptr[mat->nrows]/2);
+      for (i=0; i<mat->nrows; i++) {
+        for (j=mat->rowptr[i]; j<mat->rowptr[i+1]; j++) 
+          fprintf(fpout, " %d", mat->rowind[j]+1);
+        fprintf(fpout, "\n");
+      }
+      if (filename)
+        gk_fclose(fpout);
+      break;
+
+    case GK_CSR_FMT_BINROW:
+      if (filename == NULL)
+        gk_errexit(SIGERR, "The filename parameter cannot be NULL.\n");
+      fpout = gk_fopen(filename, "wb", "gk_csr_Write: fpout");
+
+      fwrite(&(mat->nrows), sizeof(int32_t), 1, fpout); 
+      fwrite(&(mat->ncols), sizeof(int32_t), 1, fpout); 
+      fwrite(mat->rowptr, sizeof(ssize_t), mat->nrows+1, fpout); 
+      fwrite(mat->rowind, sizeof(int32_t), mat->rowptr[mat->nrows], fpout); 
+      if (writevals)
+        fwrite(mat->rowval, sizeof(float), mat->rowptr[mat->nrows], fpout); 
+
+      gk_fclose(fpout);
+      return;
+
+      break;
+
+    case GK_CSR_FMT_BINCOL:
+      if (filename == NULL)
+        gk_errexit(SIGERR, "The filename parameter cannot be NULL.\n");
+      fpout = gk_fopen(filename, "wb", "gk_csr_Write: fpout");
+
+      fwrite(&(mat->nrows), sizeof(int32_t), 1, fpout); 
+      fwrite(&(mat->ncols), sizeof(int32_t), 1, fpout); 
+      fwrite(mat->colptr, sizeof(ssize_t), mat->ncols+1, fpout); 
+      fwrite(mat->colind, sizeof(int32_t), mat->colptr[mat->ncols], fpout); 
+      if (writevals) 
+        fwrite(mat->colval, sizeof(float), mat->colptr[mat->ncols], fpout); 
+
+      gk_fclose(fpout);
+      return;
+
+      break;
+
+    case GK_CSR_FMT_IJV:
+      if (filename == NULL)
+        gk_errexit(SIGERR, "The filename parameter cannot be NULL.\n");
+      fpout = gk_fopen(filename, "w", "gk_csr_Write: fpout");
+
+      numbering = (numbering ? 1 : 0);
+      for (i=0; i<mat->nrows; i++) {
+        for (j=mat->rowptr[i]; j<mat->rowptr[i+1]; j++) {
+          if (writevals)
+            fprintf(fpout, "%zd %d %.8f\n", i+numbering, mat->rowind[j]+numbering, mat->rowval[j]);
+          else
+            fprintf(fpout, "%zd %d\n", i+numbering, mat->rowind[j]+numbering);
+        }
+      }
+
+      gk_fclose(fpout);
+      return;
+
+      break;
+
+    case GK_CSR_FMT_BIJV:
+      if (filename == NULL)
+        gk_errexit(SIGERR, "The filename parameter cannot be NULL.\n");
+      fpout = gk_fopen(filename, "wb", "gk_csr_Write: fpout");
+
+      fwrite(&(mat->nrows), sizeof(int32_t), 1, fpout); 
+      fwrite(&(mat->ncols), sizeof(int32_t), 1, fpout); 
+      fwrite(&(mat->rowptr[mat->nrows]), sizeof(size_t), 1, fpout); 
+      fwrite(&writevals, sizeof(int32_t), 1, fpout); 
+
+      for (i=0; i<mat->nrows; i++) {
+        edge[0] = i;
+        for (j=mat->rowptr[i]; j<mat->rowptr[i+1]; j++) {
+          edge[1] = mat->rowind[j];
+          fwrite(edge, sizeof(int32_t), 2, fpout);
+          if (writevals) 
+            fwrite(&(mat->rowval[j]), sizeof(float), 1, fpout);
+        }
+      }
+
+      gk_fclose(fpout);
+      return;
+
+      break;
+
+    default:
+      if (filename)
+        fpout = gk_fopen(filename, "w", "gk_csr_Write: fpout");
+      else
+        fpout = stdout; 
+
+      if (format == GK_CSR_FMT_CLUTO) {
+        fprintf(fpout, "%d %d %zd\n", mat->nrows, mat->ncols, mat->rowptr[mat->nrows]);
+        writevals = 1;
+        numbering = 1;
+      }
+
+      for (i=0; i<mat->nrows; i++) {
+        for (j=mat->rowptr[i]; j<mat->rowptr[i+1]; j++) {
+          fprintf(fpout, " %d", mat->rowind[j]+(numbering ? 1 : 0));
+          if (writevals) 
+            fprintf(fpout, " %f", mat->rowval[j]);
+        }
+        fprintf(fpout, "\n");
+      }
+      if (filename)
+        gk_fclose(fpout);
+  }
+}
+
+
+/*************************************************************************/
+/*! Prunes certain rows/columns of the matrix. The prunning takes place 
+    by analyzing the row structure of the matrix. The prunning takes place
+    by removing rows/columns but it does not affect the numbering of the
+    remaining rows/columns.
+   
+    \param mat the matrix to be prunned,
+    \param what indicates if the rows (GK_CSR_ROW) or the columns (GK_CSR_COL)
+           of the matrix will be prunned,
+    \param minf is the minimum number of rows (columns) that a column (row) must
+           be present in order to be kept,
+    \param maxf is the maximum number of rows (columns) that a column (row) must
+          be present at in order to be kept.
+    \returns the prunned matrix consisting only of its row-based structure. 
+          The input matrix is not modified. 
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_Prune(gk_csr_t *mat, int what, int minf, int maxf)
+{
+  ssize_t i, j, nnz;
+  int nrows, ncols;
+  ssize_t *rowptr, *nrowptr;
+  int *rowind, *nrowind, *collen;
+  float *rowval, *nrowval;
+  gk_csr_t *nmat;
+
+  nmat = gk_csr_Create();
+  
+  nrows = nmat->nrows = mat->nrows;
+  ncols = nmat->ncols = mat->ncols;
+
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+  rowval = mat->rowval;
+
+  nrowptr = nmat->rowptr = gk_zmalloc(nrows+1, "gk_csr_Prune: nrowptr");
+  nrowind = nmat->rowind = gk_imalloc(rowptr[nrows], "gk_csr_Prune: nrowind");
+  nrowval = nmat->rowval = gk_fmalloc(rowptr[nrows], "gk_csr_Prune: nrowval");
+
+
+  switch (what) {
+    case GK_CSR_COL:
+      collen = gk_ismalloc(ncols, 0, "gk_csr_Prune: collen");
+
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++) {
+          ASSERT(rowind[j] < ncols);
+          collen[rowind[j]]++;
+        }
+      }
+      for (i=0; i<ncols; i++)
+        collen[i] = (collen[i] >= minf && collen[i] <= maxf ? 1 : 0);
+
+      nrowptr[0] = 0;
+      for (nnz=0, i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++) {
+          if (collen[rowind[j]]) {
+            nrowind[nnz] = rowind[j];
+            nrowval[nnz] = rowval[j];
+            nnz++;
+          }
+        }
+        nrowptr[i+1] = nnz;
+      }
+      gk_free((void **)&collen, LTERM);
+      break;
+
+    case GK_CSR_ROW:
+      nrowptr[0] = 0;
+      for (nnz=0, i=0; i<nrows; i++) {
+        if (rowptr[i+1]-rowptr[i] >= minf && rowptr[i+1]-rowptr[i] <= maxf) {
+          for (j=rowptr[i]; j<rowptr[i+1]; j++, nnz++) {
+            nrowind[nnz] = rowind[j];
+            nrowval[nnz] = rowval[j];
+          }
+        }
+        nrowptr[i+1] = nnz;
+      }
+      break;
+
+    default:
+      gk_csr_Free(&nmat);
+      gk_errexit(SIGERR, "Unknown prunning type of %d\n", what);
+      return NULL;
+  }
+
+  return nmat;
+}
+
+
+/*************************************************************************/
+/*! Eliminates certain entries from the rows/columns of the matrix. The 
+    filtering takes place by keeping only the highest weight entries whose
+    sum accounts for a certain fraction of the overall weight of the 
+    row/column.
+   
+    \param mat the matrix to be prunned,
+    \param what indicates if the rows (GK_CSR_ROW) or the columns (GK_CSR_COL)
+           of the matrix will be prunned,
+    \param norm indicates the norm that will be used to aggregate the weights
+           and possible values are 1 or 2,
+    \param fraction is the fraction of the overall norm that will be retained
+           by the kept entries.
+    \returns the filtered matrix consisting only of its row-based structure. 
+           The input matrix is not modified. 
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_LowFilter(gk_csr_t *mat, int what, int norm, float fraction)
+{
+  ssize_t i, j, nnz;
+  int nrows, ncols, ncand, maxlen=0;
+  ssize_t *rowptr, *colptr, *nrowptr;
+  int *rowind, *colind, *nrowind;
+  float *rowval, *colval, *nrowval, rsum, tsum;
+  gk_csr_t *nmat;
+  gk_fkv_t *cand;
+
+  nmat = gk_csr_Create();
+  
+  nrows = nmat->nrows = mat->nrows;
+  ncols = nmat->ncols = mat->ncols;
+
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+  rowval = mat->rowval;
+  colptr = mat->colptr;
+  colind = mat->colind;
+  colval = mat->colval;
+
+  nrowptr = nmat->rowptr = gk_zmalloc(nrows+1, "gk_csr_LowFilter: nrowptr");
+  nrowind = nmat->rowind = gk_imalloc(rowptr[nrows], "gk_csr_LowFilter: nrowind");
+  nrowval = nmat->rowval = gk_fmalloc(rowptr[nrows], "gk_csr_LowFilter: nrowval");
+
+
+  switch (what) {
+    case GK_CSR_COL:
+      if (mat->colptr == NULL) 
+        gk_errexit(SIGERR, "Cannot filter columns when column-based structure has not been created.\n");
+
+      gk_zcopy(nrows+1, rowptr, nrowptr);
+
+      for (i=0; i<ncols; i++) 
+        maxlen = gk_max(maxlen, colptr[i+1]-colptr[i]);
+
+      #pragma omp parallel private(i, j, ncand, rsum, tsum, cand)
+      {
+        cand = gk_fkvmalloc(maxlen, "gk_csr_LowFilter: cand");
+
+        #pragma omp for schedule(static)
+        for (i=0; i<ncols; i++) {
+          for (tsum=0.0, ncand=0, j=colptr[i]; j<colptr[i+1]; j++, ncand++) {
+            cand[ncand].val = colind[j];
+            cand[ncand].key = colval[j];
+            tsum += (norm == 1 ? colval[j] : colval[j]*colval[j]);
+          }
+          gk_fkvsortd(ncand, cand);
+
+          for (rsum=0.0, j=0; j<ncand && rsum<=fraction*tsum; j++) {
+            rsum += (norm == 1 ? cand[j].key : cand[j].key*cand[j].key);
+            nrowind[nrowptr[cand[j].val]] = i;
+            nrowval[nrowptr[cand[j].val]] = cand[j].key;
+            nrowptr[cand[j].val]++;
+          }
+        }
+
+        gk_free((void **)&cand, LTERM);
+      }
+
+      /* compact the nrowind/nrowval */
+      for (nnz=0, i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<nrowptr[i]; j++, nnz++) {
+          nrowind[nnz] = nrowind[j];
+          nrowval[nnz] = nrowval[j];
+        }
+        nrowptr[i] = nnz;
+      }
+      SHIFTCSR(i, nrows, nrowptr);
+
+      break;
+
+    case GK_CSR_ROW:
+      if (mat->rowptr == NULL) 
+        gk_errexit(SIGERR, "Cannot filter rows when row-based structure has not been created.\n");
+
+      for (i=0; i<nrows; i++) 
+        maxlen = gk_max(maxlen, rowptr[i+1]-rowptr[i]);
+
+      #pragma omp parallel private(i, j, ncand, rsum, tsum, cand)
+      {
+        cand = gk_fkvmalloc(maxlen, "gk_csr_LowFilter: cand");
+
+        #pragma omp for schedule(static)
+        for (i=0; i<nrows; i++) {
+          for (tsum=0.0, ncand=0, j=rowptr[i]; j<rowptr[i+1]; j++, ncand++) {
+            cand[ncand].val = rowind[j];
+            cand[ncand].key = rowval[j];
+            tsum += (norm == 1 ? rowval[j] : rowval[j]*rowval[j]);
+          }
+          gk_fkvsortd(ncand, cand);
+
+          for (rsum=0.0, j=0; j<ncand && rsum<=fraction*tsum; j++) {
+            rsum += (norm == 1 ? cand[j].key : cand[j].key*cand[j].key);
+            nrowind[rowptr[i]+j] = cand[j].val;
+            nrowval[rowptr[i]+j] = cand[j].key;
+          }
+          nrowptr[i+1] = rowptr[i]+j;
+        }
+
+        gk_free((void **)&cand, LTERM);
+      }
+
+      /* compact nrowind/nrowval */
+      nrowptr[0] = nnz = 0;
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<nrowptr[i+1]; j++, nnz++) {
+          nrowind[nnz] = nrowind[j];
+          nrowval[nnz] = nrowval[j];
+        }
+        nrowptr[i+1] = nnz;
+      }
+
+      break;
+
+    default:
+      gk_csr_Free(&nmat);
+      gk_errexit(SIGERR, "Unknown prunning type of %d\n", what);
+      return NULL;
+  }
+
+  return nmat;
+}
+
+
+/*************************************************************************/
+/*! Eliminates certain entries from the rows/columns of the matrix. The 
+    filtering takes place by keeping only the highest weight top-K entries 
+    along each row/column and those entries whose weight is greater than
+    a specified value.
+   
+    \param mat the matrix to be prunned,
+    \param what indicates if the rows (GK_CSR_ROW) or the columns (GK_CSR_COL)
+           of the matrix will be prunned,
+    \param topk is the number of the highest weight entries to keep.
+    \param keepval is the weight of a term above which will be kept. This
+           is used to select additional terms past the first topk.
+    \returns the filtered matrix consisting only of its row-based structure. 
+           The input matrix is not modified. 
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_TopKPlusFilter(gk_csr_t *mat, int what, int topk, float keepval)
+{
+  ssize_t i, j, k, nnz;
+  int nrows, ncols, ncand;
+  ssize_t *rowptr, *colptr, *nrowptr;
+  int *rowind, *colind, *nrowind;
+  float *rowval, *colval, *nrowval;
+  gk_csr_t *nmat;
+  gk_fkv_t *cand;
+
+  nmat = gk_csr_Create();
+  
+  nrows = nmat->nrows = mat->nrows;
+  ncols = nmat->ncols = mat->ncols;
+
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+  rowval = mat->rowval;
+  colptr = mat->colptr;
+  colind = mat->colind;
+  colval = mat->colval;
+
+  nrowptr = nmat->rowptr = gk_zmalloc(nrows+1, "gk_csr_LowFilter: nrowptr");
+  nrowind = nmat->rowind = gk_imalloc(rowptr[nrows], "gk_csr_LowFilter: nrowind");
+  nrowval = nmat->rowval = gk_fmalloc(rowptr[nrows], "gk_csr_LowFilter: nrowval");
+
+
+  switch (what) {
+    case GK_CSR_COL:
+      if (mat->colptr == NULL) 
+        gk_errexit(SIGERR, "Cannot filter columns when column-based structure has not been created.\n");
+
+      cand = gk_fkvmalloc(nrows, "gk_csr_LowFilter: cand");
+
+      gk_zcopy(nrows+1, rowptr, nrowptr);
+      for (i=0; i<ncols; i++) {
+        for (ncand=0, j=colptr[i]; j<colptr[i+1]; j++, ncand++) {
+          cand[ncand].val = colind[j];
+          cand[ncand].key = colval[j];
+        }
+        gk_fkvsortd(ncand, cand);
+
+        k = gk_min(topk, ncand);
+        for (j=0; j<k; j++) {
+          nrowind[nrowptr[cand[j].val]] = i;
+          nrowval[nrowptr[cand[j].val]] = cand[j].key;
+          nrowptr[cand[j].val]++;
+        }
+        for (; j<ncand; j++) {
+          if (cand[j].key < keepval) 
+            break;
+
+          nrowind[nrowptr[cand[j].val]] = i;
+          nrowval[nrowptr[cand[j].val]] = cand[j].key;
+          nrowptr[cand[j].val]++;
+        }
+      }
+
+      /* compact the nrowind/nrowval */
+      for (nnz=0, i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<nrowptr[i]; j++, nnz++) {
+          nrowind[nnz] = nrowind[j];
+          nrowval[nnz] = nrowval[j];
+        }
+        nrowptr[i] = nnz;
+      }
+      SHIFTCSR(i, nrows, nrowptr);
+
+      gk_free((void **)&cand, LTERM);
+      break;
+
+    case GK_CSR_ROW:
+      if (mat->rowptr == NULL) 
+        gk_errexit(SIGERR, "Cannot filter rows when row-based structure has not been created.\n");
+
+      cand = gk_fkvmalloc(ncols, "gk_csr_LowFilter: cand");
+
+      nrowptr[0] = 0;
+      for (nnz=0, i=0; i<nrows; i++) {
+        for (ncand=0, j=rowptr[i]; j<rowptr[i+1]; j++, ncand++) {
+          cand[ncand].val = rowind[j];
+          cand[ncand].key = rowval[j];
+        }
+        gk_fkvsortd(ncand, cand);
+
+        k = gk_min(topk, ncand);
+        for (j=0; j<k; j++, nnz++) {
+          nrowind[nnz] = cand[j].val;
+          nrowval[nnz] = cand[j].key;
+        }
+        for (; j<ncand; j++, nnz++) {
+          if (cand[j].key < keepval) 
+            break;
+
+          nrowind[nnz] = cand[j].val;
+          nrowval[nnz] = cand[j].key;
+        }
+        nrowptr[i+1] = nnz;
+      }
+
+      gk_free((void **)&cand, LTERM);
+      break;
+
+    default:
+      gk_csr_Free(&nmat);
+      gk_errexit(SIGERR, "Unknown prunning type of %d\n", what);
+      return NULL;
+  }
+
+  return nmat;
+}
+
+
+/*************************************************************************/
+/*! Eliminates certain entries from the rows/columns of the matrix. The 
+    filtering takes place by keeping only the terms whose contribution to
+    the total length of the document is greater than a user-splied multiple
+    over the average.
+
+    This routine assumes that the vectors are normalized to be unit length.
+   
+    \param mat the matrix to be prunned,
+    \param what indicates if the rows (GK_CSR_ROW) or the columns (GK_CSR_COL)
+           of the matrix will be prunned,
+    \param zscore is the multiplicative factor over the average contribution 
+           to the length of the document.
+    \returns the filtered matrix consisting only of its row-based structure. 
+           The input matrix is not modified. 
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_ZScoreFilter(gk_csr_t *mat, int what, float zscore)
+{
+  ssize_t i, j, nnz;
+  int nrows;
+  ssize_t *rowptr, *nrowptr;
+  int *rowind, *nrowind;
+  float *rowval, *nrowval, avgwgt;
+  gk_csr_t *nmat;
+
+  nmat = gk_csr_Create();
+  
+  nmat->nrows = mat->nrows;
+  nmat->ncols = mat->ncols;
+
+  nrows  = mat->nrows; 
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+  rowval = mat->rowval;
+
+  nrowptr = nmat->rowptr = gk_zmalloc(nrows+1, "gk_csr_ZScoreFilter: nrowptr");
+  nrowind = nmat->rowind = gk_imalloc(rowptr[nrows], "gk_csr_ZScoreFilter: nrowind");
+  nrowval = nmat->rowval = gk_fmalloc(rowptr[nrows], "gk_csr_ZScoreFilter: nrowval");
+
+
+  switch (what) {
+    case GK_CSR_COL:
+      gk_errexit(SIGERR, "This has not been implemented yet.\n");
+      break;
+
+    case GK_CSR_ROW:
+      if (mat->rowptr == NULL) 
+        gk_errexit(SIGERR, "Cannot filter rows when row-based structure has not been created.\n");
+
+      nrowptr[0] = 0;
+      for (nnz=0, i=0; i<nrows; i++) {
+        avgwgt = zscore/(rowptr[i+1]-rowptr[i]);
+        for (j=rowptr[i]; j<rowptr[i+1]; j++) {
+          if (rowval[j] > avgwgt) {
+            nrowind[nnz] = rowind[j];
+            nrowval[nnz] = rowval[j];
+            nnz++;
+          }
+        }
+        nrowptr[i+1] = nnz;
+      }
+      break;
+
+    default:
+      gk_csr_Free(&nmat);
+      gk_errexit(SIGERR, "Unknown prunning type of %d\n", what);
+      return NULL;
+  }
+
+  return nmat;
+}
+
+
+/*************************************************************************/
+/*! Compacts the column-space of the matrix by removing empty columns.
+    As a result of the compaction, the column numbers are renumbered. 
+    The compaction operation is done in place and only affects the row-based
+    representation of the matrix.
+    The new columns are ordered in decreasing frequency.
+   
+    \param mat the matrix whose empty columns will be removed.
+*/
+/**************************************************************************/
+void gk_csr_CompactColumns(gk_csr_t *mat)
+{
+  ssize_t i;
+  int nrows, ncols, nncols;
+  ssize_t *rowptr;
+  int *rowind, *colmap;
+  gk_ikv_t *clens;
+
+  nrows  = mat->nrows;
+  ncols  = mat->ncols;
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+
+  colmap = gk_imalloc(ncols, "gk_csr_CompactColumns: colmap");
+
+  clens = gk_ikvmalloc(ncols, "gk_csr_CompactColumns: clens");
+  for (i=0; i<ncols; i++) {
+    clens[i].key = 0;
+    clens[i].val = i;
+  }
+
+  for (i=0; i<rowptr[nrows]; i++) 
+    clens[rowind[i]].key++;
+  gk_ikvsortd(ncols, clens);
+
+  for (nncols=0, i=0; i<ncols; i++) {
+    if (clens[i].key > 0) 
+      colmap[clens[i].val] = nncols++;
+    else
+      break;
+  }
+
+  for (i=0; i<rowptr[nrows]; i++) 
+    rowind[i] = colmap[rowind[i]];
+
+  mat->ncols = nncols;
+
+  gk_free((void **)&colmap, &clens, LTERM);
+}
+
+
+/*************************************************************************/
+/*! Sorts the indices in increasing order
+    \param mat the matrix itself,
+    \param what is either GK_CSR_ROW or GK_CSR_COL indicating which set of
+           indices to sort.
+*/
+/**************************************************************************/
+void gk_csr_SortIndices(gk_csr_t *mat, int what)
+{
+  int n, nn=0;
+  ssize_t *ptr;
+  int *ind;
+  float *val;
+
+  switch (what) {
+    case GK_CSR_ROW:
+      if (!mat->rowptr)
+        gk_errexit(SIGERR, "Row-based view of the matrix does not exists.\n");
+
+      n   = mat->nrows;
+      ptr = mat->rowptr;
+      ind = mat->rowind;
+      val = mat->rowval;
+      break;
+
+    case GK_CSR_COL:
+      if (!mat->colptr)
+        gk_errexit(SIGERR, "Column-based view of the matrix does not exists.\n");
+
+      n   = mat->ncols;
+      ptr = mat->colptr;
+      ind = mat->colind;
+      val = mat->colval;
+      break;
+
+    default:
+      gk_errexit(SIGERR, "Invalid index type of %d.\n", what);
+      return;
+  }
+
+  #pragma omp parallel if (n > 100)
+  {
+    ssize_t i, j, k;
+    gk_ikv_t *cand;
+    float *tval;
+
+    #pragma omp single
+    for (i=0; i<n; i++) 
+      nn = gk_max(nn, ptr[i+1]-ptr[i]);
+  
+    cand = gk_ikvmalloc(nn, "gk_csr_SortIndices: cand");
+    tval = gk_fmalloc(nn, "gk_csr_SortIndices: tval");
+  
+    #pragma omp for schedule(static)
+    for (i=0; i<n; i++) {
+      for (k=0, j=ptr[i]; j<ptr[i+1]; j++) {
+        if (j > ptr[i] && ind[j] < ind[j-1])
+          k = 1; /* an inversion */
+        cand[j-ptr[i]].val = j-ptr[i];
+        cand[j-ptr[i]].key = ind[j];
+        tval[j-ptr[i]]     = val[j];
+      }
+      if (k) {
+        gk_ikvsorti(ptr[i+1]-ptr[i], cand);
+        for (j=ptr[i]; j<ptr[i+1]; j++) {
+          ind[j] = cand[j-ptr[i]].key;
+          val[j] = tval[cand[j-ptr[i]].val];
+        }
+      }
+    }
+
+    gk_free((void **)&cand, &tval, LTERM);
+  }
+
+}
+
+
+/*************************************************************************/
+/*! Creates a row/column index from the column/row data.
+    \param mat the matrix itself,
+    \param what is either GK_CSR_ROW or GK_CSR_COL indicating which index
+           will be created.
+*/
+/**************************************************************************/
+void gk_csr_CreateIndex(gk_csr_t *mat, int what)
+{
+  /* 'f' stands for forward, 'r' stands for reverse */
+  ssize_t i, j, k, nf, nr;
+  ssize_t *fptr, *rptr;
+  int *find, *rind;
+  float *fval, *rval;
+
+  switch (what) {
+    case GK_CSR_COL:
+      nf   = mat->nrows;
+      fptr = mat->rowptr;
+      find = mat->rowind;
+      fval = mat->rowval;
+
+      if (mat->colptr) gk_free((void **)&mat->colptr, LTERM);
+      if (mat->colind) gk_free((void **)&mat->colind, LTERM);
+      if (mat->colval) gk_free((void **)&mat->colval, LTERM);
+
+      nr   = mat->ncols;
+      rptr = mat->colptr = gk_zsmalloc(nr+1, 0, "gk_csr_CreateIndex: rptr");
+      rind = mat->colind = gk_imalloc(fptr[nf], "gk_csr_CreateIndex: rind");
+      rval = mat->colval = (fval ? gk_fmalloc(fptr[nf], "gk_csr_CreateIndex: rval") : NULL);
+      break;
+    case GK_CSR_ROW:
+      nf   = mat->ncols;
+      fptr = mat->colptr;
+      find = mat->colind;
+      fval = mat->colval;
+
+      if (mat->rowptr) gk_free((void **)&mat->rowptr, LTERM);
+      if (mat->rowind) gk_free((void **)&mat->rowind, LTERM);
+      if (mat->rowval) gk_free((void **)&mat->rowval, LTERM);
+
+      nr   = mat->nrows;
+      rptr = mat->rowptr = gk_zsmalloc(nr+1, 0, "gk_csr_CreateIndex: rptr");
+      rind = mat->rowind = gk_imalloc(fptr[nf], "gk_csr_CreateIndex: rind");
+      rval = mat->rowval = (fval ? gk_fmalloc(fptr[nf], "gk_csr_CreateIndex: rval") : NULL);
+      break;
+    default:
+      gk_errexit(SIGERR, "Invalid index type of %d.\n", what);
+      return;
+  }
+
+
+  for (i=0; i<nf; i++) {
+    for (j=fptr[i]; j<fptr[i+1]; j++)
+      rptr[find[j]]++;
+  }
+  MAKECSR(i, nr, rptr);
+  
+  if (rptr[nr] > 6*nr) {
+    for (i=0; i<nf; i++) {
+      for (j=fptr[i]; j<fptr[i+1]; j++) 
+        rind[rptr[find[j]]++] = i;
+    }
+    SHIFTCSR(i, nr, rptr);
+
+    if (fval) {
+      for (i=0; i<nf; i++) {
+        for (j=fptr[i]; j<fptr[i+1]; j++) 
+          rval[rptr[find[j]]++] = fval[j];
+      }
+      SHIFTCSR(i, nr, rptr);
+    }
+  }
+  else {
+    if (fval) {
+      for (i=0; i<nf; i++) {
+        for (j=fptr[i]; j<fptr[i+1]; j++) {
+          k = find[j];
+          rind[rptr[k]]   = i;
+          rval[rptr[k]++] = fval[j];
+        }
+      }
+    }
+    else {
+      for (i=0; i<nf; i++) {
+        for (j=fptr[i]; j<fptr[i+1]; j++) 
+          rind[rptr[find[j]]++] = i;
+      }
+    }
+    SHIFTCSR(i, nr, rptr);
+  }
+}
+
+
+/*************************************************************************/
+/*! Normalizes the rows/columns of the matrix to be unit 
+    length.
+    \param mat the matrix itself,
+    \param what indicates what will be normalized and is obtained by
+           specifying GK_CSR_ROW, GK_CSR_COL, GK_CSR_ROW|GK_CSR_COL. 
+    \param norm indicates what norm is to normalize to, 1: 1-norm, 2: 2-norm
+*/
+/**************************************************************************/
+void gk_csr_Normalize(gk_csr_t *mat, int what, int norm)
+{
+  ssize_t i, j;
+  int n;
+  ssize_t *ptr;
+  float *val, sum;
+
+
+  if (what&GK_CSR_ROW && mat->rowval) {
+    n   = mat->nrows;
+    ptr = mat->rowptr;
+    val = mat->rowval;
+
+    #pragma omp parallel for if (ptr[n] > OMPMINOPS) private(j,sum) schedule(static)
+    for (i=0; i<n; i++) {
+      sum = 0.0;
+      if (norm == 1) {
+        for (j=ptr[i]; j<ptr[i+1]; j++) 
+          sum += val[j]; /* assume val[j] > 0 */ 
+        if (sum > 0)
+          sum = 1.0/sum;
+      }
+      else if (norm == 2) {
+        for (j=ptr[i]; j<ptr[i+1]; j++) 
+          sum += val[j]*val[j];
+        if (sum > 0)
+          sum = 1.0/sqrt(sum); 
+      }
+      for (j=ptr[i]; j<ptr[i+1]; j++)
+        val[j] *= sum;
+    }
+  }
+
+  if (what&GK_CSR_COL && mat->colval) {
+    n   = mat->ncols;
+    ptr = mat->colptr;
+    val = mat->colval;
+
+    #pragma omp parallel for if (ptr[n] > OMPMINOPS) private(j,sum) schedule(static)
+    for (i=0; i<n; i++) {
+      sum = 0.0;
+      if (norm == 1) {
+        for (j=ptr[i]; j<ptr[i+1]; j++) 
+          sum += val[j]; /* assume val[j] > 0 */ 
+        if (sum > 0)
+          sum = 1.0/sum;
+      }
+      else if (norm == 2) {
+        for (j=ptr[i]; j<ptr[i+1]; j++) 
+          sum += val[j]*val[j];
+        if (sum > 0)
+          sum = 1.0/sqrt(sum); 
+      }
+      for (j=ptr[i]; j<ptr[i+1]; j++)
+        val[j] *= sum;
+    }
+  }
+
+}
+
+
+/*************************************************************************/
+/*! Applies different row scaling methods.
+    \param mat the matrix itself,
+    \param type indicates the type of row scaling. Possible values are:
+           GK_CSR_MAXTF, GK_CSR_SQRT, GK_CSR_LOG, GK_CSR_IDF, GK_CSR_MAXTF2.
+*/
+/**************************************************************************/
+void gk_csr_Scale(gk_csr_t *mat, int type)
+{
+  ssize_t i, j;
+  int nrows, ncols, nnzcols, bgfreq;
+  ssize_t *rowptr;
+  int *rowind, *collen;
+  float *rowval, *cscale, maxtf;
+  double logscale = 1.0/log(2.0);
+
+  nrows  = mat->nrows;
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+  rowval = mat->rowval;
+
+  switch (type) {
+    case GK_CSR_MAXTF: /* TF' = .5 + .5*TF/MAX(TF) */
+      #pragma omp parallel for if (rowptr[nrows] > OMPMINOPS) private(j, maxtf) schedule(static)
+      for (i=0; i<nrows; i++) {
+        maxtf = fabs(rowval[rowptr[i]]);
+        for (j=rowptr[i]; j<rowptr[i+1]; j++) 
+          maxtf = (maxtf < fabs(rowval[j]) ? fabs(rowval[j]) : maxtf);
+  
+        for (j=rowptr[i]; j<rowptr[i+1]; j++)
+          rowval[j] = .5 + .5*rowval[j]/maxtf;
+      }
+      break;
+
+    case GK_CSR_MAXTF2: /* TF' = .1 + .9*TF/MAX(TF) */
+      #pragma omp parallel for if (rowptr[nrows] > OMPMINOPS) private(j, maxtf) schedule(static)
+      for (i=0; i<nrows; i++) {
+        maxtf = fabs(rowval[rowptr[i]]);
+        for (j=rowptr[i]; j<rowptr[i+1]; j++) 
+          maxtf = (maxtf < fabs(rowval[j]) ? fabs(rowval[j]) : maxtf);
+  
+        for (j=rowptr[i]; j<rowptr[i+1]; j++)
+          rowval[j] = .1 + .9*rowval[j]/maxtf;
+      }
+      break;
+
+    case GK_CSR_SQRT: /* TF' = .1+SQRT(TF) */
+      #pragma omp parallel for if (rowptr[nrows] > OMPMINOPS) private(j) schedule(static)
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++) { 
+          if (rowval[j] != 0.0)
+            rowval[j] = .1+sign(rowval[j], sqrt(fabs(rowval[j])));
+        }
+      }
+      
+      break;
+
+    case GK_CSR_POW25: /* TF' = .1+POW(TF,.25) */
+      #pragma omp parallel for if (rowptr[nrows] > OMPMINOPS) private(j) schedule(static)
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++) { 
+          if (rowval[j] != 0.0)
+            rowval[j] = .1+sign(rowval[j], sqrt(sqrt(fabs(rowval[j]))));
+        }
+      }
+      break;
+
+    case GK_CSR_POW65: /* TF' = .1+POW(TF,.65) */
+      #pragma omp parallel for if (rowptr[nrows] > OMPMINOPS) private(j) schedule(static)
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++) { 
+          if (rowval[j] != 0.0)
+            rowval[j] = .1+sign(rowval[j], powf(fabs(rowval[j]), .65));
+        }
+      }
+      break;
+
+    case GK_CSR_POW75: /* TF' = .1+POW(TF,.75) */
+      #pragma omp parallel for if (rowptr[nrows] > OMPMINOPS) private(j) schedule(static)
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++) { 
+          if (rowval[j] != 0.0)
+            rowval[j] = .1+sign(rowval[j], powf(fabs(rowval[j]), .75));
+        }
+      }
+      break;
+
+    case GK_CSR_POW85: /* TF' = .1+POW(TF,.85) */
+      #pragma omp parallel for if (rowptr[nrows] > OMPMINOPS) private(j) schedule(static)
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++) { 
+          if (rowval[j] != 0.0)
+            rowval[j] = .1+sign(rowval[j], powf(fabs(rowval[j]), .85));
+        }
+      }
+      break;
+
+    case GK_CSR_LOG: /* TF' = 1+log_2(TF) */
+      #pragma omp parallel for if (rowptr[nrows] > OMPMINOPS) schedule(static,32)
+      for (i=0; i<rowptr[nrows]; i++) {
+        if (rowval[i] != 0.0)
+          rowval[i] = 1+(rowval[i]>0.0 ? log(rowval[i]) : -log(-rowval[i]))*logscale;
+      }
+#ifdef XXX
+      #pragma omp parallel for private(j) schedule(static)
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++) { 
+          if (rowval[j] != 0.0)
+            rowval[j] = 1+(rowval[j]>0.0 ? log(rowval[j]) : -log(-rowval[j]))*logscale;
+            //rowval[j] = 1+sign(rowval[j], log(fabs(rowval[j]))*logscale);
+        }
+      }
+#endif
+      break;
+
+    case GK_CSR_IDF: /* TF' = TF*IDF */
+      ncols  = mat->ncols;
+      cscale = gk_fmalloc(ncols, "gk_csr_Scale: cscale");
+      collen = gk_ismalloc(ncols, 0, "gk_csr_Scale: collen");
+
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++)
+          collen[rowind[j]]++;
+      }
+
+      #pragma omp parallel for if (ncols > OMPMINOPS) schedule(static)
+      for (i=0; i<ncols; i++)
+        cscale[i] = (collen[i] > 0 ? log(1.0*nrows/collen[i]) : 0.0);
+
+      #pragma omp parallel for if (rowptr[nrows] > OMPMINOPS) private(j) schedule(static)
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++)
+          rowval[j] *= cscale[rowind[j]];
+      }
+      
+      gk_free((void **)&cscale, &collen, LTERM);
+      break;
+
+    case GK_CSR_IDF2: /* TF' = TF*IDF */
+      ncols  = mat->ncols;
+      cscale = gk_fmalloc(ncols, "gk_csr_Scale: cscale");
+      collen = gk_ismalloc(ncols, 0, "gk_csr_Scale: collen");
+
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++)
+          collen[rowind[j]]++;
+      }
+
+      nnzcols = 0;
+      #pragma omp parallel for if (ncols > OMPMINOPS) schedule(static) reduction(+:nnzcols)
+      for (i=0; i<ncols; i++)
+        nnzcols += (collen[i] > 0 ? 1 : 0);
+
+      bgfreq = gk_max(10, (ssize_t)(.5*rowptr[nrows]/nnzcols));
+      printf("nnz: %zd, nnzcols: %d, bgfreq: %d\n", rowptr[nrows], nnzcols, bgfreq);
+
+      #pragma omp parallel for if (ncols > OMPMINOPS) schedule(static)
+      for (i=0; i<ncols; i++)
+        cscale[i] = (collen[i] > 0 ? log(1.0*(nrows+2*bgfreq)/(bgfreq+collen[i])) : 0.0);
+
+      #pragma omp parallel for if (rowptr[nrows] > OMPMINOPS) private(j) schedule(static)
+      for (i=0; i<nrows; i++) {
+        for (j=rowptr[i]; j<rowptr[i+1]; j++)
+          rowval[j] *= cscale[rowind[j]];
+      }
+
+      gk_free((void **)&cscale, &collen, LTERM);
+      break;
+
+    default:
+      gk_errexit(SIGERR, "Unknown scaling type of %d\n", type);
+  }
+}
+
+
+/*************************************************************************/
+/*! Computes the sums of the rows/columns
+    \param mat the matrix itself,
+    \param what is either GK_CSR_ROW or GK_CSR_COL indicating which 
+           sums to compute.
+*/
+/**************************************************************************/
+void gk_csr_ComputeSums(gk_csr_t *mat, int what)
+{
+  ssize_t i;
+  int n;
+  ssize_t *ptr;
+  float *val, *sums;
+
+  switch (what) {
+    case GK_CSR_ROW:
+      n   = mat->nrows;
+      ptr = mat->rowptr;
+      val = mat->rowval;
+
+      if (mat->rsums) 
+        gk_free((void **)&mat->rsums, LTERM);
+
+      sums = mat->rsums = gk_fsmalloc(n, 0, "gk_csr_ComputeSums: sums");
+      break;
+    case GK_CSR_COL:
+      n   = mat->ncols;
+      ptr = mat->colptr;
+      val = mat->colval;
+
+      if (mat->csums) 
+        gk_free((void **)&mat->csums, LTERM);
+
+      sums = mat->csums = gk_fsmalloc(n, 0, "gk_csr_ComputeSums: sums");
+      break;
+    default:
+      gk_errexit(SIGERR, "Invalid sum type of %d.\n", what);
+      return;
+  }
+
+  if (val) {
+    #pragma omp parallel for if (ptr[n] > OMPMINOPS) schedule(static)
+    for (i=0; i<n; i++) 
+      sums[i] = gk_fsum(ptr[i+1]-ptr[i], val+ptr[i], 1);
+  }
+  else {
+    #pragma omp parallel for if (ptr[n] > OMPMINOPS) schedule(static)
+    for (i=0; i<n; i++) 
+      sums[i] = ptr[i+1]-ptr[i];
+  }
+}
+
+
+/*************************************************************************/
+/*! Computes the norms of the rows/columns
+
+    \param mat the matrix itself,
+    \param what is either GK_CSR_ROW or GK_CSR_COL indicating which 
+           squared norms to compute.
+
+    \note If the rowval/colval arrays are NULL, the matrix is assumed
+          to be binary and the norms are computed accordingly.
+*/
+/**************************************************************************/
+void gk_csr_ComputeNorms(gk_csr_t *mat, int what)
+{
+  ssize_t i;
+  int n;
+  ssize_t *ptr;
+  float *val, *norms;
+
+  switch (what) {
+    case GK_CSR_ROW:
+      n   = mat->nrows;
+      ptr = mat->rowptr;
+      val = mat->rowval;
+
+      if (mat->rnorms) gk_free((void **)&mat->rnorms, LTERM);
+
+      norms = mat->rnorms = gk_fsmalloc(n, 0, "gk_csr_ComputeSums: norms");
+      break;
+    case GK_CSR_COL:
+      n   = mat->ncols;
+      ptr = mat->colptr;
+      val = mat->colval;
+
+      if (mat->cnorms) gk_free((void **)&mat->cnorms, LTERM);
+
+      norms = mat->cnorms = gk_fsmalloc(n, 0, "gk_csr_ComputeSums: norms");
+      break;
+    default:
+      gk_errexit(SIGERR, "Invalid norm type of %d.\n", what);
+      return;
+  }
+
+  if (val) {
+    #pragma omp parallel for if (ptr[n] > OMPMINOPS) schedule(static)
+    for (i=0; i<n; i++) 
+      norms[i] = sqrt(gk_fdot(ptr[i+1]-ptr[i], val+ptr[i], 1, val+ptr[i], 1));
+  }
+  else {
+    #pragma omp parallel for if (ptr[n] > OMPMINOPS) schedule(static)
+    for (i=0; i<n; i++) 
+      norms[i] = sqrt(ptr[i+1]-ptr[i]);
+  }
+}
+
+
+/*************************************************************************/
+/*! Computes the squared of the norms of the rows/columns
+
+    \param mat the matrix itself,
+    \param what is either GK_CSR_ROW or GK_CSR_COL indicating which 
+           squared norms to compute.
+
+    \note If the rowval/colval arrays are NULL, the matrix is assumed
+          to be binary and the norms are computed accordingly.
+*/
+/**************************************************************************/
+void gk_csr_ComputeSquaredNorms(gk_csr_t *mat, int what)
+{
+  ssize_t i;
+  int n;
+  ssize_t *ptr;
+  float *val, *norms;
+
+  switch (what) {
+    case GK_CSR_ROW:
+      n   = mat->nrows;
+      ptr = mat->rowptr;
+      val = mat->rowval;
+
+      if (mat->rnorms) gk_free((void **)&mat->rnorms, LTERM);
+
+      norms = mat->rnorms = gk_fsmalloc(n, 0, "gk_csr_ComputeSums: norms");
+      break;
+    case GK_CSR_COL:
+      n   = mat->ncols;
+      ptr = mat->colptr;
+      val = mat->colval;
+
+      if (mat->cnorms) gk_free((void **)&mat->cnorms, LTERM);
+
+      norms = mat->cnorms = gk_fsmalloc(n, 0, "gk_csr_ComputeSums: norms");
+      break;
+    default:
+      gk_errexit(SIGERR, "Invalid norm type of %d.\n", what);
+      return;
+  }
+
+  if (val) {
+    #pragma omp parallel for if (ptr[n] > OMPMINOPS) schedule(static)
+    for (i=0; i<n; i++) 
+      norms[i] = gk_fdot(ptr[i+1]-ptr[i], val+ptr[i], 1, val+ptr[i], 1);
+  }
+  else {
+    #pragma omp parallel for if (ptr[n] > OMPMINOPS) schedule(static)
+    for (i=0; i<n; i++) 
+      norms[i] = ptr[i+1]-ptr[i];
+  }
+}
+
+
+/*************************************************************************/
+/*! Returns a new matrix whose rows/columns are shuffled.
+   
+    \param mat the matrix to be shuffled,
+    \param what indicates if the rows (GK_CSR_ROW), columns (GK_CSR_COL),
+           or both (GK_CSR_ROWCOL) will be shuffled,
+    \param symmetric indicates if the same shuffling will be applied to 
+           both rows and columns. This is valid with nrows==ncols and 
+           GK_CSR_ROWCOL was specified.
+    \returns the shuffled matrix. 
+          The input matrix is not modified. 
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_Shuffle(gk_csr_t *mat, int what, int symmetric)
+{
+  ssize_t i, j;
+  int nrows, ncols;
+  ssize_t *rowptr, *nrowptr;
+  int *rowind, *nrowind;
+  int *rperm, *cperm;
+  float *rowval, *nrowval;
+  gk_csr_t *nmat;
+
+  if (what == GK_CSR_ROWCOL && symmetric && mat->nrows != mat->ncols)
+    gk_errexit(SIGERR, "The matrix is not square for a symmetric rowcol shuffling.\n");
+
+  nrows  = mat->nrows;
+  ncols  = mat->ncols;
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+  rowval = mat->rowval;
+
+  rperm = gk_imalloc(nrows, "gk_csr_Shuffle: rperm");
+  cperm = gk_imalloc(ncols, "gk_csr_Shuffle: cperm");
+
+  switch (what) {
+    case GK_CSR_ROW:
+      gk_RandomPermute(nrows, rperm, 1);
+      for (i=0; i<20; i++)
+        gk_RandomPermute(nrows, rperm, 0);
+
+      for (i=0; i<ncols; i++)
+        cperm[i] = i;
+      break;
+
+    case GK_CSR_COL:
+      gk_RandomPermute(ncols, cperm, 1);
+      for (i=0; i<20; i++)
+        gk_RandomPermute(ncols, cperm, 0);
+
+      for (i=0; i<nrows; i++)
+        rperm[i] = i;
+      break;
+
+    case GK_CSR_ROWCOL:
+      gk_RandomPermute(nrows, rperm, 1);
+      for (i=0; i<20; i++)
+        gk_RandomPermute(nrows, rperm, 0);
+
+      if (symmetric)
+        gk_icopy(nrows, rperm, cperm);
+      else {
+        gk_RandomPermute(ncols, cperm, 1);
+        for (i=0; i<20; i++)
+          gk_RandomPermute(ncols, cperm, 0);
+      }
+      break;
+
+    default:
+      gk_free((void **)&rperm, &cperm, LTERM);
+      gk_errexit(SIGERR, "Unknown shuffling type of %d\n", what);
+      return NULL;
+  }
+
+  nmat = gk_csr_Create();
+  nmat->nrows = nrows;
+  nmat->ncols = ncols;
+
+  nrowptr = nmat->rowptr = gk_zmalloc(nrows+1, "gk_csr_Shuffle: nrowptr");
+  nrowind = nmat->rowind = gk_imalloc(rowptr[nrows], "gk_csr_Shuffle: nrowind");
+  nrowval = nmat->rowval = (rowval ? gk_fmalloc(rowptr[nrows], "gk_csr_Shuffle: nrowval") : NULL) ;
+
+  for (i=0; i<nrows; i++)
+    nrowptr[rperm[i]] = rowptr[i+1]-rowptr[i];
+  MAKECSR(i, nrows, nrowptr);
+
+  for (i=0; i<nrows; i++) {
+    for (j=rowptr[i]; j<rowptr[i+1]; j++) {
+      nrowind[nrowptr[rperm[i]]] = cperm[rowind[j]];
+      if (nrowval)
+        nrowval[nrowptr[rperm[i]]] = rowval[j];
+      nrowptr[rperm[i]]++;
+    }
+  }
+  SHIFTCSR(i, nrows, nrowptr);
+
+  gk_free((void **)&rperm, &cperm, LTERM);
+
+  return nmat;
+
+}
+
+
+/*************************************************************************/
+/*! Returns the transpose of the matrix.
+   
+    \param mat the matrix to be transposed,
+    \returns the transposed matrix. 
+          The input matrix is not modified. 
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_Transpose(gk_csr_t *mat)
+{
+  int nrows, ncols;
+  ssize_t *colptr;
+  int32_t *colind;
+  float *colval;
+  gk_csr_t *nmat;
+
+  colptr = mat->colptr;
+  colind = mat->colind;
+  colval = mat->colval;
+
+  mat->colptr = NULL;
+  mat->colind = NULL;
+  mat->colval = NULL;
+
+  gk_csr_CreateIndex(mat, GK_CSR_COL);
+
+  nmat = gk_csr_Create();
+  nmat->nrows  = mat->ncols;
+  nmat->ncols  = mat->nrows;
+  nmat->rowptr = mat->colptr;
+  nmat->rowind = mat->colind;
+  nmat->rowval = mat->colval;
+
+  mat->colptr = colptr;
+  mat->colind = colind;
+  mat->colval = colval;
+
+  return nmat;
+
+}
+
+
+/*************************************************************************/
+/*! Computes the similarity between two rows/columns
+
+    \param mat the matrix itself. The routine assumes that the indices
+           are sorted in increasing order.
+    \param i1 is the first row/column,
+    \param i2 is the second row/column,
+    \param what is either GK_CSR_ROW or GK_CSR_COL indicating the type of
+           objects between the similarity will be computed,
+    \param simtype is the type of similarity and is one of GK_CSR_COS,
+           GK_CSR_JAC, GK_CSR_MIN, GK_CSR_AMIN
+    \returns the similarity between the two rows/columns.
+*/
+/**************************************************************************/
+float gk_csr_ComputeSimilarity(gk_csr_t *mat, int i1, int i2, int what, 
+          int simtype)
+{
+  int nind1, nind2;
+  int *ind1, *ind2;
+  float *val1, *val2, stat1, stat2, sim;
+
+  switch (what) {
+    case GK_CSR_ROW:
+      if (!mat->rowptr)
+        gk_errexit(SIGERR, "Row-based view of the matrix does not exists.\n");
+      nind1 = mat->rowptr[i1+1]-mat->rowptr[i1];
+      nind2 = mat->rowptr[i2+1]-mat->rowptr[i2];
+      ind1  = mat->rowind + mat->rowptr[i1];
+      ind2  = mat->rowind + mat->rowptr[i2];
+      val1  = mat->rowval + mat->rowptr[i1];
+      val2  = mat->rowval + mat->rowptr[i2];
+      break;
+
+    case GK_CSR_COL:
+      if (!mat->colptr)
+        gk_errexit(SIGERR, "Column-based view of the matrix does not exists.\n");
+      nind1 = mat->colptr[i1+1]-mat->colptr[i1];
+      nind2 = mat->colptr[i2+1]-mat->colptr[i2];
+      ind1  = mat->colind + mat->colptr[i1];
+      ind2  = mat->colind + mat->colptr[i2];
+      val1  = mat->colval + mat->colptr[i1];
+      val2  = mat->colval + mat->colptr[i2];
+      break;
+
+    default:
+      gk_errexit(SIGERR, "Invalid index type of %d.\n", what);
+      return 0.0;
+  }
+
+
+  switch (simtype) {
+    case GK_CSR_COS:
+    case GK_CSR_JAC:
+      sim = stat1 = stat2 = 0.0;
+      i1 = i2 = 0;
+      while (i1<nind1 && i2<nind2) {
+        if (i1 == nind1) {
+          stat2 += val2[i2]*val2[i2];
+          i2++;
+        }
+        else if (i2 == nind2) {
+          stat1 += val1[i1]*val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] < ind2[i2]) {
+          stat1 += val1[i1]*val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] > ind2[i2]) {
+          stat2 += val2[i2]*val2[i2];
+          i2++;
+        }
+        else {
+          sim   += val1[i1]*val2[i2];
+          stat1 += val1[i1]*val1[i1];
+          stat2 += val2[i2]*val2[i2];
+          i1++;
+          i2++;
+        }
+      }
+      if (simtype == GK_CSR_COS)
+        sim = (stat1*stat2 > 0.0 ? sim/sqrt(stat1*stat2) : 0.0);
+      else 
+        sim = (stat1+stat2-sim > 0.0 ? sim/(stat1+stat2-sim) : 0.0);
+      break;
+
+    case GK_CSR_MIN:
+      sim = stat1 = stat2 = 0.0;
+      i1 = i2 = 0;
+      while (i1<nind1 && i2<nind2) {
+        if (i1 == nind1) {
+          stat2 += val2[i2];
+          i2++;
+        }
+        else if (i2 == nind2) {
+          stat1 += val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] < ind2[i2]) {
+          stat1 += val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] > ind2[i2]) {
+          stat2 += val2[i2];
+          i2++;
+        }
+        else {
+          sim   += gk_min(val1[i1],val2[i2]);
+          stat1 += val1[i1];
+          stat2 += val2[i2];
+          i1++;
+          i2++;
+        }
+      }
+      sim = (stat1+stat2-sim > 0.0 ? sim/(stat1+stat2-sim) : 0.0);
+
+      break;
+
+    case GK_CSR_AMIN:
+      sim = stat1 = stat2 = 0.0;
+      i1 = i2 = 0;
+      while (i1<nind1 && i2<nind2) {
+        if (i1 == nind1) {
+          stat2 += val2[i2];
+          i2++;
+        }
+        else if (i2 == nind2) {
+          stat1 += val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] < ind2[i2]) {
+          stat1 += val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] > ind2[i2]) {
+          stat2 += val2[i2];
+          i2++;
+        }
+        else {
+          sim   += gk_min(val1[i1],val2[i2]);
+          stat1 += val1[i1];
+          stat2 += val2[i2];
+          i1++;
+          i2++;
+        }
+      }
+      sim = (stat1 > 0.0 ? sim/stat1 : 0.0);
+
+      break;
+
+    default:
+      gk_errexit(SIGERR, "Unknown similarity measure %d\n", simtype);
+      return -1;
+  }
+
+  return sim;
+
+}
+
+
+/*************************************************************************/
+/*! Computes the similarity between two rows/columns
+
+    \param mat_a the first matrix. The routine assumes that the indices
+           are sorted in increasing order.
+    \param mat_b the second matrix. The routine assumes that the indices
+           are sorted in increasing order.
+    \param i1 is the row/column from the first matrix (mat_a),
+    \param i2 is the row/column from the second matrix (mat_b),
+    \param what is either GK_CSR_ROW or GK_CSR_COL indicating the type of
+           objects between the similarity will be computed,
+    \param simtype is the type of similarity and is one of GK_CSR_COS,
+           GK_CSR_JAC, GK_CSR_MIN, GK_CSR_AMIN
+    \returns the similarity between the two rows/columns.
+*/
+/**************************************************************************/
+float gk_csr_ComputePairSimilarity(gk_csr_t *mat_a, gk_csr_t *mat_b, 
+          int i1, int i2, int what, int simtype)
+{
+  int nind1, nind2;
+  int *ind1, *ind2;
+  float *val1, *val2, stat1, stat2, sim;
+
+  switch (what) {
+    case GK_CSR_ROW:
+      if (!mat_a->rowptr || !mat_b->rowptr)
+        gk_errexit(SIGERR, "Row-based view of the matrix does not exists.\n");
+      nind1 = mat_a->rowptr[i1+1]-mat_a->rowptr[i1];
+      nind2 = mat_b->rowptr[i2+1]-mat_b->rowptr[i2];
+      ind1  = mat_a->rowind + mat_a->rowptr[i1];
+      ind2  = mat_b->rowind + mat_b->rowptr[i2];
+      val1  = mat_a->rowval + mat_a->rowptr[i1];
+      val2  = mat_b->rowval + mat_b->rowptr[i2];
+      break;
+
+    case GK_CSR_COL:
+      if (!mat_a->colptr || !mat_b->colptr)
+        gk_errexit(SIGERR, "Column-based view of the matrix does not exists.\n");
+      nind1 = mat_a->colptr[i1+1]-mat_a->colptr[i1];
+      nind2 = mat_b->colptr[i2+1]-mat_b->colptr[i2];
+      ind1  = mat_a->colind + mat_a->colptr[i1];
+      ind2  = mat_b->colind + mat_b->colptr[i2];
+      val1  = mat_a->colval + mat_a->colptr[i1];
+      val2  = mat_b->colval + mat_b->colptr[i2];
+      break;
+
+    default:
+      gk_errexit(SIGERR, "Invalid index type of %d.\n", what);
+      return 0.0;
+  }
+
+
+  switch (simtype) {
+    case GK_CSR_COS:
+    case GK_CSR_JAC:
+      sim = stat1 = stat2 = 0.0;
+      i1 = i2 = 0;
+      while (i1<nind1 && i2<nind2) {
+        if (i1 == nind1) {
+          stat2 += val2[i2]*val2[i2];
+          i2++;
+        }
+        else if (i2 == nind2) {
+          stat1 += val1[i1]*val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] < ind2[i2]) {
+          stat1 += val1[i1]*val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] > ind2[i2]) {
+          stat2 += val2[i2]*val2[i2];
+          i2++;
+        }
+        else {
+          sim   += val1[i1]*val2[i2];
+          stat1 += val1[i1]*val1[i1];
+          stat2 += val2[i2]*val2[i2];
+          i1++;
+          i2++;
+        }
+      }
+      if (simtype == GK_CSR_COS)
+        sim = (stat1*stat2 > 0.0 ? sim/sqrt(stat1*stat2) : 0.0);
+      else 
+        sim = (stat1+stat2-sim > 0.0 ? sim/(stat1+stat2-sim) : 0.0);
+      break;
+
+    case GK_CSR_MIN:
+      sim = stat1 = stat2 = 0.0;
+      i1 = i2 = 0;
+      while (i1<nind1 && i2<nind2) {
+        if (i1 == nind1) {
+          stat2 += val2[i2];
+          i2++;
+        }
+        else if (i2 == nind2) {
+          stat1 += val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] < ind2[i2]) {
+          stat1 += val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] > ind2[i2]) {
+          stat2 += val2[i2];
+          i2++;
+        }
+        else {
+          sim   += gk_min(val1[i1],val2[i2]);
+          stat1 += val1[i1];
+          stat2 += val2[i2];
+          i1++;
+          i2++;
+        }
+      }
+      sim = (stat1+stat2-sim > 0.0 ? sim/(stat1+stat2-sim) : 0.0);
+
+      break;
+
+    case GK_CSR_AMIN:
+      sim = stat1 = stat2 = 0.0;
+      i1 = i2 = 0;
+      while (i1<nind1 && i2<nind2) {
+        if (i1 == nind1) {
+          stat2 += val2[i2];
+          i2++;
+        }
+        else if (i2 == nind2) {
+          stat1 += val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] < ind2[i2]) {
+          stat1 += val1[i1];
+          i1++;
+        }
+        else if (ind1[i1] > ind2[i2]) {
+          stat2 += val2[i2];
+          i2++;
+        }
+        else {
+          sim   += gk_min(val1[i1],val2[i2]);
+          stat1 += val1[i1];
+          stat2 += val2[i2];
+          i1++;
+          i2++;
+        }
+      }
+      sim = (stat1 > 0.0 ? sim/stat1 : 0.0);
+
+      break;
+
+    default:
+      gk_errexit(SIGERR, "Unknown similarity measure %d\n", simtype);
+      return -1;
+  }
+
+  return sim;
+
+}
+
+/*************************************************************************/
+/*! Finds the n most similar rows (neighbors) to the query.
+
+    \param mat the matrix itself
+    \param nqterms is the number of columns in the query
+    \param qind is the list of query columns
+    \param qval is the list of correspodning query weights
+    \param simtype is the type of similarity and is one of GK_CSR_DOTP,
+           GK_CSR_COS, GK_CSR_JAC, GK_CSR_MIN, GK_CSR_AMIN. In case of 
+           GK_CSR_COS, the rows and the query are assumed to be of unit 
+           length.
+    \param nsim is the maximum number of requested most similar rows.
+           If -1 is provided, then everything is returned unsorted.
+    \param minsim is the minimum similarity of the requested most 
+           similar rows
+    \param hits is the result set. This array should be at least
+           of length nsim.
+    \param i_marker is an array of size equal to the number of rows
+           whose values are initialized to -1. If NULL is provided
+           then this array is allocated and freed internally.
+    \param i_cand is an array of size equal to the number of rows.
+           If NULL is provided then this array is allocated and freed 
+           internally.
+    \returns The number of identified most similar rows, which can be
+             smaller than the requested number of nnbrs in those cases
+             in which there are no sufficiently many neighbors.
+*/
+/**************************************************************************/
+int gk_csr_GetSimilarRows(gk_csr_t *mat, int nqterms, int *qind, 
+        float *qval, int simtype, int nsim, float minsim, gk_fkv_t *hits, 
+        int *i_marker, gk_fkv_t *i_cand)
+{
+  ssize_t i, ii, j, k;
+  int nrows, ncols, ncand;
+  ssize_t *colptr;
+  int *colind, *marker;
+  float *colval, *rnorms, mynorm, *rsums, mysum;
+  gk_fkv_t *cand;
+
+  if (nqterms == 0)
+    return 0;
+
+  nrows  = mat->nrows;
+  ncols  = mat->ncols;
+  GKASSERT((colptr = mat->colptr) != NULL);
+  GKASSERT((colind = mat->colind) != NULL);
+  GKASSERT((colval = mat->colval) != NULL);
+
+  marker = (i_marker ? i_marker : gk_ismalloc(nrows, -1, "gk_csr_SimilarRows: marker"));
+  cand   = (i_cand   ? i_cand   : gk_fkvmalloc(nrows, "gk_csr_SimilarRows: cand"));
+
+  switch (simtype) {
+    case GK_CSR_DOTP:
+    case GK_CSR_COS:
+      for (ncand=0, ii=0; ii<nqterms; ii++) {
+        i = qind[ii];
+        if (i < ncols) {
+          for (j=colptr[i]; j<colptr[i+1]; j++) {
+            k = colind[j];
+            if (marker[k] == -1) {
+              cand[ncand].val = k;
+              cand[ncand].key = 0;
+              marker[k]       = ncand++;
+            }
+            cand[marker[k]].key += colval[j]*qval[ii];
+          }
+        }
+      }
+      break;
+
+    case GK_CSR_JAC:
+      for (ncand=0, ii=0; ii<nqterms; ii++) {
+        i = qind[ii];
+        if (i < ncols) {
+          for (j=colptr[i]; j<colptr[i+1]; j++) {
+            k = colind[j];
+            if (marker[k] == -1) {
+              cand[ncand].val = k;
+              cand[ncand].key = 0;
+              marker[k]       = ncand++;
+            }
+            cand[marker[k]].key += colval[j]*qval[ii];
+          }
+        }
+      }
+
+      GKASSERT((rnorms = mat->rnorms) != NULL);
+      mynorm = gk_fdot(nqterms, qval, 1, qval, 1);
+
+      for (i=0; i<ncand; i++)
+        cand[i].key = cand[i].key/(rnorms[cand[i].val]+mynorm-cand[i].key);
+      break;
+
+    case GK_CSR_MIN:
+      for (ncand=0, ii=0; ii<nqterms; ii++) {
+        i = qind[ii];
+        if (i < ncols) {
+          for (j=colptr[i]; j<colptr[i+1]; j++) {
+            k = colind[j];
+            if (marker[k] == -1) {
+              cand[ncand].val = k;
+              cand[ncand].key = 0;
+              marker[k]       = ncand++;
+            }
+            cand[marker[k]].key += gk_min(colval[j], qval[ii]);
+          }
+        }
+      }
+
+      GKASSERT((rsums = mat->rsums) != NULL);
+      mysum = gk_fsum(nqterms, qval, 1);
+
+      for (i=0; i<ncand; i++)
+        cand[i].key = cand[i].key/(rsums[cand[i].val]+mysum-cand[i].key);
+      break;
+
+    /* Assymetric MIN  similarity */
+    case GK_CSR_AMIN:
+      for (ncand=0, ii=0; ii<nqterms; ii++) {
+        i = qind[ii];
+        if (i < ncols) {
+          for (j=colptr[i]; j<colptr[i+1]; j++) {
+            k = colind[j];
+            if (marker[k] == -1) {
+              cand[ncand].val = k;
+              cand[ncand].key = 0;
+              marker[k]       = ncand++;
+            }
+            cand[marker[k]].key += gk_min(colval[j], qval[ii]);
+          }
+        }
+      }
+
+      mysum = gk_fsum(nqterms, qval, 1);
+
+      for (i=0; i<ncand; i++)
+        cand[i].key = cand[i].key/mysum;
+      break;
+
+    default:
+      gk_errexit(SIGERR, "Unknown similarity measure %d\n", simtype);
+      return -1;
+  }
+
+  /* go and prune the hits that are bellow minsim */
+  for (j=0, i=0; i<ncand; i++) {
+    marker[cand[i].val] = -1;
+    if (cand[i].key >= minsim) 
+      cand[j++] = cand[i];
+  }
+  ncand = j;
+
+  if (nsim == -1 || nsim >= ncand) {
+    nsim = ncand;
+  }
+  else {
+    nsim = gk_min(nsim, ncand);
+    gk_dfkvkselect(ncand, nsim, cand);
+    gk_fkvsortd(nsim, cand);
+  }
+
+  gk_fkvcopy(nsim, cand, hits);
+
+  if (i_marker == NULL)
+    gk_free((void **)&marker, LTERM);
+  if (i_cand == NULL)
+    gk_free((void **)&cand, LTERM);
+
+  return nsim;
+}
+
+
+/*************************************************************************/
+/*! Returns a symmetric version of a square matrix. The symmetric version
+    is constructed by applying an A op A^T operation, where op is one of
+    GK_CSR_SYM_SUM, GK_CSR_SYM_MIN, GK_CSR_SYM_MAX, GK_CSR_SYM_AVG.
+   
+    \param mat the matrix to be symmetrized,
+    \param op indicates the operation to be performed. The possible values are
+           GK_CSR_SYM_SUM, GK_CSR_SYM_MIN, GK_CSR_SYM_MAX, and GK_CSR_SYM_AVG.
+
+    \returns the symmetrized matrix consisting only of its row-based structure. 
+          The input matrix is not modified. 
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_MakeSymmetric(gk_csr_t *mat, int op)
+{
+  ssize_t i, j, k, nnz;
+  int nrows, nadj, hasvals;
+  ssize_t *rowptr, *colptr, *nrowptr;
+  int *rowind, *colind, *nrowind, *marker, *ids;
+  float *rowval=NULL, *colval=NULL, *nrowval=NULL, *wgts=NULL;
+  gk_csr_t *nmat;
+
+  if (mat->nrows != mat->ncols) {
+    fprintf(stderr, "gk_csr_MakeSymmetric: The matrix needs to be square.\n");
+    return NULL;
+  }
+
+  hasvals = (mat->rowval != NULL);
+
+  nrows  = mat->nrows;
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+  if (hasvals)
+    rowval = mat->rowval;
+
+  /* create the column view for efficient processing */
+  colptr = gk_zsmalloc(nrows+1, 0, "colptr");
+  colind = gk_i32malloc(rowptr[nrows], "colind");
+  if (hasvals)
+    colval = gk_fmalloc(rowptr[nrows], "colval");
+
+  for (i=0; i<nrows; i++) {
+    for (j=rowptr[i]; j<rowptr[i+1]; j++) 
+      colptr[rowind[j]]++;
+  }
+  MAKECSR(i, nrows, colptr);
+
+  for (i=0; i<nrows; i++) {
+    for (j=rowptr[i]; j<rowptr[i+1]; j++) {
+      colind[colptr[rowind[j]]] = i;
+      if (hasvals)
+        colval[colptr[rowind[j]]] = rowval[j];
+      colptr[rowind[j]]++;
+    }
+  }
+  SHIFTCSR(i, nrows, colptr);
+
+
+  nmat = gk_csr_Create();
+  
+  nmat->nrows = mat->nrows;
+  nmat->ncols = mat->ncols;
+
+  nrowptr = nmat->rowptr = gk_zmalloc(nrows+1, "gk_csr_MakeSymmetric: nrowptr");
+  nrowind = nmat->rowind = gk_imalloc(2*rowptr[nrows], "gk_csr_MakeSymmetric: nrowind");
+  if (hasvals)
+    nrowval = nmat->rowval = gk_fmalloc(2*rowptr[nrows], "gk_csr_MakeSymmetric: nrowval");
+
+  marker = gk_ismalloc(nrows, -1, "marker");
+  ids    = gk_imalloc(nrows, "ids");
+  if (hasvals)
+    wgts = gk_fmalloc(nrows, "wgts");
+
+  nrowptr[0] = nnz = 0;
+  for (i=0; i<nrows; i++) {
+    nadj = 0;
+    /* out-edges */
+    for (j=rowptr[i]; j<rowptr[i+1]; j++) {
+      ids[nadj] = rowind[j]; 
+      if (hasvals)
+        wgts[nadj] = (op == GK_CSR_SYM_AVG ? 0.5*rowval[j] : rowval[j]);
+      marker[rowind[j]] = nadj++;
+    }
+
+    /* in-edges */
+    for (j=colptr[i]; j<colptr[i+1]; j++) {
+      if (marker[colind[j]] == -1) {
+        if (op != GK_CSR_SYM_MIN) {
+          ids[nadj] = colind[j]; 
+          if (hasvals) 
+            wgts[nadj] = (op == GK_CSR_SYM_AVG ? 0.5*colval[j] : colval[j]);
+          nadj++;
+        }
+      }
+      else {
+        if (hasvals) {
+          switch (op) {
+            case GK_CSR_SYM_MAX:
+              wgts[marker[colind[j]]] = gk_max(colval[j], wgts[marker[colind[j]]]);
+              break;
+            case GK_CSR_SYM_MIN:
+              wgts[marker[colind[j]]] = gk_min(colval[j], wgts[marker[colind[j]]]);
+              break;
+            case GK_CSR_SYM_SUM:
+              wgts[marker[colind[j]]] += colval[j];
+              break;
+            case GK_CSR_SYM_AVG:
+              wgts[marker[colind[j]]] = 0.5*(wgts[marker[colind[j]]] + colval[j]);
+              break;
+            default:
+              errexit("Unsupported op for MakeSymmetric!\n");
+          }
+        }
+        marker[colind[j]] = -1;
+      }
+    }
+
+    /* go over out edges again to resolve any edges that were not found in the in
+     * edges */
+    for (j=rowptr[i]; j<rowptr[i+1]; j++) {
+      if (marker[rowind[j]] != -1) {
+        if (op == GK_CSR_SYM_MIN)
+          ids[marker[rowind[j]]] = -1;
+        marker[rowind[j]] = -1;
+      }
+    }
+
+    /* put the non '-1' entries in ids[] into i's row */
+    for (j=0; j<nadj; j++) {
+      if (ids[j] != -1) {
+        nrowind[nnz] = ids[j];
+        if (hasvals)
+          nrowval[nnz] = wgts[j];
+        nnz++;
+      }
+    }
+    nrowptr[i+1] = nnz;
+  }
+
+  gk_free((void **)&colptr, &colind, &colval, &marker, &ids, &wgts, LTERM);
+
+  return nmat;
+}
+
+
+/*************************************************************************/
+/*! This function finds the connected components in a graph stored in
+    CSR format.
+
+    \param mat is the graph structure in CSR format
+    \param cptr is the ptr structure of the CSR representation of the 
+           components. The length of this vector must be mat->nrows+1.
+    \param cind is the indices structure of the CSR representation of 
+           the components. The length of this vector must be mat->nrows.
+    \param cids is an array that stores the component # of each vertex
+           of the graph. The length of this vector must be mat->nrows.
+
+    \returns the number of components that it found.
+
+    \note The cptr, cind, and cids parameters can be NULL, in which case 
+          only the number of connected components is returned.
+*/
+/*************************************************************************/
+int gk_csr_FindConnectedComponents(gk_csr_t *mat, int32_t *cptr, int32_t *cind, 
+        int32_t *cids)
+{
+  ssize_t i, ii, j, jj, k, nvtxs, first, last, ntodo, ncmps;
+  ssize_t *xadj;
+  int32_t *adjncy, *pos, *todo;
+  int32_t mustfree_ccsr=0, mustfree_where=0;
+
+  if (mat->nrows != mat->ncols) {
+    fprintf(stderr, "gk_csr_FindComponents: The matrix needs to be square.\n");
+    return -1;
+  }
+
+  nvtxs  = mat->nrows;
+  xadj   = mat->rowptr;
+  adjncy = mat->rowind;
+
+  /* Deal with NULL supplied cptr/cind vectors */
+  if (cptr == NULL) {
+    cptr = gk_i32malloc(nvtxs+1, "gk_csr_FindComponents: cptr");
+    cind = gk_i32malloc(nvtxs, "gk_csr_FindComponents: cind");
+    mustfree_ccsr = 1;
+  }
+
+  /* The list of vertices that have not been touched yet. 
+     The valid entries are from [0..ntodo). */
+  todo = gk_i32incset(nvtxs, 0, gk_i32malloc(nvtxs, "gk_csr_FindComponents: todo"));
+
+  /* For a vertex that has not been visited, pos[i] is the position in the
+     todo list that this vertex is stored. 
+     If a vertex has been visited, pos[i] = -1. */
+  pos = gk_i32incset(nvtxs, 0, gk_i32malloc(nvtxs, "gk_csr_FindComponents: pos"));
+
+
+  /* Find the connected componends */
+  ncmps = -1;
+  ntodo = nvtxs;     /* All vertices have not been visited */
+  first = last = 0;  /* Point to the first and last vertices that have been touched
+                        but not explored. 
+                        These vertices are stored in cind[first]...cind[last-1]. */
+
+  while (first < last || ntodo > 0) {
+    if (first == last) { /* Find another starting vertex */
+      cptr[++ncmps] = first;  /* Mark the end of the current CC */
+
+      /* put the first vertex in the todo list as the start of the new CC */
+      ASSERT(pos[todo[0]] != -1);
+      cind[last++] = todo[0];  
+
+      pos[todo[0]] = -1;
+      todo[0] = todo[--ntodo];
+      pos[todo[0]] = 0;
+    }
+
+    i = cind[first++];  /* Get the first visited but unexplored vertex */
+
+    for (j=xadj[i]; j<xadj[i+1]; j++) {
+      k = adjncy[j];
+      if (pos[k] != -1) {
+        cind[last++] = k;
+
+        /* Remove k from the todo list and put the last item in the todo 
+           list at the position that k was so that the todo list will be
+           consequtive. The pos[] array is updated accordingly to keep track
+           the location of the vertices in the todo[] list. */
+        todo[pos[k]] = todo[--ntodo];
+        pos[todo[pos[k]]] = pos[k];
+        pos[k] = -1;
+      }
+    }
+  }
+  cptr[++ncmps] = first;
+
+  /* see if we need to return cids */
+  if (cids != NULL) {
+    for (i=0; i<ncmps; i++) {
+      for (j=cptr[i]; j<cptr[i+1]; j++)
+        cids[cind[j]] = i;
+    }
+  }
+
+  if (mustfree_ccsr)
+    gk_free((void **)&cptr, &cind, LTERM);
+
+  gk_free((void **)&pos, &todo, LTERM);
+
+  return (int) ncmps;
+}
+
+
+/*************************************************************************/
+/*! Returns a matrix that has been reordered according to the provided
+    row/column permutation. The matrix is required to be square and the same
+    permutation is applied to both rows and columns.
+
+    \param[IN] mat is the matrix to be re-ordered.
+    \param[IN] perm is the new ordering of the rows & columns
+    \param[IN] iperm is the original ordering of the re-ordered matrix's rows & columns
+    \returns the newly created reordered matrix.
+
+    \note Either perm or iperm can be NULL but not both.
+*/
+/**************************************************************************/
+gk_csr_t *gk_csr_ReorderSymmetric(gk_csr_t *mat, int32_t *perm, int32_t *iperm)
+{
+  ssize_t j, jj;
+  ssize_t *rowptr, *nrowptr;
+  int i, k, u, v, nrows;
+  int freeperm=0, freeiperm=0;
+  int32_t *rowind, *nrowind;
+  float *rowval, *nrowval;
+  gk_csr_t *nmat;
+
+  if (mat->nrows != mat->ncols) {
+    fprintf(stderr, "gk_csr_ReorderSymmetric: The matrix needs to be square.\n");
+    return NULL;
+  }
+
+  if (perm == NULL && iperm == NULL)
+    return NULL;
+
+  nrows  = mat->nrows;
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+  rowval = mat->rowval;
+
+  nmat = gk_csr_Create();
+
+  nmat->nrows = nrows;
+  nmat->ncols = nrows;
+
+  nrowptr = nmat->rowptr = gk_zmalloc(nrows+1, "gk_csr_ReorderSymmetric: rowptr");
+  nrowind = nmat->rowind = gk_i32malloc(rowptr[nrows], "gk_csr_ReorderSymmetric: rowind");
+  nrowval = nmat->rowval = gk_fmalloc(rowptr[nrows], "gk_csr_ReorderSymmetric: rowval");
+
+  /* allocate memory for the different structures present in the matrix */
+  if (mat->rlabels)
+    nmat->rlabels = gk_i32malloc(nrows, "gk_csr_ReorderSymmetric: rlabels");
+  if (mat->rmap)
+    nmat->rmap = gk_i32malloc(nrows, "gk_csr_ReorderSymmetric: rmap");
+  if (mat->rnorms)
+    nmat->rnorms = gk_fmalloc(nrows, "gk_csr_ReorderSymmetric: rnorms");
+  if (mat->rsums)
+    nmat->rsums = gk_fmalloc(nrows, "gk_csr_ReorderSymmetric: rsums");
+  if (mat->rsizes)
+    nmat->rsizes = gk_fmalloc(nrows, "gk_csr_ReorderSymmetric: rsizes");
+  if (mat->rvols)
+    nmat->rvols = gk_fmalloc(nrows, "gk_csr_ReorderSymmetric: rvols");
+  if (mat->rwgts)
+    nmat->rwgts = gk_fmalloc(nrows, "gk_csr_ReorderSymmetric: rwgts");
+
+  if (mat->clabels)
+    nmat->clabels = gk_i32malloc(nrows, "gk_csr_ReorderSymmetric: clabels");
+  if (mat->cmap)
+    nmat->cmap = gk_i32malloc(nrows, "gk_csr_ReorderSymmetric: cmap");
+  if (mat->cnorms)
+    nmat->cnorms = gk_fmalloc(nrows, "gk_csr_ReorderSymmetric: cnorms");
+  if (mat->csums)
+    nmat->csums = gk_fmalloc(nrows, "gk_csr_ReorderSymmetric: csums");
+  if (mat->csizes)
+    nmat->csizes = gk_fmalloc(nrows, "gk_csr_ReorderSymmetric: csizes");
+  if (mat->cvols)
+    nmat->cvols = gk_fmalloc(nrows, "gk_csr_ReorderSymmetric: cvols");
+  if (mat->cwgts)
+    nmat->cwgts = gk_fmalloc(nrows, "gk_csr_ReorderSymmetric: cwgts");
+
+
+
+  /* create perm/iperm if not provided */
+  if (perm == NULL) {
+    freeperm = 1;
+    perm = gk_i32malloc(nrows, "gk_csr_ReorderSymmetric: perm"); 
+    for (i=0; i<nrows; i++)
+      perm[iperm[i]] = i;
+  }
+  if (iperm == NULL) {
+    freeiperm = 1;
+    iperm = gk_i32malloc(nrows, "gk_csr_ReorderSymmetric: iperm"); 
+    for (i=0; i<nrows; i++)
+      iperm[perm[i]] = i;
+  }
+
+  /* fill-in the information of the re-ordered matrix */
+  nrowptr[0] = jj = 0;
+  for (v=0; v<nrows; v++) {
+    u = iperm[v];
+    for (j=rowptr[u]; j<rowptr[u+1]; j++, jj++) {
+      nrowind[jj] = perm[rowind[j]];
+      nrowval[jj] = rowval[j];
+    }
+
+    if (mat->rlabels)
+      nmat->rlabels[v] = mat->rlabels[u];
+    if (mat->rmap)
+      nmat->rmap[v] = mat->rmap[u];
+    if (mat->rnorms)
+      nmat->rnorms[v] = mat->rnorms[u];
+    if (mat->rsums)
+      nmat->rsums[v] = mat->rsums[u];
+    if (mat->rsizes)
+      nmat->rsizes[v] = mat->rsizes[u];
+    if (mat->rvols)
+      nmat->rvols[v] = mat->rvols[u];
+    if (mat->rwgts)
+      nmat->rwgts[v] = mat->rwgts[u];
+
+    if (mat->clabels)
+      nmat->clabels[v] = mat->clabels[u];
+    if (mat->cmap)
+      nmat->cmap[v] = mat->cmap[u];
+    if (mat->cnorms)
+      nmat->cnorms[v] = mat->cnorms[u];
+    if (mat->csums)
+      nmat->csums[v] = mat->csums[u];
+    if (mat->csizes)
+      nmat->csizes[v] = mat->csizes[u];
+    if (mat->cvols)
+      nmat->cvols[v] = mat->cvols[u];
+    if (mat->cwgts)
+      nmat->cwgts[v] = mat->cwgts[u];
+
+    nrowptr[v+1] = jj;
+  }
+
+
+  /* free memory */
+  if (freeperm)
+    gk_free((void **)&perm, LTERM);
+  if (freeiperm)
+    gk_free((void **)&iperm, LTERM);
+
+  return nmat;
+}
+
+
+/*************************************************************************/
+/*! This function computes a permutation of the rows/columns of a symmetric
+    matrix based on a breadth-first-traversal. It can be used for re-ordering 
+    the matrix to reduce its bandwidth for better cache locality.
+
+    \param[IN]  mat is the matrix whose ordering to be computed.
+    \param[IN]  maxdegree is the maximum number of nonzeros of the rows that
+                will participate in the BFS ordering. Rows with more nonzeros
+                will be put at the front of the ordering in decreasing degree
+                order. 
+    \param[IN]  v is the starting row of the BFS. A value of -1 indicates that
+                a randomly selected row will be used.
+    \param[OUT] perm[i] stores the ID of row i in the re-ordered matrix.
+    \param[OUT] iperm[i] stores the ID of the row that corresponds to 
+                the ith vertex in the re-ordered matrix.
+
+    \note The perm or iperm (but not both) can be NULL, at which point, 
+          the corresponding arrays are not returned. Though the program
+          works fine when both are NULL, doing that is not smart.
+          The returned arrays should be freed with gk_free().
+*/
+/*************************************************************************/
+void gk_csr_ComputeBFSOrderingSymmetric(gk_csr_t *mat, int maxdegree, int v, 
+          int32_t **r_perm, int32_t **r_iperm)
+{
+  int i, k, nrows, first, last;
+  ssize_t j, *rowptr;
+  int32_t *rowind, *cot, *pos;
+
+  if (mat->nrows != mat->ncols) {
+    fprintf(stderr, "gk_csr_ComputeBFSOrderingSymmetric: The matrix needs to be square.\n");
+    return;
+  }
+  if (maxdegree < mat->nrows && v != -1) {
+    fprintf(stderr, "gk_csr_ComputeBFSOrderingSymmetric: Since maxdegree node renumbering is requested the starting row should be -1.\n");
+    return;
+  }
+  if (mat->nrows <= 0)
+    return;
+
+  nrows  = mat->nrows;
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+
+  /* This array will function like pos + touched of the CC method */
+  pos = gk_i32incset(nrows, 0, gk_i32malloc(nrows, "gk_csr_ComputeBFSOrderingSymmetric: pos"));
+
+  /* This array ([C]losed[O]pen[T]odo => cot) serves three purposes. 
+     Positions from [0...first) is the current iperm[] vector of the explored rows; 
+     Positions from [first...last) is the OPEN list (i.e., visited rows);
+     Positions from [last...nrows) is the todo list. */
+  cot = gk_i32incset(nrows, 0, gk_i32malloc(nrows, "gk_csr_ComputeBFSOrderingSymmetric: cot"));
+
+  first = last = 0;
+
+  /* deal with maxdegree handling */
+  if (maxdegree < nrows) {
+    last = nrows;
+    for (i=nrows-1; i>=0; i--) {
+      if (rowptr[i+1]-rowptr[i] < maxdegree) {
+        cot[--last] = i;
+        pos[i] = last;
+      }
+      else {
+        cot[first++] = i;
+        pos[i] = -1;
+      }
+    }
+    GKASSERT(first == last);
+
+    if (last > 0) { /* reorder them in degree decreasing order */
+      gk_ikv_t *cand = gk_ikvmalloc(first, "gk_csr_ComputeBFSOrderingSymmetric: cand");
+
+      for (i=0; i<first; i++) {
+        k = cot[i];
+        cand[i].key = (int)(rowptr[k+1]-rowptr[k]);
+        cand[i].val = k;
+      }
+
+      gk_ikvsortd(first, cand);
+      for (i=0; i<first; i++) 
+        cot[i] = cand[i].val;
+
+      gk_free((void **)&cand, LTERM);
+    }
+
+    v = cot[last + RandomInRange(nrows-last)];
+  }
+
+
+  /* swap v with the front of the todo list */
+  cot[pos[v]] = cot[last];
+  pos[cot[last]] = pos[v];
+
+  cot[last] = v;
+  pos[v] = last;
+
+
+  /* start processing the nodes */
+  while (first < nrows) {
+    if (first == last) { /* find another starting row */
+      k = cot[last];
+      GKASSERT(pos[k] != -1);
+      pos[k] = -1; /* mark node as being visited */
+      last++;
+    }
+
+    i = cot[first++];  /* the ++ advances the explored rows */
+    for (j=rowptr[i]; j<rowptr[i+1]; j++) {
+      k = rowind[j];
+      /* if a node has already been visited, its perm[] will be -1 */
+      if (pos[k] != -1) {
+        /* pos[k] is the location within iperm of where k resides (it is in the 'todo' part); 
+           It is placed in that location cot[last] (end of OPEN list) that we 
+           are about to overwrite and update pos[cot[last]] to reflect that. */
+        cot[pos[k]]    = cot[last]; /* put the head of the todo list to 
+                                       where k was in the todo list */
+        pos[cot[last]] = pos[k];    /* update perm to reflect the move */
+
+        cot[last++] = k;  /* put node at the end of the OPEN list */
+        pos[k]      = -1; /* mark node as being visited */
+      }
+    }
+  }
+
+  /* time to decide what to return */
+  if (r_perm != NULL) {
+    /* use the 'pos' array to build the perm array */
+    for (i=0; i<nrows; i++)
+      pos[cot[i]] = i;
+
+    *r_perm = pos;
+    pos = NULL;
+  }
+
+  if (r_iperm != NULL) {
+    *r_iperm = cot;
+    cot = NULL;
+  }
+
+  /* cleanup memory */
+  gk_free((void **)&pos, &cot, LTERM);
+
+}
+
+
+/*************************************************************************/
+/*! This function computes a permutation of the rows of a symmetric matrix
+    based on a best-first-traversal. It can be used for re-ordering the matrix
+    to reduce its bandwidth for better cache locality.
+
+    \param[IN]  mat is the matrix structure.
+    \param[IN]  v is the starting row of the best-first traversal.
+    \param[IN]  type indicates the criteria to use to measure the 'bestness'
+                of a row.
+    \param[OUT] perm[i] stores the ID of row i in the re-ordered matrix.
+    \param[OUT] iperm[i] stores the ID of the row that corresponds to 
+                the ith row in the re-ordered matrix.
+
+    \note The perm or iperm (but not both) can be NULL, at which point, 
+          the corresponding arrays are not returned. Though the program
+          works fine when both are NULL, doing that is not smart.
+          The returned arrays should be freed with gk_free().
+*/
+/*************************************************************************/
+void gk_csr_ComputeBestFOrderingSymmetric(gk_csr_t *mat, int v, int type, 
+          int32_t **r_perm, int32_t **r_iperm)
+{
+  ssize_t j, jj, *rowptr;
+  int i, k, u, nrows, nopen, ntodo;
+  int32_t *rowind, *perm, *degrees, *wdegrees, *sod, *level, *ot, *pos;
+  gk_i32pq_t *queue;
+
+  if (mat->nrows != mat->ncols) {
+    fprintf(stderr, "gk_csr_ComputeBestFOrderingSymmetric: The matrix needs to be square.\n");
+    return;
+  }
+  if (mat->nrows <= 0)
+    return;
+
+  nrows  = mat->nrows;
+  rowptr = mat->rowptr;
+  rowind = mat->rowind;
+
+
+  /* the degree of the vertices in the closed list */
+  degrees = gk_i32smalloc(nrows, 0, "gk_csr_ComputeBestFOrderingSymmetric: degrees");
+
+  /* the weighted degree of the vertices in the closed list for type==3 */
+  wdegrees = gk_i32smalloc(nrows, 0, "gk_csr_ComputeBestFOrderingSymmetric: wdegrees");
+
+  /* the sum of differences for type==4 */
+  sod = gk_i32smalloc(nrows, 0, "gk_csr_ComputeBestFOrderingSymmetric: sod");
+
+  /* the encountering level of a vertex type==5 */
+  level = gk_i32smalloc(nrows, 0, "gk_csr_ComputeBestFOrderingSymmetric: level");
+
+  /* The open+todo list of vertices. 
+     The vertices from [0..nopen] are the open vertices.
+     The vertices from [nopen..ntodo) are the todo vertices.
+     */
+  ot = gk_i32incset(nrows, 0, gk_i32malloc(nrows, "gk_csr_ComputeBestFOrderingSymmetric: ot"));
+
+  /* For a vertex that has not been explored, pos[i] is the position in the ot list. */
+  pos = gk_i32incset(nrows, 0, gk_i32malloc(nrows, "gk_csr_ComputeBestFOrderingSymmetric: pos"));
+
+  /* if perm[i] >= 0, then perm[i] is the order of vertex i; otherwise perm[i] == -1. */
+  perm = gk_i32smalloc(nrows, -1, "gk_csr_ComputeBestFOrderingSymmetric: perm");
+
+  /* create the queue and put the starting vertex in it */
+  queue = gk_i32pqCreate(nrows);
+  gk_i32pqInsert(queue, v, 1);
+
+  /* put v at the front of the open list */
+  pos[0] = ot[0] = v;
+  pos[v] = ot[v] = 0;
+  nopen = 1;
+  ntodo = nrows;
+
+  /* start processing the nodes */
+  for (i=0; i<nrows; i++) {
+    if (nopen == 0) { /* deal with non-connected graphs */
+      gk_i32pqInsert(queue, ot[0], 1);  
+      nopen++;
+    }
+
+    if ((v = gk_i32pqGetTop(queue)) == -1)
+      gk_errexit(SIGERR, "The priority queue got empty ahead of time [i=%d].\n", i);
+
+    if (perm[v] != -1)
+      gk_errexit(SIGERR, "The perm[%d] has already been set.\n", v);
+    perm[v] = i;
+
+    if (ot[pos[v]] != v)
+      gk_errexit(SIGERR, "Something went wrong [ot[pos[%d]]!=%d.\n", v, v);
+    if (pos[v] >= nopen)
+      gk_errexit(SIGERR, "The position of v is not in open list. pos[%d]=%d is >=%d.\n", v, pos[v], nopen);
+
+    /* remove v from the open list and re-arrange the todo part of the list */
+    ot[pos[v]]       = ot[nopen-1];
+    pos[ot[nopen-1]] = pos[v];
+    if (ntodo > nopen) {
+      ot[nopen-1]      = ot[ntodo-1];
+      pos[ot[ntodo-1]] = nopen-1;
+    }
+    nopen--;
+    ntodo--;
+
+    for (j=rowptr[v]; j<rowptr[v+1]; j++) {
+      u = rowind[j];
+      if (perm[u] == -1) {
+        /* update ot list, if u is not in the open list by putting it at the end
+           of the open list. */
+        if (degrees[u] == 0) {
+          ot[pos[u]]     = ot[nopen];
+          pos[ot[nopen]] = pos[u];
+          ot[nopen]      = u;
+          pos[u]         = nopen;
+          nopen++;
+
+          level[u] = level[v]+1;
+          gk_i32pqInsert(queue, u, 0);  
+        }
+
+
+        /* update the in-closed degree */
+        degrees[u]++;
+
+        /* update the queues based on the type */
+        switch (type) {
+          case 1: /* DFS */
+            gk_i32pqUpdate(queue, u, 1000*(i+1)+degrees[u]);
+            break;
+
+          case 2: /* Max in closed degree */
+            gk_i32pqUpdate(queue, u, degrees[u]);
+            break;
+
+          case 3: /* Sum of orders in closed list */
+            wdegrees[u] += i;
+            gk_i32pqUpdate(queue, u, wdegrees[u]);
+            break;
+
+          case 4: /* Sum of order-differences */
+            /* this is handled at the end of the loop */
+            ;
+            break;
+
+          case 5: /* BFS with in degree priority */
+            gk_i32pqUpdate(queue, u, -(1000*level[u] - degrees[u]));
+            break;
+
+          case 6: /* Hybrid of 1+2 */
+            gk_i32pqUpdate(queue, u, (i+1)*degrees[u]);
+            break;
+
+          default:
+            ;
+        }
+      }
+    }
+
+    if (type == 4) { /* update all the vertices in the open list */
+      for (j=0; j<nopen; j++) {
+        u = ot[j];
+        if (perm[u] != -1)
+          gk_errexit(SIGERR, "For i=%d, the open list contains a closed row: ot[%zd]=%d, perm[%d]=%d.\n", i, j, u, u, perm[u]);
+        sod[u] += degrees[u];
+        if (i<1000 || i%25==0)
+          gk_i32pqUpdate(queue, u, sod[u]);
+      }
+    }
+
+    /*
+    for (j=0; j<ntodo; j++) {
+      if (pos[ot[j]] != j)
+        gk_errexit(SIGERR, "pos[ot[%zd]] != %zd.\n", j, j);
+    }
+    */
+
+  }
+
+
+  /* time to decide what to return */
+  if (r_iperm != NULL) {
+    /* use the 'degrees' array to build the iperm array */
+    for (i=0; i<nrows; i++)
+      degrees[perm[i]] = i;
+
+    *r_iperm = degrees;
+    degrees = NULL;
+  }
+
+  if (r_perm != NULL) {
+    *r_perm = perm;
+    perm = NULL;
+  }
+
+
+
+
+  /* cleanup memory */
+  gk_i32pqDestroy(queue);
+  gk_free((void **)&perm, &degrees, &wdegrees, &sod, &ot, &pos, &level, LTERM);
+
+}
+

third_party/METIS/GKlib/error.c

+/*!
+\file  error.c
+\brief Various error-handling functions
+
+This file contains functions dealing with error reporting and termination
+
+\author George
+\date 1/1/2007
+\version\verbatim $Id: error.c 10711 2011-08-31 22:23:04Z karypis $ \endverbatim
+*/
+
+
+#define _GK_ERROR_C_  /* this is needed to properly declare the gk_jub* variables
+                         as an extern function in GKlib.h */
+
+#include <GKlib.h>
+
+
+/* These are the jmp_buf for the graceful exit in case of severe errors.
+   Multiple buffers are defined to allow for recursive invokation. */
+#define MAX_JBUFS 128
+__thread int gk_cur_jbufs=-1;
+__thread jmp_buf gk_jbufs[MAX_JBUFS];
+__thread jmp_buf gk_jbuf;
+
+typedef void (*gksighandler_t)(int);
+
+/* These are the holders of the old singal handlers for the trapped signals */
+static __thread gksighandler_t old_SIGMEM_handler;  /* Custom signal */
+static __thread gksighandler_t old_SIGERR_handler;  /* Custom signal */
+static __thread gksighandler_t old_SIGMEM_handlers[MAX_JBUFS];  /* Custom signal */
+static __thread gksighandler_t old_SIGERR_handlers[MAX_JBUFS];  /* Custom signal */
+
+/* The following is used to control if the gk_errexit() will actually abort or not.
+   There is always a single copy of this variable */
+static int gk_exit_on_error = 1;
+
+
+/*************************************************************************/
+/*! This function sets the gk_exit_on_error variable 
+ */
+/*************************************************************************/
+void gk_set_exit_on_error(int value)
+{
+  gk_exit_on_error = value;
+}
+
+
+
+/*************************************************************************/
+/*! This function prints an error message and exits  
+ */
+/*************************************************************************/
+void errexit(char *f_str,...)
+{
+  va_list argp;
+
+  va_start(argp, f_str);
+  vfprintf(stderr, f_str, argp);
+  va_end(argp);
+
+  if (strlen(f_str) == 0 || f_str[strlen(f_str)-1] != '\n')
+        fprintf(stderr,"\n");
+  fflush(stderr);
+
+  if (gk_exit_on_error)
+    exit(-2);
+
+  /* abort(); */
+}
+
+
+/*************************************************************************/
+/*! This function prints an error message and raises a signum signal
+ */
+/*************************************************************************/
+void gk_errexit(int signum, char *f_str,...)
+{
+  va_list argp;
+
+  va_start(argp, f_str);
+  vfprintf(stderr, f_str, argp);
+  va_end(argp);
+
+  fprintf(stderr,"\n");
+  fflush(stderr);
+
+  if (gk_exit_on_error)
+    raise(signum);
+}
+
+
+/***************************************************************************/
+/*! This function sets a number of signal handlers and sets the return point 
+    of a longjmp
+*/
+/***************************************************************************/
+int gk_sigtrap() 
+{
+  if (gk_cur_jbufs+1 >= MAX_JBUFS)
+    return 0;
+
+  gk_cur_jbufs++;
+
+  old_SIGMEM_handlers[gk_cur_jbufs]  = signal(SIGMEM,  gk_sigthrow);
+  old_SIGERR_handlers[gk_cur_jbufs]  = signal(SIGERR,  gk_sigthrow);
+
+  return 1;
+}
+  
+
+/***************************************************************************/
+/*! This function sets the handlers for the signals to their default handlers
+ */
+/***************************************************************************/
+int gk_siguntrap() 
+{
+  if (gk_cur_jbufs == -1)
+    return 0;
+
+  signal(SIGMEM,  old_SIGMEM_handlers[gk_cur_jbufs]);
+  signal(SIGERR,  old_SIGERR_handlers[gk_cur_jbufs]);
+
+  gk_cur_jbufs--;
+
+  return 1;
+}
+  
+
+/*************************************************************************/
+/*! This function is the custome signal handler, which all it does is to
+    perform a longjump to the most recent saved environment 
+ */
+/*************************************************************************/
+void gk_sigthrow(int signum)
+{
+  longjmp(gk_jbufs[gk_cur_jbufs], signum);
+}
+  
+
+/***************************************************************************
+* This function sets a number of signal handlers and sets the return point 
+* of a longjmp
+****************************************************************************/
+void gk_SetSignalHandlers() 
+{
+  old_SIGMEM_handler = signal(SIGMEM,  gk_NonLocalExit_Handler);
+  old_SIGERR_handler = signal(SIGERR,  gk_NonLocalExit_Handler);
+}
+  
+
+/***************************************************************************
+* This function sets the handlers for the signals to their default handlers
+****************************************************************************/
+void gk_UnsetSignalHandlers() 
+{
+  signal(SIGMEM,  old_SIGMEM_handler);
+  signal(SIGERR,  old_SIGERR_handler);
+}
+  
+
+/*************************************************************************
+* This function is the handler for SIGUSR1 that implements the cleaning up 
+* process prior to a non-local exit.
+**************************************************************************/
+void gk_NonLocalExit_Handler(int signum)
+{
+  longjmp(gk_jbuf, signum);
+}
+  
+
+/*************************************************************************/
+/*! \brief Thread-safe implementation of strerror() */
+/**************************************************************************/
+char *gk_strerror(int errnum)
+{
+#if defined(WIN32) || defined(__MINGW32__)
+  return strerror(errnum);
+#else 
+#ifndef SUNOS
+  static __thread char buf[1024];
+
+  strerror_r(errnum, buf, 1024);
+
+  buf[1023] = '\0';
+  return buf;
+#else
+  return strerror(errnum);
+#endif
+#endif
+}
+
+
+
+/*************************************************************************
+* This function prints a backtrace of calling functions
+**************************************************************************/
+void PrintBackTrace()
+{
+#ifdef HAVE_EXECINFO_H
+  void *array[10];
+  int i, size;
+  char **strings;
+
+  size = backtrace(array, 10);
+  strings = backtrace_symbols(array, size);
+  
+  printf("Obtained %d stack frames.\n", size);
+  for (i=0; i<size; i++) {
+    printf("%s\n", strings[i]);
+  }
+  free(strings);
+#endif
+}

third_party/METIS/GKlib/evaluate.c

+/*!
+  \file  evaluate.c
+  \brief Various routines to evaluate classification performance
+
+  \author George
+  \date 9/23/2008
+  \version\verbatim $Id: evaluate.c 13328 2012-12-31 14:57:40Z karypis $ \endverbatim
+*/
+
+#include <GKlib.h>
+
+/**********************************************************************
+ * This function computes the max accuracy score of a ranked list,
+ * given +1/-1 class list
+ **********************************************************************/
+float ComputeAccuracy(int n, gk_fkv_t *list)
+{
+  int i, P, N, TP, FN = 0;
+  float bAccuracy = 0.0;
+  float acc;
+  
+  for (P=0, i=0;i<n;i++)
+    P += (list[i].val == 1? 1 : 0);
+  N = n - P;
+  
+  TP = FN = 0;
+  
+  for(i=0; i<n; i++){
+    if (list[i].val == 1)
+      TP++; 
+    else
+      FN++;
+    
+    acc = (TP + N - FN) * 100.0/ (P + N) ;
+    if (acc > bAccuracy)
+      bAccuracy = acc;
+  }
+  
+  return bAccuracy;
+}
+
+
+/*****************************************************************************
+ * This function computes the ROC score of a ranked list, given a +1/-1 class
+ * list.
+ ******************************************************************************/
+float ComputeROCn(int n, int maxN, gk_fkv_t *list)
+{
+  int i, P, TP, FP, TPprev, FPprev, AUC;
+  float prev;
+  
+  FP = TP = FPprev = TPprev = AUC = 0;
+  prev = list[0].key -1;
+  
+  for (P=0, i=0; i<n; i++)
+    P += (list[i].val == 1 ? 1 : 0);
+  
+  for (i=0; i<n && FP < maxN; i++) {
+    if (list[i].key != prev) {
+      AUC += (TP+TPprev)*(FP-FPprev)/2;
+      prev = list[i].key;
+      FPprev = FP;
+      TPprev = TP;
+    }
+    if (list[i].val == 1) 
+      TP++;
+    else {
+      FP++;
+    }
+  }
+  AUC += (TP+TPprev)*(FP-FPprev)/2;
+
+  return (TP*FP > 0 ? (float)(1.0*AUC/(P*FP)) : 0.0);
+}
+
+
+/*****************************************************************************
+* This function computes the median rate of false positive for each positive
+* instance.
+******************************************************************************/
+float ComputeMedianRFP(int n, gk_fkv_t *list)
+{
+  int i, P, N, TP, FP;
+
+  P = N = 0;
+  for (i=0; i<n; i++) {
+    if (list[i].val == 1)
+      P++;
+    else
+      N++;
+  }
+  
+  FP = TP = 0;
+  for (i=0; i<n && TP < (P+1)/2; i++) {
+    if (list[i].val == 1) 
+      TP++;
+    else 
+      FP++;
+  }
+  
+  return 1.0*FP/N;
+}
+
+/*********************************************************
+ * Compute the mean
+ ********************************************************/
+float ComputeMean (int n, float *values)
+{
+  int i;
+  float mean = 0.0;
+
+  for(i=0; i < n; i++)
+    mean += values[i];
+  
+  return 1.0 * mean/ n;
+}
+
+/********************************************************
+ * Compute the standard deviation
+ ********************************************************/
+float ComputeStdDev(int  n, float *values)
+{
+  int i;
+  float mean = ComputeMean(n, values);
+  float stdDev = 0;
+  
+  for(i=0;i<n;i++){
+    stdDev += (values[i] - mean)* (values[i] - mean);
+  }
+  
+  return sqrt(1.0 * stdDev/n);
+}

third_party/METIS/GKlib/fkvkselect.c

+/*!
+\file  dfkvkselect.c
+\brief Sorts only the largest k values
+ 
+\date   Started 7/14/00
+\author George
+\version\verbatim $Id: fkvkselect.c 10711 2011-08-31 22:23:04Z karypis $\endverbatim
+*/
+
+
+#include <GKlib.h>
+
+/* Byte-wise swap two items of size SIZE. */
+#define QSSWAP(a, b, stmp) do { stmp = (a); (a) = (b); (b) = stmp; } while (0)
+
+
+/******************************************************************************/
+/*! This function puts the 'topk' largest values in the beginning of the array */
+/*******************************************************************************/
+int gk_dfkvkselect(size_t n, int topk, gk_fkv_t *cand)
+{
+  int i, j, lo, hi, mid;
+  gk_fkv_t stmp;
+  float pivot;
+
+  if (n <= topk)
+    return n; /* return if the array has fewer elements than we want */
+
+  for (lo=0, hi=n-1; lo < hi;) {
+    mid = lo + ((hi-lo) >> 1);
+
+    /* select the median */
+    if (cand[lo].key < cand[mid].key)
+      mid = lo;
+    if (cand[hi].key > cand[mid].key)
+      mid = hi;
+    else 
+      goto jump_over;
+    if (cand[lo].key < cand[mid].key)
+      mid = lo;
+
+jump_over:
+    QSSWAP(cand[mid], cand[hi], stmp);
+    pivot = cand[hi].key;
+
+    /* the partitioning algorithm */
+    for (i=lo-1, j=lo; j<hi; j++) {
+      if (cand[j].key >= pivot) {
+        i++;
+        QSSWAP(cand[i], cand[j], stmp);
+      }
+    }
+    i++;
+    QSSWAP(cand[i], cand[hi], stmp);
+
+
+    if (i > topk) 
+      hi = i-1;
+    else if (i < topk)
+      lo = i+1;
+    else
+      break;
+  }
+
+/*
+  if (cand[lo].key < cand[hi].key)
+    printf("Hmm Error: %d %d %d %f %f\n", i, lo, hi, cand[lo].key, cand[hi].key);
+
+
+  for (i=topk; i<n; i++) {
+    for (j=0; j<topk; j++)
+      if (cand[i].key > cand[j].key)
+        printf("Hmm Error: %d %d %f %f %d %d\n", i, j, cand[i].key, cand[j].key, lo, hi);
+  }
+*/
+
+  return topk;
+}
+
+
+/******************************************************************************/
+/*! This function puts the 'topk' smallest values in the beginning of the array */
+/*******************************************************************************/
+int gk_ifkvkselect(size_t n, int topk, gk_fkv_t *cand)
+{
+  int i, j, lo, hi, mid;
+  gk_fkv_t stmp;
+  float pivot;
+
+  if (n <= topk)
+    return n; /* return if the array has fewer elements than we want */
+
+  for (lo=0, hi=n-1; lo < hi;) {
+    mid = lo + ((hi-lo) >> 1);
+
+    /* select the median */
+    if (cand[lo].key > cand[mid].key)
+      mid = lo;
+    if (cand[hi].key < cand[mid].key)
+      mid = hi;
+    else 
+      goto jump_over;
+    if (cand[lo].key > cand[mid].key)
+      mid = lo;
+
+jump_over:
+    QSSWAP(cand[mid], cand[hi], stmp);
+    pivot = cand[hi].key;
+
+    /* the partitioning algorithm */
+    for (i=lo-1, j=lo; j<hi; j++) {
+      if (cand[j].key <= pivot) {
+        i++;
+        QSSWAP(cand[i], cand[j], stmp);
+      }
+    }
+    i++;
+    QSSWAP(cand[i], cand[hi], stmp);
+
+
+    if (i > topk) 
+      hi = i-1;
+    else if (i < topk)
+      lo = i+1;
+    else
+      break;
+  }
+
+/*
+  if (cand[lo].key > cand[hi].key)
+    printf("Hmm Error: %d %d %d %f %f\n", i, lo, hi, cand[lo].key, cand[hi].key);
+
+
+  for (i=topk; i<n; i++) {
+    for (j=0; j<topk; j++)
+      if (cand[i].key < cand[j].key)
+        printf("Hmm Error: %d %d %f %f %d %d\n", i, j, cand[i].key, cand[j].key, lo, hi);
+  }
+*/
+
+  return topk;
+}

third_party/METIS/GKlib/fs.c

+/*!
+\file  fs.c
+\brief Various file-system functions.
+
+This file contains various functions that deal with interfacing with 
+the filesystem in a portable way.
+
+\date Started 4/10/95
+\author George
+\version\verbatim $Id: fs.c 14332 2013-05-18 12:22:57Z karypis $ \endverbatim
+*/
+
+
+#include <GKlib.h>
+
+
+
+/*************************************************************************
+* This function checks if a file exists
+**************************************************************************/
+int gk_fexists(char *fname)
+{
+  struct stat status;
+
+  if (stat(fname, &status) == -1)
+    return 0;
+
+  return S_ISREG(status.st_mode);
+}
+
+
+/*************************************************************************
+* This function checks if a directory exists
+**************************************************************************/
+int gk_dexists(char *dirname)
+{
+  struct stat status;
+
+  if (stat(dirname, &status) == -1)
+    return 0;
+
+  return S_ISDIR(status.st_mode);
+}
+
+
+/*************************************************************************/
+/*! \brief Returns the size of the file in bytes
+
+This function returns the size of a file as a 64 bit integer. If there 
+were any errors in stat'ing the file, -1 is returned.
+\note That due to the -1 return code, the maximum file size is limited to
+      63 bits (which I guess is okay for now).
+*/
+/**************************************************************************/
+ssize_t gk_getfsize(char *filename)
+{
+  struct stat status;
+
+  if (stat(filename, &status) == -1) 
+    return -1;
+
+  return (size_t)(status.st_size);
+}
+
+
+/*************************************************************************/
+/*! This function gets some basic statistics about the file. 
+    \param fname is the name of the file
+    \param r_nlines is the number of lines in the file. If it is NULL,
+           this information is not returned.
+    \param r_ntokens is the number of tokens in the file. If it is NULL,
+           this information is not returned.
+    \param r_max_nlntokens is the maximum number of tokens in any line
+           in the file. If it is NULL this information is not returned.
+    \param r_nbytes is the number of bytes in the file. If it is NULL,
+           this information is not returned.
+*/
+/*************************************************************************/
+void gk_getfilestats(char *fname, size_t *r_nlines, size_t *r_ntokens, 
+        size_t *r_max_nlntokens, size_t *r_nbytes)
+{
+  size_t nlines=0, ntokens=0, max_nlntokens=0, nbytes=0, oldntokens=0, nread;
+  int intoken=0;
+  char buffer[2049], *cptr;
+  FILE *fpin;
+
+  fpin = gk_fopen(fname, "r", "gk_GetFileStats");
+
+  while (!feof(fpin)) {
+    nread = fread(buffer, sizeof(char), 2048, fpin);
+    nbytes += nread;
+
+    buffer[nread] = '\0';  /* There is space for this one */
+    for (cptr=buffer; *cptr!='\0'; cptr++) {
+      if (*cptr == '\n') {
+        nlines++;
+        ntokens += intoken;
+        intoken = 0;
+        if (max_nlntokens < ntokens-oldntokens)
+          max_nlntokens = ntokens-oldntokens;
+        oldntokens = ntokens;
+      }
+      else if (*cptr == ' ' || *cptr == '\t') {
+        ntokens += intoken;
+        intoken = 0;
+      }
+      else {
+        intoken = 1;
+      }
+    }
+  }
+  ntokens += intoken;
+  if (max_nlntokens < ntokens-oldntokens)
+    max_nlntokens = ntokens-oldntokens;
+
+  gk_fclose(fpin);
+
+  if (r_nlines != NULL)
+    *r_nlines  = nlines;
+  if (r_ntokens != NULL)
+    *r_ntokens = ntokens;
+  if (r_max_nlntokens != NULL)
+    *r_max_nlntokens = max_nlntokens;
+  if (r_nbytes != NULL)
+    *r_nbytes  = nbytes;
+}
+
+
+/*************************************************************************
+* This function takes in a potentially full path specification of a file
+* and just returns a string containing just the basename of the file.
+* The basename is derived from the actual filename by stripping the last
+* .ext part.
+**************************************************************************/
+char *gk_getbasename(char *path)
+{
+  char *startptr, *endptr;
+  char *basename;
+
+  if ((startptr = strrchr(path, '/')) == NULL) 
+    startptr = path;
+  else 
+    startptr = startptr+1;
+
+  basename = gk_strdup(startptr);
+
+  if ((endptr = strrchr(basename, '.')) != NULL) 
+    *endptr = '\0';
+
+  return basename;
+}
+
+/*************************************************************************
+* This function takes in a potentially full path specification of a file
+* and just returns a string corresponding to its file extension. The
+* extension of a file is considered to be the string right after the 
+* last '.' character.
+**************************************************************************/
+char *gk_getextname(char *path)
+{
+  char *startptr;
+
+  if ((startptr = strrchr(path, '.')) == NULL) 
+    return gk_strdup(path);
+  else 
+    return gk_strdup(startptr+1);
+}
+
+/*************************************************************************
+* This function takes in a potentially full path specification of a file
+* and just returns a string containing just the filename.
+**************************************************************************/
+char *gk_getfilename(char *path)
+{
+  char *startptr;
+
+  if ((startptr = strrchr(path, '/')) == NULL) 
+    return gk_strdup(path);
+  else 
+    return gk_strdup(startptr+1);
+}
+
+/*************************************************************************
+* This function takes in a potentially full path specification of a file
+* and extracts the directory path component if it exists, otherwise it
+* returns "./" as the path. The memory for it is dynamically allocated.
+**************************************************************************/
+char *getpathname(char *path)
+{
+  char *endptr, *tmp;
+
+  if ((endptr = strrchr(path, '/')) == NULL) {
+    return gk_strdup(".");
+  }
+  else  {
+    tmp = gk_strdup(path);
+    *(strrchr(tmp, '/')) = '\0';
+    return tmp;
+  }
+}
+
+
+
+/*************************************************************************
+* This function creates a path
+**************************************************************************/
+int gk_mkpath(char *pathname)
+{
+  char tmp[2048];
+
+  sprintf(tmp, "mkdir -p %s", pathname);
+  return system(tmp);
+}
+
+
+/*************************************************************************
+* This function deletes a directory tree and all of its contents
+**************************************************************************/
+int gk_rmpath(char *pathname)
+{
+  char tmp[2048];
+
+  sprintf(tmp, "rm -r %s", pathname);
+  return system(tmp);
+}

third_party/METIS/GKlib/getopt.c

+/*************************************************************************/
+/*! \file getopt.c
+\brief Command line parsing 
+
+This file contains a implementation of GNU's Getopt facility. The purpose
+for including it here is to ensure portability across different unix- and
+windows-based systems.
+
+\warning 
+The implementation provided here uses the \c gk_ prefix for all variables
+used by the standard Getopt facility to communicate with the program.
+So, do read the documentation here.
+
+\verbatim
+   Copyright (C) 1987,88,89,90,91,92,93,94,95,96,98,99,2000,2001
+   Free Software Foundation, Inc. This file is part of the GNU C Library.
+
+   The GNU C 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.
+
+   The GNU C 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 the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  
+\endverbatim
+*/
+/*************************************************************************/
+
+
+#include <GKlib.h>
+
+/*************************************************************************/
+/* Local function prototypes */
+/*************************************************************************/
+static void exchange (char **);
+static char *gk_getopt_initialize (int, char **, char *);
+static int gk_getopt_internal(int argc, char **argv, char *optstring, 
+        struct gk_option *longopts, int *longind, int long_only);
+
+
+
+/*************************************************************************/
+/*! \brief For communication arguments to the caller.
+
+This variable is set by getopt to point at the value of the option argument, 
+for those options that accept arguments.
+*/
+/*************************************************************************/
+char *gk_optarg;
+
+
+/*************************************************************************/
+/*! \brief Index in ARGV of the next element to be scanned. 
+
+This variable is set by getopt to the index of the next element of the argv 
+array to be processed. Once getopt has found all of the option arguments, 
+you can use this variable to determine where the remaining non-option arguments 
+begin. 
+*/
+/*************************************************************************/
+int gk_optind = 1; 
+
+
+/*************************************************************************/
+/*! \brief Controls error reporting for unrecognized options.  
+
+If the value of this variable is nonzero, then getopt prints an error 
+message to the standard error stream if it encounters an unknown option 
+character or an option with a missing required argument. This is the default 
+behavior. If you set this variable to zero, getopt does not print any messages,
+but it still returns the character ? to indicate an error.
+*/
+/*************************************************************************/
+int gk_opterr = 1;
+
+
+/*************************************************************************/
+/*! \brief Stores unknown option characters
+
+When getopt encounters an unknown option character or an option with a 
+missing required argument, it stores that option character in this 
+variable. You can use this for providing your own diagnostic messages.
+*/
+/*************************************************************************/
+int gk_optopt = '?';
+
+
+/*************************************************************************/
+/*
+Records that the getopt facility has been initialized.
+*/
+/*************************************************************************/
+int gk_getopt_initialized;
+
+
+/*************************************************************************/
+/*
+The next char to be scanned in the option-element in which the last option 
+character we returned was found.  This allows us to pick up the scan where 
+we left off.
+
+If this is zero, or a null string, it means resume the scan by advancing 
+to the next ARGV-element.  
+*/
+/*************************************************************************/
+static char *nextchar;
+
+
+/*************************************************************************/
+/*
+Value of POSIXLY_CORRECT environment variable.  
+*/
+/*************************************************************************/
+static char *posixly_correct;
+
+
+/*************************************************************************/
+/*
+Describe how to deal with options that follow non-option ARGV-elements.
+
+If the caller did not specify anything, the default is REQUIRE_ORDER if 
+the environment variable POSIXLY_CORRECT is defined, PERMUTE otherwise.
+
+REQUIRE_ORDER means don't recognize them as options; stop option processing 
+when the first non-option is seen.  This is what Unix does.  This mode of 
+operation is selected by either setting the environment variable 
+POSIXLY_CORRECT, or using `+' as the first character of the list of 
+option characters.
+
+PERMUTE is the default.  We permute the contents of ARGV as we scan, so 
+that eventually all the non-options are at the end.  This allows options
+to be given in any order, even with programs that were not written to
+expect this.
+
+RETURN_IN_ORDER is an option available to programs that were written
+to expect options and other ARGV-elements in any order and that care 
+about the ordering of the two.  We describe each non-option ARGV-element
+as if it were the argument of an option with character code 1.
+Using `-' as the first character of the list of option characters
+selects this mode of operation.
+
+The special argument `--' forces an end of option-scanning regardless
+of the value of `ordering'.  In the case of RETURN_IN_ORDER, only
+`--' can cause `getopt' to return -1 with `gk_optind' != ARGC.  
+*/
+/*************************************************************************/
+static enum
+{
+  REQUIRE_ORDER, PERMUTE, RETURN_IN_ORDER
+} ordering;
+
+
+
+/*************************************************************************/
+/* 
+Describe the part of ARGV that contains non-options that have
+been skipped.  `first_nonopt' is the index in ARGV of the first of them;
+`last_nonopt' is the index after the last of them.  
+*/
+/*************************************************************************/
+static int first_nonopt;
+static int last_nonopt;
+
+
+
+
+
+/*************************************************************************/
+/*
+Handle permutation of arguments.  
+
+Exchange two adjacent subsequences of ARGV. 
+One subsequence is elements [first_nonopt,last_nonopt)
+which contains all the non-options that have been skipped so far.
+The other is elements [last_nonopt,gk_optind), which contains all
+the options processed since those non-options were skipped.
+
+`first_nonopt' and `last_nonopt' are relocated so that they describe
+the new indices of the non-options in ARGV after they are moved.  
+*/
+/*************************************************************************/
+static void exchange (char **argv)
+{
+  int bottom = first_nonopt;
+  int middle = last_nonopt;
+  int top = gk_optind;
+  char *tem;
+
+  /* Exchange the shorter segment with the far end of the longer segment.
+     That puts the shorter segment into the right place.
+     It leaves the longer segment in the right place overall,
+     but it consists of two parts that need to be swapped next.  */
+
+  while (top > middle && middle > bottom) {
+    if (top - middle > middle - bottom) {
+      /* Bottom segment is the short one.  */
+      int len = middle - bottom;
+      register int i;
+
+      /* Swap it with the top part of the top segment.  */
+      for (i = 0; i < len; i++) {
+	tem = argv[bottom + i];
+	argv[bottom + i] = argv[top - (middle - bottom) + i];
+	argv[top - (middle - bottom) + i] = tem;
+      }
+      /* Exclude the moved bottom segment from further swapping.  */
+      top -= len;
+    }
+    else {
+      /* Top segment is the short one.  */
+      int len = top - middle;
+      register int i;
+
+      /* Swap it with the bottom part of the bottom segment.  */
+      for (i = 0; i < len; i++) {
+        tem = argv[bottom + i];
+        argv[bottom + i] = argv[middle + i];
+        argv[middle + i] = tem;
+      }
+      /* Exclude the moved top segment from further swapping.  */
+      bottom += len;
+    }
+  }
+
+  /* Update records for the slots the non-options now occupy.  */
+
+  first_nonopt += (gk_optind - last_nonopt);
+  last_nonopt = gk_optind;
+}
+
+
+
+/*************************************************************************/
+/*
+Initialize the internal data when the first call is made.  
+*/
+/*************************************************************************/
+static char *gk_getopt_initialize (int argc, char **argv, char *optstring)
+{
+  /* Start processing options with ARGV-element 1 (since ARGV-element 0
+     is the program name); the sequence of previously skipped
+     non-option ARGV-elements is empty.  */
+
+  first_nonopt = last_nonopt = gk_optind;
+
+  nextchar = NULL;
+
+  posixly_correct = getenv("POSIXLY_CORRECT");
+
+  /* Determine how to handle the ordering of options and nonoptions.  */
+  if (optstring[0] == '-') {
+    ordering = RETURN_IN_ORDER;
+    ++optstring;
+  }
+  else if (optstring[0] == '+') {
+    ordering = REQUIRE_ORDER;
+    ++optstring;
+  }
+  else if (posixly_correct != NULL)
+    ordering = REQUIRE_ORDER;
+  else
+    ordering = PERMUTE;
+
+  return optstring;
+}
+
+
+/*************************************************************************/
+/*
+   Scan elements of ARGV (whose length is ARGC) for option characters
+   given in OPTSTRING.
+
+   If an element of ARGV starts with '-', and is not exactly "-" or "--",
+   then it is an option element.  The characters of this element
+   (aside from the initial '-') are option characters.  If `getopt'
+   is called repeatedly, it returns successively each of the option characters
+   from each of the option elements.
+
+   If `getopt' finds another option character, it returns that character,
+   updating `gk_optind' and `nextchar' so that the next call to `getopt' can
+   resume the scan with the following option character or ARGV-element.
+
+   If there are no more option characters, `getopt' returns -1.
+   Then `gk_optind' is the index in ARGV of the first ARGV-element
+   that is not an option.  (The ARGV-elements have been permuted
+   so that those that are not options now come last.)
+
+   OPTSTRING is a string containing the legitimate option characters.
+   If an option character is seen that is not listed in OPTSTRING,
+   return '?' after printing an error message.  If you set `gk_opterr' to
+   zero, the error message is suppressed but we still return '?'.
+
+   If a char in OPTSTRING is followed by a colon, that means it wants an arg,
+   so the following text in the same ARGV-element, or the text of the following
+   ARGV-element, is returned in `gk_optarg'.  Two colons mean an option that
+   wants an optional arg; if there is text in the current ARGV-element,
+   it is returned in `gk_optarg', otherwise `gk_optarg' is set to zero.
+
+   If OPTSTRING starts with `-' or `+', it requests different methods of
+   handling the non-option ARGV-elements.
+   See the comments about RETURN_IN_ORDER and REQUIRE_ORDER, above.
+
+   Long-named options begin with `--' instead of `-'.
+   Their names may be abbreviated as long as the abbreviation is unique
+   or is an exact match for some defined option.  If they have an
+   argument, it follows the option name in the same ARGV-element, separated
+   from the option name by a `=', or else the in next ARGV-element.
+   When `getopt' finds a long-named option, it returns 0 if that option's
+   `flag' field is nonzero, the value of the option's `val' field
+   if the `flag' field is zero.
+
+   LONGOPTS is a vector of `struct gk_option' terminated by an
+   element containing a name which is zero.
+
+   LONGIND returns the index in LONGOPT of the long-named option found.
+   It is only valid when a long-named option has been found by the most
+   recent call.
+
+   If LONG_ONLY is nonzero, '-' as well as '--' can introduce
+   long-named options.  
+*/
+/*************************************************************************/
+static int gk_getopt_internal(int argc, char **argv, char *optstring, 
+        struct gk_option *longopts, int *longind, int long_only)
+{
+  int print_errors = gk_opterr;
+  if (optstring[0] == ':')
+    print_errors = 0;
+
+  if (argc < 1)
+    return -1;
+
+  gk_optarg = NULL;
+
+  if (gk_optind == 0 || !gk_getopt_initialized) {
+    if (gk_optind == 0)
+      gk_optind = 1;	/* Don't scan ARGV[0], the program name.  */
+
+    optstring = gk_getopt_initialize (argc, argv, optstring);
+    gk_getopt_initialized = 1;
+  }
+
+  /* Test whether ARGV[gk_optind] points to a non-option argument.
+     Either it does not have option syntax, or there is an environment flag
+     from the shell indicating it is not an option.  The later information
+     is only used when the used in the GNU libc.  */
+# define NONOPTION_P (argv[gk_optind][0] != '-' || argv[gk_optind][1] == '\0')
+
+  if (nextchar == NULL || *nextchar == '\0') {
+    /* Advance to the next ARGV-element.  */
+
+    /* Give FIRST_NONOPT & LAST_NONOPT rational values if OPTIND has been
+       moved back by the user (who may also have changed the arguments).  */
+    if (last_nonopt > gk_optind)
+      last_nonopt = gk_optind;
+    if (first_nonopt > gk_optind)
+      first_nonopt = gk_optind;
+
+    if (ordering == PERMUTE) {
+      /* If we have just processed some options following some non-options,
+	 exchange them so that the options come first.  */
+
+      if (first_nonopt != last_nonopt && last_nonopt != gk_optind)
+	exchange ((char **) argv);
+      else if (last_nonopt != gk_optind)
+	first_nonopt = gk_optind;
+
+      /* Skip any additional non-options
+	 and extend the range of non-options previously skipped.  */
+
+      while (gk_optind < argc && NONOPTION_P)
+        gk_optind++;
+
+      last_nonopt = gk_optind;
+    }
+
+    /* The special ARGV-element `--' means premature end of options.
+       Skip it like a null option,
+       then exchange with previous non-options as if it were an option,
+       then skip everything else like a non-option.  */
+
+    if (gk_optind != argc && !strcmp (argv[gk_optind], "--")) {
+      gk_optind++;
+
+      if (first_nonopt != last_nonopt && last_nonopt != gk_optind)
+        exchange ((char **) argv);
+      else if (first_nonopt == last_nonopt)
+        first_nonopt = gk_optind;
+      last_nonopt = argc;
+
+      gk_optind = argc;
+    }
+
+    /* If we have done all the ARGV-elements, stop the scan
+       and back over any non-options that we skipped and permuted.  */
+
+    if (gk_optind == argc) {
+      /* Set the next-arg-index to point at the non-options
+	 that we previously skipped, so the caller will digest them.  */
+      if (first_nonopt != last_nonopt)
+	gk_optind = first_nonopt;
+      return -1;
+    }
+
+    /* If we have come to a non-option and did not permute it,
+       either stop the scan or describe it to the caller and pass it by.  */
+
+    if (NONOPTION_P) {
+      if (ordering == REQUIRE_ORDER)
+	return -1;
+      gk_optarg = argv[gk_optind++];
+      return 1;
+    }
+
+    /* We have found another option-ARGV-element.
+       Skip the initial punctuation.  */
+
+    nextchar = (argv[gk_optind] + 1 + (longopts != NULL && argv[gk_optind][1] == '-'));
+  }
+
+  /* Decode the current option-ARGV-element.  */
+
+  /* Check whether the ARGV-element is a long option.
+
+     If long_only and the ARGV-element has the form "-f", where f is
+     a valid short option, don't consider it an abbreviated form of
+     a long option that starts with f.  Otherwise there would be no
+     way to give the -f short option.
+
+     On the other hand, if there's a long option "fubar" and
+     the ARGV-element is "-fu", do consider that an abbreviation of
+     the long option, just like "--fu", and not "-f" with arg "u".
+
+     This distinction seems to be the most useful approach.  */
+
+  if (longopts != NULL && (argv[gk_optind][1] == '-' || (long_only && (argv[gk_optind][2] || !strchr(optstring, argv[gk_optind][1]))))) {
+    char *nameend;
+    struct gk_option *p;
+    struct gk_option *pfound = NULL;
+    int exact = 0;
+    int ambig = 0;
+    int indfound = -1;
+    int option_index;
+
+    for (nameend = nextchar; *nameend && *nameend != '='; nameend++)
+      /* Do nothing.  */ ;
+
+    /* Test all long options for either exact match or abbreviated matches.  */
+    for (p = longopts, option_index = 0; p->name; p++, option_index++) {
+      if (!strncmp (p->name, nextchar, nameend - nextchar)) {
+        if ((unsigned int) (nameend - nextchar) == (unsigned int) strlen (p->name)) {
+	  /* Exact match found.  */
+	  pfound = p;
+	  indfound = option_index;
+	  exact = 1;
+	  break;
+	}
+	else if (pfound == NULL) {
+          /* First nonexact match found.  */
+	  pfound = p;
+	  indfound = option_index;
+	}
+	else if (long_only || pfound->has_arg != p->has_arg || pfound->flag != p->flag || pfound->val != p->val)
+	  /* Second or later nonexact match found.  */
+	  ambig = 1;
+      }
+    }
+
+    if (ambig && !exact) {
+      if (print_errors)
+        fprintf(stderr, "%s: option `%s' is ambiguous\n", argv[0], argv[gk_optind]);
+
+      nextchar += strlen (nextchar);
+      gk_optind++;
+      gk_optopt = 0;
+      return '?';
+    }
+
+    if (pfound != NULL) {
+      option_index = indfound;
+      gk_optind++;
+      if (*nameend) {
+	/* Don't test has_arg with >, because some C compilers don't allow it to be used on enums.  */
+	if (pfound->has_arg)
+	  gk_optarg = nameend + 1;
+	else {
+	  if (print_errors) {
+	    if (argv[gk_optind - 1][1] == '-')
+	      /* --option */
+	      fprintf(stderr, "%s: option `--%s' doesn't allow an argument\n", argv[0], pfound->name);
+	    else
+	      /* +option or -option */
+	      fprintf(stderr, "%s: option `%c%s' doesn't allow an argument\n", argv[0], argv[gk_optind - 1][0], pfound->name);
+	  }
+
+	  nextchar += strlen (nextchar);
+
+	  gk_optopt = pfound->val;
+	  return '?';
+	}
+      }
+      else if (pfound->has_arg == 1) {
+	if (gk_optind < argc)
+	  gk_optarg = argv[gk_optind++];
+	else {
+	  if (print_errors)
+	    fprintf(stderr, "%s: option `%s' requires an argument\n", argv[0], argv[gk_optind - 1]);
+	  nextchar += strlen (nextchar);
+	  gk_optopt = pfound->val;
+	  return optstring[0] == ':' ? ':' : '?';
+	}
+      }
+      nextchar += strlen (nextchar);
+      if (longind != NULL)
+        *longind = option_index;
+      if (pfound->flag) {
+	*(pfound->flag) = pfound->val;
+	return 0;
+      }
+      return pfound->val;
+    }
+
+    /* Can't find it as a long option.  If this is not getopt_long_only,
+       or the option starts with '--' or is not a valid short
+        option, then it's an error. Otherwise interpret it as a short option.  */
+    if (!long_only || argv[gk_optind][1] == '-' || strchr(optstring, *nextchar) == NULL) {
+      if (print_errors) {
+	if (argv[gk_optind][1] == '-')
+	  /* --option */
+	  fprintf(stderr, "%s: unrecognized option `--%s'\n", argv[0], nextchar);
+	else
+	  /* +option or -option */
+	  fprintf(stderr, "%s: unrecognized option `%c%s'\n", argv[0], argv[gk_optind][0], nextchar);
+      }
+      nextchar = (char *) "";
+      gk_optind++;
+      gk_optopt = 0;
+      return '?';
+    }
+  }
+
+  /* Look at and handle the next short option-character.  */
+  {
+    char c = *nextchar++;
+    char *temp = strchr(optstring, c);
+
+    /* Increment `gk_optind' when we start to process its last character.  */
+    if (*nextchar == '\0')
+      ++gk_optind;
+
+    if (temp == NULL || c == ':') {
+      if (print_errors) {
+        if (posixly_correct)
+	  /* 1003.2 specifies the format of this message.  */
+	  fprintf(stderr, "%s: illegal option -- %c\n", argv[0], c);
+	else
+	  fprintf(stderr, "%s: invalid option -- %c\n", argv[0], c);
+      }
+      gk_optopt = c;
+      return '?';
+    }
+
+    /* Convenience. Treat POSIX -W foo same as long option --foo */
+    if (temp[0] == 'W' && temp[1] == ';') {
+      char *nameend;
+      struct gk_option *p;
+      struct gk_option *pfound = NULL;
+      int exact = 0;
+      int ambig = 0;
+      int indfound = 0;
+      int option_index;
+
+      /* This is an option that requires an argument.  */
+      if (*nextchar != '\0') {
+	gk_optarg = nextchar;
+	/* If we end this ARGV-element by taking the rest as an arg,
+	   we must advance to the next element now.  */
+	gk_optind++;
+      }
+      else if (gk_optind == argc) {
+	if (print_errors) {
+	  /* 1003.2 specifies the format of this message.  */
+	  fprintf(stderr, "%s: option requires an argument -- %c\n", argv[0], c);
+	}
+	gk_optopt = c;
+	if (optstring[0] == ':')
+	  c = ':';
+	else
+	  c = '?';
+	return c;
+      }
+      else
+	/* We already incremented `gk_optind' once; increment it again when taking next ARGV-elt as argument.  */
+	gk_optarg = argv[gk_optind++];
+
+      /* gk_optarg is now the argument, see if it's in the table of longopts.  */
+
+      for (nextchar = nameend = gk_optarg; *nameend && *nameend != '='; nameend++)
+	/* Do nothing.  */ ;
+
+      /* Test all long options for either exact match or abbreviated matches.  */
+      for (p = longopts, option_index = 0; p->name; p++, option_index++) {
+	if (!strncmp (p->name, nextchar, nameend - nextchar)) {
+	  if ((unsigned int) (nameend - nextchar) == strlen (p->name)) {
+	    /* Exact match found.  */
+	    pfound = p;
+	    indfound = option_index;
+	    exact = 1;
+	    break;
+	  }
+	  else if (pfound == NULL) {
+	    /* First nonexact match found.  */
+	    pfound = p;
+	    indfound = option_index;
+	  }
+	  else
+	    /* Second or later nonexact match found.  */
+	    ambig = 1;
+	}
+      }
+      if (ambig && !exact) {
+	if (print_errors)
+	  fprintf(stderr, "%s: option `-W %s' is ambiguous\n", argv[0], argv[gk_optind]);
+	nextchar += strlen (nextchar);
+	gk_optind++;
+	return '?';
+      }
+      if (pfound != NULL) {
+	option_index = indfound;
+	if (*nameend) {
+	  /* Don't test has_arg with >, because some C compilers don't allow it to be used on enums.  */
+	  if (pfound->has_arg)
+	    gk_optarg = nameend + 1;
+	  else {
+	    if (print_errors)
+	      fprintf(stderr, "%s: option `-W %s' doesn't allow an argument\n", argv[0], pfound->name);
+
+	    nextchar += strlen (nextchar);
+	    return '?';
+	  }
+	}
+	else if (pfound->has_arg == 1) {
+	  if (gk_optind < argc)
+	    gk_optarg = argv[gk_optind++];
+	  else {
+	    if (print_errors)
+	      fprintf(stderr, "%s: option `%s' requires an argument\n", argv[0], argv[gk_optind - 1]);
+	    nextchar += strlen (nextchar);
+	    return optstring[0] == ':' ? ':' : '?';
+	  }
+        }
+	nextchar += strlen (nextchar);
+	if (longind != NULL)
+	  *longind = option_index;
+	if (pfound->flag) {
+	  *(pfound->flag) = pfound->val;
+	  return 0;
+	}
+	return pfound->val;
+      }
+      nextchar = NULL;
+      return 'W';	/* Let the application handle it.   */
+    }
+
+    if (temp[1] == ':') {
+      if (temp[2] == ':') {
+	/* This is an option that accepts an argument optionally.  */
+	if (*nextchar != '\0') {
+  	  gk_optarg = nextchar;
+	  gk_optind++;
+	}
+	else
+	  gk_optarg = NULL;
+	nextchar = NULL;
+      }
+      else {
+	/* This is an option that requires an argument.  */
+	if (*nextchar != '\0') {
+	  gk_optarg = nextchar;
+	  /* If we end this ARGV-element by taking the rest as an arg, we must advance to the next element now.  */
+	  gk_optind++;
+	}
+	else if (gk_optind == argc) {
+	  if (print_errors) {
+	    /* 1003.2 specifies the format of this message.  */
+	    fprintf(stderr, "%s: option requires an argument -- %c\n", argv[0], c);
+	  }
+	  gk_optopt = c;
+	  if (optstring[0] == ':')
+	    c = ':';
+	  else
+	    c = '?';
+	}
+	else
+	  /* We already incremented `gk_optind' once; increment it again when taking next ARGV-elt as argument.  */
+	  gk_optarg = argv[gk_optind++];
+	nextchar = NULL;
+      }
+    }
+    return c;
+  }
+}
+
+
+
+/*************************************************************************/
+/*! \brief Parse command-line arguments
+
+The gk_getopt() function gets the next option argument from the argument 
+list specified by the \c argv and \c argc arguments. Normally these values 
+come directly from the arguments received by main().
+
+\param argc is the number of command line arguments passed to main().
+\param argv is an array of strings storing the above command line 
+       arguments.
+\param options is a string that specifies the option characters that 
+       are valid for this program. An option character in this string 
+       can be followed by a colon (`:') to indicate that it takes a 
+       required argument. If an option character is followed by two 
+       colons (`::'), its argument is optional; this is a GNU extension.
+
+\return  
+It returns the option character for the next command line option. When no 
+more option arguments are available, it returns -1. There may still be 
+more non-option arguments; you must compare the external variable 
+#gk_optind against the \c argc parameter to check this.
+
+\return  
+If the option has an argument, gk_getopt() returns the argument by storing 
+it in the variable #gk_optarg. You don't ordinarily need to copy the 
+#gk_optarg string, since it is a pointer into the original \c argv array, 
+not into a static area that might be overwritten.
+
+\return  
+If gk_getopt() finds an option character in \c argv that was not included 
+in options, or a missing option argument, it returns `?' and sets the 
+external variable #gk_optopt to the actual option character. 
+If the first character of options is a colon (`:'), then gk_getopt() 
+returns `:' instead of `?' to indicate a missing option argument. 
+In addition, if the external variable #gk_opterr is nonzero (which is 
+the default), gk_getopt() prints an error message.  This variable is 
+set by gk_getopt() to point at the value of the option argument, 
+for those options that accept arguments.
+
+
+gk_getopt() has three ways to deal with options that follow non-options 
+\c argv elements. The special argument <tt>`--'</tt> forces in all cases 
+the end of option scanning.
+  - The default is to permute the contents of \c argv while scanning it 
+    so that eventually all the non-options are at the end. This allows 
+    options to be given in any order, even with programs that were not 
+    written to expect this.
+  - If the options argument string begins with a hyphen (`-'), this is 
+    treated specially. It permits arguments that are not options to be 
+    returned as if they were associated with option character `\\1'.
+  - POSIX demands the following behavior: The first non-option stops 
+    option processing. This mode is selected by either setting the 
+    environment variable POSIXLY_CORRECT or beginning the options
+    argument string with a plus sign (`+'). 
+
+*/
+/*************************************************************************/
+int gk_getopt(int argc, char **argv, char *options)
+{
+  return gk_getopt_internal(argc, argv, options, NULL, NULL, 0);
+}
+
+
+/*************************************************************************/
+/*! \brief Parse command-line arguments with long options
+
+This function accepts GNU-style long options as well as single-character 
+options. 
+
+\param argc is the number of command line arguments passed to main().
+\param argv is an array of strings storing the above command line 
+       arguments.
+\param options describes the short options to accept, just as it does 
+       in gk_getopt(). 
+\param long_options describes the long options to accept. See the 
+       defintion of ::gk_option for more information.
+\param opt_index this is a returned variable.  For any long option, 
+       gk_getopt_long() tells you the index in the array \c long_options 
+       of the options definition, by storing it into <tt>*opt_index</tt>. 
+       You can get the name of the option with <tt>longopts[*opt_index].name</tt>. 
+       So you can distinguish among long options either by the values 
+       in their val fields or by their indices. You can also distinguish 
+       in this way among long options that set flags.
+
+
+\return
+When gk_getopt_long() encounters a short option, it does the same thing 
+that gk_getopt() would do: it returns the character code for the option, 
+and stores the options argument (if it has one) in #gk_optarg.
+
+\return
+When gk_getopt_long() encounters a long option, it takes actions based 
+on the flag and val fields of the definition of that option.
+
+\return
+If flag is a null pointer, then gk_getopt_long() returns the contents 
+of val to indicate which option it found. You should arrange distinct 
+values in the val field for options with different meanings, so you 
+can decode these values after gk_getopt_long() returns. If the long 
+option is equivalent to a short option, you can use the short option's 
+character code in val.
+
+\return
+If flag is not a null pointer, that means this option should just set 
+a flag in the program. The flag is a variable of type int that you 
+define. Put the address of the flag in the flag field. Put in the 
+val field the value you would like this option to store in the flag. 
+In this case, gk_getopt_long() returns 0.
+
+\return
+When a long option has an argument, gk_getopt_long() puts the argument 
+value in the variable #gk_optarg before returning. When the option has 
+no argument, the value in #gk_optarg is a null pointer. This is
+how you can tell whether an optional argument was supplied.
+
+\return
+When gk_getopt_long() has no more options to handle, it returns -1, 
+and leaves in the variable #gk_optind the index in argv of the next 
+remaining argument. 
+*/
+/*************************************************************************/
+int gk_getopt_long( int argc, char **argv, char *options, 
+       struct gk_option *long_options, int *opt_index)
+{
+  return gk_getopt_internal (argc, argv, options, long_options, opt_index, 0);
+}
+
+
+
+/*************************************************************************/
+/*! \brief Parse command-line arguments with only long options
+
+Like gk_getopt_long(), but '-' as well as '--' can indicate a long option.
+If an option that starts with '-' (not '--') doesn't match a long option,
+but does match a short option, it is parsed as a short option instead.  
+*/
+/*************************************************************************/
+int gk_getopt_long_only(int argc, char **argv, char *options, 
+       struct gk_option *long_options, int *opt_index)
+{
+  return gk_getopt_internal(argc, argv, options, long_options, opt_index, 1);
+}
+

third_party/METIS/GKlib/gk_arch.h

+/*!
+\file gk_arch.h
+\brief This file contains various architecture-specific declerations
+
+\date   Started 3/27/2007
+\author George
+\version\verbatim $Id: gk_arch.h 21637 2018-01-03 22:37:24Z karypis $ \endverbatim
+*/
+
+#ifndef _GK_ARCH_H_
+#define _GK_ARCH_H_
+
+/*************************************************************************
+* Architecture-specific differences in header files
+**************************************************************************/
+#ifdef LINUX
+#if !defined(__USE_XOPEN)
+#define __USE_XOPEN
+#endif
+#if !defined(_XOPEN_SOURCE)
+#define _XOPEN_SOURCE 600
+#endif
+#if !defined(__USE_XOPEN2K)
+#define __USE_XOPEN2K
+#endif
+#endif
+
+
+#ifdef HAVE_EXECINFO_H
+#include <execinfo.h>
+#endif
+
+
+#ifdef __MSC__ 
+  #include "ms_stdint.h"
+  #include "ms_inttypes.h"
+  #include "ms_stat.h"
+#else
+#ifndef SUNOS
+  #include <stdint.h>
+#endif
+  #include <inttypes.h>
+  #include <sys/types.h>
+  #include <sys/resource.h>
+  #include <sys/time.h>
+#endif
+
+
+/*************************************************************************
+* Architecture-specific modifications
+**************************************************************************/
+#ifdef WIN32
+typedef ptrdiff_t ssize_t;
+#endif
+
+
+#ifdef SUNOS
+#define PTRDIFF_MAX  INT64_MAX
+#endif
+
+/* MSC does not have INFINITY defined */
+#ifndef INFINITY
+#define INFINITY FLT_MAX
+#endif
+
+#endif

third_party/METIS/GKlib/gk_defs.h

+/*!
+\file gk_defs.h
+\brief This file contains various constants definitions
+
+\date   Started 3/27/2007
+\author George
+\version\verbatim $Id: gk_defs.h 21012 2017-05-19 03:05:23Z karypis $ \endverbatim
+*/
+
+#ifndef _GK_DEFS_H_
+#define _GK_DEFS_H_
+
+
+#define LTERM                   (void **) 0     /* List terminator for GKfree() */
+
+/* mopt_t types */
+#define GK_MOPT_MARK            1
+#define GK_MOPT_CORE            2
+#define GK_MOPT_HEAP            3
+
+#define HTABLE_EMPTY            -1
+#define HTABLE_DELETED          -2
+#define HTABLE_FIRST             1
+#define HTABLE_NEXT              2
+
+/* pdb corruption bit switches */
+#define CRP_ALTLOCS    1
+#define CRP_MISSINGCA  2
+#define CRP_MISSINGBB  4
+#define CRP_MULTICHAIN 8
+#define CRP_MULTICA    16
+#define CRP_MULTIBB    32
+
+#define MAXLINELEN 300000
+
+/* GKlib signals to standard signal mapping */
+#define SIGMEM  SIGABRT
+#define SIGERR  SIGTERM
+
+
+/* CSR-related defines */
+#define GK_CSR_ROW      1
+#define GK_CSR_COL      2
+#define GK_CSR_ROWCOL   3
+
+#define GK_CSR_MAXTF    1
+#define GK_CSR_SQRT     2
+#define GK_CSR_POW25    3
+#define GK_CSR_POW65    4
+#define GK_CSR_POW75    5
+#define GK_CSR_POW85    6
+#define GK_CSR_LOG      7
+#define GK_CSR_IDF      8
+#define GK_CSR_IDF2     9
+#define GK_CSR_MAXTF2   10
+
+#define GK_CSR_DOTP     1
+#define GK_CSR_COS      2
+#define GK_CSR_JAC      3
+#define GK_CSR_MIN      4
+#define GK_CSR_AMIN     5
+
+#define GK_CSR_FMT_AUTO         2
+#define GK_CSR_FMT_CLUTO        1
+#define GK_CSR_FMT_CSR          2
+#define GK_CSR_FMT_METIS        3
+#define GK_CSR_FMT_BINROW       4
+#define GK_CSR_FMT_BINCOL       5
+#define GK_CSR_FMT_IJV          6
+#define GK_CSR_FMT_BIJV         7
+
+#define GK_CSR_SYM_SUM          1
+#define GK_CSR_SYM_MIN          2
+#define GK_CSR_SYM_MAX          3
+#define GK_CSR_SYM_AVG          4
+
+
+#define GK_GRAPH_FMT_METIS      1
+#define GK_GRAPH_FMT_IJV        2
+
+#endif

third_party/METIS/GKlib/gk_externs.h

+/*!
+\file gk_externs.h
+\brief This file contains definitions of external variables created by GKlib
+
+\date   Started 3/27/2007
+\author George
+\version\verbatim $Id: gk_externs.h 10711 2011-08-31 22:23:04Z karypis $ \endverbatim
+*/
+
+#ifndef _GK_EXTERNS_H_
+#define _GK_EXTERNS_H_
+
+
+/*************************************************************************
+* Extern variable definition. Hopefully, the __thread makes them thread-safe.
+**************************************************************************/
+#ifndef _GK_ERROR_C_
+/* declared in error.c */
+extern __thread int gk_cur_jbufs;
+extern __thread jmp_buf gk_jbufs[];
+extern __thread jmp_buf gk_jbuf;
+
+#endif
+
+#endif

third_party/METIS/GKlib/gk_getopt.h

+/*!
+\file gk_getopt.h
+\brief This file contains GNU's externs/structs/prototypes
+
+\date   Started 3/27/2007
+\author George
+\version\verbatim $Id: gk_getopt.h 10711 2011-08-31 22:23:04Z karypis $ \endverbatim
+*/
+
+#ifndef _GK_GETOPT_H_
+#define _GK_GETOPT_H_
+
+
+/* Externals from getopt.c */
+extern char *gk_optarg;
+extern int gk_optind;
+extern int gk_opterr;
+extern int gk_optopt;
+
+
+/*! \brief The structure that stores the information about the command-line options 
+
+This structure describes a single long option name for the sake of 
+gk_getopt_long(). The argument <tt>long_options</tt> must be an array 
+of these structures, one for each long option. Terminate the array with 
+an element containing all zeros.
+*/
+struct gk_option {
+  char *name;       /*!< This field is the name of the option. */
+  int has_arg;      /*!< This field says whether the option takes an argument. 
+                         It is an integer, and there are three legitimate values: 
+                         no_argument, required_argument and optional_argument. 
+                         */
+  int *flag;        /*!< See the discussion on ::gk_option#val */
+  int val;          /*!< These fields control how to report or act on the option 
+                         when it occurs. 
+                         
+                         If flag is a null pointer, then the val is a value which 
+                         identifies this option. Often these values are chosen 
+                         to uniquely identify particular long options.
+
+                         If flag is not a null pointer, it should be the address 
+                         of an int variable which is the flag for this option. 
+                         The value in val is the value to store in the flag to 
+                         indicate that the option was seen. */
+};
+
+/* Names for the values of the `has_arg' field of `struct gk_option'.  */
+#define no_argument		0
+#define required_argument	1
+#define optional_argument	2
+
+
+/* Function prototypes */
+extern int gk_getopt(int __argc, char **__argv, char *__shortopts);
+extern int gk_getopt_long(int __argc, char **__argv, char *__shortopts,
+              struct gk_option *__longopts, int *__longind);
+extern int gk_getopt_long_only (int __argc, char **__argv,
+              char *__shortopts, struct gk_option *__longopts, int *__longind);
+
+
+
+#endif
+