///////////////////////////////////////////////////////////////////////
//matrix.h
#ifndef MATRIX_H
#define MATRIX_H
#include "vector.h"

class Matrix{
  //class members
  int nRows, nCols;
  void pointerCheck() const;
  void pointerCheck(int i) const;

 protected:
  double** pA;

 public:
  //constructors and destructor
  Matrix(int nrows = 0, int ncols = 0, double a = 0.);
  Matrix(int nrows, int ncols, const double* const* pa);
  Matrix(const Matrix& A);
  virtual ~Matrix();

  //overloaded operators
  Matrix& operator=(const Matrix& B);
  Matrix operator+(const Matrix& B) const;
  Matrix& operator+=(const Matrix& B);
  Matrix operator-(const Matrix& B) const;
  Matrix& operator-=(const Matrix& B);

  Matrix operator*(const Matrix& B) const; 
  Matrix operator*(double b) const;
  Vector operator*(const Vector& v) const;

  const double* operator[](int i) const{
    return pA[i];
  }

  //matrix operations
  Matrix trans() const;

  //solution of a linear system
  virtual Vector backward(const Vector& RHS) const;
  virtual Vector forward(const Vector& RHS) const;
  Vector gauss(const Vector& RHS) const;
  virtual Vector cholesky(const Vector& RHS) const;
  //banded systems
  Vector bandBack(const Vector& RHS, int bWidth) const;
  Vector bandFor(const Vector& RHS, int bWidth) const;
  Vector bandChol(const Vector& RHS, Matrix& G, int bWidth) const;

  void QR(Matrix& Q, Matrix& R) const;

  //computation of spectra 
  Vector eigen(int Nvals, Vector& v, Matrix& U, int itMax = 38, 
               double delta = .00000001) const;
  Vector SVD(int Nvals, Vector& v, Matrix& U, Matrix& V, 
             int itMax = 38, double delta = .00000001) const;

  //utilities and accessors
  void printMatrix() const;
  Vector matToVec(int j) const;//jth column to a vector
  void vecToMat(int j, const Vector& v);//vector into jth column
  int getnRows() const{
    return nRows;
    }
  int getnCols() const{
    return nCols;
  }

  friend Matrix operator*(const Vector& v1, const Vector& v2);
};

class pdBand : public Matrix{
  //derived class member
  int bW;

 public:

  pdBand(int nrows = 0, int bw = 0, double a = 0.);
  pdBand(int nrows, int bw, const double* const* pa);
  pdBand(const pdBand& A);
 
  virtual ~pdBand(){};

  pdBand& operator=(const pdBand& B);
  pdBand operator+(const pdBand& B) const;
  Matrix operator+(const Matrix& B) const;


  pdBand trans() const{return pdBand(*this);}

  virtual Vector backward(const Vector& RHS) const;
  virtual Vector forward(const Vector& RHS) const;
  virtual Vector cholesky(const Vector& RHS) const;

  int getBW() const{return bW;}
  Matrix bandToFull() const;
};

Matrix operator*(const Vector& v1, const Vector& v2);

Matrix operator*(double, const Matrix& A);

#endif
///////////////////////////////////////////////////////////////////////
//matrix.cpp
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <limits>
#include "matrix.h"

using std::cout; using std::endl;

//Constructors/destructors:

//diagonal matrix constructor
Matrix:: Matrix(int nrows, int ncols, double a){
  nRows = nrows;
  nCols = ncols;
  //set pA to null pointer in default case
  if(nRows == 0 || nCols == 0){
    pA = 0;
    return;
  }
  pA = new(std::nothrow) double*[nRows];
  pointerCheck();

  for(int i = 0; i < nRows; i++){
    pA[i] = new(std::nothrow) double[nCols];
    pointerCheck(i);
    for(int j = 0; j < nCols; j++){
      if(i == j)
        pA[i][j] = a;
      else
        pA[i][j] = 0.;
    }
  }
}
    
//constructor from a pointer-to-pointer 
Matrix:: Matrix(int nrows, int ncols, const double* const* pa){
  nRows = nrows;
  nCols = ncols;
  pA = new(std::nothrow) double*[nRows];
  pointerCheck();
      
  for(int i = 0; i < nRows; i++){
    pA[i] = new(std::nothrow) double[nCols];
    pointerCheck(i);
  }

  for(int i = 0; i < nRows; i++)
    for(int j = 0; j < nCols; j++)
      pA[i][j] = pa[i][j];
}

