packmat.cc.txt

压缩文件中是Error Correction Coding - Mathematical Methods and Algorithms(Wiley 2005)作者:(Todd K. Moon )的配

#include "PackMat.h"

using namespace std;

int PackMat::bitsper = 8*sizeof(PKSIZE);
PKSIZE  PackMat::masker = (1<<bitsper)-1;

int PackMat::shifter;				// amount to shift (e.g., >> 5 to /32)
int *PackMat::wtct = NULL;
int PackMat::PkperCt = sizeof(PKSIZE)/sizeof(CTSIZE);
StackMode PackMat::sidestackmode = Top;
StackMode PackMat::topstackmode = Left;

static PackMat U;				// the upper matrix for the LU decomposition
static int numind;				// number of linearly independent rows
static PackMat L;				// the lower matrix for the LU decomposition
static int *luindex= NULL;
static char printspace[5] = " ";

void
PackMat::initPackMat(void )
{
   int nbits = 8*sizeof(CTSIZE);
   int i,j;

   if(wtct) {					// if already allocated
	  delete[] wtct;
   }
   int N = (1<< (8*sizeof(CTSIZE)));// number of different weights to consider
   int wt;
   wtct = new int[N];
   for(i = 0; i < N; i++) {		// for each of these, find the weight
	  wt = 0;
	  for(j = 0; j < nbits; j++) {
		 if(i & (1<<j)) {
			wt++;
		 }
	  }
	  wtct[i] = wt;
   }
   shifter = 0;
   j = bitsper;
   while(j != 1) {
	  j /= 2;
	  shifter++;
   }
}


PackMat::PackMat(void)
{
   numcols = 0;
   numrows = 0;
   numpackcols = 0;
   numpackrows = 0;
   numpackcols1 = numpackcols2 = 0;
   numpackrows1 = numpackrows2 = 0;
   fracbits = 0;
   packmode = HORIZPACK;
   Matrix = 0;
}


PackMat::PackMat(PackMat &D)
{
   int i,j;
   setsize(*this,D.numrows,D.numcols,D.packmode);
   callocmat(*this);
   for(i = 0; i < D.numpackrows; i++) {
	  for(j = 0; j < D.numpackcols; j++) {
		 Matrix[i][j] = D.Matrix[i][j];
	  }
   }
}

inline void
PackMat::setsize(PackMat &A,int rows, int cols, PackMode pkmode)
{
   int i,j;
   A.numcols = cols;
   A.numrows = rows;
   A.packmode = pkmode;
   if(A.packmode == HORIZPACK) {
	  A.numpackrows = A.numpackrows1 = A.numrows;
	  A.numpackrows2 = 0;

	  A.numpackcols1 = A.numpackcols = A.numcols/bitsper;
	  A.numpackcols2 = 0;
	  A.fracbits = A.numcols%bitsper;
	  if(A.fracbits) {
		 A.numpackcols++;			// increment to hold the fractional int
		 A.numpackcols2 = 1;
	  }
   }
   else if(A.packmode == VERTPACK) {
	  A.numpackcols1 = A.numpackcols = A.numcols;
	  A.numpackcols2 = 0;

	  A.numpackrows1 = A.numpackrows = A.numrows/bitsper;
	  A.numpackrows2 = 0;
	  A.fracbits = A.numrows%bitsper;
	  if(A.fracbits) {
		 A.numpackrows++;
		 A.numpackrows2 = 1;
	  }
   }
   else {
	  cerr<<"Error: Invalid packmode"<<endl;
	  exit(-1);
   }
}

inline void
PackMat::callocmat(PackMat &A)
{
   int i,j;
   CALLOCMATRIX(A.Matrix,PKSIZE,A.numpackrows,A.numpackcols);
   // clear out the matrix
   for(i = 0; i < A.numpackrows; i++) {
	  for(j = 0; j < A.numpackcols; j++) {
		 A.Matrix[i][j] = 0;
	  }
   }
}

