/*
    This file is part of the FElt finite element analysis package.
    Copyright (C) 1993-2000 Jason I. Gobat and Darren C. Atkinson

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/

/************************************************************************
 * File:	matrixfunc.c						*
 *									*
 * Description:	This file contains the function definitions for the	*
 *		matrix-related intrinsic functions.			*
 ************************************************************************/

# include <math.h>
# include <time.h>
# include "debug.h"
# include "error.h"
# include "coerce.h"
# include "execute.h"
# include "matrixfunc.h"
# include "our-stdlib.h"

# ifdef DOS
# define srand48 srand
# define drand48 rand
# endif

/************************************************************************
 * Function:	chol_func						*
 *									*
 * Description:	Pops and computes the cholesky factorization of the	*
 *		descriptor on the top of the stack and places the	*
 *		result on the stack.  The following types are legal:	*
 *									*
 *		chol (double) -> double (scalar square root)		*
 *		chol (matrix) -> matrix (cholesky factorization)	*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int chol_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    double	value;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    int		status;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);


    status = 0;
    type_error = F_False;


    switch (D_Type (arg)) {
    case T_Double:
	if ((value = *D_Double (arg)) < 0.0) {
	    MathException ("cholesky factorization of negative number");
	    status = 1;
	} else {
	    CreateData (result, arg, NULL, T_Double);
	    *D_Double (result) = sqrt (value);
	}
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	if ((status = CholeskyFactorMatrix (b, a)))
	    MatrixError ("chol", a, NULL, status, F_True);
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("chol", arg, NULL, NULL, F_True);


    RecycleData (arg);
    d_printf ("chol ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	cols_func						*
 *									*
 * Description:	Pops the descriptor on the top of the stack and pushes	*
 *		its number of columns on the stack.  The following	*
 *		types are legal:					*
 *									*
 *		cols (double) -> double (always 1)			*
 *		cols (matrix) -> double (number of columns)		*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int cols_func (n)
    int n;
{
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;


    switch (D_Type (arg)) {
    case T_Double:
	D_Type    (result) = T_Double;
	D_Temp    (result) = F_False;
	D_Trapped (result) = F_False;
	D_Double  (result) = dbllit (1);
	break;


    case T_Matrix:
	CreateData (result, NULL, NULL, T_Double);
	*D_Double (result) = Mcols (D_Matrix (arg));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("cols", arg, NULL, NULL, F_True);


    RecycleData (arg);
    d_printf ("cols ans =\n");
    d_PrintData (result);

    return type_error == F_True;
}


/************************************************************************
 * Function:	compact_func						*
 *									*
 * Description:	Pops and compacts the descriptor on the top of the	*
 *		stack and pushes the result on the stack.  The		*
 *		following types are legal:				*
 *									*
 *		compact (double) -> double (identity operation)		*
 *		compact (matrix) -> matrix (matrix compaction)		*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int compact_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    int		status;


    result = top ( );
    temp = *result;
    arg = &temp;


    status = 0;
    type_error = F_False;
    arg = CoerceData (deref (arg), T_Double);


    switch (D_Type (arg)) {
    case T_Double:
	CreateData (result, arg, NULL, T_Double);
	*D_Double (result) = *D_Double (arg);
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	if (IsCompact (a)) {
	    CreateData (result, NULL, NULL, T_Matrix, Mrows (a), Mcols (a));
	    CopyMatrix (D_Matrix (result), a);
	} else if (IsSymmetricMatrix (a)) {
	    b = MakeCompactFromFull (a);
	    D_Type    (result) = T_Matrix;
	    D_Temp    (result) = F_True;
	    D_Trapped (result) = F_False;
	    D_Matrix  (result) = b;
	} else {
	    status = M_NOTSYMMETRIC;
	    MatrixError ("compact", a, NULL, status, F_True);
	}
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error)
	TypeError ("compact", arg, NULL, NULL, F_True);


    RecycleData (arg);
    d_printf ("compact ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	det_func						*
 *									*
 * Description:	Pops and computes the determinant of the descriptor on	*
 *		the top of the stack and places the result on the	*
 *		stack.  The following types are legal:			*
 *									*
 *		det (double) -> double (identity operation)		*
 *		det (matrix) -> matrix (matrix determinant)		*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int det_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    Matrix	p;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    int		status;
    int		info;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;
    status = 0;


    switch (D_Type (arg)) {
    case T_Double:
	CreateData (result, arg, NULL, T_Double);
	*D_Double (result) = *D_Double (arg);
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	p = CreateColumnVector (Mrows (a));
	b = D_Writable (arg) ? a : CreateFullMatrix (Mrows (a), Mcols (a));

	if ((status = LUFactorMatrix (b, a, p, &info)))
	    MatrixError ("det", a, NULL, status, F_True);
	else if (info) {
	    MathException ("singular matrix in det() function");
	    status = 1;
	} else {
	    CreateData (result, NULL, NULL, T_Double);
	    if ((status = DeterminantMatrix (D_Double (result), b, p)))
		MatrixError ("det", b, NULL, status, F_True);
	}

	DestroyMatrix (p);
	if (a != b)
	    DestroyMatrix (b);
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("det", arg, NULL, NULL, F_True);


    RecycleData (arg);
    d_printf ("det ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	eig_func						*
 *									*
 * Description:	Pops and computes the eigenvalue decomposition of the	*
 *		descriptor on the top of the stack and places the	*
 *		result on the stack.  The result is the vector of	*
 *		eigenvalues.  The matrix of eigenvectors of a symmetric	*
 *		matrix may be obtained by specifying it as the		*
 *		remaining parameter.  The following types are legal:	*
 *									*
 *		eig (double, ...) -> double (identity operation)	*
 *		eig (matrix, ...) -> matrix (eigenvalue decomposition)	*
 *									*
 *		An attempt is first made to coerce the first argument	*
 *		to a double value.					*
 ************************************************************************/