//copy constructor
Matrix::Matrix(const Matrix& A){
    nRows = A.nRows;
    nCols = A.nCols;
    pA = new(std::nothrow) double*[nRows];
    pointerCheck();
 
    for(int i = 0; i < nRows; i++){
      pA[i] = new(std::nothrow) double[nCols];
      pointerCheck(i);
    }

    for(int i = 0; i < nRows; i++)
      for(int j = 0; j < nCols; j++)
        pA[i][j] = A.pA[i][j];
}    
      
//destructor
Matrix::~Matrix(){   
  if(pA != 0){
    for(int i = 0; i < nRows; i++)
      delete[] pA[i];
    delete[] pA;
  }
}
    
//Overloaded operators:

Matrix& Matrix::operator=(const Matrix& A){

  if(this == &A) //avoid self assignment
    return *this;
  //if dimensions match we can just overwrite; otherwise....
  if((nRows != A.nRows)||(nCols != A.nCols)){
    
    if(pA != 0){//check first before releasing memory
      for(int i = 0; i < nRows; i++)//release current memory
        delete[] pA[i];
      delete[] pA;
    }
    
    //define/redefine object's members
    nRows = A.nRows; 
    nCols = A.nCols;
    pA = new(std::nothrow) double*[nRows];
    pointerCheck();      
 
    for(int i = 0; i < nRows; i++){
      pA[i] = new(std::nothrow) double[nCols];
      pointerCheck(i);
    }
  }
  for(int i = 0; i < nRows; i++)
    for(int j = 0; j < nCols; j++)
      pA[i][j] = A.pA[i][j];
    
  return *this;
}    

Matrix Matrix::operator+(const Matrix& B) const{
  if(nRows != B.nRows || nCols != B.nCols){
    cout << "Bad row and/or column dimensions in +!" << endl;
    exit(1);
  }
  Matrix temp(nRows, nCols);
  for(int i = 0; i < nRows; i++)
    for(int j = 0; j < nCols; j++)
      temp.pA[i][j] = this->pA[i][j] + B.pA[i][j];
  
  return temp;
}

Matrix& Matrix::operator+=(const Matrix& B){
  if(nRows != B.nRows || nCols != B.nCols){
    cout << "Bad row and/or column dimensions in +=!" << endl;
    exit(1);
  }
  for(int i = 0; i < nRows; i++)
    for(int j = 0; j < nCols; j++)
      pA[i][j] += B.pA[i][j];

  return *this;
}

Matrix Matrix::operator-(const Matrix& B) const{
  if(nRows != B.nRows || nCols != B.nCols){
    cout << "Bad row and/or column dimensions in -!" << endl;
    exit(1);
  }
  Matrix temp(nRows, nCols);
  for(int i = 0; i < nRows; i++)
    for(int j = 0; j < nCols; j++)
      temp.pA[i][j] = this->pA[i][j] - B.pA[i][j];

  return temp;
}

Matrix& Matrix::operator-=(const Matrix& B){
  if(nRows != B.nRows || nCols != B.nCols){
    cout << "Bad row and/or column dimensions in -=!" << endl;
    exit(1);
  }
  for(int i = 0; i < nRows; i++)
    for(int j = 0; j < nCols; j++)
      pA[i][j] -= B.pA[i][j];

  return *this;
}

Matrix Matrix::operator*(const Matrix& B) const{
  if(nCols != B.nRows){
    cout << "Bad row and column dimensions in *!" << endl;
    exit(1);
  }
  Matrix C(nRows, B.nCols);//matrix of all 0s
  for(int i = 0; i < nRows; i++)
    for(int j = 0; j < B.nCols; j++)
      for(int k = 0; k < B.nRows; k++)
        C.pA[i][j] += pA[i][k]*B.pA[k][j];
  
  return C;
}

Matrix Matrix::operator*(double b) const{

  Matrix C(nRows, nCols);//matrix of all 0s
  for(int i = 0; i < nRows; i++)
    for(int j = 0; j < nCols; j++)
      C.pA[i][j] = b*pA[i][j];

  return C;
}