void
PackMat::resize(PackMat &A,int rows, int cols, PackMode pkmode)
// set A to the right size, but keep all the old stuff that was there before
{
   PKSIZE **NewMatrix;
   int oldnumpackrows, oldnumpackcols;
   int i,j;

   oldnumpackrows = A.numpackrows;
   oldnumpackcols = A.numpackcols;
   setsize(A,rows,cols,pkmode);

   CALLOCMATRIX(NewMatrix,PKSIZE,A.numpackrows,A.numpackcols);
   // clear out the matrix
   for(i = 0; i < oldnumpackrows; i++) {
	  for(j = 0; j < oldnumpackcols; j++) {
		 NewMatrix[i][j] = A.Matrix[i][j];
	  }
   }
   if(A.Matrix) FREEMATRIX(A.Matrix);
   A.Matrix = NewMatrix;
}



PackMat::PackMat(int rows, int cols, PackMode pkmode)
{
   int i,j;

   double temp1, temp2;

   setsize(*this,rows, cols, pkmode);
   callocmat(*this);
}

int
PackMat::ReadFile(char *fname,PackMode rowcolpack)
{
   ifstream Densefile(fname);
   if(!Densefile) {
	  cerr<<"Error: Unable to open input file "<<fname<<endl;
	  return 0;
   }
   PackMat_is(Densefile,rowcolpack);
   Densefile.close();
   return 1;
}

PackMat::PackMat(char *fname,PackMode rowcolpack)
{
   Matrix = 0;
   ifstream Densefile(fname);
   if(!Densefile) {
	  cerr<<"Error: Unable to open input file"<<fname<<endl;
	  exit(-1);
   }
   PackMat_is(Densefile,rowcolpack);
   Densefile.close();
}

PackMat::PackMat(string desc,PackMode rowcolpack)
{

   Matrix = 0;
   istringstream Densefile(desc);
   PackMat_is(Densefile,rowcolpack);
}


void
PackMat::fill(string desc)
{

   istringstream Densefile(desc);
   fill(Densefile);
}

PackMat::PackMat(istream &Densefile,PackMode rowcolpack)
{
   Matrix = 0;
   PackMat_is(Densefile,rowcolpack);
}

void
PackMat::PackMat_is(istream &Densefile,PackMode rowcolpack)
{

   Densefile>>numrows;
   Densefile>>numcols;
   setsize(*this,numrows, numcols,rowcolpack);
   if(Matrix) FREEMATRIX(Matrix);
   callocmat(*this);

   fill(Densefile);
}

void
PackMat::fill(istream &Densefile)
{
   int i, j, bit, k1;
   PKSIZE blob;

   if(packmode == HORIZPACK) {
	  for(i=0; i<numrows; i++) {
		 for(j=0;j< numpackcols1;j++) {
			blob = 0;
			for(k1 = 0; k1 < bitsper; k1++) {
			   Densefile >> bit;
//cout << bit << " ";
				blob = (bit << k1) | blob;
			}
			Matrix[i][j] = blob;
		 }
		 if(fracbits) {
			j = numpackcols1;
			blob = 0;
			for(k1 = 0; k1 < fracbits; k1++) {
			   Densefile >> bit;
//cout << bit << " ";
			   blob = (bit << k1) | blob;
			}
			Matrix[i][j] = blob;
		 }
	  }
   }
   else if(packmode == VERTPACK) {
	  for(i=0;i< numpackrows1;i++) {
		 for(j = 0; j < numcols; j++) Matrix[i][j] = 0;
		 for(k1 = 0; k1 < bitsper; k1++) {
			for(j = 0; j < numcols; j++) {
			   Densefile >> bit;
			   Matrix[i][j] = (bit << k1) | Matrix[i][j];
			}
		 }
	  }
	  if(fracbits) {
		 i = numpackrows1;
		 for(j = 0; j < numcols; j++) Matrix[i][j] = 0;
		 for(k1 = 0; k1 < fracbits; k1++) {
			for(j = 0; j < numcols; j++) {
			   Densefile >> bit;
			   Matrix[i][j] = (bit << k1) | Matrix[i][j];
			}
		 }
	  }
   }
}


