Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/parser.cpp

    r042f82 r29812d  
    88
    99#include "helpers.hpp"
     10#include "memoryallocator.hpp"
    1011#include "parser.hpp"
    1112
    … …  
    5657MatrixContainer::MatrixContainer() {
    5758  Indices = NULL;
    58   Header = (char *) Malloc(sizeof(char)*1023, "MatrixContainer::MatrixContainer: *Header");
    59   Matrix = (double ***) Malloc(sizeof(double **)*(1), "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
    60   RowCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *RowCounter");
     59  Header = Malloc<char*>(1, "MatrixContainer::MatrixContainer: **Header");
     60  Matrix = Malloc<double**>(1, "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
     61  RowCounter = Malloc<int>(1, "MatrixContainer::MatrixContainer: *RowCounter");
     62  ColumnCounter = Malloc<int>(1, "MatrixContainer::MatrixContainer: *ColumnCounter");
     63  Header[0] = NULL;
    6164  Matrix[0] = NULL;
    6265  RowCounter[0] = -1;
    6366  MatrixCounter = 0;
    64   ColumnCounter = 0;
     67  ColumnCounter[0] = -1;
    6568};
    6669
    … …  
    7073  if (Matrix != NULL) {
    7174    for(int i=MatrixCounter;i--;) {
    72       if (RowCounter != NULL) {
     75      if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
    7376          for(int j=RowCounter[i]+1;j--;)
    74             Free((void **)&Matrix[i][j], "MatrixContainer::~MatrixContainer: *Matrix[][]");
    75         Free((void **)&Matrix[i], "MatrixContainer::~MatrixContainer: **Matrix[]");
     77            Free(&Matrix[i][j]);
     78        Free(&Matrix[i]);
    7679      }
    7780    }
    78     if ((RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
     81    if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
    7982      for(int j=RowCounter[MatrixCounter]+1;j--;)
    80         Free((void **)&Matrix[MatrixCounter][j], "MatrixContainer::~MatrixContainer: *Matrix[MatrixCounter][]");
     83        Free(&Matrix[MatrixCounter][j]);
    8184    if (MatrixCounter != 0)
    82       Free((void **)&Matrix[MatrixCounter], "MatrixContainer::~MatrixContainer: **Matrix[MatrixCounter]");
    83     Free((void **)&Matrix, "MatrixContainer::~MatrixContainer: ***Matrix");
     85      Free(&Matrix[MatrixCounter]);
     86    Free(&Matrix);
    8487  }
    8588  if (Indices != NULL)
    8689    for(int i=MatrixCounter+1;i--;) {
    87       Free((void **)&Indices[i], "MatrixContainer::~MatrixContainer: *Indices[]");
    88     }
    89   Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices");
    90 
    91   Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header");
    92   Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
    93 };
    94 
     90      Free(&Indices[i]);
     91    }
     92  Free(&Indices);
     93 
     94  if (Header != NULL)
     95    for(int i=MatrixCounter+1;i--;)
     96      Free(&Header[i]);
     97  Free(&Header);
     98  Free(&RowCounter);
     99  Free(&ColumnCounter);
     100};
     101
     102/** Either copies index matrix from another MatrixContainer or initialises with trivial mapping if NULL.
     103 * This either copies the index matrix or just maps 1 to 1, 2 to 2 and so on for all fragments.
     104 * \param *Matrix pointer to other MatrixContainer
     105 * \return true - copy/initialisation sucessful, false - dimension false for copying
     106 */
     107bool MatrixContainer::InitialiseIndices(class MatrixContainer *Matrix)
     108{
     109  cout << "Initialising indices";
     110  if (Matrix == NULL) {
     111    cout << " with trivial mapping." << endl;
     112    Indices = Malloc<int*>(MatrixCounter + 1, "MatrixContainer::InitialiseIndices: **Indices");
     113    for(int i=MatrixCounter+1;i--;) {
     114      Indices[i] = Malloc<int>(RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
     115      for(int j=RowCounter[i];j--;)
     116        Indices[i][j] = j;
     117    }
     118  } else {
     119    cout << " from other MatrixContainer." << endl;
     120    if (MatrixCounter != Matrix->MatrixCounter)
     121      return false;
     122    Indices = Malloc<int*>(MatrixCounter + 1, "MatrixContainer::InitialiseIndices: **Indices");
     123    for(int i=MatrixCounter+1;i--;) {
     124      if (RowCounter[i] != Matrix->RowCounter[i])
     125        return false;
     126      Indices[i] = Malloc<int>(Matrix->RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
     127      for(int j=Matrix->RowCounter[i];j--;) {
     128        Indices[i][j] = Matrix->Indices[i][j];
     129        //cout << Indices[i][j] << "\t";
     130      }
     131      //cout << endl;
     132    }
     133  }
     134  return true;
     135};
    95136
    96137/** Parsing a number of matrices.
    … …  
    121162  }
    122163
    123   // skip some initial lines
     164  // parse header
     165  Header[MatrixNr] = Malloc<char>(1024, "MatrixContainer::ParseMatrix: *Header[]");
    124166  for (int m=skiplines+1;m--;)
    125     input.getline(Header, 1023);
    126 
     167    input.getline(Header[MatrixNr], 1023);
     168 
    127169  // scan header for number of columns
    128   line.str(Header);
     170  line.str(Header[MatrixNr]);
    129171  for(int k=skipcolumns;k--;)
    130     line >> Header;
     172    line >> Header[MatrixNr];
    131173  //cout << line.str() << endl;
    132   ColumnCounter=0;
     174  ColumnCounter[MatrixNr]=0;
    133175  while ( getline(line,token, '\t') ) {
    134176    if (token.length() > 0)
    135       ColumnCounter++;
     177      ColumnCounter[MatrixNr]++;
    136178  }
    137179  //cout << line.str() << endl;
    138   //cout << "ColumnCounter: " << ColumnCounter << "." << endl;
    139   if (ColumnCounter == 0)
    140     cerr << "ColumnCounter: " << ColumnCounter << " from file " << name << ", this is probably an error!" << endl;
    141 
     180  //cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
     181  if (ColumnCounter[MatrixNr] == 0)
     182    cerr << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
     183 
    142184  // scan rest for number of rows/lines
    143185  RowCounter[MatrixNr]=-1;    // counts one line too much
    … …  
    155197
    156198  // allocate matrix if it's not zero dimension in one direction
    157   if ((ColumnCounter > 0) && (RowCounter[MatrixNr] > -1)) {
    158     Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    159 
     199  if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
     200    Matrix[MatrixNr] = Malloc<double*>(RowCounter[MatrixNr] + 1, "MatrixContainer::ParseMatrix: **Matrix[]");
     201 
    160202    // parse in each entry for this matrix
    161203    input.clear();
    162204    input.seekg(ios::beg);
    163205    for (int m=skiplines+1;m--;)
    164       input.getline(Header, 1023);    // skip header
    165     line.str(Header);
     206      input.getline(Header[MatrixNr], 1023);    // skip header
     207    line.str(Header[MatrixNr]);
    166208    for(int k=skipcolumns;k--;)  // skip columns in header too
    167209      line >> filename;
    168     strncpy(Header, line.str().c_str(), 1023);
     210    strncpy(Header[MatrixNr], line.str().c_str(), 1023); 
    169211    for(int j=0;j<RowCounter[MatrixNr];j++) {
    170       Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     212      Matrix[MatrixNr][j] = Malloc<double>(ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[][]");
    171213      input.getline(filename, 1023);
    172214      stringstream lines(filename);
    … …  
    174216      for(int k=skipcolumns;k--;)
    175217        lines >> filename;
    176       for(int k=0;(k<ColumnCounter) && (!lines.eof());k++) {
     218      for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
    177219        lines >> Matrix[MatrixNr][j][k];
    178220        //cout << " " << setprecision(2) << Matrix[MatrixNr][j][k];;
    179221      }
    180222      //cout << endl;
    181       Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
    182       for(int j=ColumnCounter;j--;)
     223      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = Malloc<double>(ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[RowCounter[MatrixNr]][]");
     224      for(int j=ColumnCounter[MatrixNr];j--;)
    183225        Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
    184226    }
    185227  } else {
    186     cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl;
     228    cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl;
    187229  }
    188230  input.close();
    … …  
    233275
    234276  cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
    235   Matrix = (double ***) ReAlloc(Matrix, sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
    236   RowCounter = (int *) ReAlloc(RowCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
     277  Header = ReAlloc<char*>(Header, MatrixCounter + 1, "MatrixContainer::ParseFragmentMatrix: **Header"); // one more each for the total molecule
     278  Matrix = ReAlloc<double**>(Matrix, MatrixCounter + 1, "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
     279  RowCounter = ReAlloc<int>(RowCounter, MatrixCounter + 1, "MatrixContainer::ParseFragmentMatrix: *RowCounter");
     280  ColumnCounter = ReAlloc<int>(ColumnCounter, MatrixCounter + 1, "MatrixContainer::ParseFragmentMatrix: *ColumnCounter");
    237281  for(int i=MatrixCounter+1;i--;) {
    238282    Matrix[i] = NULL;
     283    Header[i] = NULL;
    239284    RowCounter[i] = -1;
     285    ColumnCounter[i] = -1;
    240286  }
    241287  for(int i=0; i < MatrixCounter;i++) {
    … …  
    246292    if (!ParseMatrix(file.str().c_str(), skiplines, skipcolumns, i))
    247293      return false;
    248     Free((void **)&FragmentNumber, "MatrixContainer::ParseFragmentMatrix: *FragmentNumber");
     294    Free(&FragmentNumber);
    249295  }
    250296  return true;
    … …  
    252298
    253299/** Allocates and resets the memory for a number \a MCounter of matrices.
    254  * \param *GivenHeader Header line
     300 * \param **GivenHeader Header line for each matrix
    255301 * \param MCounter number of matrices
    256302 * \param *RCounter number of rows for each matrix
    257  * \param CCounter number of columns for all matrices
    258  * \return Allocation successful
    259  */
    260 bool MatrixContainer::AllocateMatrix(const char *GivenHeader, int MCounter, int *RCounter, int CCounter)
    261 {
    262   Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseFragmentMatrix: *EnergyHeader");
    263   strncpy(Header, GivenHeader, 1023);
    264 
     303 * \param *CCounter number of columns for each matrix
     304 * \return Allocation successful
     305 */
     306bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
     307{
    265308  MatrixCounter = MCounter;
    266   ColumnCounter = CCounter;
    267   Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
    268   RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
     309  Header = Malloc<char*>(MatrixCounter + 1, "MatrixContainer::AllocateMatrix: *Header");
     310  Matrix = Malloc<double**>(MatrixCounter + 1, "MatrixContainer::AllocateMatrix: ***Matrix"); // one more each for the total molecule
     311  RowCounter = Malloc<int>(MatrixCounter + 1, "MatrixContainer::AllocateMatrix: *RowCounter");
     312  ColumnCounter = Malloc<int>(MatrixCounter + 1, "MatrixContainer::AllocateMatrix: *ColumnCounter");
    269313  for(int i=MatrixCounter+1;i--;) {
     314    Header[i] = Malloc<char>(1024, "MatrixContainer::AllocateMatrix: *Header[i]");
     315    strncpy(Header[i], GivenHeader[i], 1023);
    270316    RowCounter[i] = RCounter[i];
    271     Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     317    ColumnCounter[i] = CCounter[i];
     318    Matrix[i] = Malloc<double*>(RowCounter[i] + 1, "MatrixContainer::AllocateMatrix: **Matrix[]");
    272319    for(int j=RowCounter[i]+1;j--;) {
    273       Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    274       for(int k=ColumnCounter;k--;)
     320      Matrix[i][j] = Malloc<double>(ColumnCounter[i], "MatrixContainer::AllocateMatrix: *Matrix[][]");
     321      for(int k=ColumnCounter[i];k--;)
    275322        Matrix[i][j][k] = 0.;
    276323    }
    … …  
    286333  for(int i=MatrixCounter+1;i--;)
    287334    for(int j=RowCounter[i]+1;j--;)
    288       for(int k=ColumnCounter;k--;)
     335      for(int k=ColumnCounter[i];k--;)
    289336        Matrix[i][j][k] = 0.;
    290337   return true;
    … …  
    299346  for(int i=MatrixCounter+1;i--;)
    300347    for(int j=RowCounter[i]+1;j--;)
    301       for(int k=ColumnCounter;k--;)
     348      for(int k=ColumnCounter[i];k--;)
    302349        if (fabs(Matrix[i][j][k]) > max)
    303350          max = fabs(Matrix[i][j][k]);
    … …  
    315362  for(int i=MatrixCounter+1;i--;)
    316363    for(int j=RowCounter[i]+1;j--;)
    317       for(int k=ColumnCounter;k--;)
     364      for(int k=ColumnCounter[i];k--;)
    318365        if (fabs(Matrix[i][j][k]) < min)
    319366          min = fabs(Matrix[i][j][k]);
    … …  
    331378{
    332379  for(int j=RowCounter[MatrixCounter]+1;j--;)
    333     for(int k=skipcolumns;k<ColumnCounter;k++)
     380    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
    334381      Matrix[MatrixCounter][j][k] = value;
    335382   return true;
    … …  
    344391{
    345392  for(int j=RowCounter[MatrixCounter]+1;j--;)
    346     for(int k=skipcolumns;k<ColumnCounter;k++)
     393    for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
    347394      Matrix[MatrixCounter][j][k] = values[j][k];
    348395   return true;
    349396};
    350397
    351 /** Sums the energy with each factor and put into last element of \a ***Energies.
     398/** Sums the entries with each factor and put into last element of \a ***Matrix.
    352399 * Sums over "E"-terms to create the "F"-terms
    353  * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
     400 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
    354401 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    355402 * \param Order bond order
    … …  
    384431              }
    385432              if (Order == SubOrder) { // equal order is always copy from Energies
    386                 for(int l=ColumnCounter;l--;) // then adds/subtract each column
     433                for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column
    387434                  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
    388435              } else {
    389                 for(int l=ColumnCounter;l--;)
     436                for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;)
    390437                  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
    391438              }
    392439            }
    393             //if ((ColumnCounter>1) && (RowCounter[0]-1 >= 1))
     440            //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
    394441             //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
    395442          }
    … …  
    420467    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    421468    line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
    422     Free((void **)&FragmentNumber, "*FragmentNumber");
     469    Free(&FragmentNumber);
    423470    output.open(line.str().c_str(), ios::out);
    424471    if (output == NULL) {
    … …  
    426473      return false;
    427474    }
    428     output << Header << endl;
     475    output << Header[i] << endl;
    429476    for(int j=0;j<RowCounter[i];j++) {
    430       for(int k=0;k<ColumnCounter;k++)
     477      for(int k=0;k<ColumnCounter[i];k++)
    431478        output << scientific << Matrix[i][j][k] << "\t";
    432479      output << endl;
    … …  
    455502    return false;
    456503  }
    457   output << Header << endl;
     504  output << Header[MatrixCounter] << endl;
    458505  for(int j=0;j<RowCounter[MatrixCounter];j++) {
    459     for(int k=0;k<ColumnCounter;k++)
     506    for(int k=0;k<ColumnCounter[MatrixCounter];k++)
    460507      output << scientific << Matrix[MatrixCounter][j][k] << "\t";
    461508    output << endl;
    … …  
    474521{
    475522  cout << "Parsing energy indices." << endl;
    476   Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "EnergyMatrix::ParseIndices: **Indices");
     523  Indices = Malloc<int*>(MatrixCounter + 1, "EnergyMatrix::ParseIndices: **Indices");
    477524  for(int i=MatrixCounter+1;i--;) {
    478     Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "EnergyMatrix::ParseIndices: *Indices[]");
     525    Indices[i] = Malloc<int>(RowCounter[i], "EnergyMatrix::ParseIndices: *Indices[]");
    479526    for(int j=RowCounter[i];j--;)
    480527      Indices[i][j] = j;
    … …  
    485532/** Sums the energy with each factor and put into last element of \a EnergyMatrix::Matrix.
    486533 * Sums over the "F"-terms in ANOVA decomposition.
    487  * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
     534 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
    488535 * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
    489536 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    … …  
    498545    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
    499546      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
    500         for(int k=ColumnCounter;k--;)
     547        for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
    501548          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
    502549  else
    503550    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
    504551      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
    505         for(int k=ColumnCounter;k--;)
     552        for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
    506553          Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
    507554  return true;
    … …  
    524571    // count maximum of columns
    525572    RowCounter[MatrixCounter] = 0;
    526     for(int j=0; j < MatrixCounter;j++) // (energy matrix might be bigger than number of atoms in terms of rows)
     573    ColumnCounter[MatrixCounter] = 0;
     574    for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
    527575      if (RowCounter[j] > RowCounter[MatrixCounter])
    528576        RowCounter[MatrixCounter] = RowCounter[j];
     577      if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
     578        ColumnCounter[MatrixCounter] = ColumnCounter[j];
     579    }
    529580    // allocate last plus one matrix
    530     cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
    531     Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     581    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
     582    Matrix[MatrixCounter] = Malloc<double*>(RowCounter[MatrixCounter] + 1, "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    532583    for(int j=0;j<=RowCounter[MatrixCounter];j++)
    533       Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    534 
     584      Matrix[MatrixCounter][j] = Malloc<double>(ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     585   
    535586    // try independently to parse global energysuffix file if present
    536587    strncpy(filename, name, 1023);
    … …  
    556607
    557608  cout << "Parsing force indices for " << MatrixCounter << " matrices." << endl;
    558   Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "ForceMatrix::ParseIndices: **Indices");
     609  Indices = Malloc<int*>(MatrixCounter + 1, "ForceMatrix::ParseIndices: **Indices");
    559610  line << name << FRAGMENTPREFIX << FORCESFILE;
    560611  input.open(line.str().c_str(), ios::in);
    … …  
    569620    line.str(filename);
    570621    // parse the values
    571     Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
     622    Indices[i] = Malloc<int>(RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
    572623    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    573624    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
    574     Free((void **)&FragmentNumber, "ForceMatrix::ParseIndices: *FragmentNumber");
     625    Free(&FragmentNumber);
    575626    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
    576627      line >> Indices[i][j];
    … …  
    579630    //cout << endl;
    580631  }
    581   Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
     632  Indices[MatrixCounter] = Malloc<int>(RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
    582633  for(int j=RowCounter[MatrixCounter];j--;) {
    583634    Indices[MatrixCounter][j] = j;
    … …  
    589640
    590641/** Sums the forces and puts into last element of \a ForceMatrix::Matrix.
    591  * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
     642 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
    592643 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    593644 * \param Order bond order
    … …  
    609660      if (j != -1) {
    610661        //if (j == 0) cout << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
    611         for(int k=2;k<ColumnCounter;k++)
     662        for(int k=2;k<ColumnCounter[MatrixCounter];k++)
    612663          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
    613664      }
    … …  
    656707    input.close();
    657708
     709    ColumnCounter[MatrixCounter] = 0;
     710    for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
     711      if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
     712        ColumnCounter[MatrixCounter] = ColumnCounter[j];
     713    }
     714 
    658715    // allocate last plus one matrix
    659     cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter << " columns." << endl;
    660     Matrix[MatrixCounter] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixCounter]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     716    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
     717    Matrix[MatrixCounter] = Malloc<double*>(RowCounter[MatrixCounter] + 1, "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    661718    for(int j=0;j<=RowCounter[MatrixCounter];j++)
    662       Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     719      Matrix[MatrixCounter][j] = Malloc<double>(ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     720
     721    // try independently to parse global forcesuffix file if present
     722    strncpy(filename, name, 1023);
     723    strncat(filename, prefix, 1023-strlen(filename));
     724    strncat(filename, suffix.c_str(), 1023-strlen(filename));
     725    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
     726  }
     727 
     728
     729  return status;
     730};
     731
     732// ======================================= CLASS HessianMatrix =============================
     733
     734/** Parsing force Indices of each fragment
     735 * \param *name directory with \a ForcesFile
     736 * \return parsing successful
     737 */
     738bool HessianMatrix::ParseIndices(char *name)
     739{
     740  ifstream input;
     741  char *FragmentNumber = NULL;
     742  char filename[1023];
     743  stringstream line;
     744 
     745  cout << "Parsing hessian indices for " << MatrixCounter << " matrices." << endl;
     746  Indices = Malloc<int*>(MatrixCounter + 1, "HessianMatrix::ParseIndices: **Indices");
     747  line << name << FRAGMENTPREFIX << FORCESFILE;
     748  input.open(line.str().c_str(), ios::in);
     749  //cout << "Opening " << line.str() << " ... "  << input << endl;
     750  if (input == NULL) {
     751    cout << endl << "Unable to open " << line.str() << ", is the directory correct?" << endl;
     752    return false;
     753  }
     754  for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
     755    // get the number of atoms for this fragment
     756    input.getline(filename, 1023);
     757    line.str(filename);
     758    // parse the values
     759    Indices[i] = Malloc<int>(RowCounter[i], "HessianMatrix::ParseIndices: *Indices[]");
     760    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
     761    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
     762    Free(&FragmentNumber);
     763    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
     764      line >> Indices[i][j];
     765      //cout << " " << Indices[i][j];
     766    }
     767    //cout << endl;
     768  }
     769  Indices[MatrixCounter] = Malloc<int>(RowCounter[MatrixCounter], "HessianMatrix::ParseIndices: *Indices[]");
     770  for(int j=RowCounter[MatrixCounter];j--;) {
     771    Indices[MatrixCounter][j] = j;
     772  }
     773  input.close();
     774  return true;
     775};
     776
     777
     778/** Sums the hessian entries and puts into last element of \a HessianMatrix::Matrix.
     779 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     780 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
     781 * \param Order bond order
     782 *  \param sign +1 or -1
     783 * \return true if summing was successful
     784 */
     785bool HessianMatrix::SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
     786{
     787  int FragmentNr;
     788  // sum forces
     789  for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
     790    FragmentNr = KeySet.OrderSet[Order][i];
     791    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
     792      int j = Indices[ FragmentNr ][l];
     793      if (j > RowCounter[MatrixCounter]) {
     794        cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
     795        return false;
     796      }
     797      if (j != -1) {
     798        for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
     799          int k = Indices[ FragmentNr ][m];
     800          if (k > ColumnCounter[MatrixCounter]) {
     801            cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
     802            return false;
     803          }
     804          if (k != -1) {
     805            //cout << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
     806            Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][m];
     807          }
     808        }
     809      }
     810    }
     811  }
     812  return true;
     813};
     814
     815/** Constructor for class HessianMatrix.
     816 */
     817HessianMatrix::HessianMatrix() : MatrixContainer()
     818{
     819   IsSymmetric = true;
     820}
     821
     822/** Sums the hessian entries with each factor and put into last element of \a ***Matrix.
     823 * Sums over "E"-terms to create the "F"-terms
     824 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     825 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
     826 * \param Order bond order
     827 * \return true if summing was successful
     828 */
     829bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order)
     830{
     831  // go through each order
     832  for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
     833    //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     834    // then go per order through each suborder and pick together all the terms that contain this fragment
     835    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
     836      for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
     837        if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
     838          //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
     839          // if the fragment's indices are all in the current fragment
     840          for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
     841            int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
     842            //cout << "Current row index is " << k << "/" << m << "." << endl;
     843            if (m != -1) { // if it's not an added hydrogen
     844              for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
     845                //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
     846                if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
     847                  m = l;
     848                  break; 
     849                }
     850              }
     851              //cout << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
     852              if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
     853                cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
     854                return false;
     855              }
     856             
     857              for(int l=0;l<ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l++) {
     858                int n = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][l];
     859                //cout << "Current column index is " << l << "/" << n << "." << endl;
     860                if (n != -1) { // if it's not an added hydrogen
     861                  for (int p=0;p<ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
     862                    //cout << "Comparing " << n << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
     863                    if (n == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p]) {
     864                      n = p;
     865                      break; 
     866                    }
     867                  }
     868                  //cout << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
     869                  if (n > ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
     870                    cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
     871                    return false;
     872                  }
     873                  if (Order == SubOrder) { // equal order is always copy from Energies
     874                    //cout << "Adding " << MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
     875                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     876                  } else {
     877                    //cout << "Subtracting " << Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
     878                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     879                  }
     880                }
     881              }
     882            }
     883            //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
     884             //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
     885          }
     886        } else {
     887          //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     888        }
     889      }
     890    }
     891   //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
     892  }
     893 
     894  return true;
     895};
     896
     897/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
     898 * \param *name directory with files
     899 * \param *prefix prefix of each matrix file
     900 * \param *suffix suffix of each matrix file
     901 * \param skiplines number of inital lines to skip
     902 * \param skiplines number of inital columns to skip
     903 * \return parsing successful
     904 */
     905bool HessianMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
     906{
     907  char filename[1023];
     908  ifstream input;
     909  stringstream file;
     910  int nr;
     911  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);
     912
     913  if (status) {
     914    // count number of atoms for last plus one matrix
     915    file << name << FRAGMENTPREFIX << KEYSETFILE;
     916    input.open(file.str().c_str(), ios::in);
     917    if (input == NULL) {
     918      cout << endl << "Unable to open " << file.str() << ", is the directory correct?" << endl;
     919      return false;
     920    }
     921    RowCounter[MatrixCounter] = 0;
     922    ColumnCounter[MatrixCounter] = 0;
     923    while (!input.eof()) {
     924      input.getline(filename, 1023);
     925      stringstream zeile(filename);
     926      while (!zeile.eof()) {
     927        zeile >> nr;
     928        //cout << "Current index: " << nr << "." << endl;
     929        if (nr > RowCounter[MatrixCounter]) {
     930          RowCounter[MatrixCounter] = nr;
     931          ColumnCounter[MatrixCounter] = nr;
     932        }
     933      }
     934    }
     935    RowCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
     936    ColumnCounter[MatrixCounter]++;    // nr start at 0, count starts at 1
     937    input.close();
     938 
     939    // allocate last plus one matrix
     940    cout << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl;
     941    Matrix[MatrixCounter] = Malloc<double*>(RowCounter[MatrixCounter] + 1, "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     942    for(int j=0;j<=RowCounter[MatrixCounter];j++)
     943      Matrix[MatrixCounter][j] = Malloc<double>(ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    663944
    664945    // try independently to parse global forcesuffix file if present
    … …  
    690971KeySetsContainer::~KeySetsContainer() {
    691972  for(int i=FragmentCounter;i--;)
    692     Free((void **)&KeySets[i], "KeySetsContainer::~KeySetsContainer: *KeySets[]");
     973    Free(&KeySets[i]);
    693974  for(int i=Order;i--;)
    694     Free((void **)&OrderSet[i], "KeySetsContainer::~KeySetsContainer: *OrderSet[]");
    695   Free((void **)&KeySets, "KeySetsContainer::~KeySetsContainer: **KeySets");
    696   Free((void **)&OrderSet, "KeySetsContainer::~KeySetsContainer: **OrderSet");
    697   Free((void **)&AtomCounter, "KeySetsContainer::~KeySetsContainer: *AtomCounter");
    698   Free((void **)&FragmentsPerOrder, "KeySetsContainer::~KeySetsContainer: *FragmentsPerOrder");
     975    Free(&OrderSet[i]);
     976  Free(&KeySets);
     977  Free(&OrderSet);
     978  Free(&AtomCounter);
     979  Free(&FragmentsPerOrder);
    699980};
    700981
    … …  
    713994  FragmentCounter = FCounter;
    714995  cout << "Parsing key sets." << endl;
    715   KeySets = (int **) Malloc(sizeof(int *)*FragmentCounter, "KeySetsContainer::ParseKeySets: **KeySets");
     996  KeySets = Malloc<int*>(FragmentCounter, "KeySetsContainer::ParseKeySets: **KeySets");
    716997  for(int i=FragmentCounter;i--;)
    717998    KeySets[i] = NULL;
    … …  
    7231004  }
    7241005
    725   AtomCounter = (int *) Malloc(sizeof(int)*FragmentCounter, "KeySetsContainer::ParseKeySets: *RowCounter");
     1006  AtomCounter = Malloc<int>(FragmentCounter, "KeySetsContainer::ParseKeySets: *RowCounter");
    7261007  for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
    7271008    stringstream line;
    7281009    AtomCounter[i] = ACounter[i];
    7291010    // parse the values
    730     KeySets[i] = (int *) Malloc(sizeof(int)*AtomCounter[i], "KeySetsContainer::ParseKeySets: *KeySets[]");
     1011    KeySets[i] = Malloc<int>(AtomCounter[i], "KeySetsContainer::ParseKeySets: *KeySets[]");
    7311012    for(int j=AtomCounter[i];j--;)
    7321013      KeySets[i][j] = -1;
    7331014    FragmentNumber = FixedDigitNumber(FragmentCounter, i);
    7341015    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
    735     Free((void **)&FragmentNumber, "KeySetsContainer::ParseKeySets: *FragmentNumber");
     1016    Free(&FragmentNumber);
    7361017    input.getline(filename, 1023);
    7371018    line.str(filename);
    … …  
    7671048
    7681049  // scan through all to determine fragments per order
    769   FragmentsPerOrder = (int *) Malloc(sizeof(int)*Order, "KeySetsContainer::ParseManyBodyTerms: *FragmentsPerOrder");
     1050  FragmentsPerOrder = Malloc<int>(Order, "KeySetsContainer::ParseManyBodyTerms: *FragmentsPerOrder");
    7701051  for(int i=Order;i--;)
    7711052    FragmentsPerOrder[i] = 0;
    … …  
    7811062
    7821063  // scan through all to gather indices to each order set
    783   OrderSet = (int **) Malloc(sizeof(int *)*Order, "KeySetsContainer::ParseManyBodyTerms: **OrderSet");
     1064  OrderSet = Malloc<int*>(Order, "KeySetsContainer::ParseManyBodyTerms: **OrderSet");
    7841065  for(int i=Order;i--;)
    785     OrderSet[i] = (int *) Malloc(sizeof(int)*FragmentsPerOrder[i], "KeySetsContainer::ParseManyBodyTermsKeySetsContainer::ParseManyBodyTerms: *OrderSet[]");
     1066    OrderSet[i] = Malloc<int>(FragmentsPerOrder[i], "KeySetsContainer::ParseManyBodyTermsKeySetsContainer::ParseManyBodyTerms: *OrderSet[]");
    7861067  for(int i=Order;i--;)
    7871068    FragmentsPerOrder[i] = 0;
Note: See TracChangeset for help on using the changeset viewer.