{
		if(ptm->mat[i][j] != elem)
		  {
		    tmp->mat[k][0]=i;
		    tmp->mat[k][1]=j;
                    k++;
		  }
	      }
	  }
	if(k==0)   /* No entries different from the element */
	{
	  retmat = mmake(1,1);
	  retmat->row=0;
	  retmat->col=0;
	}          /* Copy coordinates to a [ 2 | # of instances] matrix */
	else
	  {
	  retmat = mmake(k,2);
	  for(i=0; i < k; i++)
	  {
	    retmat->mat[i][0] = tmp->mat[i][0];
	    retmat->mat[i][1] = tmp->mat[i][1];
	  }
	}
	mfree(tmp);
	return retmat;
}



/* 
 * MDIAG :
 * -------
 * This function initializes a matrix to all zeros except for the
 * main diagonal elements, which will be assigned the specified value.
 * Inputs : ptm - Pointer to matrix that will be initialized.
 *          diag - Diagonal element value.
 */
void mdiag( matrix *ptm, double diag )
{
	register int i;

#if RUNCHK

	if (ptm->row != ptm->col) \
	merror("matrix not quadratic error in mdiag!");

#endif

	minit(ptm);
	for ( i=0; i < ptm->row; i++ ) ptm->mat[i][i] = diag;
}


/*
 * MUNIT :
 * -------
 * This function initializes a matrix to the unit or identity matrix.
 * Input : ptm - Pointer to the matrix.
 */
void munit( matrix *ptm )
{
	register int i;

#if RUNCHK

	if (ptm->row != ptm->col)
		merror("matrix not quadratic error in mdiag!");

#endif

	minit(ptm);
	for ( i=0; i < ptm->row; i++ ) ptm->mat[i][i] = 1.0;
}


/*
 * MSET :
 * ------
 * This function copies the content of a matrix to another matrix.
 * ( A := B ) == mset(a,b); Assignment of matrix.
 * Inputs : ptm1, ptm2 - Pointers to matrices.
 */
void mset( matrix *ptm1, matrix *ptm2 )
{
	register int bytes;

#if RUNCHK

	if ( !( (ptm1->row == ptm2->row) && (ptm1->col == ptm2->col) ) ) \
	merror("Dimension mismatch error in mset!");

#endif

	bytes = 8*(ptm1->row)*(ptm1->col);
	memcpy( ptm1->mat[0], ptm2->mat[0], bytes );
}


/*
 * MADD
 * ----
 *
 * Matrix addition
 *
 * Input:  *ptm1 - Pointer to result matrix
 *                 (can be equal to *ptm1 and/or *ptm2)
 *         *ptm2 - Pointer to first argument matrix
 *         *ptm3 - Pointer to second argument matrix
 */
void madd( matrix *ptm1, matrix *ptm2, matrix *ptm3 )
{
	register int i,j;

#if RUNCHK

	if (!((ptm1->row == ptm2->row) && (ptm2->row == ptm3->row))) \
	merror("Dimension mismatch error in madd!");
	if (!((ptm1->col == ptm2->col) && (ptm2->col == ptm3->col))) \
	merror("Dimension mismatch error in madd!");

#endif

	/* Add the two matrices element by element */
	for ( i=0; i < ptm2->row; i++ )
	{
		for ( j=0; j < ptm2->col; j++ )
		{
			ptm1->mat[i][j] = ptm2->mat[i][j] + ptm3->mat[i][j];
		}
	}
}


/*
 * MSUB
 * ----
 *
 * Matrix subtraction
 *
 * Input:  *ptm1 - Pointer to result matrix
 *                 (can be equal to *ptm1 and/or *ptm2)
 *         *ptm2 - Pointer to first argument matrix
 *         *ptm3 - Pointer to second argument matrix
 */