int eig_func (n)
    int n;
{
    Matrix	a;
    Matrix	l;
    Matrix	v;
    descriptor *arg1;
    descriptor *arg2;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    int		status;


    arg2 = pop ( );
    result = top ( );
    temp = *result;
    arg1 = &temp;


    status = 0;
    type_error = F_False;
    arg1 = CoerceData (deref (arg1), T_Double);


    if (!assignable (arg2)) {
	RecycleData (arg2);
	arg2 = NULL;
    } else
	arg2 = deref (arg2);


    switch (D_Type (arg1)) {
    case T_Double:
	CreateData (result, arg1, NULL, T_Double);
	*D_Double (result) = *D_Double (arg1);

	if (arg2)
	    status = AssignObject (arg2, T_Double, F_False, dbllit (1));
	break;


    case T_Matrix:
	a = D_Matrix (arg1);
	CreateData (result, NULL, NULL, T_Matrix, Mrows (a), 1);
	l = D_Matrix (result);

	if (IsSymmetricMatrix (a)) {
	    v = CreateFullMatrix (Mrows (a), Mcols (a));

	    if ((status = SymmetricMatrixEigenModes (a, l, v, 0.0))) {
		MatrixError ("eig", a, NULL, status, F_True);
		DestroyMatrix (v);

	    } else if (!arg2)
		DestroyMatrix (v);
	    else if ((status = AssignObject (arg2, T_Matrix, F_True, v)))
		DestroyMatrix (v);

	} else if ((status = GeneralMatrixEigenModes (a, l, 0.0, 0.0)))
	    MatrixError ("eig", a, NULL, status, F_True);
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("eig", arg1, NULL, NULL, F_True);


    RecycleData (arg1);
    d_printf ("eig ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	eye_func						*
 *									*
 * Description:	Pops two descriptors from the top of the stack and	*
 *		creates and pushes on the stack an identity matrix.	*
 *		The first descriptor specifies the number of rows and	*
 *		the second descriptor specifies the number of columns.	*
 *		The following types are legal:				*
 *									*
 *		eye (double, null)   -> matrix (matrix creation)	*
 *		eye (double, double) -> matrix (matrix creation)	*
 *									*
 *		An attempt is first made to coerce both arguments to	*
 *		double values.  The default value for the second	*
 *		argument is the value of the first argument.  If the	*
 *		resulting matrix consists of only a single scalar then	*
 *		the result will be a double scalar.			*
 ************************************************************************/

int eye_func (n)
    int n;
{
    Matrix	a;
    double	val1;
    double	val2;
    descriptor *arg1;
    descriptor *arg2;
    descriptor *result;
    descriptor *old;
    descriptor	temp;
    int		type_error;
    int		status;


    arg2 = pop ( );
    result = top ( );
    temp = *result;
    arg1 = &temp;


    arg1 = CoerceData (deref (arg1), T_Double);
    arg2 = CoerceData (deref (arg2), T_Double);


    status = 0;
    type_error = F_False;


    if (D_Type (arg2) == T_Null) {
	RecycleData (arg2);
	arg2 = arg1;
	old = NULL;
    } else
	old = arg2;


    if (D_Type (arg1) == T_Double && D_Type (arg2) == T_Double) {
	val1 = *D_Double (arg1);
	val2 = *D_Double (arg2);
	if (val1 < 1) {
	    rterror ("non-positive number of rows in eye() function");
	    status = 1;
	} else if (val2 < 1) {
	    rterror ("non-positive number of columns in eye() function");
	    status = 1;
	} else if ((int) val1 == 1 && (int) val2 == 1) {
	    D_Type    (result) = T_Double;
	    D_Temp    (result) = F_False;
	    D_Trapped (result) = F_False;
	    D_Double  (result) = dbllit (1);
	} else {
	    CreateData (result, NULL, NULL, T_Matrix, (int) val1, (int) val2);
	    a = D_Matrix (result);
	    if ((status = IdentityMatrix (a)))
		MatrixError ("eye", a, NULL, status, F_True);
	}
    } else
	type_error = F_True;


    if (type_error == F_True)
	TypeError ("eye", arg1, old, NULL, F_True);


    RecycleData (arg1);
    RecycleData (old);
    d_printf ("eye ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	inv_func						*
 *									*
 * Description:	Pops and computes the inverse of the descriptor on the	*
 *		top of the stack and places the result on the stack.	*
 *		The following types are legal:				*
 *									*
 *		inv (double) -> double (scalar inversion)		*
 *		inv (matrix) -> matrix (matrix inversion)		*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int inv_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    Matrix	p;
    Matrix	c;
    double	value;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    int		status;
    int		info;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;
    status = 0;


    switch (D_Type (arg)) {
    case T_Double:
	if ((value = *D_Double (arg)) == 0.0) {
	    MathException ("inverse of zero in inv() function");
	    status = 1;
	} else {
	    CreateData (result, arg, NULL, T_Double);
	    *D_Double (result) = 1.0 / *D_Double (arg);
	}
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	c = CreateFullMatrix (Mrows (a), Mcols (a));
	p = CreateColumnVector (Mrows (a));

	if ((status = LUFactorMatrix (c, a, p, &info)))
	    MatrixError ("inv", a, NULL, status, F_True);
	else if (info)
	    MathException ("singular matrix in inv() function"), status = 1;
	else if ((status = InvertMatrix (b, c, p)))
	    MatrixError ("inv", a, NULL, status, F_True);

	DestroyMatrix (c);
	DestroyMatrix (p);
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("inv", arg, NULL, NULL, F_True);


    RecycleData (arg);
    d_printf ("inv ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	lu_func							*
 *									*
 * Description:	Pops and computes the LU factorization of the		*
 *		descriptor on the top of the stack and places the	*
 *		result on the stack.  The result is a row permuted	*
 *		superposition of L and U, with the diagonal of L not	*
 *		being stored since it is always ones.  The individual	*
 *		L, U, and permutation matrices can be retrieved by	*
 *		specifying them as the remaining parameters.  The	*
 *		following types are legal:				*
 *									*
 *		lu (double, ...) -> double (identity operation)		*
 *		lu (matrix, ...) -> matrix (LU factorization)		*
 *									*
 *		An attempt is first made to coerce the first argument	*
 *		to a double value.					*
 ************************************************************************/

int lu_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    Matrix	p;
    Matrix	L;
    Matrix	U;
    Matrix	P;
    descriptor *arg1;
    descriptor *arg2;
    descriptor *arg3;
    descriptor *arg4;
    descriptor *result;
    descriptor	scratch;
    descriptor	temp;
    int		type_error;
    int		status;
    int		info;


    arg4 = pop ( );
    arg3 = pop ( );
    arg2 = pop ( );
    result = top ( );
    temp = *result;
    arg1 = &temp;


    status = 0;
    type_error = F_False;
    arg1 = CoerceData (deref (arg1), T_Double);


    if (!assignable (arg2)) {
	RecycleData (arg2);
	arg2 = NULL;
    } else
	arg2 = deref (arg2);

    if (!assignable (arg3)) {
	RecycleData (arg3);
	arg3 = NULL;
    } else
	arg3 = deref (arg3);

    if (!assignable (arg4)) {
	RecycleData (arg4);
	arg4 = NULL;
    } else
	arg4 = deref (arg4);


    switch (D_Type (arg1)) {
    case T_Double:
	CreateData (result, arg1, NULL, T_Double);
	*D_Double (result) = *D_Double (arg1);
	if (arg2)
	    status = AssignObject (arg2, T_Double, F_False, dbllit (1));

	if (arg3) {
	    CreateData (&scratch, NULL, NULL, T_Double);
	    *D_Double (&scratch) = *D_Double (arg1);
	    status = AssignObject (arg3, T_Double, F_True, D_Double (&scratch));
	    if (status)
		FreeData (&scratch);
	}

	if (arg4)
	    status = AssignObject (arg4, T_Double, F_False, dbllit (1));
	break;


    case T_Matrix:
	a = D_Matrix (arg1);
	CreateData (result, arg1, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	p = CreateColumnVector (Mrows (a));

	if ((status = LUFactorMatrix (b, a, p, &info)))
	    MatrixError ("lu", a, NULL, status, F_True);
	else if (info)
	    MathException ("singular matrix in lu() function"), status = 1;
	else if (arg2 || arg3 || arg4) {
	    L = CreateFullMatrix (Mrows (a), Mcols (a));
	    U = CreateFullMatrix (Mrows (a), Mcols (a));
	    P = CreateFullMatrix (Mrows (a), Mcols (a));

	    if ((status = FormLUPMatrices (L, U, P, b, p)))
		MatrixError ("lu", a, NULL, status, F_True);

	    if (!arg2 || (status = AssignObject (arg2, T_Matrix, F_True, L)))
		DestroyMatrix (L);

	    if (!arg3 || (status = AssignObject (arg3, T_Matrix, F_True, U)))
		DestroyMatrix (U);

	    if (!arg4 || (status = AssignObject (arg4, T_Matrix, F_True, P)))
		DestroyMatrix (P);
	}

	DestroyMatrix (p);
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("lu", arg1, NULL, NULL, F_True);


    RecycleData (arg1);
    d_printf ("lu ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	norm_func						*
 *									*
 * Description:	Pops two descriptors from the top of the stack and	*
 *		places the norm of the first on the stack.  The second	*
 *		descriptor specifies the type of norm.  The following	*
 *		types are legal:					*
 *									*
 *		norm (double, string) -> double (absolute value)	*
 *		norm (matrix, string) -> matrix (matrix or vector norm)	*
 *									*
 *		An attempt is first made to coerce the first argument	*
 *		to a double value and the second argument to a string	*
 *		value.  The value of the result is based on the values	*
 *		of the arguments as follows:				*
 *									*
 *			null		"1"	"2"	"fro"		*
 *		scalar	fabs		fabs	fabs	fabs		*
 *		vector	2-norm		1-norm	2-norm	frobenius-norm	*
 *		matrix	frobenius-norm	1-norm	error	frobenius-norm	*
 *									*
 *		The frobenius-norm of a vector is identical to the	*
 *		2-norm of the vector.					*
 ************************************************************************/

int norm_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    char       *string;
    descriptor *arg1;
    descriptor *arg2;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    int		status;
    int		norm;
    unsigned	i;


    arg2 = pop ( );
    result = top ( );
    temp = *result;
    arg1 = &temp;


    arg1 = CoerceData (deref (arg1), T_Double);
    arg2 = CoerceData (deref (arg2), T_String);


    norm = 0;
    status = 0;
    type_error = F_False;


    if (D_Type (arg2) != T_Null && D_Type (arg2) != T_String)
	type_error = F_True;
    else if (D_Type (arg2) == T_String) {
	string = *D_String (arg2);
	if (!strcmp (string, "1"))
	    norm = 1;
	else if (!strcmp (string, "2"))
	    norm = 2;
	else if (!strcmp (string, "inf"))
	    norm = 'i';
	else if (!strcmp (string, "fro"))
	    norm = 'f';
	else {
	    rterror ("illegal norm type in norm() function: %s", string);
	    status = 1;
	}
    }


    if (type_error == F_False && status == 0)
	switch (D_Type (arg1)) {
	case T_Double:
	    CreateData (result, arg1, NULL, T_Double);
	    *D_Double (result) = fabs (*D_Double (arg1));
	    break;


	case T_Matrix:
	    a = D_Matrix (arg1);
	    CreateData (result, NULL, NULL, T_Double);
	    if (Mrows (a) == 1) {
		b = CreateColumnVector (Mcols (a));
		for (i = 1; i <= Mcols (a); i ++)
		    sdata (b, i, 1) = mdata (a, 1, i);

		if (norm == 1) {
		    if ((status = PNormVector (D_Double (result), b, "1")))
			MatrixError ("norm", a, NULL, status, F_True);
		} else if (norm == 2 || norm == 0) {
		    if ((status = PNormVector (D_Double (result), b, "2")))
			MatrixError ("norm", a, NULL, status, F_True);
		} else if (norm == 'i') {
		    if ((status = PNormVector (D_Double (result), b, "inf")))
			MatrixError ("norm", a, NULL, status, F_True);
		} else if (norm == 'f') {
		    if ((status = FrobeniusNormMatrix (D_Double (result), b)))
			MatrixError ("norm", a, NULL, status, F_True);
		}

		DestroyMatrix (b);

	    } else if (Mcols (a) == 1) {
		if (norm == 1) {
		    if ((status = PNormVector (D_Double (result), a, "1")))
			MatrixError ("norm", a, NULL, status, F_True);
		} else if (norm == 2 || norm == 0) {
		    if ((status = PNormVector (D_Double (result), a, "2")))
			MatrixError ("norm", a, NULL, status, F_True);
		} else if (norm == 'i') {
		    if ((status = PNormVector (D_Double (result), a, "inf")))
			MatrixError ("norm", a, NULL, status, F_True);
		} else if (norm == 'f') {
		    if ((status = FrobeniusNormMatrix (D_Double (result), a)))
			MatrixError ("norm", a, NULL, status, F_True);
		}

	    } else {
		if (norm == 2) {
		    rterror ("2-norm of matrix undefined in norm() function");
		    status = 1;
		} else if (norm == 1) {
		    if ((status = PNormMatrix (D_Double (result), a, "1")))
			MatrixError ("norm", a, NULL, status, F_True);
		} else if (norm == 'i') {
		    if ((status = PNormMatrix (D_Double (result), a, "inf")))
			MatrixError ("norm", a, NULL, status, F_True);
		} else if (norm == 'f' || norm == 0) {
		    if ((status = FrobeniusNormMatrix (D_Double (result), a)))
			MatrixError ("norm", a, NULL, status, F_True);
		}
	    }
	    break;


	default:
	    type_error = F_True;
	    break;
	}


    if (type_error == F_True)
	TypeError ("norm", arg1, arg2, NULL, F_True);


    RecycleData (arg1);
    RecycleData (arg2);
    d_printf ("norm ans=\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	ones_func						*
 *									*
 * Description:	Pops two descriptors from the top of the stack and	*
 *		creates and pushes on the stack a matrix of ones.  The	*
 *		first descriptor specifies the number of rows and the	*
 *		second descriptor specifies the number of columns.  The	*
 *		following types are legal:				*
 *									*
 *		ones (double, null)   -> matrix (matrix creation)	*
 *		ones (double, double) -> matrix (matrix creation)	*
 *									*
 *		An attempt is first made to coerce the arguments to	*
 *		double values.  The default value for the second	*
 *		argument is the value of the first argument.  If the	*
 *		resulting matrix consists of only a single scalar then	*
 *		the result will be a double scalar.			*
 ************************************************************************/

int ones_func (n)
    int n;
{
    Matrix	a;
    double	val1;
    double	val2;
    descriptor *arg1;
    descriptor *arg2;
    descriptor *result;
    descriptor *old;
    descriptor	temp;
    int		type_error;
    int		status;


    arg2 = pop ( );
    result = top ( );
    temp = *result;
    arg1 = &temp;


    arg1 = CoerceData (deref (arg1), T_Double);
    arg2 = CoerceData (deref (arg2), T_Double);


    status = 0;
    type_error = F_False;


    if (D_Type (arg2) == T_Null) {
	RecycleData (arg2);
	arg2 = arg1;
	old = NULL;
    } else
	old = arg2;


    if (D_Type (arg1) == T_Double && D_Type (arg2) == T_Double) {
	val1 = *D_Double (arg1);
	val2 = *D_Double (arg2);
	if (val1 < 1) {
	    rterror ("non-positive number of rows in ones() function");
	    status = 1;
	} else if (val2 < 1) {
	    rterror ("non-positive number of columns in ones() function");
	    status = 1;
	} else if ((int) val1 == 1 && (int) val2 == 1) {
	    D_Type    (result) = T_Double;
	    D_Temp    (result) = F_False;
	    D_Trapped (result) = F_False;
	    D_Double  (result) = dbllit (0);
	} else {
	    CreateData (result, NULL, NULL, T_Matrix, (int) val1, (int) val2);
	    a = D_Matrix (result);
	    if ((status = ScaleMatrix (a, a, 0.0, 1.0)))
		MatrixError ("ones", a, NULL, status, F_True);
	}
    } else
	type_error = F_True;


    if (type_error == F_True)
	TypeError ("ones", arg1, old, NULL, F_True);


    RecycleData (arg1);
    RecycleData (old);
    d_printf ("ones ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	qr_func							*
 *									*
 * Description:	Pops and computes the QR factorization of the		*
 *		descriptor on the top of the stack and places the	*
 *		result on the stack.  The result is the R matrix.  The	*
 *		Q and R matrices may be obtained by specifying them as	*
 *		the remaining parameters.  The following types are	*
 *		legal:							*
 *									*
 *		qr (double, ...) -> double (identity operation)		*
 *		qr (matrix, ...) -> matrix (QR factorization)		*
 *									*
 *		An attempt is first made to coerce the first argument	*
 *		to a double value.					*
 ************************************************************************/

int qr_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    Matrix	q;
    Matrix	r;
    descriptor *arg1;
    descriptor *arg2;
    descriptor *arg3;
    descriptor *result;
    descriptor	scratch;
    descriptor	temp;
    int		type_error;
    int		status;
    unsigned	i;
    unsigned	j;
    unsigned	cols;
    unsigned	rows;


    arg3 = pop ( );
    arg2 = pop ( );
    result = top ( );
    temp = *result;
    arg1 = &temp;


    status = 0;
    type_error = F_False;
    arg1 = CoerceData (deref (arg1), T_Double);


    if (!assignable (arg2)) {
	RecycleData (arg2);
	arg2 = NULL;
    } else
	arg2 = deref (arg2);

    if (!assignable (arg3)) {
	RecycleData (arg3);
	arg3 = NULL;
    } else
	arg3 = deref (arg3);


    switch (D_Type (arg1)) {
    case T_Double:
	CreateData (result, arg1, NULL, T_Double);
	*D_Double (result) = *D_Double (arg1);

	if (arg2)
	    status = AssignObject (arg2, T_Double, F_False, dbllit (1));

	if (arg3) {
	    CreateData (&scratch, NULL, NULL, T_Double);
	    *D_Double (&scratch) = *D_Double (arg1);
	    status = AssignObject (arg3, T_Double, F_True, D_Double (&scratch));
	    if (status)
		FreeData (&scratch);
	}
	break;


    case T_Matrix:
	a = D_Matrix (arg1);
	rows = Mrows (a);
	cols = Mcols (a);
	CreateData (result, arg1, NULL, T_Matrix, rows, cols);
	b = D_Matrix (result);

	q = CreateFullMatrix (rows, cols);
	r = CreateFullMatrix (rows, cols);

	if ((status = QRFactorMatrix (q, r, a))) {
	    MatrixError ("qr", a, NULL, status, F_True);
	    DestroyMatrix (q);
	    DestroyMatrix (r);

	} else {
	    for (i = 2; i <= rows; i ++)
		for (j = 1; j < i; j ++)
		    sdata (r, i, j) = 0;

	    CopyMatrix (b, r);

	    if (!arg2 || (status = AssignObject (arg2, T_Matrix, F_True, q)))
		DestroyMatrix (q);

	    if (!arg3 || (status = AssignObject (arg3, T_Matrix, F_True, r)))
		DestroyMatrix (r);
	}
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("qr", arg1, NULL, NULL, F_True);


    RecycleData (arg1);
    d_printf ("qr ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	rand_func						*
 *									*
 * Description:	Pops three descriptors from the top of the stack and	*
 *		creates and pushes on the stack a random matrix.  The	*
 *		first descriptor specifies the number of rows and the	*
 *		second descriptor specifies the number of columns.  The	*
 *		third descriptor is the seed for the random number	*
 *		generator.  The following types are legal:		*
 *									*
 *		rand (null, null, null)	      -> scalar (random double)	*
 *		rand (double, null, null)     -> matrix (random matrix)	*
 *		rand (double, double, null)   -> matrix (random matrix)	*
 *		rand (double, double, double) -> matrix (random matrix)	*
 *									*
 *		An attempt is first made to coerce all arguments to	*
 *		double values.  The default value for the second	*
 *		argument is the value of the first argument.  If the	*
 *		third argument is null then the random generator is not	*
 *		seeded.  If the resulting matrix consists of only a	*
 *		single scalar then the result will be a double scalar.	*
 ************************************************************************/

int rand_func (n)
    int n;
{
    Matrix	a;
    double	val1;
    double	val2;
    double	val3;
    descriptor *arg1;
    descriptor *arg2;
    descriptor *arg3;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    int		status;
    static int  seeded;


    arg3 = pop ( );
    arg2 = pop ( );
    result = top ( );
    temp = *result;
    arg1 = &temp;


    arg1 = CoerceData (deref (arg1), T_Double);
    arg2 = CoerceData (deref (arg2), T_Double);
    arg3 = CoerceData (deref (arg3), T_Double);


    status = 0;
    val1 = val2 = val3 = 0;
    type_error = F_False;


    switch (D_Type (arg1)) {
    case T_Null:
	val1 = 1;
	arg1 = NULL;
	break;


    case T_Double:
	val1 = *D_Double (arg1);
	break;


    default:
	type_error = F_True;
	break;
    }

    switch (D_Type (arg2)) {
    case T_Null:
	val2 = val1;
	arg2 = NULL;
	break;


    case T_Double:
	val2 = *D_Double (arg2);
	break;


    default:
	type_error = F_True;
	break;
    }

    switch (D_Type (arg3)) {
    case T_Null:
	val3 = 0;
	arg3 = NULL;
	break;


    case T_Double:
	val3 = *D_Double (arg3);
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_False)
	if (val1 < 1) {
	    rterror ("non-positive number of rows in rand() function");
	    status = 1;
	} else if (val2 < 1) {
	    rterror ("non-positive number of columns in rand() function");
	    status = 1;
	} else if ((int) val1 == 1 && (int) val2 == 1) {
	    if (val3)
		srand48 ((int) val3);
	    else if (!seeded) {
		seeded = 1;
		srand48 (time (NULL));
	    }
	    CreateData (result, arg1, arg2, T_Double);
	    *D_Double (result) = drand48 ( );
	} else {
	    CreateData (result, NULL, NULL, T_Matrix, (int) val1, (int) val2);
	    a = D_Matrix (result);
	    if ((status = RandomMatrix (a, (int) val3)))
		MatrixError ("rand", a, NULL, status, F_True);
	}


    if (type_error == F_True)
	TypeError ("rand", arg1, arg2, arg3, F_True);


    RecycleData (arg1);
    RecycleData (arg2);
    RecycleData (arg3);
    d_printf ("rand ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	rows_func						*
 *									*
 * Description:	Pops the descriptor on the top of the stack and pushes	*
 *		its number of rows on the stack.  The following types	*
 *		are legal:						*
 *									*
 *		rows (double) -> double (always 1)			*
 *		rows (matrix) -> double (number of rows)		*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int rows_func (n)
    int n;
{
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;


    switch (D_Type (arg)) {
    case T_Double:
	D_Type    (result) = T_Double;
	D_Temp    (result) = F_False;
	D_Trapped (result) = F_False;
	D_Double  (result) = dbllit (1);
	break;


    case T_Matrix:
	CreateData (result, NULL, NULL, T_Double);
	*D_Double (result) = Mrows (D_Matrix (arg));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("rows", arg, NULL, NULL, F_True);


    RecycleData (arg);
    d_printf ("rows ans =\n");
    d_PrintData (result);

    return type_error == F_True;
}


/************************************************************************
 * Function:	zeros_func						*
 *									*
 * Description:	Pops two descriptors from the top of the stack and	*
 *		creates and pushes on the stack a zero matrix.  The	*
 *		first descriptor specifies the number of rows and the	*
 *		second descriptor specified the number of columns.  The	*
 *		following types are legal:				*
 *									*
 *		zeros (double, null)   -> matrix (matrix creation)	*
 *		zeros (double, double) -> matrix (matrix creation)	*
 *									*
 *		An attempt is first made to coerce the arguments to	*
 *		double values.   The default value for the second	*
 *		argument is the value of the first argument.  If the	*
 *		resulting matrix consists of only a single scalar then	*
 *		the result will be a double scalar.			*
 ************************************************************************/

int zeros_func (n)
    int n;
{
    Matrix	a;
    double	val1;
    double	val2;
    descriptor *arg1;
    descriptor *arg2;
    descriptor *result;
    descriptor *old;
    descriptor	temp;
    int		type_error;
    int		status;


    arg2 = pop ( );
    result = top ( );
    temp = *result;
    arg1 = &temp;


    arg1 = CoerceData (deref (arg1), T_Double);
    arg2 = CoerceData (deref (arg2), T_Double);


    status = 0;
    type_error = F_False;


    if (D_Type (arg2) == T_Null) {
	RecycleData (arg2);
	arg2 = arg1;
	old = NULL;
    } else
	old = arg2;


    if (D_Type (arg1) == T_Double && D_Type (arg2) == T_Double) {
	val1 = *D_Double (arg1);
	val2 = *D_Double (arg2);
	if (val1 < 1) {
	    rterror ("non-positive number of rows in zeros() function");
	    status = 1;
	} else if (val2 < 1) {
	    rterror ("non-positive number of columns in zeros() function");
	    status = 1;
	} else if ((int) val1 == 1 && (int) val2 == 1) {
	    D_Type    (result) = T_Double;
	    D_Temp    (result) = F_False;
	    D_Trapped (result) = F_False;
	    D_Double  (result) = dbllit (0);
	} else {
	    CreateData (result, NULL, NULL, T_Matrix, (int) val1, (int) val2);
	    a = D_Matrix (result);
	    if ((status = ZeroMatrix (a)))
		MatrixError ("zeros", a, NULL, status, F_True);
	}
    } else
	type_error = F_True;


    if (type_error == F_True)
	TypeError ("zeros", arg1, old, NULL, F_True);


    RecycleData (arg1);
    RecycleData (old);
    d_printf ("zeros ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}