/*
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: arithmetic.c *
* *
* Description: This file contains the function definitions for the *
* arithmetic virtual machine instructions. *
************************************************************************/
# include <math.h>
# include "debug.h"
# include "coerce.h"
# include "execute.h"
# include "arithmetic.h"
/************************************************************************
* Function: add_op *
* *
* Description: Pops and adds the two descriptors on the top of the *
* stack and places the result on the stack. The *
* following types are legal for addition: *
* *
* double + double -> double (scalar addition) *
* double + matrix -> matrix (matrix scaling) *
* matrix + double -> matrix (matrix scaling) *
* matrix + matrix -> matrix (matrix addition) *
* *
* An attempt is first made to coerce both operands to *
* double values. *
************************************************************************/
int add_op ( )
{
Matrix a;
Matrix b;
Matrix c;
descriptor *left;
descriptor *right;
descriptor *result;
descriptor temp;
int type_error;
int status;
right = pop ( );
result = top ( );
temp = *result;
left = &temp;
left = CoerceData (deref (left), T_Double);
right = CoerceData (deref (right), T_Double);
status = 0;
type_error = F_False;
switch (D_Type (left)) {
case T_Double:
switch (D_Type (right)) {
case T_Double:
CreateData (result, left, right, T_Double);
*D_Double (result) = *D_Double (left) + *D_Double (right);
break;
case T_Matrix:
a = D_Matrix (right);
CreateData (result, left, right, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
if ((status = ScaleMatrix (b, a, 1.0, *D_Double (left))))
MatrixError ("+", a, NULL, status, F_False);
break;
default:
type_error = F_True;
break;
}
break;
case T_Matrix:
switch (D_Type (right)) {
case T_Double:
a = D_Matrix (left);
CreateData (result, left, right, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
if ((status = ScaleMatrix (b, a, 1.0, *D_Double (right))))
MatrixError ("+", a, NULL, status, F_False);
break;
case T_Matrix:
a = D_Matrix (left);
b = D_Matrix (right);
CreateData (result, left, right, T_Matrix, Mrows (a), Mcols (a));
c = D_Matrix (result);
if ((status = AddMatrices (c, a, b)))
MatrixError ("+", a, b, status, F_False);
break;
default:
type_error = F_True;
break;
}
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("+", left, right, NULL, F_False);
RecycleData (left);
RecycleData (right);
d_printf ("add ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: bkslv_op *
* *
* Description: Pops and left divides the two descriptors on the top of *
* the stack and places the result on the stack. The *
* following types are legal for left division: *
* *
* double \ double -> double (scalar left division) *
* double \ matrix -> matrix (matrix scaling) *
* matrix \ double -> matrix (matrix inversion + scaling) *
* matrix \ matrix -> matrix (matrix backsolve) *
* *
* An attempt is first made to coerce the operands to *
* double values. *
************************************************************************/
int bkslv_op ( )
{
Matrix a;
Matrix b;
Matrix c;
Matrix l;
Matrix p;
Matrix v;
double lvalue;
descriptor *left;
descriptor *right;
descriptor *result;
descriptor temp;
int type_error;
int status;
int info;
unsigned i;
unsigned j;
unsigned cols;
unsigned rows;
right = pop ( );
result = top ( );
temp = *result;
left = &temp;
left = CoerceData (deref (left), T_Double);
right = CoerceData (deref (right), T_Double);
status = 0;
type_error = F_False;
switch (D_Type (left)) {
case T_Double:
switch (D_Type (right)) {
case T_Double:
if ((lvalue = *D_Double (left)) == 0) {
MathException ("left division by zero");
status = 1;
} else {
CreateData (result, left, right, T_Double);
*D_Double (result) = *D_Double (right) / lvalue;
}
break;
case T_Matrix:
if ((lvalue = *D_Double (left)) == 0) {
MathException ("left division by zero");
status = 1;
} else {
a = D_Matrix (right);
CreateData (result, NULL, right, T_Matrix, Mrows(a), Mcols(a));
b = D_Matrix (result);
if ((status = ScaleMatrix (b, a, 1.0 / lvalue, 0.0)))
MatrixError ("\\", NULL, a, status, F_False);
}
break;
default:
type_error = F_True;
break;
}
break;
case T_Matrix:
switch (D_Type (right)) {
case T_Double:
a = D_Matrix (left);
CreateData (result, left, NULL, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
p = CreateColumnVector (Mrows (a));
if ((status = LUFactorMatrix (b, a, p, &info)))
MatrixError ("\\", a, NULL, status, F_False);
else if (info)
MathException ("singular matrix in left division"), status = 1;
else if ((status = InvertMatrix (b, b, p)))
MatrixError ("\\", a, NULL, status, F_False);
else if ((status = ScaleMatrix (b, b, *D_Double (right), 0.0)))
MatrixError ("\\", a, NULL, status, F_False);
DestroyMatrix (p);
break;
case T_Matrix:
a = D_Matrix (left);
b = D_Matrix (right);
CreateData (result, right, NULL, T_Matrix, Mrows (b), Mcols (b));
c = D_Matrix (result);
rows = Mrows (a);
cols = Mcols (a);
if (IsCompact (a) && rows == cols) {
l = D_Temp (left) ? a : CreateCopyMatrix (a);
CroutFactorMatrix (l);
if (Mcols (b) == 1) {
CopyMatrix (c, b);
if ((status = CroutBackSolveMatrix (l, c)))
MatrixError ("\\", a, b, status, F_False);
} else {
v = CreateColumnVector (Mrows (b));
for (j = 1; j <= Mcols (b); j ++) {
for (i = 1; i <= Mrows (b); i ++)
sdata (v, i, 1) = mdata (b, i, j);
if ((status = CroutBackSolveMatrix (l, v))) {
MatrixError ("\\", a, b, status, F_False);
j = Mcols (b);
break;
}
for (i = 1; i <= Mrows (b); i ++)
sdata (c, i, j) = mdata (v, i, 1);
}
DestroyMatrix (v);
}
} else {
p = CreateColumnVector (rows);
l = D_Writable (left) ? a : CreateFullMatrix (rows, cols);
if ((status = LUFactorMatrix (l, a, p, &info)))
MatrixError ("\\", a, b, status, F_False);
else if (info) {
MathException ("singular matrix in left division");
status = 1;
} else if (Mcols (b) != 1) {
v = CreateColumnVector (Mrows (b));
for (j = 1; j <= Mcols (b); j ++) {
for (i = 1; i <= Mrows (b); i ++)
sdata (v, i, 1) = mdata (b, i, j);
if ((status = LUBackSolveMatrix (v, l, v, p))) {
MatrixError ("\\", a, b, status, F_False);
j = Mcols (b);
break;
}
for (i = 1; i <= Mrows (b); i ++)
sdata (c, i, j) = mdata (v, i, 1);
}
DestroyMatrix (v);
} else if ((status = LUBackSolveMatrix (c, l, b, p)))
MatrixError ("\\", a, b, status, F_False);
DestroyMatrix (p);
}
if (a != l)
DestroyMatrix (l);
break;
default:
type_error = F_True;
break;
}
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("\\", left, right, NULL, F_False);
RecycleData (left);
RecycleData (right);
d_printf ("bkslv ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: div_op *
* *
* Description: Pops and right divides the two descriptors on the top *
* of the stack and places the result on the stack. The *
* following types are legal for right division: *
* *
* double / double -> double (scalar right division) *
* double / matrix -> matrix (matrix inversion + scaling) *
* matrix / double -> matrix (matrix scaling) *
* matrix / matrix -> matrix (matrix transposed backsolve) *
* *
* An attempt is first made to coerce the operands to *
* double values. *
************************************************************************/
int div_op ( )
{
Matrix a;
Matrix b;
Matrix c;
Matrix l;
Matrix p;
Matrix v;
Matrix at;
Matrix bt;
double rvalue;
descriptor *left;
descriptor *right;
descriptor *result;
descriptor temp;
int type_error;
int status;
int info;
unsigned i;
unsigned j;
right = pop ( );
result = top ( );
temp = *result;
left = &temp;
left = CoerceData (deref (left), T_Double);
right = CoerceData (deref (right), T_Double);
status = 0;
type_error = F_False;
switch (D_Type (left)) {
case T_Double:
switch (D_Type (right)) {
case T_Double:
if ((rvalue = *D_Double (right)) == 0) {
MathException ("right division by zero");
status = 1;
} else {
CreateData (result, left, right, T_Double);
*D_Double (result) = *D_Double (left) / rvalue;
}
break;
case T_Matrix:
a = D_Matrix (right);
CreateData (result, right, NULL, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
p = CreateColumnVector (Mrows (a));
if ((status = LUFactorMatrix (b, a, p, &info)))
MatrixError ("/", NULL, a, status, F_False);
else if (info)
MathException ("singular matrix in right division"), status = 1;
else if ((status = InvertMatrix (b, b, p)))
MatrixError ("/", NULL, a, status, F_False);
else if ((status = ScaleMatrix (b, b, *D_Double (left), 0.0)))
MatrixError ("/", NULL, a, status, F_False);
DestroyMatrix (p);
break;
default:
type_error = F_True;
break;
}
break;
case T_Matrix:
switch (D_Type (right)) {
case T_Double:
if ((rvalue = *D_Double (right)) == 0) {
MathException ("right division by zero");
status = 1;
} else {
a = D_Matrix (left);
CreateData (result, left, NULL, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
if ((status = ScaleMatrix (b, a, 1.0 / rvalue, 0.0)))
MatrixError ("/", a, NULL, status, F_False);
}
break;
case T_Matrix:
a = D_Matrix (left);
b = D_Matrix (right);
CreateData (result, left, NULL, T_Matrix, Mrows (a), Mcols (a));
c = D_Matrix (result);
at = CreateFullMatrix (Mcols (a), Mrows (a));
bt = CreateFullMatrix (Mcols (b), Mrows (b));
TransposeMatrix (at, a);
TransposeMatrix (bt, b);
p = CreateColumnVector (Mrows (bt));
l = bt;
if ((status = LUFactorMatrix (l, bt, p, &info)))
MatrixError ("/", a, b, status, F_False);
else if (info) {
MathException ("singular matrix in right division");
status = 1;
} else if (Mcols (at) != 1) {
v = CreateColumnVector (Mrows (at));
for (j = 1; j <= Mcols (at); j ++) {
for (i = 1; i <= Mrows (at); i ++)
sdata (v, i, 1) = mdata (at, i, j);
if ((status = LUBackSolveMatrix (v, l, v, p))) {
MatrixError ("/", a, b, status, F_False);
j = Mcols (at);
break;
}
for (i = 1; i <= Mrows (at); i ++)
sdata (c, j, i) = mdata (v, i, 1);
}
DestroyMatrix (v);
} else if ((status = LUBackSolveMatrix (at, l, at, p)))
MatrixError ("/", a, b, status, F_False);
else
TransposeMatrix (c, at);
DestroyMatrix (at);
DestroyMatrix (bt);
DestroyMatrix (p);
break;
default:
type_error = F_True;
break;
}
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("/", left, right, NULL, F_False);
RecycleData (left);
RecycleData (right);
d_printf ("div ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: mod_op *
* *
* Description: Pops and modulos the two descriptors on the top of the *
* stack and places the result on the stack. The *
* following types are legal for modulo: *
* *
* double % double -> double (scalar modulo) *
* double % matrix -> matrix (matrix element-wise modulo) *
* matrix % double -> matrix (matrix element-wise modulo) *
* matrix % matrix -> matrix (matrix element-wise modulo) *
* *
* An attempt is first made to coerce both operands to *
* double values. *
************************************************************************/
int mod_op ( )
{
Matrix a;
Matrix b;
Matrix c;
double lvalue;
double rvalue;
descriptor *left;
descriptor *right;
descriptor *result;
descriptor temp;
int type_error;
int status;
unsigned i;
unsigned j;
right = pop ( );
result = top ( );
temp = *result;
left = &temp;
left = CoerceData (deref (left), T_Double);
right = CoerceData (deref (right), T_Double);
status = 0;
type_error = F_False;
switch (D_Type (left)) {
case T_Double:
switch (D_Type (right)) {
case T_Double:
CreateData (result, left, right, T_Double);
if ((rvalue = *D_Double (right)) == 0) {
MathException ("remainder by zero");
status = 1;
} else
*D_Double (result) = fmod (*D_Double (left), rvalue);
break;
case T_Matrix:
a = D_Matrix (right);
CreateData (result, left, right, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
lvalue = *D_Double (left);
for (i = 1; i <= Mrows (a); i ++)
for (j = 1; j <= Mcols (a); j ++)
sdata (b, i, j) = fmod (lvalue, mdata (a, i, j));
break;
default:
type_error = F_True;
break;
}
break;
case T_Matrix:
switch (D_Type (right)) {
case T_Double:
a = D_Matrix (left);
CreateData (result, left, right, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
rvalue = *D_Double (right);
for (i = 1; i <= Mrows (a); i ++)
for (j = 1; j <= Mcols (a); j ++)
sdata (b, i, j) = fmod (mdata (a, i, j), rvalue);
break;
case T_Matrix:
a = D_Matrix (left);
b = D_Matrix (right);
CreateData (result, left, right, T_Matrix, Mrows (a), Mcols (a));
c = D_Matrix (result);
if (((status = ModMatrices (c, a, b))))
MatrixError ("%", a, b, status, F_False);
break;
default:
type_error = F_True;
break;
}
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("%", left, right, NULL, F_False);
RecycleData (left);
RecycleData (right);
d_printf ("mod ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: mul_op *
* *
* Description: Pops and multiplies the two descriptors on the top of *
* the stack and places the result on the stack. The *
* following types are legal for multiplication: *
* *
* double * double -> double (scalar multiplication) *
* double * matrix -> matrix (matrix scaling) *
* matrix * double -> matrix (matrix scaling) *
* matrix * matrix -> matrix (matrix multiplication) *
* *
* An attempt is first made to coerce both operands to *
* double values. If the resulting matrix consists of *
* only a single scalar then the result will be a scalar. *
************************************************************************/
int mul_op ( )
{
Matrix a;
Matrix b;
Matrix c;
descriptor *left;
descriptor *right;
descriptor *result;
descriptor temp;
int type_error;
int status;
right = pop ( );
result = top ( );
temp = *result;
left = &temp;
left = CoerceData (deref (left), T_Double);
right = CoerceData (deref (right), T_Double);
status = 0;
type_error = F_False;
switch (D_Type (left)) {
case T_Double:
switch (D_Type (right)) {
case T_Double:
CreateData (result, left, right, T_Double);
*D_Double (result) = *D_Double (left) * *D_Double (right);
break;
case T_Matrix:
a = D_Matrix (right);
CreateData (result, left, right, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
if ((status = ScaleMatrix (b, a, *D_Double (left), 0.0)))
MatrixError ("*", a, NULL, status, F_False);
break;
default:
type_error = F_True;
break;
}
break;
case T_Matrix:
switch (D_Type (right)) {
case T_Double:
a = D_Matrix (left);
CreateData (result, left, right, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
if ((status = ScaleMatrix (b, a, *D_Double (right), 0.0)))
MatrixError ("*", a, NULL, status, F_False);
break;
case T_Matrix:
a = D_Matrix (left);
b = D_Matrix (right);
if (Mrows (a) == 1 && Mcols (b) == 1) {
CreateData (result, NULL, NULL, T_Double);
if ((status = DotBProduct (D_Double (result), a, b)))
MatrixError ("*", a, b, status, F_False);
} else {
CreateData (result, NULL, NULL, T_Matrix, Mrows (a), Mcols (b));
c = D_Matrix (result);
if ((status = MultiplyMatrices (c, a, b)))
MatrixError ("*", a, b, status, F_False);
}
break;
default:
type_error = F_True;
break;
}
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("*", left, right, NULL, F_False);
RecycleData (left);
RecycleData (right);
d_printf ("mul ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: neg_op *
* *
* Description: Pops and negates the descriptor on the top of the stack *
* and places the result on the stack. The following *
* types are legal for negation: *
* *
* - scalar -> scalar (scalar arithmetic negation) *
* - matrix -> matrix (matrix element-wise negation) *
* *
* An attempt is first made to coerce the operand to a *
* double value. *
************************************************************************/
int neg_op ( )
{
Matrix a;
Matrix b;
descriptor *op;
descriptor *result;
descriptor temp;
int type_error;
int status;
result = top ( );
temp = *result;
op = &temp;
op = CoerceData (deref (op), T_Double);
status = 0;
type_error = F_False;
switch (D_Type (op)) {
case T_Double:
CreateData (result, op, NULL, T_Double);
*D_Double (result) = - *D_Double (op);
break;
case T_Matrix:
a = D_Matrix (op);
CreateData (result, op, NULL, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
if ((status = ScaleMatrix (b, a, -1.0, 0.0)))
MatrixError ("-", NULL, a, status, F_False);
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("-", NULL, op, NULL, F_False);
RecycleData (op);
d_printf ("neg ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: not_op *
* *
* Description: Pops and logically negates the descriptor on the top of *
* the stack and places the result on the stack. The *
* following types are legal for logical negation: *
* *
* ! scalar -> scalar (scalar logical negation) *
* ! matrix -> matrix (matrix element-wise negation) *
* *
* An attempt is first made to coerce the operand to a *
* double value. *
************************************************************************/
int not_op ( )
{
Matrix a;
Matrix b;
descriptor *op;
descriptor *result;
descriptor temp;
int type_error;
unsigned i;
unsigned j;
result = top ( );
temp = *result;
op = &temp;
op = CoerceData (deref (op), T_Double);
type_error = F_False;
switch (D_Type (op)) {
case T_Double:
CreateData (result, op, NULL, T_Double);
*D_Double (result) = ! *D_Double (op);
break;
case T_Matrix:
a = D_Matrix (op);
CreateData (result, op, NULL, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
for (i = 1; i <= Mrows (a); i ++)
for (j = 1; j <= Mcols (a); j ++)
sdata (b, i, j) = ! mdata (a, i, j);
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("!", NULL, op, NULL, F_False);
RecycleData (op);
d_printf ("not ans =\n");
d_PrintData (result);
return type_error == F_True;
}
/************************************************************************
* Function: plus_op *
* *
* Description: Pops and computes the positive of the descriptor on the *
* top of the stack and places the result on the stack. *
* The following types are legal: *
* *
* + double -> double (identity operation) *
* + matrix -> matrix (identity operation) *
* *
* An attempt is first made to coerce the operand to a *
* double value. *
************************************************************************/
int plus_op ( )
{
Matrix a;
Matrix b;
descriptor *op;
descriptor *result;
descriptor temp;
int type_error;
int status;
result = top ( );
temp = *result;
op = &temp;
op = CoerceData (deref (op), T_Double);
status = 0;
type_error = F_False;
switch (D_Type (op)) {
case T_Double:
CreateData (result, op, NULL, T_Double);
*D_Double (result) = *D_Double (op);
break;
case T_Matrix:
a = D_Matrix (op);
CreateData (result, op, NULL, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
if ((status = CopyMatrix (b, a)))
MatrixError ("+", NULL, a, status, F_False);
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("+", NULL, op, NULL, F_False);
RecycleData (op);
d_printf ("plus ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: pow_op *
* *
* Description: Pops and exponentiates the two descriptors on the top *
* of the stack and places the result on the stack. The *
* following types are legal for exponentiation: *
* *
* double ^ double -> double (scalar exponentation) *
* *
* An attempt is first made to coerce the operands to *
* double values. *
************************************************************************/
int pow_op ( )
{
double lvalue;
double rvalue;
descriptor *left;
descriptor *right;
descriptor *result;
descriptor temp;
int type_error;
int status;
right = pop ( );
result = top ( );
temp = *result;
left = &temp;
left = CoerceData (deref (left), T_Double);
right = CoerceData (deref (right), T_Double);
status = 0;
type_error = F_False;
switch (D_Type (left)) {
case T_Double:
switch (D_Type (right)) {
case T_Double:
lvalue = *D_Double (left);
rvalue = *D_Double (right);
if (lvalue < 0 && rvalue != (int) rvalue) {
MathException ("illegal base and exponent");
status = 1;
} else {
CreateData (result, left, right, T_Double);
*D_Double (result) = pow (lvalue, rvalue);
}
break;
default:
type_error = F_True;
break;
}
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("^", left, right, NULL, F_False);
RecycleData (left);
RecycleData (right);
d_printf ("pow ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: sub_op *
* *
* Description: Pops and subtracts the two descriptors on the top of *
* the stack and places the result on the stack. The *
* following types are legal for subtraction: *
* *
* double - double -> double (scalar subtraction) *
* double - matrix -> matrix (matrix scaling) *
* matrix - double -> matrix (matrix scaling) *
* matrix - matrix -> matrix (matrix subtraction) *
* *
* An attempt is first made to coerce both operands to *
* double values. *
************************************************************************/
int sub_op ( )
{
Matrix a;
Matrix b;
Matrix c;
descriptor *left;
descriptor *right;
descriptor *result;
descriptor temp;
int type_error;
int status;
right = pop ( );
result = top ( );
temp = *result;
left = &temp;
left = CoerceData (deref (left), T_Double);
right = CoerceData (deref (right), T_Double);
status = 0;
type_error = F_False;
switch (D_Type (left)) {
case T_Double:
switch (D_Type (right)) {
case T_Double:
CreateData (result, left, right, T_Double);
*D_Double (result) = *D_Double (left) - *D_Double (right);
break;
case T_Matrix:
a = D_Matrix (right);
CreateData (result, left, right, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
if ((status = ScaleMatrix (b, a, -1.0, *D_Double (left))))
MatrixError ("-", a, NULL, status, F_False);
break;
default:
type_error = F_True;
break;
}
break;
case T_Matrix:
switch (D_Type (right)) {
case T_Double:
a = D_Matrix (left);
CreateData (result, left, right, T_Matrix, Mrows (a), Mcols (a));
b = D_Matrix (result);
if ((status = ScaleMatrix (b, a, 1.0, -*D_Double (right))))
MatrixError ("-", a, NULL, status, F_False);
break;
case T_Matrix:
a = D_Matrix (left);
b = D_Matrix (right);
CreateData (result, left, right, T_Matrix, Mrows (a), Mcols (a));
c = D_Matrix (result);
if ((status = SubtractMatrices (c, a, b)))
MatrixError ("-", a, b, status, F_False);
break;
default:
type_error = F_True;
break;
}
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("-", left, right, NULL, F_False);
RecycleData (left);
RecycleData (right);
d_printf ("sub ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
/************************************************************************
* Function: trans_op *
* *
* Description: Pops and transposes the descriptor on the top of the *
* stack and places the result on the stack. The *
* following types are legal for transposition: *
* *
* double ' -> double (identity operation) *
* matrix ' -> matrix (matrix transposition) *
* *
* An attempt is first made to coerce the operand to a *
* double value. *
************************************************************************/
int trans_op ( )
{
Matrix a;
Matrix b;
descriptor *op;
descriptor *result;
descriptor temp;
int type_error;
int status;
result = top ( );
temp = *result;
op = &temp;
op = CoerceData (deref (op), T_Double);
status = 0;
type_error = F_False;
switch (D_Type (op)) {
case T_Double:
CreateData (result, op, NULL, T_Double);
*D_Double (result) = *D_Double (op);
break;
case T_Matrix:
a = D_Matrix (op);
CreateData (result, NULL, NULL, T_Matrix, Mcols (a), Mrows (a));
b = D_Matrix (result);
if ((status = TransposeMatrix (b, a)))
MatrixError ("'", a, NULL, status, F_False);
break;
default:
type_error = F_True;
break;
}
if (type_error == F_True)
TypeError ("'", op, NULL, NULL, F_False);
RecycleData (op);
d_printf ("trans ans =\n");
d_PrintData (result);
return type_error == F_True || status != 0;
}
syntax highlighted by Code2HTML, v. 0.9.1