PackMat::PackMat(PackVec &vec, PackMode rowcolpack)
{
   int i;
   if(rowcolpack == HORIZPACK) { // pack horizontally; make row vector
	  setsize(*this,1, vec.numrows,HORIZPACK);
	  callocmat(*this);
	  for(i = 0; i < numpackcols; i++) {
		 Matrix[0][i] = vec.Vec[i];
	  }
   }
   else {						// pack vertically; make a column vector
	  setsize(*this,vec.numrows,1,VERTPACK);
	  callocmat(*this);
	  for(i = 0; i < numpackrows; i++) {
		 Matrix[i][0] = vec.Vec[i];
	  }
   }
}


void
PackMat::Size(int rows, int cols, PackMode pkmode)
// resize, but do not attempt to preserve the contents
{
   int i,j;

   checkandset(*this,pkmode,rows,cols);
   for(i = 0; i < numpackcols; i++) {
	  for(j = 0; j < numpackrows; j++) {
		 Matrix[i][j] = 0;
	  }
   }
}

void
PackMat::Size(PackMat &A)
// resize to the same size as A, but do not attempt to preserve the contents
{
   Size(A.numrows,A.numcols,A.packmode);
}

PackMat &PackMat::operator=(PackMat &DM2)
{
   int i,j;

   checkandset(*this,DM2.packmode,DM2.numrows, DM2.numcols);
   
   for(i=0;i<numpackrows;i++)
	  for(j=0;j<numpackcols;j++)
		 Matrix[i][j] = DM2.Matrix[i][j];

   return *this;
}

PackMat::~PackMat(void)
{
   if(Matrix) FREEMATRIX(Matrix);
}

int PackMat::GetVal(int row, int col)
{
   int i,j,k;
   PKSIZE value;
   int col1, row1;
   int vb;
   int bitcol, bitrow;

   if(row >= numrows || col >= numcols) {
	  cerr<<"Error: Index exceeds matrix dimensions"<<endl;
	  exit(-1);
   }

   value = 0;

   if(packmode == HORIZPACK) {
	  col1 = col >> shifter;	// determine the integer number
	  value = Matrix[row][col1];

	  bitcol = col % bitsper;
	  vb = (value >> bitcol) & 1;
   }
   else if(packmode == VERTPACK) {
	  row1 = row >> shifter;
	  value = Matrix[row1][col];
	  bitrow = row % bitsper;
	  vb =  (value >> bitrow) & 1;
   }
   return vb;
}

void PackMat::SetVal(int row, int col, int val)
{
   int i,j,k, col1, row1;
   int bitcol, bitrow;
   PKSIZE value;
   
   if(row >= numrows || col >= numcols) {
	  cerr<<"Error: Index exceeds matrix dimensions"<<endl;
	  exit(-1);
   }
   
//    if(val != 0 && val != 1) {
// 	  cerr<<"Error: Illegal input value"<<endl;
// 	  cout<<"val = "<<val<<endl;
// 	  exit(-1);
//    }
   
   if(packmode == HORIZPACK) {
	  col1 = col >> shifter;
	  bitcol = col % bitsper;
	  if(!val) {					// set value to 0
		 Matrix[row][col1] &= ~(1<<bitcol);
	  }
	  else {						// set value to 1
		 Matrix[row][col1] |= (1<<bitcol);
	  }
   }
   
   else if(packmode == VERTPACK) {
	  row1 = row >> shifter;
	  bitrow = row % bitsper;
	  if(!val) {				// set value to 0
		 Matrix[row1][col] &= ~(1<<bitrow);
	  }
	  else {
		 Matrix[row1][col] |= (1<<bitrow);
	  }
   }
}



void PackMat::deletecolumn(int colno)
{
   int i,j;
   if(colno < 0 || colno >= numcols) return;
   if(colno == numcols-1) {
	  for(i = 0; i < numrows; i++) {
		 SetVal(i,colno,0);
	  }
   }
   else {
	  for(j = colno; j < numcols-1; j++) {
		 for(i = 0; i < numrows; i++) {
			SetVal(i,j,GetVal(i,j+1));
		 }
	  }
   }
   setsize(*this,numrows,numcols-1,packmode);
}


int PackMat::getrowweight(int rowno)
{
   int i;
   int wt = 0;
   if(rowno < 0 || rowno >= numrows) return 0;

   if(packmode == HORIZPACK) {	// the easy way
	  for(i = 0; i < numpackcols; i++)
		 wt += getWeight(Matrix[rowno][i]);
	  return wt;
   }
   else {
	  for(i = 0; i < numcols; i++) {
		 wt += GetVal(rowno,i);
	  }
   }
   return wt;
}

