/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Emc/EmcCalib/EmcBhaCalib/EmcBhaCalib-00-00-34/src/EmcLSSMatrix.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Description:
//      Class EmcLSSMatrix - Implementation of a Large Sparse Symmetric Matrix
//
// Environment:
//      Software developed for the BESIII Detector at the BEPCII.
//
// Author List:
//      Chunxiu Liu
//
// Copyright Information:
//      Copyright (C) 2005           IHEP
//
//------------------------------------------------------------------------

//-----------------------
// This Class's Header --
//-----------------------
#include "EmcBhaCalib/EmcLSSMatrix.h"

//-------------
// C Headers --
//-------------
extern "C" {
}

//---------------
// C++ Headers --
//---------------

#include <fstream>

using namespace std;

//              ----------------------------------------
//              -- Public Function Member Definitions --
//              ----------------------------------------

//----------------
// Constructors --
//----------------
EmcLSSMatrix::EmcLSSMatrix()
  : _size(0), _nrrows(0), _nrcol(0), _matrix(0), _columns(0), _rows(0),
    _nothing(0),_verb(false)  
{
}

EmcLSSMatrix::EmcLSSMatrix( int rows, int nonzero_col)
  : _size(rows * nonzero_col),
    _nrrows(rows),
    _nrcol(nonzero_col)
{
  _matrix = new double[_size];
  _rows = new int[_size];
  _columns = new int[_size];
  for (int i=0; i<_size;i++) 
    {
      _matrix[i] =0.;
      _rows[i] = -1;
      _columns[i] = -1;
    }
  _nothing = new double;
  _verb = false;

  _MsgFlag =0;
}

EmcLSSMatrix::EmcLSSMatrix( int rows, int nonzero_col, int MsgFlag)
  : _size(rows * nonzero_col),
    _nrrows(rows),
    _nrcol(nonzero_col)
{
  _matrix = new double[_size];
  _rows = new int[_size];
  _columns = new int[_size];
  for (int i=0; i<_size;i++) 
    {
      _matrix[i] =0.;
      _rows[i] = -1;
      _columns[i] = -1;
    }
  _nothing = new double;
  _verb = false;

  _MsgFlag = MsgFlag;
}

EmcLSSMatrix::EmcLSSMatrix( const EmcLSSMatrix &m1 ) 
  : _size(m1._size),
    _nrrows(m1._nrrows),
    _nrcol(m1._nrcol)

{
  _matrix = new double[_size]; 
  _rows = new int[_size];
  _columns = new int[_size];
  _nothing = m1._nothing;
  _verb = m1._verb;

  for ( long int i=0;i<_size;i++) 
    {
      _matrix[i] = m1._matrix[i];
      _rows[i] = m1._rows[i];
      _columns[i] = m1._columns[i];
    };
}

//--------------
// Destructor --
//--------------
EmcLSSMatrix::~EmcLSSMatrix() {
  if ( 0 != _matrix) 
    {
      delete [] _matrix;
      _matrix = 0;
    }
  if ( 0 != _rows) 
    {
      delete [] _rows;
      _rows = 0;
    }
  if ( 0 != _columns) 
    {
      delete [] _columns;
      _columns = 0;
    }
  if ( 0 != _nothing) 
    {
      delete _nothing;
      _nothing = 0;
    }
}

