Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/parser.cpp

    r29812d r042f82  
    88
    99#include "helpers.hpp"
    10 #include "memoryallocator.hpp"
    1110#include "parser.hpp"
    1211
    … …  
    5756MatrixContainer::MatrixContainer() {
    5857  Indices = NULL;
    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;
     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");
    6461  Matrix[0] = NULL;
    6562  RowCounter[0] = -1;
    6663  MatrixCounter = 0;
    67   ColumnCounter[0] = -1;
     64  ColumnCounter = 0;
    6865};
    6966
    … …  
    7370  if (Matrix != NULL) {
    7471    for(int i=MatrixCounter;i--;) {
    75       if ((ColumnCounter != NULL) && (RowCounter != NULL)) {
     72      if (RowCounter != NULL) {
    7673          for(int j=RowCounter[i]+1;j--;)
    77             Free(&Matrix[i][j]);
    78         Free(&Matrix[i]);
     74            Free((void **)&Matrix[i][j], "MatrixContainer::~MatrixContainer: *Matrix[][]");
     75        Free((void **)&Matrix[i], "MatrixContainer::~MatrixContainer: **Matrix[]");
    7976      }
    8077    }
    81     if ((ColumnCounter != NULL) && (RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
     78    if ((RowCounter != NULL) && (Matrix[MatrixCounter] != NULL))
    8279      for(int j=RowCounter[MatrixCounter]+1;j--;)
    83         Free(&Matrix[MatrixCounter][j]);
     80        Free((void **)&Matrix[MatrixCounter][j], "MatrixContainer::~MatrixContainer: *Matrix[MatrixCounter][]");
    8481    if (MatrixCounter != 0)
    85       Free(&Matrix[MatrixCounter]);
    86     Free(&Matrix);
     82      Free((void **)&Matrix[MatrixCounter], "MatrixContainer::~MatrixContainer: **Matrix[MatrixCounter]");
     83    Free((void **)&Matrix, "MatrixContainer::~MatrixContainer: ***Matrix");
    8784  }
    8885  if (Indices != NULL)
    8986    for(int i=MatrixCounter+1;i--;) {
    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  */
    107 bool 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 };
     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
    13695
    13796/** Parsing a number of matrices.
    … …  
    162121  }
    163122
    164   // parse header
    165   Header[MatrixNr] = Malloc<char>(1024, "MatrixContainer::ParseMatrix: *Header[]");
     123  // skip some initial lines
    166124  for (int m=skiplines+1;m--;)
    167     input.getline(Header[MatrixNr], 1023);
    168  
     125    input.getline(Header, 1023);
     126
    169127  // scan header for number of columns
    170   line.str(Header[MatrixNr]);
     128  line.str(Header);
    171129  for(int k=skipcolumns;k--;)
    172     line >> Header[MatrixNr];
     130    line >> Header;
    173131  //cout << line.str() << endl;
    174   ColumnCounter[MatrixNr]=0;
     132  ColumnCounter=0;
    175133  while ( getline(line,token, '\t') ) {
    176134    if (token.length() > 0)
    177       ColumnCounter[MatrixNr]++;
     135      ColumnCounter++;
    178136  }
    179137  //cout << line.str() << endl;
    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  
     138  //cout << "ColumnCounter: " << ColumnCounter << "." << endl;
     139  if (ColumnCounter == 0)
     140    cerr << "ColumnCounter: " << ColumnCounter << " from file " << name << ", this is probably an error!" << endl;
     141
    184142  // scan rest for number of rows/lines
    185143  RowCounter[MatrixNr]=-1;    // counts one line too much
    … …  
    197155
    198156  // allocate matrix if it's not zero dimension in one direction
    199   if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
    200     Matrix[MatrixNr] = Malloc<double*>(RowCounter[MatrixNr] + 1, "MatrixContainer::ParseMatrix: **Matrix[]");
    201  
     157  if ((ColumnCounter > 0) && (RowCounter[MatrixNr] > -1)) {
     158    Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     159
    202160    // parse in each entry for this matrix
    203161    input.clear();
    204162    input.seekg(ios::beg);
    205163    for (int m=skiplines+1;m--;)
    206       input.getline(Header[MatrixNr], 1023);    // skip header
    207     line.str(Header[MatrixNr]);
     164      input.getline(Header, 1023);    // skip header
     165    line.str(Header);
    208166    for(int k=skipcolumns;k--;)  // skip columns in header too
    209167      line >> filename;
    210     strncpy(Header[MatrixNr], line.str().c_str(), 1023); 
     168    strncpy(Header, line.str().c_str(), 1023);
    211169    for(int j=0;j<RowCounter[MatrixNr];j++) {
    212       Matrix[MatrixNr][j] = Malloc<double>(ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[][]");
     170      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    213171      input.getline(filename, 1023);
    214172      stringstream lines(filename);
    … …  
    216174      for(int k=skipcolumns;k--;)
    217175        lines >> filename;
    218       for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
     176      for(int k=0;(k<ColumnCounter) && (!lines.eof());k++) {
    219177        lines >> Matrix[MatrixNr][j][k];
    220178        //cout << " " << setprecision(2) << Matrix[MatrixNr][j][k];;
    221179      }
    222180      //cout << endl;
    223       Matrix[MatrixNr][ RowCounter[MatrixNr] ] = Malloc<double>(ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[RowCounter[MatrixNr]][]");
    224       for(int j=ColumnCounter[MatrixNr];j--;)
     181      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
     182      for(int j=ColumnCounter;j--;)
    225183        Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
    226184    }
    227185  } else {
    228     cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl;
     186    cerr << "ERROR: Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl;
    229187  }
    230188  input.close();
    … …  
    275233
    276234  cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
    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");
     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");
    281237  for(int i=MatrixCounter+1;i--;) {
    282238    Matrix[i] = NULL;
    283     Header[i] = NULL;
    284239    RowCounter[i] = -1;
    285     ColumnCounter[i] = -1;
    286240  }
    287241  for(int i=0; i < MatrixCounter;i++) {
    … …  
    292246    if (!ParseMatrix(file.str().c_str(), skiplines, skipcolumns, i))
    293247      return false;
    294     Free(&FragmentNumber);
     248    Free((void **)&FragmentNumber, "MatrixContainer::ParseFragmentMatrix: *FragmentNumber");
    295249  }
    296250  return true;
    … …  
    298252
    299253/** Allocates and resets the memory for a number \a MCounter of matrices.
    300  * \param **GivenHeader Header line for each matrix
     254 * \param *GivenHeader Header line
    301255 * \param MCounter number of matrices
    302256 * \param *RCounter number of rows for each matrix
    303  * \param *CCounter number of columns for each matrix
    304  * \return Allocation successful
    305  */
    306 bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
    307 {
     257 * \param CCounter number of columns for all matrices
     258 * \return Allocation successful
     259 */
     260bool 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
    308265  MatrixCounter = MCounter;
    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");
     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");
    313269  for(int i=MatrixCounter+1;i--;) {
    314     Header[i] = Malloc<char>(1024, "MatrixContainer::AllocateMatrix: *Header[i]");
    315     strncpy(Header[i], GivenHeader[i], 1023);
    316270    RowCounter[i] = RCounter[i];
    317     ColumnCounter[i] = CCounter[i];
    318     Matrix[i] = Malloc<double*>(RowCounter[i] + 1, "MatrixContainer::AllocateMatrix: **Matrix[]");
     271    Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
    319272    for(int j=RowCounter[i]+1;j--;) {
    320       Matrix[i][j] = Malloc<double>(ColumnCounter[i], "MatrixContainer::AllocateMatrix: *Matrix[][]");
    321       for(int k=ColumnCounter[i];k--;)
     273      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     274      for(int k=ColumnCounter;k--;)
    322275        Matrix[i][j][k] = 0.;
    323276    }
    … …  
    333286  for(int i=MatrixCounter+1;i--;)
    334287    for(int j=RowCounter[i]+1;j--;)
    335       for(int k=ColumnCounter[i];k--;)
     288      for(int k=ColumnCounter;k--;)
    336289        Matrix[i][j][k] = 0.;
    337290   return true;
    … …  
    346299  for(int i=MatrixCounter+1;i--;)
    347300    for(int j=RowCounter[i]+1;j--;)
    348       for(int k=ColumnCounter[i];k--;)
     301      for(int k=ColumnCounter;k--;)
    349302        if (fabs(Matrix[i][j][k]) > max)
    350303          max = fabs(Matrix[i][j][k]);
    … …  
    362315  for(int i=MatrixCounter+1;i--;)
    363316    for(int j=RowCounter[i]+1;j--;)
    364       for(int k=ColumnCounter[i];k--;)
     317      for(int k=ColumnCounter;k--;)
    365318        if (fabs(Matrix[i][j][k]) < min)
    366319          min = fabs(Matrix[i][j][k]);
    … …  
    378331{
    379332  for(int j=RowCounter[MatrixCounter]+1;j--;)
    380     for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
     333    for(int k=skipcolumns;k<ColumnCounter;k++)
    381334      Matrix[MatrixCounter][j][k] = value;
    382335   return true;
    … …  
    391344{
    392345  for(int j=RowCounter[MatrixCounter]+1;j--;)
    393     for(int k=skipcolumns;k<ColumnCounter[MatrixCounter];k++)
     346    for(int k=skipcolumns;k<ColumnCounter;k++)
    394347      Matrix[MatrixCounter][j][k] = values[j][k];
    395348   return true;
    396349};
    397350
    398 /** Sums the entries with each factor and put into last element of \a ***Matrix.
     351/** Sums the energy with each factor and put into last element of \a ***Energies.
    399352 * Sums over "E"-terms to create the "F"-terms
    400  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     353 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
    401354 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    402355 * \param Order bond order
    … …  
    431384              }
    432385              if (Order == SubOrder) { // equal order is always copy from Energies
    433                 for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;) // then adds/subtract each column
     386                for(int l=ColumnCounter;l--;) // then adds/subtract each column
    434387                  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
    435388              } else {
    436                 for(int l=ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l--;)
     389                for(int l=ColumnCounter;l--;)
    437390                  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][l] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
    438391              }
    439392            }
    440             //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
     393            //if ((ColumnCounter>1) && (RowCounter[0]-1 >= 1))
    441394             //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
    442395          }
    … …  
    467420    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    468421    line << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix;
    469     Free(&FragmentNumber);
     422    Free((void **)&FragmentNumber, "*FragmentNumber");
    470423    output.open(line.str().c_str(), ios::out);
    471424    if (output == NULL) {
    … …  
    473426      return false;
    474427    }
    475     output << Header[i] << endl;
     428    output << Header << endl;
    476429    for(int j=0;j<RowCounter[i];j++) {
    477       for(int k=0;k<ColumnCounter[i];k++)
     430      for(int k=0;k<ColumnCounter;k++)
    478431        output << scientific << Matrix[i][j][k] << "\t";
    479432      output << endl;
    … …  
    502455    return false;
    503456  }
    504   output << Header[MatrixCounter] << endl;
     457  output << Header << endl;
    505458  for(int j=0;j<RowCounter[MatrixCounter];j++) {
    506     for(int k=0;k<ColumnCounter[MatrixCounter];k++)
     459    for(int k=0;k<ColumnCounter;k++)
    507460      output << scientific << Matrix[MatrixCounter][j][k] << "\t";
    508461    output << endl;
    … …  
    521474{
    522475  cout << "Parsing energy indices." << endl;
    523   Indices = Malloc<int*>(MatrixCounter + 1, "EnergyMatrix::ParseIndices: **Indices");
     476  Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "EnergyMatrix::ParseIndices: **Indices");
    524477  for(int i=MatrixCounter+1;i--;) {
    525     Indices[i] = Malloc<int>(RowCounter[i], "EnergyMatrix::ParseIndices: *Indices[]");
     478    Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "EnergyMatrix::ParseIndices: *Indices[]");
    526479    for(int j=RowCounter[i];j--;)
    527480      Indices[i][j] = j;
    … …  
    532485/** Sums the energy with each factor and put into last element of \a EnergyMatrix::Matrix.
    533486 * Sums over the "F"-terms in ANOVA decomposition.
    534  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     487 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
    535488 * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
    536489 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    … …  
    545498    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
    546499      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
    547         for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
     500        for(int k=ColumnCounter;k--;)
    548501          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
    549502  else
    550503    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
    551504      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
    552         for(int k=ColumnCounter[ KeySet.OrderSet[Order][i] ];k--;)
     505        for(int k=ColumnCounter;k--;)
    553506          Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
    554507  return true;
    … …  
    571524    // count maximum of columns
    572525    RowCounter[MatrixCounter] = 0;
    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)
     526    for(int j=0; j < MatrixCounter;j++) // (energy matrix might be bigger than number of atoms in terms of rows)
    575527      if (RowCounter[j] > RowCounter[MatrixCounter])
    576528        RowCounter[MatrixCounter] = RowCounter[j];
    577       if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
    578         ColumnCounter[MatrixCounter] = ColumnCounter[j];
    579     }
    580529    // allocate last plus one 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[]");
     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[]");
    583532    for(int j=0;j<=RowCounter[MatrixCounter];j++)
    584       Matrix[MatrixCounter][j] = Malloc<double>(ColumnCounter[MatrixCounter], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    585    
     533      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     534
    586535    // try independently to parse global energysuffix file if present
    587536    strncpy(filename, name, 1023);
    … …  
    607556
    608557  cout << "Parsing force indices for " << MatrixCounter << " matrices." << endl;
    609   Indices = Malloc<int*>(MatrixCounter + 1, "ForceMatrix::ParseIndices: **Indices");
     558  Indices = (int **) Malloc(sizeof(int *)*(MatrixCounter+1), "ForceMatrix::ParseIndices: **Indices");
    610559  line << name << FRAGMENTPREFIX << FORCESFILE;
    611560  input.open(line.str().c_str(), ios::in);
    … …  
    620569    line.str(filename);
    621570    // parse the values
    622     Indices[i] = Malloc<int>(RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
     571    Indices[i] = (int *) Malloc(sizeof(int)*RowCounter[i], "ForceMatrix::ParseIndices: *Indices[]");
    623572    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    624573    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
    625     Free(&FragmentNumber);
     574    Free((void **)&FragmentNumber, "ForceMatrix::ParseIndices: *FragmentNumber");
    626575    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
    627576      line >> Indices[i][j];
    … …  
    630579    //cout << endl;
    631580  }
    632   Indices[MatrixCounter] = Malloc<int>(RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
     581  Indices[MatrixCounter] = (int *) Malloc(sizeof(int)*RowCounter[MatrixCounter], "ForceMatrix::ParseIndices: *Indices[]");
    633582  for(int j=RowCounter[MatrixCounter];j--;) {
    634583    Indices[MatrixCounter][j] = j;
    … …  
    640589
    641590/** Sums the forces and puts into last element of \a ForceMatrix::Matrix.
    642  * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     591 * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
    643592 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
    644593 * \param Order bond order
    … …  
    660609      if (j != -1) {
    661610        //if (j == 0) cout << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
    662         for(int k=2;k<ColumnCounter[MatrixCounter];k++)
     611        for(int k=2;k<ColumnCounter;k++)
    663612          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
    664613      }
    … …  
    707656    input.close();
    708657
    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  
    715658    // allocate last plus one 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[]");
     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[]");
    718661    for(int j=0;j<=RowCounter[MatrixCounter];j++)
    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  */
    738 bool 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  */
    785 bool 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  */
    817 HessianMatrix::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  */
    829 bool 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  */
    905 bool 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[][]");
     662      Matrix[MatrixCounter][j] = (double *) Malloc(sizeof(double)*ColumnCounter, "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
    944663
    945664    // try independently to parse global forcesuffix file if present
    … …  
    971690KeySetsContainer::~KeySetsContainer() {
    972691  for(int i=FragmentCounter;i--;)
    973     Free(&KeySets[i]);
     692    Free((void **)&KeySets[i], "KeySetsContainer::~KeySetsContainer: *KeySets[]");
    974693  for(int i=Order;i--;)
    975     Free(&OrderSet[i]);
    976   Free(&KeySets);
    977   Free(&OrderSet);
    978   Free(&AtomCounter);
    979   Free(&FragmentsPerOrder);
     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");
    980699};
    981700
    … …  
    994713  FragmentCounter = FCounter;
    995714  cout << "Parsing key sets." << endl;
    996   KeySets = Malloc<int*>(FragmentCounter, "KeySetsContainer::ParseKeySets: **KeySets");
     715  KeySets = (int **) Malloc(sizeof(int *)*FragmentCounter, "KeySetsContainer::ParseKeySets: **KeySets");
    997716  for(int i=FragmentCounter;i--;)
    998717    KeySets[i] = NULL;
    … …  
    1004723  }
    1005724
    1006   AtomCounter = Malloc<int>(FragmentCounter, "KeySetsContainer::ParseKeySets: *RowCounter");
     725  AtomCounter = (int *) Malloc(sizeof(int)*FragmentCounter, "KeySetsContainer::ParseKeySets: *RowCounter");
    1007726  for(int i=0;(i<FragmentCounter) && (!input.eof());i++) {
    1008727    stringstream line;
    1009728    AtomCounter[i] = ACounter[i];
    1010729    // parse the values
    1011     KeySets[i] = Malloc<int>(AtomCounter[i], "KeySetsContainer::ParseKeySets: *KeySets[]");
     730    KeySets[i] = (int *) Malloc(sizeof(int)*AtomCounter[i], "KeySetsContainer::ParseKeySets: *KeySets[]");
    1012731    for(int j=AtomCounter[i];j--;)
    1013732      KeySets[i][j] = -1;
    1014733    FragmentNumber = FixedDigitNumber(FragmentCounter, i);
    1015734    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
    1016     Free(&FragmentNumber);
     735    Free((void **)&FragmentNumber, "KeySetsContainer::ParseKeySets: *FragmentNumber");
    1017736    input.getline(filename, 1023);
    1018737    line.str(filename);
    … …  
    1048767
    1049768  // scan through all to determine fragments per order
    1050   FragmentsPerOrder = Malloc<int>(Order, "KeySetsContainer::ParseManyBodyTerms: *FragmentsPerOrder");
     769  FragmentsPerOrder = (int *) Malloc(sizeof(int)*Order, "KeySetsContainer::ParseManyBodyTerms: *FragmentsPerOrder");
    1051770  for(int i=Order;i--;)
    1052771    FragmentsPerOrder[i] = 0;
    … …  
    1062781
    1063782  // scan through all to gather indices to each order set
    1064   OrderSet = Malloc<int*>(Order, "KeySetsContainer::ParseManyBodyTerms: **OrderSet");
     783  OrderSet = (int **) Malloc(sizeof(int *)*Order, "KeySetsContainer::ParseManyBodyTerms: **OrderSet");
    1065784  for(int i=Order;i--;)
    1066     OrderSet[i] = Malloc<int>(FragmentsPerOrder[i], "KeySetsContainer::ParseManyBodyTermsKeySetsContainer::ParseManyBodyTerms: *OrderSet[]");
     785    OrderSet[i] = (int *) Malloc(sizeof(int)*FragmentsPerOrder[i], "KeySetsContainer::ParseManyBodyTermsKeySetsContainer::ParseManyBodyTerms: *OrderSet[]");
    1067786  for(int i=Order;i--;)
    1068787    FragmentsPerOrder[i] = 0;
Note: See TracChangeset for help on using the changeset viewer.