// utils.cc
//  provides functions for 
//  integer and double vectors and matrices
//  memory allocation and clean-up
//
//  DJCM Thu 1/11/07    Free Software, released into the public domain

#include <iostream>
using namespace std;

#include "utils.h"

// allocates memory for an array of integers, say
//          b[low]..b[high]
// returns zero if failure
int *ivector ( int low, int high ) {
  int N = high-low+1 ;
  int *a ; 
  if ( N <= 0 ) {
    cout << "Illegal range in ivector: "
	 << low << ", " << high << endl ;
    return 0 ;
  }
  a = new int[N] ; // allocates the memory for a[0]..a[N-1]
  if(!a) {
    cout << "Memory allocation failed\n" ;
    return 0 ; 
  }
  return (a-low) ; // offset the pointer by low.
  //                  the user uses b[low]..b[high]
}

// allocates memory for a pointer to an array of pointers
// and calls ivector to allocate memory for a matrix
//          m[low1][low2] .. m[high1][high2].
// returns 0 if something goes wrong.
int **imatrix ( int low1, int high1, int low2, int high2 ) {
  int N1 = high1-low1+1 ;
  int **m ; 
  if ( N1 <= 0 ) {
    cout << "Illegal range in imatrix: "
	 << low1 << ", " << high1 << endl ;
    return 0 ;
  }
  m = new int*[N1] ; // allocates the space for m[0]..m[N-1]
  // Note m is zero-offset here
  if(!m) {
    cout << "Memory allocation failed\n" ;
    return 0 ; 
  }
  for ( int n=0 ; n < N1 ; n ++ ) {
    m[n] = ivector(low2,high2) ; // allocate the nth row of matrix
    if(!m[n]) return 0 ; 
  }
  return (m-low1) ; // offset the pointer by low,
  //                   so that the user uses m[low][]..m[high][]
}

// allocates memory for an array of doubles, say
//          b[low]..b[high]
// returns zero if failure
double *dvector ( int low, int high ) {
  int N = high-low+1 ;
  double *a ; 
  if ( N <= 0 ) {
    cout << "Illegal range in dvector: "
	 << low << ", " << high << endl ;
    return 0 ;
  }
  a = new double[N] ; // allocates the memory for a[0]..a[N-1]
  if(!a) {
    cout << "Memory allocation failed\n" ;
    return 0 ; 
  }
  return (a-low) ; // offset the pointer by low.
  //                  the user uses b[low]..b[high]
}

// allocates memory for a pointer to an array of pointers
// and calls ivector to allocate memory for a matrix
//          m[low1][low2] .. m[high1][high2].
// returns 0 if something goes wrong.
double **dmatrix ( int low1, int high1, int low2, int high2 ) {
  int N1 = high1-low1+1 ;
  double **m ; 
  if ( N1 <= 0 ) {
    cout << "Illegal range in imatrix: "
	 << low1 << ", " << high1 << endl ;
    return 0 ;
  }
  m = new double*[N1] ; // allocates the space for m[0]..m[N-1]
  // Note m is zero-offset here
  if(!m) {
    cout << "Memory allocation failed\n" ;
    return 0 ; 
  }
  for ( int n=0 ; n < N1 ; n ++ ) {
    m[n] = dvector(low2,high2) ; // allocate the nth row of matrix
    if(!m[n]) return 0 ; 
  }
  return (m-low1) ; // offset the pointer by low,
  //                   so that the user uses m[low][]..m[high][]
}

void freeivector ( int *a , int low ) {
  delete [] &(a[low]) ; // The '[]' indicate that what's
} //                       being freed is an array

int freeimatrix ( int **m , int low1, int high1, int low2, int high2 ) {
  int N1 = high1-low1+1 ;
  if ( N1 <= 0 ) {
    cout << "Illegal range in freeimatrix: "
	 << low1 << ", " << high1 << endl ;
    return 0 ;
  }
  // Note m is offset by low1 here
  for ( int n=low1 ; n <= high1 ; n ++ ) {
    delete [] &(m[n][low2]) ; // free the nth row of the matrix
  }
  delete [] &(m[low1]) ; 
  return 1 ; // 1 indicates success
}

void freedvector ( double *a , int low ) {
  delete [] &(a[low]) ; // The '[]' indicate that what's
} //                       being freed is an array

int freedmatrix ( double **m , int low1, int high1, int low2, int high2 ) {
  int N1 = high1-low1+1 ;
  if ( N1 <= 0 ) {
    cout << "Illegal range in freedmatrix: "
	 << low1 << ", " << high1 << endl ;
    return 0 ;
  }
  // Note m is offset by low1 here
  for ( int n=low1 ; n <= high1 ; n ++ ) {
    delete [] &(m[n][low2]) ; // free the nth row of the matrix
  }
  delete [] &(m[low1]) ; 
  return 1 ; // 1 indicates success
}

double innerProduct ( double *a , double *b , int low , int high ) {
  double  answer = 0.0;
  for ( int n = low ; n <= high ; n ++ ) {
    answer += a[n] * b[n] ; 
  }
  return answer ;
}

// sets vector to constant
void   constantdvector( double *y , int lo, int hi, double c  ) {
  for ( int n = lo ; n <= hi ; n ++ ) {
    y[n] = c;
  }  
}

// a += epsilon * b
void incrementdvector ( double *a , double epsilon, double *b , int low , int high ) {
  for ( int n = low ; n <= high ; n ++ ) {
    a[n] += epsilon * b[n] ; 
  }
}

void printMatrix( double **M , int low1, int high1, int low2, int high2 ) {
  for ( int m = low1 ; m <= high1 ; m ++ ) {
    for ( int n = low2 ; n <= high2 ; n ++ ) {
      cout << M[m][n]  ;
      if ( n == high2 ) cout << endl ;
      else cout << "\t" ; 
    }
  }
}
// print matrix transposed
void printMatrixT( double **M , int lo1, int hi1, int lo2, int hi2 ) {
  for ( int n = lo2 ; n <= hi2 ; n ++ ) {
    for ( int m = lo1 ; m <= hi1 ; m ++ ) {
      cout << M[m][n]  ;
      if ( m == hi1 ) cout << endl ;
      else cout << "\t" ; 
    }
  }
}

void printVector( double * b , int lo , int hi , int style=1 ) {
  for ( int n = lo ; n <= hi ; n ++ ) {
    //    cout << "b[" << n << "] = " ;
    cout << b[n] ;
    if(style) {
      if ( n == hi ) 
	cout << endl;
      else
	cout << "\t" ;
    } else {
      cout << endl;
    }
  }
}

void fprintVector( double * b , int lo , int hi , ostream &fout, int style  ) {
  for ( int n = lo ; n <= hi ; n ++ ) {
    //    cout << "b[" << n << "] = " ;
    fout << b[n] ;
    if(style) {
      if ( n == hi ) 
	fout << endl;
      else
	fout << "\t" ;
    } else {
      fout << endl;
    }
  }
}