//-------------
// Operators --
//-------------
double&
EmcLSSMatrix::operator()(int row, int col) {

  long int found = find( row, col );
  
  
  if ( -1 == found ) 
    {
      if (_MsgFlag <= 5) {
        std::cout << "EmcLSSMatrix:: ERROR "
                  << "EmcLSSMatrix: matrix element not found !!" << endl
                  << "EmcLSSMatrix: Return zero !" << endl;
      }
      _nothing = 0;
      return *_nothing;
    }
  
  return (*(_matrix+found));
}


    
//-------------
// Selectors --
//-------------
long int
EmcLSSMatrix::find( int row, int col) 
{  
  int smaller,larger;
  
 
  if ( col > row ) 
    {
      smaller = row;
      larger = col; 
    }
  else 
    {
      smaller = col;
      larger = row;
    }
  
  int* col_p;
  int* row_p;
  double* matr_p;
  double* matr_row_max;

  if (larger < 0 || larger > _nrrows || smaller < 0 || smaller > _nrrows ) 
    {
      if (_MsgFlag <= 5) {
        std::cout <<"EmcLSSMatrix::ERROR"    
                  << "!!! ERROR in bound check of EmcLSSMatrix !!!"
                  << "!!! Return zero !!!" 
                  << endl;
      }
      matr_p = 0;
    }
  else 
    {     
      col_p = (_columns + (larger * _nrcol));
      row_p = (_rows + (larger * _nrcol));
      matr_p = (_matrix + (larger * _nrcol));
      matr_row_max = (matr_p + _nrcol);
      
      while ( matr_p < matr_row_max ) {
        
        if (_MsgFlag <= 1) {
          std::cout <<"EmcLSSMatrix::VERBOSE" 
                    << "C: " << larger << " " << smaller << " "
                    << col_p << " " << *col_p << " " 
                    << matr_p << " " << *matr_p << " "
                    << (_matrix+(matr_p-_matrix)) << " "
                    << *(_matrix+(matr_p-_matrix)) << endl;
        }
        /*
          cout<< "C: " << larger << " " << smaller << " "
          << col_p << " " << *col_p << " " 
          << matr_p << " " << *matr_p << " "
          << (_matrix+(matr_p-_matrix)) << " "
          << *(_matrix+(matr_p-_matrix)) <<endl;
        */
        if ( matr_p == (matr_row_max-1) ) 
          {
            if (_MsgFlag <= 4) {
              std::cout <<"EmcLSSMatrix::WARNING " 
                        << "!! WARNING: Reached maximum number of columns "
                        << "in LSSMatrix when searching for row " 
                        << larger << " column " << smaller << " !!" 
                        << endl
                        << "!! Return zero pointer !! " << endl;
            }
            matr_p = 0;
            break;
          }

        //if row does already exist
        if ( *col_p == smaller ) 
          {
            break;
          }

        //if at the end of the list, use this as new element
        if ( (*matr_p == 0.) ) 
          {
            *col_p = smaller;
            *row_p = larger;
            //  if (*row_p == 1616 ) {
            //    nun=(row_p-(_rows + (larger * _nrcol)));
            //if (_MsgFlag <= 2) {
            //  std::cout <<"EmcLSSMatrix::DEBUG " 
            //          << nun << " " << larger << " " << smaller  << " "
            //          << *row_p << " " << *col_p << " " << *matr_p << endl;
            //}
            //  }
            break;
          }
      
        matr_p++;
        col_p++;
        row_p++;
        
      }
    }
  
  long int diff = matr_p-_matrix;

  if (matr_p == 0) 
    {
      diff = -1;
    }
  
  return (diff);
}

double*
EmcLSSMatrix::matrix( const int& rowind ) const 
{
  double* here=0;
  
  if (rowind < _nrrows) 
    {
      here = (_matrix+(rowind*_nrcol));
    }
  else 
    {
      if (_MsgFlag <= 5) {
        std::cout <<"EmcLSSMatrix::ERROR " 
                  << "EmcLSSMatrix::matrix: Error "
                  << "- larger row index than existing requested !" 
                  << endl;
            }
      here = 0;
    }
  return here;
}

int*
EmcLSSMatrix::column( const int& rowind ) const 
{
  int* here=0;
   
  if (rowind < _nrrows) 
    {
      here =  (_columns+(rowind*_nrcol));
    }
  else 
    {
      
      if (_MsgFlag <= 5) {
        std::cout <<"EmcLSSMatrix::ERROR " 
                  << "EmcLSSMatrix.column: Error "
                  << "- larger row index than existing requested !" 
                  << endl;
      }
      here = 0;
    }
  return here;
}

