Guidelines for Developers

Altering the Program

To develop the program further, the following rules should be obeyed. Any rules concerning the PRCS archive program only hold for the members of the Heidelberg group.

Check out the latest code version from the PRCS archive to your own directory, or merge your personal version with the latest one. Update the PRCS keywords by using prcs rekey for each project.
Make additions so as to change the existing code as little as possible. When changing the program, always keep in mind that others also need to use and understand your new code. This means for instance:
- Split larger problems into smaller ones and write an individual subroutine for each.
- Try to be as general as possible.
- Use self-explanatory variable and routine names.
- Comments in the code are wellcome.
Update the revision number of any sub-project or analyse program changed. See Version and Revision Numbers
Describe your changes in the corresponding documentation files that have been checked out to your local directory. It is good (and efficient!) practice to document changes "on the fly", i.e. parallel to the programmings. Major changes should be reported also in the Changes of the Code documentation.
Test the new code with the Automatic Program Test. Add any new test cases that are necessary for your new options. If changes had to be made that set the Automatic Program Test out of order (e.g. a new output format), immediately update the test files. Check your code on different platforms.
If everything works correctly, again merge your code and documentation with the latest version, as it might have changed during your work. Then check in code and documentation to the PRCS archive. Enter an expressive version log.
Members of the Heidelberg group are to inform their colleagues by email about their alterations. Others are invited to send their changed files to the Heidelberg group for being included in future versions of the MCTDH package.

Debugging the Program

In order to debug the code set the option -d of the compile shell script when compiling the program. By this, the compiler produces additional symbol table information for full symbolic debugging instead of optimising the code. Note that a "d" is added to the name of the executable, e.g. mctdh82d. A debugger (e.g dbx or gdb) can then aid you in detecting a bug. It may also be helpful to compile and run your code on another platform, since different compilers often perform different checks of the code and produce differently rigorous compiler and run-time warning messages.

Known bugs that have not been fixed yet are compiled in the To Be Done documentation.

Memory Allocation

One of the more involving issues of Fortran programming (and MCTDH is mostly written in this language) is memory allocation. In Fortran 77, which is still used in high-performance computing because it produces code of outstanding speed, memory can only be allocated statically.
This is largely responsible for Fortran 77's speed advantage - the CPU does not need to check whether the address of an array has changed since the last access, as pointers can not be accessed by the programmer. Such arrays can then be held in the CPU cache without incurring any cache misses or pipeline stalls due to speculative execution on the wrong data. This advantage is intrinsic to the Fortran 77 language and hence even future C compilers with better optimisations will not reach the speed of a Fortran 77 program compiled with the same optimisations. This means that Fortran 90 with its possibility of dynamic memory allocation will also be slower than comparable Fortran 77 implementations.
The problem for complex software such as MCTDH is, of course, that the program does not know in advance how big the job the users sends its way is going to be.

"Quasistatic" Memory Allocation in MCTDH

This problem is circumvented in MCTDH while retaining the speed advantage afforded by Fortran 77: the memory needed is allocated at the start of the program in a small C routine which then calls the main part of the program, passing the pointers to the buffer arrays to the Fortran 77 sections of the code. Hence the pointers to the buffer arrays are set only once and are never changed later. The box below describes this process in more detail and explains how you can use it to add new modules to MCTDH.

Besides the buffer arrays for the wave function, which we do not want to touch, there are four general purpose buffers for complex, real, integer and logical data, called mc, mr, mi, and ml, respectively.

Each of these consists of smaller blocks of data, one each for different components of the program. If you would like to add a new code section with its own data block that is accessed throughout the running time of MCTDH (see below), you need to insert an additional block. The following description is based on the MCTDH application only but allocation, where used, proceeds in a similar fashion in the other programs contained within the MCTDH program package.

You can see examples of how this is done in the subroutine memprop in source/propwf/memprop.F. At the end of the appropriate code section, add the appropriate code to increase the size of the array type needed. If you need an array of complex numbers, for example, your code would look something like this:
```
      my_c_index = mcdim
      c_nums_needed = find_out(how_much_I_need)
      mcdim = my_c_index + c_nums_needed + 1
```
The variable my_c_index will help your new code section to find your data block and has to be inserted into the corresponding common block. The 1 added at the end is important; you can use the resulting additional field to check whether your code has written over the end of its data block (an access violation).
MCTDH then allocates mcdim memory cells for complex numbers in CALLPROP in the file source/propwf/fp_alloc.c and passes them on to runprop in source/propwf/runprop.F. Have a look at these files now to see what happens.

