--- SRC/MatOps.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/MatOps.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,11 +1,5 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h> 
-#include "globheads.h"
 #include "protos.h"
 
-
 int diag_scal( vbsptr vbmat ){
 /*----------------------------------------------------------------------------
  * Diagonal scaling:
--- SRC/PQ.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/PQ.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,8 +1,3 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include "globheads.h"
 #include "protos.h"
 
 #define  ALPHA  0.00001
@@ -126,217 +121,7 @@
 |-----end-of-indsetPQ--------------------------------------------------
 |--------------------------------------------------------------------*/
 
-int add2is(int *last, int nod, int *iord, int *riord)
-{
-/*----------------------------------------------------------------------
-|   adds element nod to independent set
-|---------------------------------------------------------------------*/
-   (*last)++;
-   iord[nod] = *last;
-   riord[*last] = nod;
-   return 0;
-}
-/*---------------------------------------------------------------------
-|---- end of add2is ---------------------------------------------------
-|--------------------------------------------------------------------*/
-int add2com(int *nback, int nod, int *iord, int *riord) 
-{
-/*----------------------------------------------------------------------
-|   adds element nod to independent set
-|---------------------------------------------------------------------*/
-  iord[nod] = *nback;
-  riord[*nback] = nod;
-  (*nback)--;
-  return 0;
-}
-/*---------------------------------------------------------------------
-|---- end of add2com --------------------------------------------------
-|--------------------------------------------------------------------*/
-int indsetC(csptr mat, int bsize, int *iord, int *nnod, double tol) 
-{
-/*--------------------------------------------------------------------- 
-| greedy algorithm for independent set ordering -- 
-|----------------------------------------------------------------------
-|     Input parameters:
-|     -----------------
-|     (mat)  =  matrix in SpaFmt format
-|     
-|     bsize  =  integer (input) the target size of each block.
-|               each block is of size >= bsize. 
-|
-|     w      =  weight factors for the selection of the elements in the
-|               independent set. If w(i) is small i will be left for the
-|               vertex cover set. 
-|
-|     tol    =  a tolerance for excluding a row from independent set.
-|
-|     Output parameters:
-|     ------------------ 
-|     iord   = permutation array corresponding to the independent set 
-|     ordering.  Row number i will become row number iord[i] in 
-|     permuted matrix.
-|     
-|     nnod   = (output) number of elements in the independent set. 
-|     
-|----------------------------------------------------------------------- 
-|     the algorithm searches nodes in lexicographic order and groups
-|     the (BSIZE-1) nearest nodes of the current to form a block of
-|     size BSIZE. The current algorithm does not use values of the matrix.
-|---------------------------------------------------------------------*/ 
-/*   local variables   */
-   int nod, jcount, lastlev, begin, last0, last, nback, mid,
-     j1, j2, jcol, inod, jnod, j, k, jcount0, begin0, *rowj;
-   int prog, n=mat->n, *riord;
-   double *w;
-   csptr matT,gmat;  
-
-/*-----------------------------------------------------------------------*/
-   riord = (int *) Malloc(n*sizeof(int), "indsetC:1" );
-   w     = (double *) Malloc(n*sizeof(double), "indsetC:2" );
-   matT  = (csptr) Malloc(sizeof(SparMat), "indsetC:3" );
-/*  	 call weights to compute the weights for  input matrix.. */
-   setupCS(matT, mat->n,1);
-   SparTran(mat, matT, 1, 0);
-   SparTran(matT, mat, 1, 1); 
-   weightsC(mat, w); 
-/*---------------------------------------------------------------------- 
-| scan all nodes first to eliminate those not satisfying DD criterion 
-+----------------------------------------------------------------------*/
-   nback = n-1; 
-   nod = 0;
-   for(j=0; j<n; j++)
-     iord[j] = -1; 
-   for(j=0; j<n; j++) {
-     if (w[j] < tol) {
-       add2com(&nback, j, iord, riord);
-       nod++;
-     }
-   }
-   last = -1;
-   for (nod=0; nod<n; nod++) {	
-     while (iord[nod] != -1)   {
-       if (++nod >= mat->n) goto label50;
-     }
-/*-------------------- initialize level-set - contains nod (only)*/
-     add2is(&last, nod, iord, riord);
-     begin   = last;
-     begin0  = begin; 
-     lastlev = begin;
-     jcount  = 1;
-/*----------------------------------------------------------------------
-|     put all the nearest neighbor nodes of the current node into
-|     the block until the number is BSIZE.
-|---------------------------------------------------------------------*/
-     prog = 1;
-     while (jcount < bsize && prog) {
-/*--------------------   traverse all the current level-set   */
-       last0 = last;
-       jcount0 = jcount;
-       for (inod=begin; inod<=last0; inod++) {
-	 jnod = riord[inod]; 
-/*--------------------   This assumes A is not symmetric.   */
-	 gmat = mat; 
-	 for (k=0; k<2; k++) {
-	   rowj = gmat->ja[jnod];
-	   for (j=0; j<gmat->nzcount[jnod]; j++) {
-	     jcol = rowj[j];
-	     if (iord[jcol] == -1 ) {	
-	       add2is(&last, jcol, iord, riord);
-	       jcount++;
-	     }
-	   }
-	   gmat = matT; 	
-	 }
-       }
-       prog = jcount > jcount0 ? 1 : 0;
-       lastlev = begin;
-       begin = last0+1;
-     }
-/*-----------------------------------------------------------------------
-| the neighbors of elements of last level go to the complement   
-| gmat loop over original matrix and its transpose 
-+-----------------------------------------------------------------------*/ 
-     gmat = mat; 
-     for (k=0; k<2; k++) {
-       for (inod=lastlev; inod<=last; inod++)  {	
-	 jnod = riord[inod]; 
-	 rowj = gmat->ja[jnod];
-	 for (j=0; j<gmat->nzcount[jnod]; j++){	
-	   jcol = rowj[j];
-	   if (iord[jcol] == -1) 
-	     add2com(&nback, jcol, iord, riord);
-	 }
-       }
-       gmat = matT; 	
-     }
-/*   reverse ordering for this level   */
-     mid = (begin0+last) / 2;
-     for (inod=begin0; inod<=mid; inod++) {
-       j = last - inod + begin0;
-       jnod = riord[inod];
-       riord[inod] = riord[j];
-       riord[j] = jnod;
-     }
-   }
-/*--------------------------------------------------
-|  end-main-loop
-|-------------------------------------------------*/
-/*-------------------- relabel nodes of vertex cover   */
-label50:
-   *nnod = last;
-   j1 = *nnod;
-   for (j2=*nnod+1; j2<n; j2++) { 
-     if (iord[riord[j2]] > -1) {
-       if (++j1 != j2) {
-	 j = riord[j2];
-	 riord[j2] = riord[j1];
-	 riord[j1] = j;
-       }
-     }
-   }
-/*-------------------- obtain reverse permutation array   */
-   for (j=0; j<n; j++)
-     iord[riord[j]] = j;
-   (*nnod)++;
-   cleanCS(matT); 
-   free(riord);
-   free(w);
-   return 0;
-}
-/*---------------------------------------------------------------------
-|-----end-of-indsetC---------------------------------------------------
-|--------------------------------------------------------------------*/
-int weightsC(csptr mat, double *w)       
-{
-/*---------------------------------------------------------------------
-|     defines weights based on diagonal dominance ratios
-|--------------------------------------------------------------------*/
-   int irow, k, n=mat->n, *kj, kz;
-   double tdia, wmax=0.0, tnorm, *kr;
-   for (irow=0; irow<n; irow++) {
-      kz = mat->nzcount[irow];
-      kr = mat->ma[irow];
-      kj = mat->ja[irow];
-      tnorm = 0.0;
-      tdia = 0.0;
-      for (k=0; k<kz; k++) {
-	 if (kj[k] == irow) tdia = fabs(kr[k]); 
-	 tnorm += fabs(kr[k]);
-      }
-      if (tnorm > 0.0)
-	 tnorm =  tdia / tnorm;
-      w[irow] = tnorm;
-      if (tnorm > wmax) wmax = tnorm;
-   }
-   for (irow=0; irow<n; irow++)
-      w[irow] = w[irow]/wmax;
-   return 0;
-}
-/*---------------------------------------------------------------------
-|---- end of weightsC -------------------------------------------------
-|--------------------------------------------------------------------*/
-
-int preSel(csptr mat, int *icor, int *jcor, int job, double tol, int *count) 
+int preSel(csptr mat, int *icor, int *jcor, int job, double tol, int *count)
 {
 /*---------------------------------------------------------------------
 | does a preselection of possible diagonal entries. will return a list
--- SRC/arms2.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/arms2.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,9 +1,5 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
 #define  PERMTOL  0.99   /*  0 --> no permutation 0.01 to 0.1 good  */
-#include "globheads.h"
+
 #include "protos.h" 
 
 /*-------------------- end protos */
--- SRC/auxill.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/auxill.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,8 +1,3 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h> 
-#include <math.h>
-#include "globheads.h"
 #include "protos.h"
 #include "ios.h"
 
--- SRC/fgmr.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/fgmr.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,8 +1,3 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include "globheads.h"
 #include "protos.h"
 
 #define  epsmac  1.0e-16
--- SRC/iluk.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/iluk.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,19 +1,7 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include "globheads.h"
 #include "protos.h"
 
-#ifndef min
-#define min(a,b) (((a)>(b))?(b):(a))
-#endif
-#ifndef max
-#define max(a,b) (((a)>(b))?(a):(b))
-#endif
-
 /*--------------------protos */
-int lofC( int lofM, csptr csmat, iluptr lu, FILE *fp ); 
+static int lofC( int lofM, csptr csmat, iluptr lu, FILE *fp ); 
 /*--------------------end protos */
 
 int ilukC( int lofM, csptr csmat, iluptr lu, FILE *fp )
@@ -149,7 +137,7 @@
     return 0;
 }
 
-int lofC( int lofM, csptr csmat, iluptr lu, FILE *fp )
+static int lofC( int lofM, csptr csmat, iluptr lu, FILE *fp )
 {
 /*--------------------------------------------------------------------
  * symbolic ilu factorization to calculate structure of ilu matrix
--- SRC/ilut.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/ilut.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,7 +1,3 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include "globheads.h"
 #include "protos.h"
 
 /*-------------------- end protos*/
--- SRC/ilutc.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/ilutc.c	2013-10-08 23:09:17.000000000 -0400
@@ -5,10 +5,6 @@
  *                                                                       *
  * Report bugs / send comments to: saad@cs.umn.edu, nli@cs.umn.edu       *
  *-----------------------------------------------------------------------*/
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include "globheads.h"
 #include "protos.h"
 
 #ifndef MAXFLOAT
@@ -23,13 +19,12 @@
 static int Lnnz, *Lfirst, *Llist, *Lid, Unnz, *Ufirst, *Ulist, *Uid;
 static double Mnorm, *wL, *wU, *w, *D;
 static csptr L;
-static csptr U; 
+static csptr U;
 
 /*-------------------- protos */
-int update_diagonals(iluptr lu, int i ); 
-int comp(const void *fst, const void *snd ); 
-int std_drop(int lfil, int i, double tolL, double tolU, double toldiag) ;
-int lumsolC(double *y, double *x, iluptr lu );
+static int update_diagonals(iluptr lu, int i ); 
+static int comp(const void *fst, const void *snd ); 
+static int std_drop(int lfil, int i, double tolL, double tolU, double toldiag) ;
 /*-------------------- end protos */
 
 int ilutc(iluptr mt, iluptr lu, int lfil, double tol, int drop, FILE *fp )
@@ -465,7 +460,7 @@
   return 0;
 }
 
-int update_diagonals(iluptr lu, int i )
+static int update_diagonals(iluptr lu, int i )
 {
 /*---------------------------------------------------------------------
  * update diagonals D_{i+1,...,n}
@@ -513,7 +508,7 @@
   return 0;
 }
 
-int comp(const void *fst, const void *snd )
+static int comp(const void *fst, const void *snd )
 {
 /*-------------------------------------------------------------------
  * compares two integers
@@ -525,7 +520,7 @@
   return 0;
 }
 
-int std_drop(int lfil, int i, double tolL, double tolU, double toldiag) 
+static int std_drop(int lfil, int i, double tolL, double tolU, double toldiag) 
 {
 /*---------------------------------------------------------------------
  * Standard Dual drop-off strategy 
--- SRC/ilutpC.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/ilutpC.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,8 +1,3 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include "globheads.h"
 #include "protos.h"
 
 #ifdef ILUTIME
--- SRC/indsetC.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/indsetC.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,12 +1,5 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include "globheads.h"
 #include "protos.h"
 
-void *Malloc( int, char * );
-
 int add2is(int *last, int nod, int *iord, int *riord)
 {
 /*----------------------------------------------------------------------
--- SRC/misc.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/misc.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,11 +1,6 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include "globheads.h"
 #include "protos.h"
 
-int dumpCooMat(csptr A, int nglob, int,  FILE *ft);
+static int dumpCooMat(csptr A, int nglob, int,  FILE *ft);
 
 int qsplitC(double *a, int *ind, int n, int Ncut)
 {
@@ -535,7 +530,7 @@
   return(0);
 }
 
-int dumpCooMat(csptr A, int shiftR, int shiftC, FILE *ft){
+static int dumpCooMat(csptr A, int shiftR, int shiftC, FILE *ft){
   /*-------------------- dump matrix in coo format
     -------------------- matlab style */
   int n, i, k, nzi;
@@ -633,8 +628,3 @@
    qsortR1I(wa, cor1, left, last-1);
    qsortR1I(wa, cor1, last+1, right);
 }
-
-
-
-  
-  
--- SRC/piluNEW.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/piluNEW.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,8 +1,3 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include "globheads.h"
 #include "protos.h"
 
 int pilu(p4ptr amat, csptr B, csptr C, double *droptol, 
--- SRC/setblks.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/setblks.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,6 +1,3 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include "globheads.h"
 #include "protos.h"
 
 typedef struct __KeyType
--- SRC/sets.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/sets.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,9 +1,4 @@
-#include <stdio.h>
-#include <stdlib.h>
 #include <stdarg.h>
-#include <string.h>
-#include <math.h>
-#include "globheads.h"
 #include "protos.h"
 
 void errexit( char *f_str, ... ){
--- SRC/svdInvC.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/svdInvC.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,11 +1,6 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <stdarg.h>
-#include "globheads.h"
 #include "protos.h"
 
 #define TOL 1.e-17
-#define max(a,b) (((a)>(b))?(a):(b))
 
 int invGauss(int nn, double *A) {
   /* *-------------------- inversion by svd
--- SRC/vbiluk.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/vbiluk.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,36 +1,11 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-
-#include "globheads.h"
-
-#ifndef min
-#define min(a,b) (((a)>(b))?(b):(a))
-#endif
-#ifndef max
-#define max(a,b) (((a)>(b))?(a):(b))
-#endif
+#include "protos.h"
+
 #define SVD 1
-#define dgemm dgemm_
 
 /*-------------------- protos */
-void *Malloc(int nbytes, char *msg); 
-int vblusolC(double *y, double *x, vbiluptr lu); 
-int invGauss(int nn, double *A); 
-int invSVD(int nn, double *A) ;
-void dgemm(char*, char*, int*, int*, int*, double*, double*, int*, 
-	   double*, int*, double*, double*, int*) ; 
-int setupVBMat(vbsptr vbmat, int n, int *nB);
-int mallocVBRow(vbiluptr lu, int nrow); 
-void zrmC(int m, int n, BData data); 
-int lofC( int lofM, vbsptr vbmat, vbiluptr lu, FILE *fp ); 
-int setupVBILU(vbiluptr lu, int n, int *bsz);
-void copyBData(int m, int n, BData dst, BData src, int isig);
-int lofC( int lofM, vbsptr vbmat, vbiluptr lu, FILE *fp ); 
+static int lofC( int lofM, vbsptr vbmat, vbiluptr lu, FILE *fp ); 
 /*-------------------- END of protos */
 
-
 int vbilukC( int lofM, vbsptr vbmat, vbiluptr lu, FILE *fp )
 {
 /*----------------------------------------------------------------------------
@@ -188,7 +163,7 @@
     return 0;
 }
 
-int lofC( int lofM, vbsptr vbmat, vbiluptr lu, FILE *fp )
+static int lofC( int lofM, vbsptr vbmat, vbiluptr lu, FILE *fp )
 {
 /*--------------------------------------------------------------------
  * symbolic ilu factorization to calculate structure of ilu matrix
--- SRC/vbilut.c.orig	2013-10-08 23:09:05.000000000 -0400
+++ SRC/vbilut.c	2013-10-08 23:09:17.000000000 -0400
@@ -1,32 +1,4 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-
-#include "globheads.h"
-
-#ifndef min
-#define min(a,b) (((a)>(b))?(b):(a))
-#endif
-#ifndef max
-#define max(a,b) (((a)>(b))?(a):(b))
-#endif
-
-#define qsplit qsplit_ 
-#define gauss gauss_
-#define bxinv bxinv_
-#define dgemm dgemm_
-/*-------------------- protos */
-void *Malloc(int nbytes, char *msg); 
-void zrmC(int m, int n, BData data); 
-void copyBData(int m, int n, BData dst, BData src, int isig);
-void dgemm(char*, char*, int*, int*, int*, double*, double*, int*, 
-	   double*, int*, double*, double*, int*) ; 
-int vblusolC(double *y, double *x, vbiluptr lu); 
-void gauss (int *, double*, int*); 
-void bxinv (int*, int*, double*,double*,double*);
-int setupVBILU(vbiluptr lu, int n, int *bsz);
-void qsplit(double*, int*, int*, int*) ;
-/*-------------------- END protos */
+#include "protos.h"
 
 int vbilutC( vbsptr vbmat, vbiluptr lu, int lfil, double tol,
              BData *w, FILE *fp )