void PackMat::SetMode(PackMode newmode,PackMat &DM2)
{

   int i,j,k;
   int numr, numc;
   PKSIZE **NewMat;

   if(packmode == newmode) { 	// keep the mode
	  DM2 = *this;
	  return;
   }

   CALLOCMATRIX(NewMat,PKSIZE,numrows,numcols);

   for(i=0;i<numrows;i++) {
	  for(j=0;j<numcols;j++) {
		 k = GetVal(i,j);
		 setval(NewMat,i,j,k,newmode);
	  }
   }
   setsize(DM2,numrows,numcols,newmode);
   if(DM2.Matrix) FREEMATRIX(DM2.Matrix);
   DM2.Matrix = NewMat;
}

void PackMat::Submatrix(PackMat &submat,int row1, int row2, int col1, int col2,
						int subrow, int subcol,PackMode pkmode)
// extract the rectangular matrix, placing it starting at (subrow,subcol)
// in the submatrix, packing in the pkmode order
// if out of range in any direction, fill in with zeros
{
   int i,i1, j,j1, nrows, ncols;
   if(row1 < 0) row1 = 0;
   if(row1 >= numrows) row1 = numrows-1;
   if(col1 < 0) col1 = 0;
   if(col2 >= numcols) col2 = numcols-1;
   nrows = (row2-row1+1) + subrow;
   ncols = (col2-col1+1) + subcol;
   checkandset2(submat,pkmode,nrows,ncols);
   for(i=0, i1=subrow; i < (row2-row1+1); i++,i1++) {
	  for(j = 0, j1 = subcol; j < (col2-col1+1); j++,j1++) {
		 submat.SetVal(i1,j1,GetVal(i+row1,j+col1));
	  }
   }
}

void PackMat::Subvector(PackVec &subvec,int row1, int row2,int col1, int col2,
						int subrow) 
// extract the row or vector and place it (starting at subrow)
// in the vector subvec
{
   int i,i1, j,j1, nrows, ncols;
   if(row1 < 0) row1 = 0;
   if(row1 >= numrows) row1 = numrows-1;
   if(col1 < 0) col1 = 0;
   if(col2 >= numcols) col2 = numcols-1;

   if(row1==row2) {				// extract a row
	  ncols = (col2-col1+1) + subrow;
	  subvec.checkandset2(subvec,ncols);
	  for(i=0, i1=subrow; i < (col2-col1+1); i++,i1++) {
		 subvec.SetVal(i1,GetVal(row1,i+col1));
	  }
   }
   else if(col1==col2) {		// extract a column
	  nrows = (row2-row1+1) + subrow;
	  subvec.checkandset2(subvec,nrows);
	  for(i=0, i1=subrow; i < (row2-row1+1); i++,i1++) {
		 subvec.SetVal(i1,GetVal(i+row1,col1));
	  }
   }
   else {
	  cerr << "Either row or column must be constant\n";
	  exit(-1);
   }
}


void PackMat::Mult(PackMat &DM, PackMat &Prod)
{
// This function tests to make sure things are the right size
// and then calls the multiplication function

   if(packmode == HORIZPACK && DM.packmode == VERTPACK) {
	  if(numcols != DM.numrows) {  // correct size to multiply
		 cerr << "Error: Matrices not conformable in product\n";
		 exit(-1);
	  }
	  checkandset(Prod,HORIZPACK,numrows,DM.numcols);
	  Multnocheck(DM,Prod);
   }
   else {
	  cerr << "Multiplication not implemented for these matrix modes\n";
	  exit(-1);
   }
}