runprop then calls its subroutines with pointers to the arrays needed. This often looks something like this:

      call rdrst(time,psi,linwf,chkdvr,chkgrd,chkpsi,chkprp,
     +           mc(mcintpsi),lanczord,mr(mrworkr),iflag,mi(mijindx))

As seen in the code fragment above, the data buffers are often passed using a particular index for the data block needed in the subroutine. This desingn decision was made early on in the development of MCTDH. It has the following advantages and disadvantages:

The subroutines do not have to add the pointer to their data block at each access, i.e. rather than accessing the array's elements as
mc(my_c_index), mc(my_c_index+1), mc(my_c_index+2), ...
The elements are indexed as in mc(1), mc(2), mc(3), .... This is convenient for the programmer and slightly increases the performance of the code.
One disadvantage is felt when your new subroutine is called from several different subroutines, each of which got passed a different pointer to mc (mr, mi, ml). If there are only two or three of these, your subroutine may receive an extra argument indicating where it was called from and compensate accordingly. However, for functions called in the application of the Hamiltonian, for example, this is not practicable, as these are called from all over the place.
Another disadvantage, and this is the most important one, is waste of memory. A lot of the time, memory is only needed as buffer for certain operations in a very localised section of code (e.g. a subroutine). Globally allocating this buffer memory for every such subroutine results in large sections of memory that are unused most of the time.
It is therefore desirable to take these working buffers from the same section of memory. This can be done using the techniques given below.

Work Arrays for Temporary Data Storage in MCTDH

The subroutine memprop.F already contains provisions for temporary storage in the form of subarrays of the mc, mr, mi, and ml arrays.

These sections of the corresponding array are pointed to by the mcworkc, mrworkr, miworki, and mlworkl indices. The size of these sections of is given by the variables workcdim, workrdim, workidim, workldim.

If you want to use one of these work arrays, you will have to make sure that they fit the amount of data you want to write into them. This can be done by consulting the corresponding work?dim or setting its size prior to memory allocation.

The latter is demonstrated for example in subroutine wkdiagham in source/propwf/diagham.F - the variable workrdim is set to the maximum of its previous value and the value needed by this module using the construct
workrdim = max(workrdim,lanczorder**2+4*lanczorder-1)

You have to write a corresponding subroutine for your code and make sure that it is called by including a call statement into the appropriate subroutine (here propwf/zeigprop.F).

Note: The work arrays are highly volatile, i.e. if you should call any subroutine from a module that uses work memory, it is your responsibility to make sure that these do not overwrite the work memory.
The same is true if you change existing subroutines of functions to overwrite work memory - you have to make sure that these are not called from within sections of the code that rely on the work memory not to be overwritten.

Dynamic Memory Allocation in MCTDH

The problems indicated by the last points in the section "Quasistatic" Memory Allocation and immediately above are alleviated by using true dynamic memory allocation. For this purpose, MCTDH contains a small C library found in lib/utilities/dynaccess.c. This allows dynamic allocation and acces to dynamically allocated data from Fortran 77 sections of the code. Note, however, that this indirect access of data is very slow compared to native Fortran 77 data accesses. (Every time a dynamic data set is accessed, the pointer to the data, the index to the particular element and the address of the variable to write the array element to are put on the stack, a C function is called which has to dereference the pointer, add the index and transfer the array element to the variable provided.)

This means that dynamic data allocation should be used if one or several of the following points are satisfied:

For temporary storage. An example for this use can be found in the outgpop1 subroutine in source/propdo/d1outgpop.F. The help array of complex numbers is used to store the wavefunction psi while it is transformed and used to temporarily store data. At the end of the subroutine, the original contents of psi are read back from the help array which is consequently deallocated.
Generally for blockwise access. A (fixed length) Fortran 77 array (cache) may be used for blockwise access to the dynamically allocated data. This way, the number of indirect accesses via the dynaccess C library can be minimised, while keeping the advantage of dynamical access.
For rarely used operations or out of necessity. For subroutines which are not called in the inner loop of MCTDH, such as input and output operations which are typically called only once at the start or the end of processing, speed is not an issue and hence full advantage can be taken of the liberties afforded by dynamical data allocation.
On the other hand, necessity may force the programmer to use dynamical data in the inner loop because the above offset problem shift problem occurs and cannot be cured with justifiable programming effort. An example for such a use of the dynaccess library can be found in source/opfuncs/linint.F.

