Variable size multi-dimensional arrays

Multi-dimensional arrays crop up in a wide variety of different applications in scientific programming. Moreover, it is very common for the sizes of the arrays employed in a scientific code to vary from run to run, depending on the particular values of the input parameters. For instance, in a standard fluid code the array sizes depend on the number of grid points, which, in turn, depends on the requested accuracy. Thus, a crucial test of the suitability of a programming language for scientific purposes is its ability to deal with variable size, multi-dimensional arrays in a convenient manner. Unfortunately, C fails this test rather badly, since variable size matrix declarations of the form

void function(a, m, n) 
{
  int m, n;
  double a[m][n];
  . . .

are illegal (unlike in FORTRAN 77, where variable size matrix declarations are perfectly acceptable). Indeed, this unfortunate (and quite unnecessary!) omission in the C language definition is often put forward as a reason why C should not be used for scientific purposes. Fortunately, we can get around the absence of variable size multi-dimensional arrays in C by making use of a freely available C++ package called the Blitz++ library--see http://www.oonumerics.org/blitz/.

The program listed below illustrates the use of the Blitz++ library. The program adds together two matrices whose dimensions and elements are input by the user, and then prints out the result.

/* addmatrix.c */
/*
  Program to add two variable dimension matrices input by user
*/

#include <stdio.h>
#include <stdlib.h>
#include <blitz/array.h>

using namespace blitz;

/* Function prototypes */
void readin(Array<double,2>);
void writeout(Array<double,2>);
void addmatrices(Array<double,2>, Array<double,2>, Array<double,2>);

int main() 
{

  int n, m;
  
  /* Input number of rows and columns */
  printf("\nPlease input number of rows, n, and number of columns, m: ");
  scanf("%d %d", &n, &m);

  /* Check that n, m are positive integers */
  if (n <= 0 || m <= 0) 
    {
      printf("\nError: invalid values for n and/or m\n");
      exit(1);
    }

  /* Array declarations */
  Array<double,2> A(n, m), B(n, m), C(n, m);
     
  /* Read in elements of A, row by row */
  printf("\nReading in elements of A:\n");
  readin(A);

  /* Read in elements of B, row by row */
  printf("\nReading in elements of B:\n");
  readin(B);

  /* Write out elements of A, row by row */
  printf("\nWriting out elements of A:\n");
  writeout(A);

  /* Write out elements of B, row by row */
  printf("\n\nWriting out elements of B:\n");
  writeout(B);

  /* Add matrices A and B */
  addmatrices(A, B, C);

  /* Write out matrix C = A + B, row by row */
  printf("\n\nWriting out elements of C = A + B:\n");
  writeout(C);
  printf("\n");

  return 0;
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

/*
  Read in elements of matrix M, row by row
*/  

void readin(Array<double,2> M) 
{
  int n = M.extent(0);
  int m = M.extent(1);

  for (int i = 0; i < n; i++) 
    {
      printf("\nRow %d: ", i + 1);
      for (int j = 0; j < m; j++) scanf("%lf", &M(i, j)); 
    }
  return;
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
/*
  Write out elements of matrix M, row by row
*/  

void writeout(Array<double,2> M) 
{ 
  int n = M.extent(0);
  int m = M.extent(1);

  for (int i = 0; i < n; i++)
    {
      printf("\nRow %d: ", i + 1);
      for (int j = 0; j < m; j++) printf("%7.2f ", M(i, j)); 
    }
  return;
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
/*
    Add matrices M and N and store result in matrix P
*/

void addmatrices(Array<double,2> M, Array<double,2> N, Array<double,2> P) 
{ 
  int n = M.extent(0);
  int m = M.extent(1);
  
  for (int i = 0; i < n; i++) 
    for (int j = 0; j < m; j++) 
      P(i, j) = M(i, j) + N(i, j); 

   return;
}

Typical output from the program looks like

Please input number of rows, n, and number of columns, m: 3 3

Reading in elements of A:

Row 1: 1 4 5

Row 2: 6 7 8

Row 3: 3 6 0

Reading in elements of B:

Row 1: 1 6 0

Row 2: 4 7 8

Row 3: 4 8 8

Writing out elements of A:

Row 1:    1.00    4.00    5.00 
Row 2:    6.00    7.00    8.00 
Row 3:    3.00    6.00    0.00 

Writing out elements of B:

Row 1:    1.00    6.00    0.00 
Row 2:    4.00    7.00    8.00 
Row 3:    4.00    8.00    8.00 

Writing out elements of C = A + B:

Row 1:    2.00   10.00    5.00 
Row 2:   10.00   14.00   16.00 
Row 3:    7.00   14.00    8.00

The header file for the Blitz++ library is called blitz/array.h. Moreover, any program file which makes use of Blitz++ must include the cryptic line using namespace blitz; before the first call, otherwise compilation errors will ensue. The declaration Array<double,2> A(n, m) declares a two-dimensional n by m array of doubles. The generalization to an array of integers, or a higher dimension array, is fairly obvious. Note that the i, j element of matrix A is refered to simply as A(i,j). The function call M.extent(0) returns the size of array M in its first dimension. Likewise, M.extent(1) returns the size of array M in its second dimension, etc. More details of the operation of Blitz++ can be found by reading the extensive documentation which accompanies this library. Unfortunately, the Blitz++ library slows down compilation considerably since it makes use of some very advanced templating features of the C++ language.