src/gromacs/topology/block.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
   7  * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "block.h"
  41
  42 #include <cstdio>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/utility/listoflists.h"
  47 #include "gromacs/utility/smalloc.h"
  48 #include "gromacs/utility/txtdump.h"
  49
  50 void gmx::RangePartitioning::setAllBlocksSizeOne(int numBlocksToSet)
  51 {
  52     if (!allBlocksHaveSizeOne())
  53     {
  54         clear();
  55     }
  56     if (numBlocksToSet < numBlocks())
  57     {
  58         index_.resize(numBlocksToSet + 1);
  59     }
  60     else if (numBlocksToSet > numBlocks())
  61     {
  62         for (int b = numBlocks(); b < numBlocksToSet; b++)
  63         {
  64             appendBlock(1);
  65         }
  66     }
  67 }
  68
  69 void init_block(t_block* block)
  70 {
  71     block->nr           = 0;
  72     block->nalloc_index = 1;
  73     snew(block->index, block->nalloc_index);
  74     block->index[0] = 0;
  75 }
  76
  77 void init_block_null(t_block* block)
  78 {
  79     block->nr           = 0;
  80     block->nalloc_index = 0;
  81     block->index        = nullptr;
  82 }
  83
  84 void init_blocka(t_blocka* block)
  85 {
  86     block->nr           = 0;
  87     block->nra          = 0;
  88     block->nalloc_index = 1;
  89     snew(block->index, block->nalloc_index);
  90     block->index[0] = 0;
  91     block->nalloc_a = 0;
  92     block->a        = nullptr;
  93 }
  94
  95 void init_blocka_null(t_blocka* block)
  96 {
  97     block->nr           = 0;
  98     block->nra          = 0;
  99     block->nalloc_index = 0;
 100     block->index        = nullptr;
 101     block->nalloc_a     = 0;
 102     block->a            = nullptr;
 103 }
 104
 105 t_blocka* new_blocka()
 106 {
 107     t_blocka* block = nullptr;
 108
 109     snew(block, 1);
 110     snew(block->index, 1);
 111
 112     return block;
 113 }
 114
 115 void done_block(t_block* block)
 116 {
 117     block->nr = 0;
 118     sfree(block->index);
 119     block->index        = nullptr;
 120     block->nalloc_index = 0;
 121 }
 122
 123 void done_blocka(t_blocka* block)
 124 {
 125     block->nr  = 0;
 126     block->nra = 0;
 127     sfree(block->index);
 128     sfree(block->a);
 129     block->index        = nullptr;
 130     block->a            = nullptr;
 131     block->nalloc_index = 0;
 132     block->nalloc_a     = 0;
 133 }
 134
 135 void stupid_fill_block(t_block* grp, int natom, gmx_bool bOneIndexGroup)
 136 {
 137     if (bOneIndexGroup)
 138     {
 139         grp->nalloc_index = 2;
 140         srenew(grp->index, grp->nalloc_index);
 141         grp->index[0] = 0;
 142         grp->index[1] = natom;
 143         grp->nr       = 1;
 144     }
 145     else
 146     {
 147         grp->nalloc_index = natom + 1;
 148         srenew(grp->index, grp->nalloc_index);
 149         for (int i = 0; i <= natom; ++i)
 150         {
 151             grp->index[i] = i;
 152         }
 153         grp->nr = natom;
 154     }
 155 }
 156
 157 void stupid_fill_blocka(t_blocka* grp, int natom)
 158 {
 159     grp->nalloc_a = natom;
 160     snew(grp->a, grp->nalloc_a);
 161     for (int i = 0; i < natom; ++i)
 162     {
 163         grp->a[i] = i;
 164     }
 165     grp->nra = natom;
 166
 167     grp->nalloc_index = natom + 1;
 168     snew(grp->index, grp->nalloc_index);
 169     for (int i = 0; i <= natom; ++i)
 170     {
 171         grp->index[i] = i;
 172     }
 173     grp->nr = natom;
 174 }
 175
 176 void copy_blocka(const t_blocka* src, t_blocka* dest)
 177 {
 178     dest->nr = src->nr;
 179     /* Workaround for inconsistent handling of nalloc_index in
 180      * other parts of the code. Often nalloc_index and nalloc_a
 181      * are not set.
 182      */
 183     dest->nalloc_index = std::max(src->nalloc_index, dest->nr + 1);
 184     snew(dest->index, dest->nalloc_index);
 185     for (int i = 0; i < dest->nr + 1; ++i)
 186     {
 187         dest->index[i] = src->index[i];
 188     }
 189     dest->nra = src->nra;
 190     /* See above. */
 191     dest->nalloc_a = std::max(src->nalloc_a, dest->nra);
 192     snew(dest->a, dest->nalloc_a);
 193     for (int i = 0; i < dest->nra; ++i)
 194     {
 195         dest->a[i] = src->a[i];
 196     }
 197 }
 198
 199 static int pr_block_title(FILE* fp, int indent, const char* title, const t_block* block)
 200 {
 201     if (available(fp, block, indent, title))
 202     {
 203         indent = pr_title(fp, indent, title);
 204         pr_indent(fp, indent);
 205         fprintf(fp, "nr=%d\n", block->nr);
 206     }
 207     return indent;
 208 }
 209
 210 static int pr_blocka_title(FILE* fp, int indent, const char* title, const t_blocka* block)
 211 {
 212     if (available(fp, block, indent, title))
 213     {
 214         indent = pr_title(fp, indent, title);
 215         pr_indent(fp, indent);
 216         fprintf(fp, "nr=%d\n", block->nr);
 217         pr_indent(fp, indent);
 218         fprintf(fp, "nra=%d\n", block->nra);
 219     }
 220     return indent;
 221 }
 222
 223 static int pr_listoflists_title(FILE* fp, int indent, const char* title, const gmx::ListOfLists<int>* lists)
 224 {
 225     if (available(fp, lists, indent, title))
 226     {
 227         indent = pr_title(fp, indent, title);
 228         pr_indent(fp, indent);
 229         fprintf(fp, "numLists=%zu\n", lists->size());
 230         pr_indent(fp, indent);
 231         fprintf(fp, "numElements=%d\n", lists->numElements());
 232     }
 233     return indent;
 234 }
 235
 236 static void low_pr_blocka(FILE* fp, int indent, const char* title, const t_blocka* block, gmx_bool bShowNumbers)
 237 {
 238     if (available(fp, block, indent, title))
 239     {
 240         indent = pr_blocka_title(fp, indent, title, block);
 241         for (int i = 0; i <= block->nr; i++)
 242         {
 243             pr_indent(fp, indent + INDENT);
 244             fprintf(fp, "%s->index[%d]=%d\n", title, bShowNumbers ? i : -1, block->index[i]);
 245         }
 246         for (int i = 0; i < block->nra; i++)
 247         {
 248             pr_indent(fp, indent + INDENT);
 249             fprintf(fp, "%s->a[%d]=%d\n", title, bShowNumbers ? i : -1, block->a[i]);
 250         }
 251     }
 252 }
 253
 254 void pr_block(FILE* fp, int indent, const char* title, const t_block* block, gmx_bool bShowNumbers)
 255 {
 256     if (available(fp, block, indent, title))
 257     {
 258         indent    = pr_block_title(fp, indent, title, block);
 259         int start = 0;
 260         if (block->index[start] != 0)
 261         {
 262             fprintf(fp, "block->index[%d] should be 0\n", start);
 263         }
 264         else
 265         {
 266             for (int i = 0; i < block->nr; i++)
 267             {
 268                 int end = block->index[i + 1];
 269                 pr_indent(fp, indent);
 270                 if (end <= start)
 271                 {
 272                     fprintf(fp, "%s[%d]={}\n", title, i);
 273                 }
 274                 else
 275                 {
 276                     fprintf(fp,
 277                             "%s[%d]={%d..%d}\n",
 278                             title,
 279                             bShowNumbers ? i : -1,
 280                             bShowNumbers ? start : -1,
 281                             bShowNumbers ? end - 1 : -1);
 282                 }
 283                 start = end;
 284             }
 285         }
 286     }
 287 }
 288
 289 void pr_blocka(FILE* fp, int indent, const char* title, const t_blocka* block, gmx_bool bShowNumbers)
 290 {
 291     bool ok = false;
 292
 293     if (available(fp, block, indent, title))
 294     {
 295         indent    = pr_blocka_title(fp, indent, title, block);
 296         int start = 0;
 297         int end   = start;
 298         ok        = (block->index[start] == 0);
 299         if (!ok)
 300         {
 301             fprintf(fp, "block->index[%d] should be 0\n", start);
 302         }
 303         else
 304         {
 305             for (int i = 0; i < block->nr; i++)
 306             {
 307                 end      = block->index[i + 1];
 308                 int size = pr_indent(fp, indent);
 309                 if (end <= start)
 310                 {
 311                     size += fprintf(fp, "%s[%d]={", title, i);
 312                 }
 313                 else
 314                 {
 315                     size += fprintf(fp,
 316                                     "%s[%d][%d..%d]={",
 317                                     title,
 318                                     bShowNumbers ? i : -1,
 319                                     bShowNumbers ? start : -1,
 320                                     bShowNumbers ? end - 1 : -1);
 321                 }
 322                 for (int j = start; j < end; j++)
 323                 {
 324                     if (j > start)
 325                     {
 326                         size += fprintf(fp, ", ");
 327                     }
 328                     if ((size) > (USE_WIDTH))
 329                     {
 330                         fprintf(fp, "\n");
 331                         size = pr_indent(fp, indent + INDENT);
 332                     }
 333                     size += fprintf(fp, "%d", block->a[j]);
 334                 }
 335                 fprintf(fp, "}\n");
 336                 start = end;
 337             }
 338         }
 339         if ((end != block->nra) || (!ok))
 340         {
 341             pr_indent(fp, indent);
 342             fprintf(fp, "tables inconsistent, dumping complete tables:\n");
 343             low_pr_blocka(fp, indent, title, block, bShowNumbers);
 344         }
 345     }
 346 }
 347
 348 void pr_listoflists(FILE* fp, int indent, const char* title, const gmx::ListOfLists<int>* lists, gmx_bool bShowNumbers)
 349 {
 350     if (available(fp, lists, indent, title))
 351     {
 352         indent = pr_listoflists_title(fp, indent, title, lists);
 353         for (gmx::index i = 0; i < lists->ssize(); i++)
 354         {
 355             int                      size = pr_indent(fp, indent);
 356             gmx::ArrayRef<const int> list = (*lists)[i];
 357             if (list.empty())
 358             {
 359                 size += fprintf(fp, "%s[%d]={", title, int(i));
 360             }
 361             else
 362             {
 363                 size += fprintf(fp, "%s[%d][num=%zu]={", title, bShowNumbers ? int(i) : -1, list.size());
 364             }
 365             bool isFirst = true;
 366             for (const int j : list)
 367             {
 368                 if (!isFirst)
 369                 {
 370                     size += fprintf(fp, ", ");
 371                 }
 372                 if ((size) > (USE_WIDTH))
 373                 {
 374                     fprintf(fp, "\n");
 375                     size = pr_indent(fp, indent + INDENT);
 376                 }
 377                 size += fprintf(fp, "%d", j);
 378                 isFirst = false;
 379             }
 380             fprintf(fp, "}\n");
 381         }
 382     }
 383 }
 384
 385 void copy_block(const t_block* src, t_block* dst)
 386 {
 387     dst->nr = src->nr;
 388     /* Workaround for inconsistent handling of nalloc_index in
 389      * other parts of the code. Often nalloc_index and nalloc_a
 390      * are not set.
 391      */
 392     dst->nalloc_index = std::max(src->nalloc_index, dst->nr + 1);
 393     snew(dst->index, dst->nalloc_index);
 394     for (int i = 0; i < dst->nr + 1; ++i)
 395     {
 396         dst->index[i] = src->index[i];
 397     }
 398 }