/* 
 This program demonstrates a technique for allocating 
 memory to a 2-dimensional matrix which allows you to 
 use double subscripting when the pointer to the 
 matrix is passed to a function.  
*/

main()
{
  double **x,*xnow;
  long i,j;
  long m=3,n=5;
  double findelt(double **,long,long),findelt1(double *,long,long,long);

/*
    As an example, a 3 x 5 matrix is allocated.  Notice that
    by allocating all the memory for x in a single malloc call 
    and then setting the row pointers, a dereferenced version
    of the matrix (either *x or x[0]) can be passed to routines
    which expect the matrix to be stored as a vector.  One
    caveat:  FORTRAN stores it's matrices by columns, while this
    scheme forces the matrix to be stored by rows.  This needs
    to be considered if any FORTRAN subroutines or functions are
    used with matrices.                                    

    The first allocation only provides memory for the row pointers,
    *not* for the actual data values which will be stored in x. 
*/

  x = (double**)malloc((unsigned)(m * sizeof(double*)));
  
/*
    The allocation to xnow is for the entire m by n array - then the
    previous allocated row pointers can be set to beginnings of each
    of the m rows, by incrementing the xnow pointer by n elements
*/

  xnow = (double*)malloc((unsigned)(m * n * sizeof(double)));
  for(i=0;i<3;i++,xnow += n)x[i] = xnow;

/*
    To test this technique, 10 * i + j is put into the i,j th
    element of x, and then retreived via two different functions 
*/

    for(i=0;i<m;i++)for(j=0;j<n;j++)x[i][j] = (double)(10 * i + j);

    printf("1,3 element from findelt = %g\n",findelt(x,1,3));
    printf("2,4 element from findelt1 = %g\n",findelt1(*x,2,4,n));

/* 
    It's left as an exercise to write a function which would 
    return the properly allocated memory, i.e. the function 
    would be passed the dimension of the matrix, and it would 
    allocate the contiguous memory and return a double ** 
    set up as described above 
*/

}

/*  
    There's no need to pass the dimensions of the matrix to a routine
    which is being passed a double**, since neither the compiler or the 
    programmer has to do any subscripting computations. 
*/
    
double findelt(mat,i,j)
 double **mat;
 long i,j;
{
  return(mat[i][j]);
}

/*
    Notice that when you use single subscripts to store matrices,
    you must always know the number of columns (or rows, if the 
    matrix is stored by columns) in order to properly find the 
    element
*/

double findelt1(mat,i,j,nc)
 double *mat;
 long i,j,nc;

{
  return(mat[i * nc + j]);
}