Vector Matrix::operator*(const Vector& v) const{
  double temp;
  Vector c(nRows);//matrix of all 0s
  for(int i = 0; i < nRows; i++){
    temp = 0;
    for(int j = 0; j < nCols; j++)
      temp += pA[i][j]*v[j];
    
    c.pA[i]=temp;
  }

  return c;
}

//matrix operations

Matrix Matrix::trans() const{
  Matrix B(nCols, nRows, 0.);
  for(int i = 0; i < nCols; i++)
    for(int j = 0; j < nRows; j++)
      B.pA[i][j] = pA[j][i];
  return B;
}


//solution of a linear system

//back-solve an upper triangular system 
Vector Matrix::backward(const Vector& RHS) const{
  double temp;//temporary storage 

  Vector b(nRows, 0.);//solution vector

  //initialize recursion on last row of all right-hand-sides
  if(pA[nRows - 1][nRows - 1] != 0)
    b.pA[nRows - 1] = RHS.pA[nRows - 1]/pA[nRows - 1][nRows - 1];
  else{
    cout << "Singular system!" << endl;
    exit(1);
  }
  
  //now work through the remaining rows
  for(int i = (nRows - 2); i >=0 ; i--){
    if(pA[i][i] != 0){
      temp = RHS.pA[i];
      for(int k = (i + 1); k < nRows; k++)
        temp -= b.pA[k]*pA[i][k]; 
      b.pA[i] = temp/pA[i][i];
    }
    else{
      cout << "Singular system!" << endl;
      exit(1);
    }
  }
  return b;
}

//forward-solve a lower triangular system 
Vector Matrix::forward(const Vector& RHS) const{
  double temp;//temporary storage

  Vector b(nRows, 0.);//solution matrix

  //initialize recursion on first row of all right-hand-sides
  if(pA[0][0] != 0)
    b.pA[0] = RHS.pA[0]/pA[0][0];
  else{
    cout << "Singular system!" << endl;
    exit(1);
    }
  
  //now work through the remaining rows
  for(int i = 1; i < nRows; i++){
    if(pA[i][i] != 0){
      temp = RHS.pA[i];
      for(int k = 0; k < i; k++)
        temp -= b.pA[k]*pA[i][k]; 
      b.pA[i] = temp/pA[i][i];
    }
    else{
      cout << "Singular system!" << endl;
      exit(1);
    }
  }
  return b;
}

Vector Matrix::gauss(const Vector& RHS) const{
  double multiplier = 0;
  Matrix G(nRows, nRows, pA);
  Vector h(nRows, RHS.pA);
  
  for(int j = 0; j < nCols; j++){//column to be swept
    
    for(int jj = j + 1; jj < nCols; jj++){//current operation column

      if(G[j][j] == 0){
        cout << "Oops! Division by 0!" << endl;
        exit(1);
      }

      for(int i = j + 1; i < nRows; i++){//work down rows
        multiplier = G[i][j]/G[j][j];
        G.pA[i][jj] = G.pA[i][jj] - multiplier*G.pA[j][jj];
      }
    }
    
    for(int i = j + 1; i < nRows; i++){//do the same to the RHS
      multiplier = G[i][j]/G[j][j];
      h.pA[i] = h.pA[i] - multiplier*h.pA[j];
    }
  }
  //now backsolve
  Vector b = G.backward(h);
  return b;
}

Vector Matrix::cholesky(const Vector& RHS) const{
  double temp = 0;
  Matrix G(nRows, nCols, 0.); 
  for(int j = 0; j < nCols; j++){//proceed by columns
    if(pA[j][j] == 0){
      cout << "Singular system!" << endl;
      exit(1);
    }
    temp = pA[j][j];//starting value for diagonal element recursion
    
    if(j > 0)
      for(int k = 0; k < j; k++)
        temp -= G[j][k]*G[j][k];
    
    G.pA[j][j] = sqrt(temp);
    
    for(int i = (j + 1); i < nRows; i++){//now do the rest
      temp = pA[j][i];
      for(int k = 0; k < j; k++)
        temp -= G[j][k]*G[i][k];
      G.pA[i][j] = temp/G[j][j];
    }
  }
  cout << "Cholesky factor:" << endl;
  G.printMatrix();
  Vector h = G.forward(RHS);
  Vector b = G.trans().backward(h);
  return b;
}    