void msub( matrix *ptm1, matrix *ptm2, matrix *ptm3 )
{
	register int i,j;

#if RUNCHK

	if (!((ptm1->row == ptm2->row) && (ptm2->row == ptm3->row))) \
	merror("Dimension mismatch error in msub!");
	if (!((ptm1->col == ptm2->col) && (ptm2->col == ptm3->col))) \
	merror("Dimension mismatch error in msub!");

#endif

	/* Subtract the two matrices element by element */
	for ( i=0; i < ptm2->row; i++ )
	{
		for ( j=0; j < ptm2->col; j++ )
		{
			ptm1->mat[i][j] = ptm2->mat[i][j] - ptm3->mat[i][j];
		}
	}
}


/*
 * MMUL
 * ----
 *
 * Matrix multiplication
 *
 * Input:  *ptm1 - Pointer to result matrix (Not equal to *ptm1 or *ptm2) 
 *         *ptm2 - Pointer to first argument matrix
 *         *ptm3 - Pointer to second argument matrix
 *
 */
void mmul( matrix *ptm1, matrix *ptm2, matrix *ptm3 )
{
	register int i,j,k;

#if RUNCHK

	if ((ptm1==ptm2) || (ptm1==ptm3)) \
	merror("Memory conflict error in mmul!");

	if ( !( (ptm2->col == ptm3->row) && \
	(ptm2->row == ptm1->row) && (ptm3->col == ptm1->col) ) ) \
	merror("Dimension mismatch error in mmul!");

#endif

	for ( i=0; i < ptm2->row; i++ )
	{
		for ( j=0; j < ptm3->col; j++ )
		{
			ptm1->mat[i][j] = 0.0;
			for ( k=0; k < ptm2->col; k++ )
			{
				ptm1->mat[i][j] += ptm2->mat[i][k] * ptm3->mat[k][j];
			}
		}
	}
}


/*
 * SMUL
 * ----
 *
 * Multiply matrix with scalar
 *
 * Input:  *ptm1  - Pointer to result matrix
 *         *ptm2  - Pointer to argument matrix
 *         factor - Scalar factor to be multiplied with argument matrix  
 */
void smul( matrix *ptm1, matrix *ptm2, double factor)
{
	register int i,j;

#if RUNCHK

	if (!((ptm1->row==ptm2->row) && (ptm1->col==ptm2->col))) \
	merror("Dimension mismatch error in smul!");

#endif

	for( i=0; i < ptm2->row; i++ )
	{
		for( j=0; j < ptm2->col; j++ )
		{
			ptm1->mat[i][j] = factor*ptm2->mat[i][j];
		}
	}
}


/*
 * TRACE :
 * -------
 * Calculates the trace (the sum of the main diagonal elements) of a matrix.
 * Input  : ptm   - Pointer to the matrix.
 * Output : trace - Trace of the matrix.
 */
double trace( matrix* ptm )
{
	register double trace = 0.0;
	register int i;

#if RUNCHK

	/* Check whether matrix is square */
	if (ptm->row != ptm->col) \
	merror("Matrix not square error in trace!");

#endif

	/* Calculate sum of diagonal elements */
	for ( i=0; i < ptm->row; i++)
	{
		trace += ptm->mat[i][i];
	}

	/* Return with trace */
	return trace;
}


/* SPROD :
 * -------
 * This function calculates the scalar-product of two vectors of equal
 * length. The vectors may be either row- or column- vectors or any 
 * combination of those. (The usual mmul routine can be used to calculate
 * scalar-products, but it requires a row- and a column-vector.)
 * Inputs : ptv1 - Pointer to first factor vector.
 *          ptv2 - Pointer to second factor vector.
 * Output : prod - Scalar product of the two vectors.
 */
double sprod( matrix* ptv1, matrix* ptv2 )
{
	register double prod = 0.0;
	register int i, elements;

#if RUNCHK

	/* Check whether the arguments are vectors.   */
	if (!(((ptv1->row==1)||(ptv1->col==1))&&((ptv2->row==1)||(ptv2->col==1)))) \
	merror("One of the inputs is not a vector error in sprod!");

	/* Check whether vectors has the same length. */
	if ((ptv1->row + ptv1->col) != (ptv2->row + ptv2->col)) \
	merror("Dimension mismatch error in sprod!");

#endif
	/* Calculate scalar product of the vectors.   */

	elements = ptv1->row + ptv1->col - 1;

	for ( i=0; i<elements; i++ ) {
		prod += (vget(ptv1,i))*(vget(ptv2,i));
	}

	return prod;
}


