#include <stdio.h>

#ifndef vms
# include <stdlib.h>	/* for malloc() and calloc() */
#endif

#define X(i,j) x[i][j]
#define Y(i,j) y[i][j]
#define Z(i,j) z[i][j]

void matmul();

/*
 * matcreate: allocates space for m-by-n matrix, returns pointer if
 *  successful, NULL if not.  The matrix elements are initialized to zero.
 *  The matrix elements are contiguous, since space for all of the rows
 *  is obtained with one calloc() call, thus the matrix can be passed
 *  to f77 or pascal routines. a[-1] (as well as a[0]) is set to point to the
 *  storage so that rows may be interchanged by swapping row pointers, and
 *  matfree() can still reference the storage through a[-1] even if a[0] was
 *  changed.
 */
double **matcreate( m, n)
int m;	/* number of rows */
int n;	/* number of columns */
{
    double **a;		/* temporary pointer */
    int i;		/* row index */

    if( m < 1 || n < 1)
	return NULL;
    /*
     * first, get space for m+1 pointers to double
     */
    if( (a = (double **) malloc( (m+1) * sizeof(double *)) ) == NULL)
	return NULL;
    /*
     * get space for all of the rows, m*n doubles
     */
    if( (a[0] = a[1] = (double *) calloc( m*n, sizeof(double)) ) == NULL) {
	free(a);
	return NULL;
    }
    /*
     * set the pointers a[1]..a[m-1] to point to their row
     */
    ++a;
    for( i = 1; i < m; i++)
	a[i] = a[i-1] + n;

    return a;
}

void matfree(a)
double **a;
{
    if( a) {
	free(a[-1]);	/* free the matrix data */
	free(a-1);	/* free the row pointers */
    }
}

main()
{
    double **x, **y, **z;
    int i, j;

    x = matcreate(100,100);
    y = matcreate(100,100);
    z = matcreate(100,100);

    if( !x || !y || !z) {
	printf("matcreate() failed");
	exit(1);
    }

    for( i = 0; i < 100; ++i)
	for( j = 0; j < 100; ++j)
	    Y(i,j) = Z(i,j) = 1;

    matmul( x, y, z, 100, 100, 100);

    exit(0);
}