 /*
  * Khoros: $Id$
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id$";
#endif

 /*
  * $Log$
  */

/*
 *   Copyright, 1991, The Regents of the University of California.
 *   This software was produced under a U.S. Government contract
 *   (W-7405-ENG-36) by the Los Alamos National Laboratory, which is
 *   operated by the University of California for the U.S. Department
 *   of Energy.  The U.S. Government is licensed to use, reproduce,
 *   and distribute this software.  Neither the Government nor the
 *   University makes any warranty, express or implied, or assumes
 *   any liability responsibility for the use of this software.
 */


#include <stdio.h>


/*
 *  SM_get_matrix ()  Allocate memory for a symmetric matrix.  The
 *                    matrix returned is compact (only the upper
 *                    portion of the matrix is actually allocated)
 *                    and needs to be treated with care!
 *
 *  NOTES:  This technique is useful for programming in C, simply
 *          because arrays are actually pointers to pointers to.....
 *          and for this reason, this code should not be implemented
 *          like this in say a FORTRAN program if speed of program
 *          execution is more important than size requirements.
 *
 *  WARNING:  BECAUSE OF THE WAY THE POINTERS ARE ARRANGED FOR THE
 *            RETURNED MATRIX, ONLY THE UPPER POSITIONS SHOULD BE
 *            ACCESSED AND THE LOWER PORTIONS SHOULD *NEVER* BE 
 *            ACCESSED !!!!!!!!!  THUS, THE COLUMN INDEX SHOULD
 *            *ALWAYS* BE GREATER THAN OR EQUAL TO THE ROW INDEX
 *            WHEN ACCESSING THE MATRIX.
 *
 *  Written by:  Patrick M. Kelly
 *
 *  Date:  10/3/90
 */

float **SM_get_matrix ( dim )

int dim;					/* The dimension of the matrix */

{
	int row;				/* Loop control */
	int place;			/* Position in contiguous memory */
	int size;				/* Size of allocated block */
	float *contig;			/* Matrix values */
	float **row_indices;	/* Row indices into contig */

	/*
	 *  STEP I:  Allocate the memory for ALL of the matrix values 
	 *           contiguously.  The number of values stored will 
	 *           be:
	 *
	 *               dim  +  (dim - 1)  +  (dim - 2)  +  ....  +  1
	 *           OR
	 *               (dim + 1) * dim / 2
	 */
	size = ( ( dim + 1 ) * dim ) / 2;
	contig = (float *) malloc ( (unsigned) ( sizeof (float) * size ) );
	if ( contig == NULL ) {
		fprintf (stderr,"SM_get_matrix ()  MEMORY ALLOCATION FAILURE\n");
		exit ( -1 );
	}

	/*
	 *  STEP II:  Allocate memory for row indices.
	 */
	row_indices = (float **) malloc ((unsigned) (sizeof (float *) * dim ));
	if ( row_indices == NULL ) {
		fprintf (stderr,"SM_get_matrix ()  MEMORY ALLOCATION FAILURE\n");
		exit ( -1 );
	}

	/*
	 *  STEP III:  Set up the row pointers into the contiguous memory
	 *             block.
	 */
	for ( row = 0 ; row < dim ; row ++ ) {
		place = ( row * ( 2 * dim - row - 1 ) ) / 2 ;
		row_indices [row] = & ( contig [ place ] ) ;
	}

	return ( row_indices );
}