void PackMat::Multnocheck(PackMat &DM, PackMat &Prod)
{
// This function does not testing --- it assumes that
// everything is the right size and shape for multiplication

   int i,j,k,nj,j1,k1;
   int prodwt;

   nj = DM.numcols/bitsper;
   for(i = 0; i < numrows; i++) {
// cout << i << " " << flush;
	  // do all the ones that have the full 32 bits
	  j1 = 0;
	  for(j = 0; j < nj; j++) {
		 Prod.Matrix[i][j] = 0;
		 for(k1 = 0; k1 < bitsper; k1++) {
			prodwt = 0;
			for(k = 0; k < numpackcols;k ++) {
			   prodwt += getWeight(Matrix[i][k] & DM.Matrix[k][j1]);
			}
			Prod.Matrix[i][j] = ((prodwt % 2)<<k1) + Prod.Matrix[i][j];
			j1++;			   
		 }
	  }
	  // do the last fractbits
	  if(Prod.fracbits) {
		 j = nj;
		 Prod.Matrix[i][j] = 0;
		 for(k1 = 0; k1 < Prod.fracbits; k1++) {
			prodwt = 0;
			for(k = 0; k < numpackcols;k ++) {
			   prodwt += getWeight(Matrix[i][k] & DM.Matrix[k][j1]);
			}
			Prod.Matrix[i][j] = ((prodwt % 2)<<k1) + Prod.Matrix[i][j]; 
			j1++;
		 }
	  }
   }
}


void PackMat::Mult(PackVec &DM, PackVec &Prod)
// Multiply a matrix times a vector
{
// This function tests to make sure things are the right size
// and then calls the multiplication function

   if(packmode == HORIZPACK) {
	  if(numcols != numrows) {  // correct size to multiply
		 cerr << "Error: Matrices not conformable in product\n";
		 exit(-1);
	  }
	  Prod.checkandset(Prod,numrows);
	  Multnocheck(DM,Prod);
   }
   else {
	  cerr << "Multiplication not implemented for these matrix modes\n";
	  exit(-1);
   }
}


void PackMat::Multnocheck(PackVec &DM, PackVec &Prod)
{
// This function does not testing --- it assumes that
// everything is the right size and shape for multiplication

   int i,k,ibig,ilit;
   int prodwt;

   for(i = 0; i < DM.numpackrows; i++) {Prod.Vec[i] = 0;}
   for(i = 0; i < numrows; i++) {
	  ibig = i/bitsper;
	  ilit = i % bitsper;
	  prodwt = 0;
	  for(k = 0; k < numpackcols;k ++) {
		 prodwt += getWeight(Matrix[i][k] & DM.Vec[k]);
	  }
	  Prod.Vec[ibig] = ((prodwt % 2)<<ilit) + Prod.Vec[ibig];
   }
}


inline void
PackMat::checkandset(PackMat &B,PackMode bpackmode,int brows, int bcols)
// set up exactly brows and bcols in B, if it is not already there
{
   if(!(B.packmode == bpackmode && B.numrows == brows && B.numcols == bcols)) {
	  if(B.Matrix) {FREEMATRIX(B.Matrix);}
	  setsize(B,brows, bcols, bpackmode);
	  callocmat(B);
   }	  
}

inline void
PackMat::checkandset2(PackMat &B,PackMode bpackmode,int brows, int bcols)
// make sure there is at least brows and bcols space,
// but possibly more
{
   if(B.packmode != bpackmode || brows > B.numrows || bcols > B.numcols) {
	  // wrong mode, or not enough room already
	  if(B.Matrix) {FREEMATRIX(B.Matrix);}
	  setsize(B,MAX(brows,B.numrows), MAX(bcols,B.numcols), bpackmode);
	  callocmat(B);
   }	  
}

void
PackMat::checkandresize2(PackMat &B,PackMode bpackmode,int brows, int bcols)
// make sure there is at least brows and bcols space,
// but possibly more
// If necessary, resize, but keep the old stuff
{
   if(B.packmode != bpackmode || brows > B.numrows || bcols > B.numcols) {
	  // wrong mode, or not enough room already
//	  if(B.Matrix) {FREEMATRIX(B.Matrix);}
	  resize(B,MAX(brows,B.numrows), MAX(bcols,B.numcols), bpackmode);
   }	  
}




void PackMat::transpose(PackMat &T,PackMode newpackmode)
{
   int numc, numr;
   PKSIZE **NewMat;

   numc = numcols;
   numr = numrows;

   CALLOCMATRIX(NewMat,PKSIZE,numcols, numrows);