//back-solve a banded  upper triangular system 
Vector Matrix::bandBack(const Vector& RHS, int bWidth) const{
  double temp;//temporary storage
  int up;

  Vector b(nRows, 0.);//solution vector

  //initialize recursion on last row of right-hand-side
  if(pA[nRows-1][nRows-1] != 0)
        b.pA[nRows - 1] = RHS.pA[nRows - 1]/pA[nRows - 1][nRows - 1];
  else{
    cout << "Singular system!" << endl;
    exit(1);
  }
  
  //now work through the remaining rows
  for(int i = (nRows - 2); i >= 0; i--){
    if(pA[i][i] != 0){
      temp = RHS.pA[i];
      up = std::min(nCols, i + bWidth + 1);
      for(int k = (i + 1); k < up; k++)
        temp -= b.pA[k]*pA[i][k]; 
      b.pA[i] = temp/pA[i][i];
    }
    else{
      cout << "Singular system!" << endl;
      exit(1);
    }
  }
  return b;
}

//forward-solve a banded lower triangular system 
Vector Matrix::bandFor(const Vector& RHS, int bWidth) const{
  double temp;//temporary storage
  int low;

  Vector b(nRows, 0.);//solution vector

  //initialize recursion on first row of right-hand-side
  if(pA[0][0] != 0){
    b.pA[0] = RHS.pA[0]/pA[0][0];
  }
  else{
    cout << "Singular system!" << endl;
    exit(1);
  }
  
  //now work through the remaining rows
  for(int i = 1; i < nRows; i++){
    if(pA[i][i] != 0){
      temp = RHS.pA[i];
      low = std::max(0, i - bWidth);
      for(int k = low; k < i; k++)
        temp -= b.pA[k]*pA[i][k]; 
      b.pA[i] = temp/pA[i][i];
    }
    else{
      cout << "Singular system!" << endl;
      exit(1);
    }
  }
  return b;
}

//solve a banded system via Cholesky factorization
Vector Matrix::bandChol(const Vector& RHS, Matrix& G, 
			int bWidth) const{
  double temp = 0;
  int low, up;

  for(int j = 0; j < nCols; j++){
    if(pA[j][j] == 0){
      cout << "Singular system!" << endl;
      exit(1);
    }
    temp = pA[j][j];
    
    low = std::max(0, j - bWidth);
      
    for(int k = low; k < j; k++)
      temp -= G[j][k]*G[j][k];
      
    G.pA[j][j] = sqrt(temp);

    up = std::min(j + bWidth, nRows - 1);
      
    for(int i = (j + 1); i <= up; i++){
      temp = pA[j][i];
      for(int k = low; k < j; k++)
        temp -= G[j][k]*G[i][k];
      G.pA[i][j] = temp/G[j][j];
    }
  }
  Vector h = G.bandFor(RHS, bWidth);
  Vector b = G.trans().bandBack(h, bWidth);
  return b;
}    

void Matrix::QR(Matrix& Q, Matrix&R) const{
  if(nRows < nCols){
    cout << "The number of rows is less than the number of columns!" << endl;
    exit(1);
  }

  //make work matrix copy of *this
  Matrix Work = Matrix(nRows, nCols, this->pA);
  Vector temp;
  Vector q;
  for(int i = 0; i < nCols; i++){
    q = Work.matToVec(i);
    R.pA[i][i] = sqrt(q.dotProd(q));
    q = q*(1/R.pA[i][i]);
    Q.vecToMat(i, q);
    temp = Work.trans()*q;
    for(int j = (i + 1); j < nCols; j++)
      R.pA[i][j] = temp[j];
    Work -= q*temp;//outer product update
  }
  cout << "The Q matrix" << endl;
  Q.printMatrix();
  cout << "The R matrix" << endl;
  R.printMatrix();
  cout << "The product QR" << endl;
  (Q*R).printMatrix();
}

//eigenvalues, eigenvectors and SVD

