// partition.h

#include <iostream.h>

#ifndef _partition_h
#define _partition_h

#include "matrix.h"
#include "tarray1.h"

typedef double Real;

#define SP << " " <<

void error(const char*);

class Precision
    {
public:
    int ndigits;
    double data;
    Precision(double d,int nd) : data(d), ndigits(nd) { }
    };
inline ostream& operator << (ostream &s, Precision p) 
    {
    int oldp = cout.precision(p.ndigits);
    s << p.data;
    s.precision(oldp);
    return s;
    }

class Partition1D
    {
    void init2(int ind1, int ind2) 
	{
	lastind->ReDimension(0,2);
	lastind[0] = 0;
	lastind[1] = ind1;
	lastind[2] = ind2;
	}
public:
    IntArray1 lastind;
    int Length() { return lastind(lastind->Upper()); }
    int NumParts() { return lastind.Size() - 1; }

    Partition1D(int l1, int l2)		// The size of each part
    	{ init2(l1,l1+l2); }

    Partition1D(const VectorRef& v, int l1) 
    	{ init2(l1,v.Length()); }

    Partition1D(const VectorRef& v1, const VectorRef& v2) 
    	{ init2(v1.Length(),v1.Length()+v2.Length()); }

    Partition1D(const IntArray1& inds);

    VectorRef Sub(const VectorRef& v, int part) const
	{ return v.SubVector(lastind(part-1)+1,lastind(part)); }

    VectorRef FirstPart(const VectorRef& v) const
	{ return Sub(v,1); }

    VectorRef SecondPart(const VectorRef& v) const
	{ return Sub(v,2); }

    VectorRef LastPart(const VectorRef& v) const
	{ return Sub(v,lastind->Upper()); }
    };

class Partition2D
    {
public:
    Partition1D left;
    Partition1D right;

    Partition2D(int l1, int l2, int r1, int r2) 
    	: left(l1,l2), right(r1,r2) {}

    Partition2D(const MatrixRef& M, int l1, int r1) 
    	: left(l1,M.Nrows()-l1), right(r1,M.Ncols()-r1) {}

    Partition2D(const MatrixRef& M1, const MatrixRef& M2) 
    	: left(M1.Nrows(),M2.Nrows()), right(M1.Ncols(),M2.Ncols()) {}

    Partition2D(const IntArray1& l,const IntArray1& r)
	: left(l), right(r) {}

    MatrixRef Sub(const MatrixRef& m, int l, int r)
	{
	return m.SubMatrix(left.lastind(l-1)+1,left.lastind(l),
			  right.lastind(r-1)+1,right.lastind(r));
	}

    int Nrows() { return left.Length(); }
    int Ncols() { return right.Length(); }

    MatrixRef TopLeft(const MatrixRef& m)
	{ return Sub(m,1,1); }

    MatrixRef TopRight(const MatrixRef& m)
	{ return Sub(m,1,right.NumParts()); }

    MatrixRef BottomRight(const MatrixRef& m)
	{ return Sub(m,left.NumParts(),right.NumParts()); }

    MatrixRef BottomLeft(const MatrixRef& m)
	{ return Sub(m,left.NumParts(),1); }
    };

#endif