int
EmcLSSMatrix::num_filled_cols( const int row ) const 
{
  double * search_i = _matrix + ( row * _nrcol );
  double * max_i = search_i + _nrcol;
  int nonZeroCol=0;
  
  while ( (*search_i != 0.) && (search_i < max_i) ) 
    {
      nonZeroCol++;
      search_i++;
    }
  return nonZeroCol;
} 

int
EmcLSSMatrix::num_filled_rows( const int col ) const 
{
  int nonZeroRows = 0;
  for ( long int i=0; i<_size; i++ )
    { 
      if ( (_matrix[i] != 0.) && (_columns[i] == col) ) 
        {
          nonZeroRows++;  
        }
    }

  return nonZeroRows;
}


long int
EmcLSSMatrix::num_nonZeros() 
{
  int* col_p = _columns;
  double* ele_p = _matrix;
  double* mat_max_p = (_matrix + _size);
  long int nrele = 0;
  
  while ( ele_p < mat_max_p ) 
    {
      if ( *ele_p != 0.) nrele++;
      col_p++;
      ele_p++;
    }
  
  return nrele;
}

//-------------
// Modifiers --
//-------------
void
EmcLSSMatrix::reset() 
{
  for (int i=0; i<_size;i++) 
    {
      _matrix[i] =0.;
      _rows[i] = -1;
      _columns[i] = -1;
    }

}

bool
EmcLSSMatrix::reduce_Matrix(int* xRef_list) 
{
  bool successful = true;

  //delete all zero elements in matrix
  //save only non zero elements

  long int _newIndx = 0;
  
  
  for ( long int _arrayIndx = 0; 
        _arrayIndx < _size; _arrayIndx++) 
    {
      
      //add 1 to matrix indices because SLAP wants indices 1..N,
      //but Xtals counting in geometry starts with 0
      
      if ( _matrix[_arrayIndx] > 0. ) 
        {
          
          if ( (xRef_list[(_rows[_arrayIndx])]) >= 0
               && (xRef_list[(_columns[_arrayIndx])]) >= 0 ) 
            {         
              _matrix[_newIndx] = _matrix[_arrayIndx];
              _rows[_newIndx] = ((xRef_list[(_rows[_arrayIndx])])+1);
              _columns[_newIndx] = ((xRef_list[(_columns[_arrayIndx])])+1);       
              _newIndx++;
            } 
          else 
            {

              //int indxtal;
              if (xRef_list[(_rows[_arrayIndx])] < 0 ) 
                {
                  if (_MsgFlag <= 5) {
                    std::cout <<"EmcLSSMatrix::ERROR " 
                              << "EmcLSSMatrix:  Xtal index "
                              << _rows[_arrayIndx]
                              << " appears in matrix, " 
                              << "but not in vector !!! "
                              << _rows[_arrayIndx] << " "
                              << _columns[_arrayIndx]
                              << endl;
                  }
                } 
              else 
                {
                  
                  if (_MsgFlag <= 5) {
                    std::cout <<"EmcLSSMatrix::ERROR " 
                              << "EmcLSSMatrix:  Xtal index "
                              << _columns[_arrayIndx]
                              << " appears in matrix, " 
                              << "but not in vector !!! "
                              << _rows[_arrayIndx] << " "
                              << _columns[_arrayIndx]
                              << endl;
                  }
                }
              successful=false;
            }
        }
    }

  if (_verb) 
    {
      if (_MsgFlag <= 2) {
        std::cout <<"EmcLSSMatrix::DEBUG " 
                  << "Reduced LSSMatrix !!! Number of non zeros: " 
                  << _newIndx << endl;
      }
    }
  
  return successful;
}