/* SPROD2 :
 * --------
 * This function calculates the scalar-product of two vectors of equal
 * length. The vectors may be either row- or column- vectors or any 
 * combination of those. (The usual mmul routine can be used to calculate
 * scalar-products, but it requires a row- and a column-vector.) In this
 * function, a starting and ending index can be specified for each vector.
 * In which case, only the specified part of the vectors will be used for
 * calculating the scalar-product.
 *
 * Inputs : ptv1 - Pointer to first factor vector.
 *          ptv2 - Pointer to second factor vector.
 *          begX - Beginning index of vectors 1 and 2.
 *          endX - Ending index of vectors 1 and 2.
 * Output : prod - Scalar product of the two vectors.
 */
double sprod2( matrix* ptv1, matrix* ptv2, int beg1, int end1, int beg2, int end2 )
{
	register double prod = 0.0;
	register int i, elements;

#if RUNCHK

	/* Check whether the arguments are vectors.   */
	if (!(((ptv1->row==1)||(ptv1->col==1))&&((ptv2->row==1)||(ptv2->col==1)))) \
	merror("One of the inputs is not a vector error in sprod!");

	/* Check whether vectors has the same length. */
	if ((end1 - beg1) != (end2 - beg2)) \
	merror("Dimension mismatch error in sprod!");

#else
	end2 = end2;  /* Avoid a warning. */
	
#endif
	/* Calculate scalar product of the vectors.   */

	elements = end1 - beg1 + 1;

	for ( i=0; i<elements; i++ ) {
		prod += (vget(ptv1,i+beg1))*(vget(ptv2,i+beg2));
	}

	return prod;
}

	
/*
 * MPUT :
 * ------
 * This function assigns a specified value to a specified element in
 * a matrix.
 * Inputs : ptm     - Pointer to the matrix.
 *          value   - Value to assign in matrix.
 *          row_pos - The row in which the element is to be assigned.
 *          col_pos - The column in which the element is to be assigned.
 */
void mput( matrix *ptm, int row_pos, int col_pos, double value )
{

#if RUNCHK

	/* Check if indices are inside matrix bounds. */
	if ( !((row_pos >= 0) && (row_pos < ptm->row)) ) \
	merror("Index out of range error in mput!");
	if ( !((col_pos >= 0) && (col_pos < ptm->col)) ) \
	merror("Index out of range error in mput!");

#endif

	/* Assign the value to the element */
	ptm->mat[row_pos][col_pos] = value;
}

	
/*
 * VPUT :
 * ------
 * This function assigns a specified value to a specified element in
 * a vector.
 * Inputs : ptv   - Pointer to the vector.
 *          value - Value to assign in matrix.
 *          pos   - The position in which the element is to be assigned.
 */
void vput( matrix *ptv, int pos, double value )
{

#if RUNCHK

	/* Check if index is inside vector bounds. */
	if ( !((pos>=0) && (pos <= ptv->row + ptv->col - 2)) ) \
	merror("Index out of range error in vput!");

#endif

	/* Assign the value to the element */
	if (ptv->row == 1)
	{
		ptv->mat[0][pos] = value; /* Row vector */
	}
	else if (ptv->col == 1)
	{
		ptv->mat[pos][0] = value; /* Column vector */
	}
	else merror("Input is not a vector error in vput!");

}



/*
 * MGET :
 * ------
 * This function returns a pointer to a specified element in a matrix.
 * Inputs : ptm     - Pointer to the matrix.
 *          row_pos - row of element.
 *          col_pos - column of element.
 */