Vector Matrix::eigen(int Nvals, Vector& v, Matrix& U, int itMax, 
		     double delta) const{
  double change, temp;
  int niter;
  
  Matrix ACopy(nRows, nCols, pA);//work copy of A
  Vector vCopy;//work copy of v
  
  Vector lambda(Nvals, 0.);
  for(int i = 0; i < Nvals; i++){
    change = std::numeric_limits<double>::infinity();
    niter = 0;
    vCopy = v;
    lambda.pA[i] = 0;
    
    while((change > delta)&&(niter < itMax)){
      vCopy = ACopy*vCopy;
      vCopy = (1./sqrt(vCopy.dotProd(vCopy)))*vCopy;
      lambda.pA[i] = vCopy.dotProd(ACopy*vCopy);
      if(niter == 0)
        temp = lambda.pA[i];
      else if(temp != 0.){
        change = fabs(1. - lambda.pA[i]/temp);
        temp = lambda.pA[i];
      }
      niter++;
    }
    U.vecToMat(i, vCopy);
    v -= (vCopy.dotProd(vCopy))*vCopy; 
    ACopy -= lambda.pA[i]*(vCopy*vCopy);
  }
  return lambda;
}
    
Vector Matrix::SVD(int Nvals, Vector& v, Matrix& U, Matrix& V, 
		   int itMax, double delta) const{
  Vector temp;
  Matrix ACopy(nRows, nCols, pA); 
  
  if(nRows > nCols){//work with A^TA
    temp = (ACopy.trans()*ACopy).eigen(Nvals, v, V, itMax, delta);
    Matrix Temp(Nvals, Nvals, 0.);
    for(int i = 0; i < Nvals; i++)
      Temp.pA[i][i] = 1./sqrt(temp[i]);
    U = ACopy*V*Temp;
  }
  else{//work with AA^T
    temp = (ACopy*ACopy.trans()).eigen(Nvals, v, U, itMax, delta);
    Matrix Temp(Nvals, Nvals, 0.);
    for(int i = 0; i < Nvals; i++)
      Temp.pA[i][i] = 1./sqrt(temp[i]);
    V = ACopy.trans()*U*Temp;
  }
  Vector lambda(Nvals, 0.);
  for(int i = 0; i < Nvals; i++)
    lambda.pA[i] = sqrt(temp[i]);
  return lambda;
}

// Utilities:

void Matrix::printMatrix() const{
  for(int i = 0; i < nRows; i++){
    for(int j = 0; j < nCols; j++)
      cout << " " << pA[i][j] << "  ";
    cout << endl;
  }
}

Vector Matrix::matToVec(int j) const{
  Vector v(nRows, 0.);
  for(int i = 0; i < nRows; i++)
    v.pA[i] = pA[i][j];
  return v;
}

void Matrix::vecToMat(int j, const Vector& v){
  for(int i = 0; i < nRows; i++)
    pA[i][j] = v.pA[i];
}

void Matrix::pointerCheck() const{
  if(pA == 0){
    cout << "Memory allocation for pA failed" << endl;
    exit(1);
  }
}

void Matrix::pointerCheck(int i) const{
  if(pA[i] == 0){
    cout << "Memory allocation for pA[" << i << "] failed" << endl;
    exit(1);
  }
}

//Overloaded multiplication operators
Matrix operator*(const Vector& v1, const Vector& v2){
  Matrix B(v1.getnRows(), v2.getnRows());
  for(int i = 0; i < v1.getnRows(); i++)
    for(int j = 0; j < v2.getnRows(); j++)
      B.pA[i][j] = v1[i]*v2[j];
  return B;
}

Matrix operator*(double b, const Matrix& A){
  return A*b;
}

//pdBand methods
 
pdBand::pdBand(int r, int bw, double a) 
  : Matrix(r, bw + 1, a) {
  bW = bw;
}

pdBand::pdBand(int r, int bw, const double* const* pa) 
  : Matrix(r, bw + 1, pa) {
  bW = bw;
}

pdBand::pdBand(const pdBand& A) 
  : Matrix(A) {
  bW = A.bW;
}

pdBand& pdBand::operator=(const pdBand& A){

  if(this == &A)//avoid self assignment
    return *this;
  this->Matrix::operator=(A);
  bW = A.bW;

  return *this;
}
 
Matrix pdBand::operator+(const Matrix& B) const{
  Matrix A = this->bandToFull();
  return (A + B);
}
  
pdBand pdBand::operator+(const pdBand& B) const{
  pdBand temp(getnRows(), B.bW);
  temp.Matrix::operator=(this->Matrix::operator+(B));

  return temp;
}
  