In the following, a brief example is given for using blockwise access to a dynamical storage array. Note that pointers are stored in arrays consisting of two integer numbers for portability. As such, they may be passed between Fortran 77 subroutines and functions.
dynaccess also contains routines to access array elements of type character, short int, integer, long int, real, double precision, single complex, double complex and string individually. When necessary, you may read a detailed description of the functions provided by dynaccess. source/opfuncs/linint.F, for example, uses some of these functions.

Block example:

Simple temporary storage example. In this example, some data in the MCTDH array mctdh_array is overwritten and written back to it from temporary storage temp after the subroutine finishes.

integer
            temp(2)

holds the pointer to the array.

      complex*16   mctdh_array(dim)
      integer      start,index,needed,temp(2),tempsize

c ___ allocate temp array of length tempsize
      tempsize = needed*complexsize()
      call alloc(temp,tempsize)
c ___ read contents of mctdh_array from element start into buffer temp
      index = (start - 1)*complexsize()
      call setblock(mctdh_array, temp, index, tempsize )
c ___ now mess around with mctdh_array a little
      mctdh_array(start) = (0.0, 1.0)
      mctdh_array(start+needed/2-1) = (0.5, 0.5)
      mctdh_array(start+needed-1) = (1.0, 0.0)
c ___ now we're done, get back original contents and deallocate memory
      call getblock(mctdh_array, temp, index, tempsize )
      call dealloc(temp, tempsize)

More involved example with cached access. In this example, data in a dynamical array is accessed sequentially by blockwise operations for speed. Individual elements in the data blocks are accessed via the auxiliary Fortran 77 array aux.

      integer auxlen
      parameter (auxlen = 16)
      complex*16   aux(auxlen)
      integer      index,needed,pointer(2),size
      integer      auxsize,ii,jj

c ___ allocate temp array of length size
      size = needed*complexsize()
      call alloc(pointer,size)
c ___ write something into pointer
      auxsize = auxlen*complexsize()
      index = 0
      do ii=1,needed,auxlen
        do jj=1,auxlen
          aux(jj) = cmplx(ii+jj, 1.0-(ii+jj))
        end do
        call setblock(aux, pointer, index, auxlen )
        index = index + auxlen
      end do
c ___ now pass around pointer array a little
      call my_routine1(pointer, size)
      call my_routine2(pointer, size)
c ___ now we're done, deallocate memory
      call dealloc(pointer, size)

In the subroutines, the data is accessed in just the same way with getblock.

Link: Detailed description of the functions provided by dynaccess.

Frequently Made Mistakes

Below you find a list of things that are frequently forgotten and likely to lead to problems. Again any items concerning the PRCS archive program are only relevant for the members of the Heidelberg group.

Follow the rules given in Altering the Program.
Analyse programs must have a program statement at the beginning (required by compile).
If one adds/removes include files to/from a routine, it has to be checked that the corresponding dependencies in the Makefile are still correct. Run $MCTDH_DIR/install/mkdepend to generate a new list of dependencies.
Code has to be merged (prcs merge) before it is checked in (prcs checkin).
When populating the archive (prcs populate), explicitly name the files to be added (to avoid unintentional population of files).
When renaming files, this must also be noted in the corresponding *.prj file (see the PRCS documentation for details).
The length of comment lines should not exceed 72 characters.
Real numbers in the code, input, and operator should be explicitly double precision, e.g. 1.0d0 rather than 1.0 (AIX machines otherwise interpret them as single precision).
Include files must be included also in the main program (to ensure that common blocks are saved throughout the program).
Use the mygetarg and myiargc routines rather than getarg and iargc (the former ones are platform-independent).

Add any new mistakes that brought you in trouble here and please send us an email, so we can extend this list for the benefit of all developers!

Version and Revision Numbers.

The MCTDH program project version number is listed in the include file versions.inc. This number has the form V.rxxx, where V is the version number, r the release, and xxx the revision number. Thus 8.1001 denotes the first revision of version 8 release 1.

This version number is important, as it not only allows the code to be identified, but it also allows the read-write and data file structures to be changed, while retaining backward compatibility.

If the project is changed in such a way that it should be differentiated from previous revisions, then the revision number should be incremented by 1. Such changes are the removal of bugs, minor changes in file formats, major re-writing of code etc. A more major change may then lead to a change in the release, or even version, number.

This information is written, e.g., to the log file and can be used to help narrow down problems if they arise after revisions.