double mget( matrix *ptm, int row_pos, int col_pos )
{

#if RUNCHK

	/* Check whether indices is inside matrix bounds. */
	if ( !((row_pos >= 0) && (row_pos < ptm->row)) ) \
	merror("Index out of range error in mget!");
	if ( !((col_pos >= 0) && (col_pos < ptm->col)) ) \
	merror("Index out of range error in mget!");

#endif

	/* Get address of specified element */
	return (ptm->mat[row_pos][col_pos]);
}


/*
 * VGET :
 * ------
 * Gets an element of either a column- or row-vector.
 * Inputs : ptv - Pointer to the vector.
 *          pos - Position of element.
 * Output : pte - Pointer to the specified element.
 */
double vget( matrix* ptv, int pos )
{
	register double *pte;

#if RUNCHK

	/* Check if index is inside matrix bounds. */
	if ( pos > ptv->row + ptv->col - 2 ) \
	merror("Index out of range error in vget!");

#endif

	if ( ptv->row == 1 )
	{
		pte = &(ptv->mat[0][pos]);
	}
	else if ( ptv->col == 1 )
	{
		pte = &(ptv->mat[pos][0]);
	}
	else merror("Input is not a vector in vget!");

	return *pte;
}


/*
 * ADDCOLS :
 * ---------
 * Create a new matrix, by putting two matrices next to one another
 *
 *                            [    :    ]
 *                       M1 = [ M2 : M3 ]
 *                            [    :    ]
 *
 * OBS: The number of rows in M1, M2 and M3 must be the same, and the number 
 *      of columns in M1 must be equal to the sum of rows in M2 and M3.
 *
 * Inputs:   *ptm1       : Pointer to result matrix (can be equal to *ptm2)
 *           *ptm2       : Pointer to the matrix to be placed left
 *           *ptm3       : Pointer to the matrix to be placed right
 */
void addcols(matrix *ptm1, matrix *ptm2, matrix *ptm3)
{
	register int i, j;

#if RUNCHK

	/* Check whether M1, M2 and M3 have the same number of rows        */
	if ( !( (ptm2->row == ptm3->row) && (ptm2->row == ptm1->row) ) ) \
	merror("Dimension mismatch in addcols!");
	
	/* Check that there aren't too many columns in M1                  */
	if ( ptm1->col != (ptm2->col + ptm3->col) ) \
	merror("Number of columns out of range in addcols!");

#endif

	/* Copy M2 to the left side of M1 */
	for( j = 0; j < ptm2->col; j++ )
	{
		for( i = 0; i < ptm1->row; i++ ) 
		{
			ptm1->mat[i][j] = ptm2->mat[i][j];
		}
	}

	/* Copy M3 to the right side of M1 */
	for( j = 0; j < ptm3->col; j++)
	{
		for( i = 0; i < ptm1->row; i++ )
		{
			ptm1->mat[i][(j+ptm2->col)] = ptm3->mat[i][j];
		}
	}

}

/*
 * ADDROWS :
 * ---------
 * Create a new matrix, by putting two matrices on top of one another
 *
 *                            [ M2 ]
 *                       M1 = [----]
 *                            [ M3 ]
 *
 * OBS: The number of columns in M2 and M3 must be the same, and the size of 
 *      M1 must fit M2 and M3.
 *
 * Inputs:   *ptm1       : Pointer to result matrix 
 *           *ptm2       : Pointer to matrix to be placed in the top
 *           *ptm3       : Pointer to matrix to be placed in the bottom
 */
void addrows( matrix *ptm1, matrix *ptm2, matrix *ptm3 )
{
	register int i, j;

#if RUNCHK
       
	/* Check whether M1, M2 and M3 have the same number of columns     */
	if (!((ptm2->col == ptm3->col) && (ptm1->col == ptm2->col))) \
	merror("Dimension mismatch in addrows!");

	/* Check whether M1 has the right dimensions to fit M2 & M3        */
	if (ptm1->row != (ptm2->row + ptm3->row)) \
	merror("Dimension mismatch in addrows!");

#endif