/*
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;
}
syntax highlighted by Code2HTML, v. 0.9.1