void
EmcLSSMatrix::print_NonZeros() 
{
  int* col_p = _columns;
  int* row_p = _rows;
  double* ele_p = _matrix;
  double* mat_max_p = (_matrix + _size);
  long int wo = 0;
  long int nrele = 0;


  if (_MsgFlag <= 2) {
    std::cout <<"EmcLSSMatrix::DEBUG "  
              <<"List of nonzero Matrix elements (Matrix size: " 
              << _size << " ) : " << endl;
  }
  
  while ( ele_p < mat_max_p ) 
    {
      if ( *ele_p != 0. ) 
        {
          nrele++;
          if (_MsgFlag <= 2) {
            std::cout <<"EmcLSSMatrix::DEBUG " 
                      << "nr: " << nrele
                      << " M( " <<  *row_p << " ,  " << *col_p
                      << " ):  " << *ele_p 
                      << "   array index: " << wo << endl;
          }
          /*  
              cout<< "nr: " << nrele
              << " M( " <<  *row_p << " ,  " << *col_p
              << " ):  " << *ele_p 
              << "   array index: " << wo << endl;
          */
        }
      wo++;
      col_p++;
      row_p++;
      ele_p++;
    }
  
}


void
EmcLSSMatrix::print_row(int rownr) 
{
  int* col_p =  _columns+(_nrcol*rownr);
  int* row_p = _rows+(_nrcol*rownr);
  double* ele_p = _matrix+(_nrcol*rownr);
  double* mat_max_p = (ele_p + _nrcol);
  long int wo = 0;
  long int nrele = 0;
  
  if (_MsgFlag <= 2) {
    std::cout <<"EmcLSSMatrix::DEBUG " 
              <<"row length: " << num_filled_cols(rownr) << endl;
  }
  while ( ele_p < mat_max_p ) 
    {
      if ( *ele_p != 0. ) 
        {
          nrele++;
          if (_MsgFlag <= 2) {
            std::cout <<"EmcLSSMatrix::DEBUG " 
                      << "nr: " << nrele
                      << " M( " <<  *row_p << " ,  " << *col_p
                      << " ):  " << *ele_p 
                      << "   array index: " << wo << endl;
          }       
          /*
            cout<< "nr: " << nrele
            << " M( " <<  *row_p << " ,  " << *col_p
            << " ):  " << *ele_p 
            << "   array index: " << wo << endl;
          */
        }
      wo++;
      col_p++;
      row_p++;
      ele_p++;
    }
  
}

void
EmcLSSMatrix::writeOut( ostream& Out) 
{
  int* col_p = _columns;
  int* row_p =  _rows;
  double* ele_p = _matrix;
  double* mat_max_p = (_matrix + _size);
  //  long int nrele = 0;

  long int nonz = num_nonZeros();
  Out << nonz << " ";
  
  while ( ele_p < mat_max_p ) 
    {
      if ( *ele_p != 0.) 
        {
          Out<< *row_p << " " 
             << *col_p << " " 
             << *ele_p << " ";
        }
      col_p++;
      row_p++;
      ele_p++;
    }
  
  if (_MsgFlag <= 2) {
    std::cout <<"EmcLSSMatrix::DEBUG " 
              << "Wrote " << nonz 
              << " non zero matrix elements to file !!!" << endl;
  }
  
}


void
EmcLSSMatrix::readIn( istream& In) 
{  
  //  long int nonz = num_nonZeros();
  long int nonz = 0;
  In >> nonz;
  
  cout<<"nonz="<<nonz<<endl;
  
  int theRow;
  int theCol;
  long int index;
  double theEle;
  

  for (long int i=0; i<nonz; i++ ) 
    {
      In >> theRow >> theCol >> theEle;
      index = find (theRow,theCol);
      /*
      cout<<"index = "<<index
          <<"row="<<theRow
          <<"col="<<theCol<<endl;
      */
      if ( -1 != index ) 
        {
          _matrix[index] += theEle;
          if (_verb) 
            {
              if ( i < 50 || i > (nonz-10) )
                {
                  if (_MsgFlag <= 2) {
                    std::cout <<"EmcLSSMatrix::DEBUG " 
                              << "M: " << _rows[index] << " " << _columns[index]
                              << " " << _matrix[index] << " " << index << endl;
                  }   
                }
            }
        }
    }
  
  if (_verb) 
    {
      if (_MsgFlag <= 2) {
        std::cout <<"EmcLSSMatrix::DEBUG " 
                  << "Have read in " << nonz 
                  << " non zero matrix elements from file !!!" << endl;
            }
    }
}

---