//back-solve a banded  upper triangular system 
Vector pdBand::backward(const Vector& RHS) const{
  double temp;//temporary storage
  int up;
  int r = this->getnRows();
  Vector b(r, 0.);//solution vector

  //initialize recursion on last row of right-hand-side
  if(pA[r - 1][0] != 0)
    b.pA[r - 1] = RHS.pA[r - 1]/pA[r - 1][0];
  else{
    cout << "Singular system!" << endl;
    exit(1);
  }
  
  //now work through the remaining rows
  for(int i = (r - 2); i >= 0; i--){
    if(pA[i][0] != 0){
      temp = RHS.pA[i];
      up = std::min(r, bW + 1);
      for(int k = 1; k < up; k++)
        temp -= b.pA[i+k]*pA[i][k]; 
      b.pA[i] = temp/pA[i][0];
    }
    else{
      cout << "Singular system!" << endl;
      exit(1);
    }
  }
  return b;
}

//forward-solve a banded lower triangular system 
Vector pdBand::forward(const Vector& RHS) const{
  double temp;//temporary storage
  int low;
  int r = this->getnRows();
  Vector b(r, 0.);//solution vector

  //initialize recursion on first row of right-hand-side
  if(pA[0][0] != 0){
    b.pA[0] = RHS.pA[0]/pA[0][0];
  }
  else{
    cout << "Singular system!" << endl;
    exit(1);
  }
  
  //now work through the remaining rows
  for(int i = 1; i < r; i++){
    if(pA[i][0] != 0){
      temp = RHS.pA[i];
      low = std::max(0, i - bW);
      for(int k = low; k < i; k++)
        temp -= b.pA[k]*pA[k][i - k]; 
      b.pA[i] = temp/pA[i][0];
    }
    else{
      cout << "Singular system!" << endl;
      exit(1);
    }
  }
  return b;
}

//solve a banded system via Cholesky factorization
Vector pdBand::cholesky(const Vector& RHS) const{
  double temp = 0;
  int low, up;
  int r = this->getnRows();
  pdBand G(r, bW);

  for(int i = 0; i < r; i++){
    if(pA[i][0] == 0){
      cout << "Singular system!" << endl;
      exit(1);
    }
    temp = pA[i][0];
    low = std::max(0, i - bW);
      
    for(int k = low; k < i; k++)
      temp -= G[k][i - k]*G[k][i - k];
      
    G.pA[i][0] = sqrt(temp);

    up = std::min(i + bW, r - 1);
    
    for(int j = (i + 1); j <= up; j++){
      temp = pA[i][j - i];
      low = std::max(0, j - bW);
      for(int k = low; k < i; k++)
        temp -= G[k][i - k]*G[k][j - k];
      G.pA[i][j - i] = temp/G[i][0];
    }
  }
  cout << "Cholesky factor in band storage" << endl;
  G.printMatrix();
  Vector h = G.forward(RHS);
  Vector b = G.backward(h);
  std::cout << "Solution vector" << endl;
  return b;
}    

Matrix pdBand::bandToFull() const{
  double** pTemp = new(std::nothrow) double*[getnRows()];
  if(pTemp == 0){
    cout << "Memory allocation for pA failed" << endl;
    exit(1);
  }
  for(int i = 0; i < getnRows(); i++){
    pTemp[i] = new(std::nothrow) double[getnRows()];
    if(pA[i] == 0){
      cout << "Memory allocation for pA failed" << endl;
      exit(1);
    }
    for(int j = i; j < getnRows(); j++){
      if(j <= std::min(getnRows(), i + bW))
         pTemp[i][j] = pA[i][j - i];
      else pTemp[i][j] = 0.;
    }
  }
  for(int i = 1; i < getnRows(); i ++)
    for(int j = 0; j < i; j++)
      pTemp[i][j] = pTemp[j][i];
  Matrix A(getnRows(), getnRows(), pTemp);
  for(int i = 0; i < getnRows(); i++)
    delete[] pTemp[i];
  delete[] pTemp;
  return A;
}
///////////////////////////////////////////////////////////////////////
//vector.h
#ifndef VECTOR_H
#define VECTOR_H

class Vector{
  // class members
  double* pA; 
  int nRows;
  void pointerCheck() const;

 public:

  //constructors and destructor
  
  Vector(int nrows, double const* pa);//constructor 
  Vector(int nrows = 0, double b = 0.);//default constructor 
  Vector(const Vector& v);
  ~Vector();
  
  //overloaded operators
  Vector& operator=(const Vector& v);
  Vector operator+(const Vector& v) const;
  Vector& operator+=(const Vector& v);
  Vector operator-(const Vector& v) const;
  Vector& operator-=(const Vector& v);
  Vector operator*(double b) const;
  double operator[](int i) const{
    return pA[i];
  }

  //utilities and accessors
  double dotProd(const Vector& v) const;
  void printVec() const;
  int getnRows() const{return nRows;}
  
  //friends
  friend class Matrix;
  friend class pdBand;
  
};

Vector operator*(double b, const Vector& v);
#endif
///////////////////////////////////////////////////////////////////////
//vector.cpp
#include "vector.h"
#include <iostream>
#include <cstdlib>

using std::cout; using std::endl;

Vector::Vector(int nrows, double const* pa){
  nRows = nrows;
  pA = new(std::nothrow) double[nRows];
  pointerCheck();
      
  for(int i = 0; i < nRows; i++)
    pA[i] = pa[i];
}

Vector::Vector(int nrows, double b){
  nRows = nrows;
  if(nRows == 0){
    pA = 0;
    return;
  }
  pA = new(std::nothrow) double[nRows];
  pointerCheck();
      
  for(int i = 0; i < nRows; i++)
    pA[i] = b;
}

Vector::Vector(const Vector& v){
  nRows = v.nRows;
  pA = new(std::nothrow) double[nRows];
  pointerCheck();
      
  for(int i = 0; i < nRows; i++)
    pA[i] = v.pA[i];
}    

Vector::~Vector(){
  delete[] pA;
}

double Vector::dotProd(const Vector& v) const{
  double dot = 0;
  for(int i = 0; i < nRows; i++)
    dot += v.pA[i]*pA[i];
  return dot;
}

Vector& Vector::operator=(const Vector& v){

  if(this == &v)//avoid self assignment
    return *this;

  //define/redefine object's members
  if(nRows != v.nRows){

    if(pA != 0)
      delete[] pA;

    nRows = v.nRows; 
    pA = new(std::nothrow) double[nRows];
    pointerCheck();
  }

  for(int i = 0; i < nRows; i++)
    pA[i] = v.pA[i];

  return *this;
}    

Vector Vector::operator+(const Vector& v) const{
  if(nRows != v.nRows){
    cout << "Bad row dimensions in +!" << endl;
    exit(1);
  }
  Vector temp(nRows);
  for(int i = 0; i < nRows; i++)
    temp.pA[i] = this->pA[i] + v.pA[i];
  
  return temp;
}

Vector& Vector::operator+=(const Vector& v){
  if(nRows != v.nRows){
    cout << "Bad row dimensions in +=!" << endl;
    exit(1);
    }
  for(int i = 0; i < nRows; i++)
    pA[i] += v.pA[i];

  return *this;
}

Vector Vector::operator-(const Vector& v) const{
  if(nRows != v.nRows){
    cout << "Bad row dimensions in -!" << endl;
    exit(1);
    }
  Vector temp(nRows);
  for(int i = 0; i < nRows; i++)
    temp.pA[i] = this->pA[i] - v.pA[i];

  return temp;
}

Vector& Vector::operator-=(const Vector& v){
  if(nRows != v.nRows){
    cout << "Bad row dimensions in -=!" << endl;
    exit(1);
  }
  for(int i = 0; i < nRows; i++)
    pA[i] -= v.pA[i];

  return *this;
}

Vector Vector::operator*(double b) const{
  Vector temp(nRows);
  for(int i = 0; i < nRows; i++)
    temp.pA[i] = b*pA[i];

  return temp;
}

void Vector::pointerCheck() const{
  if(pA == 0){
    cout << "Memory allocation for pA failed" << endl;
    exit(1);
  }
}

void Vector::printVec() const{
  for(int i = 0; i < nRows; i++)
    cout << pA[i] << " " << endl;
}

Vector operator*(double b, const Vector& v){
  Vector B = v;
  return B*b;
}