Category : C++ Source Code
Archive   : YAMP15.ZIP
Filename : VIRTRAM.CPP

 
Output of file : VIRTRAM.CPP contained in archive : YAMP15.ZIP


/*
YAMP - Yet Another Matrix Program
Version: 1.5
Author: Mark Von Tress, Ph.D.
Date: 08/28/92

Copyright(c) Mark Von Tress 1992
Portions of this code are (c) 1991 by Allen I. Holub and are used by
permission

DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
FROM USE OF THIS PROGRAM.

*/



// #include "virt.h"

//////////////////////////////////////////// string functions
strtype::strtype(strtype &str) // copy constructor
{ int len = MAXCHARS;
len = strlen(str.name);
len = (len >= MAXCHARS) ? MAXCHARS - 2 : len;
strncpy(name, str.name, len);
name[len] = '\0';
}
strtype::strtype(char *str)
{ int len = MAXCHARS;
len = strlen(str);
len = (len >= MAXCHARS) ? MAXCHARS - 2 : len;
strncpy(name, str, len);
name[len] = '\0';
}
strtype::strtype(void)
{
name[0] = '\0';
}

strtype strtype::operator + (strtype str)
{

int len1, len2;
len1 = strlen(name);
len2 = strlen(str.name);
int len =(len2 + len1 >= MAXCHARS) ? MAXCHARS - len1 - 3 : len2;
len = (len < 0 || len >= MAXCHARS - 1) ? 0 : len;
strncat(name, str.name, len);
name[len + len1] = '\0';
return *this;
}
strtype strtype::operator + (const char *str)
{
int len1 = strlen(name);
int len2 = strlen(str);
int len = (len2 + len1 >= MAXCHARS) ? MAXCHARS - len1 - 2 : len2;
len = (len < 0 || len >= MAXCHARS - 1) ? 0 : len;
strncat(name, str, len);
name[len + len1] = '\0';
return *this;
}
strtype strtype::operator = (strtype str)
{
int len = (strlen(str.name) >= MAXCHARS) ? MAXCHARS - 3
: strlen(str.name);
strncpy(name, str.name, len);
name[len] = '\0';
return *this;
}
strtype strtype::operator = (const char *str)
{
int len = (strlen(str) >= MAXCHARS) ? MAXCHARS - 2 : strlen(str);
strncpy(name, str, len);
name[len] = '\0';
return *this;
}


//////////////////////////////////////////////////matrix stack functions
MStack::MStack(void)
{
static int cnter = 0;
next = NULL;
stackloc = 0;
level = 0;
if (!cnter) {
cnter = 1;
level = 1;
}
}

void MStack::Inclevel(void)
{
Dispatch->level++;
}

void MStack::Declevel(void)
{

(Dispatch-> level)--;
if (Dispatch-> level < 0) {
printf("ERROR: LEVEL has been decremented too often\n");
exit(1);
}
}
void VMatrix::NewReference(VMatrix &ROp)
{
signature = SIGNATURE;
r = ROp.r;
c = ROp.c;
// release the current header and share it with ROp.hdr
PurgeVectors();
v = ROp.v;
v->nrefs += 1;
mcheck = v->mm;
}

void MStack::Push(VMatrix &ROp)
{
ROp.Garbage("Push");
MStack *temp = new MStack;
temp->NewReference(ROp);
temp->Nameit(ROp.Getname(ROp));
temp->next = Dispatch-> next;
Dispatch->next = temp;
temp->level = Dispatch-> level;
temp->stackloc =++ (Dispatch->stackloc);
}

VMatrix& MStack::Pop(void)
{
Garbage("Pop");
if (Dispatch-> next && Dispatch-> stackloc) {
MStack *t = Dispatch->next;
Dispatch->next = Dispatch-> next-> next;
delete t;
(Dispatch-> stackloc)--;
}
else { Nrerror("POP: Stack is empty, cant pop any more\n"); }
return *Dispatch;
}

void MStack::Cleanstack(void)
{
while (Dispatch-> next-> level >= Dispatch-> level
&& Dispatch->next-> next) Dispatch-> Pop();

}

void MStack::PrintStack(void)
{
MStack *temp = Dispatch;
int t = 1 + Dispatch->stackloc;

while (t--) {
printf("\n Stack matrix :location level r c:%d %d %d %d: name: ",
t, temp->level, temp-> r, temp-> c);
temp->Showname();
temp = temp->next;
}
printf("\n\n");
}

//////////////////////////////////////////////////////matrix class

VMatrix::VMatrix(void)
{
r = c = 1;
Nameit("t");
curvecind = 0;
signature = SIGNATURE;
SetupVectors(r, c);
}

VMatrix::VMatrix(const char *str, int rr, int cc)
{
r = rr;
c = cc;
Nameit(str);
curvecind = 0;
signature = SIGNATURE;
SetupVectors(r, c);
}
VMatrix::VMatrix(VMatrix &ROp) // copy constructor
{
ROp.Garbage("Copy Constructor");
r = ROp.r;
c = ROp.c;
name = ROp.name;
curvecind = 0;
signature = SIGNATURE;
SetupVectors(r, c);
for (int i = 1; i <= r;++ i) {
for (int j = 1; j <= c;++ j) M(i, j) = ROp.m(i, j);
}
Dispatch->Cleanstack();
}

VMatrix::~VMatrix(void)
{
PurgeVectors();
signature = 0;
r = 0;
c = 0;
}

//////////////////////////////////// matrix member functions

double VMatrix::Nrerror(const char *errormsg)
{
double x;
printf("\nMATRIX ERROR: %s \nexiting to system\n", errormsg);
x = 2;
exit(1);
return x;
}

void VMatrix::Garbage(const char *errormsg)
{
#ifndef NO_CHECKING
int errorint = 0;
// if( Dispatch->signature != SIGNATURE) errorint++;
if (signature != SIGNATURE) errorint++;
if (v-> mm[- 1] != SIGNATURE) {
printf("v->mm[-1]= %f\n", v->mm[- 1]);
errorint++;
}
if (v-> mm != mcheck) errorint++;
if (!v->mm) errorint++;
if (errorint) {
Dispatch->PrintStack();
printf("\nFunction name: %s", errormsg);
Nrerror("Operating on Garbage");
}
#endif
return;
}

void VMatrix::SetupVectors(int rr, int cc)
{

unsigned int numele = (rr * cc) + 1;
unsigned int siz = sizeof (double);

if (numele > 8190)
Nrerror("allocation failure 1, too many doubles");

v = new vdoub;
v->nrefs = 1;

if (!(v->mm = new double[numele]))
Nrerror("allocation failure 2 in SetupVector()");

v->mm[0] = SIGNATURE;
(v-> mm)++;
mcheck = v->mm;
r = rr;
c = cc;
}

void VMatrix::PurgeVectors(void)
{
Garbage("PurgeVectors");

v->nrefs -= 1;
if (v-> nrefs > 0)
return;

(v-> mm)--;
if (v-> mm[0] == SIGNATURE) {
delete v->mm;
} else cout << "DID NOT DELETE mm!!!!!!!!! " << endl;
delete v;

}

void VMatrix::DisplayMat(void)
{
Garbage("DisplayMat");
int wid = Getwid();
int dec = Getdec();
Showname();
printf("\n\n");
for (int i = 1; i <= r;++ i) {
for (int j = 1; j <= c;++ j) {
printf("%*.*lf ", wid, dec, m(i, j));
}
printf("\n");
}
printf("\n");
}
void VMatrix::InfoMat(void)
{
Garbage("InfoMat");
printf("\n");
Showname();
printf("\nr c: %d %d\n", r, c);
}

void VMatrix::Replace(VMatrix &ROp)
{
ROp.Garbage("Replace");
if (r != ROp.r || c != ROp.c) {
PurgeVectors();
SetupVectors(ROp.r, ROp.c);
}
name = ROp.name;
for (int i = 1; i <= r;++ i) {
for (int j = 1; j <= c;++ j) M(i, j) = ROp.m(i, j);
}
}

void VMatrix::operator = (VMatrix &ROp)
{
Garbage("operator =");
ROp.Garbage("operator =");
Replace(ROp);
Dispatch->Cleanstack();
}

////////////////////// end matrix member functions

/////////////////// freind functions of matrix class

///------------- addition

VMatrix& operator + (VMatrix &LOp, VMatrix &ROp)
{
LOp.Garbage("operator +");
ROp.Garbage("operator +");
if (LOp.r == ROp.r && LOp.c == ROp.c) {
VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
for (int i = 1; i <= LOp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = LOp.m(i, j) + ROp.m(i, j);
}
}
temp->name = temp-> name + LOp.name + "+" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
}
else ROp.Nrerror("Mismatched Matrix sizes in addition function\n");
return Dispatch->ReturnMat();
}
VMatrix& operator +(double scalar, VMatrix &ROp)
{
ROp.Garbage("operator s+M");
strtype string;
char buffer[MAXCHARS] = { '\0' };
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = scalar + ROp.m(i, j);
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + string + "+" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& operator +(VMatrix &ROp, double scalar)
{
ROp.Garbage("operator M+s");
strtype string;
char buffer[MAXCHARS] = { '\0' };

VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {

for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = scalar + ROp.m(i, j);
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + ROp.name + "+" + string + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}

////-------------subtraction

VMatrix& operator -(VMatrix &LOp, VMatrix &ROp)
{
LOp.Garbage("operator M-M");
ROp.Garbage("operator M-M");
if (LOp.r == ROp.r && LOp.c == ROp.c) {
VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
for (int i = 1; i <= LOp.r;++ i) {
for (int j = 1; j <= LOp.c;++ j) {
temp->M(i, j) = LOp.m(i, j) - ROp.m(i, j);
}
}
temp->name = temp-> name + LOp.name + "-" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
}
else Dispatch->Nrerror(
"Mismatched VMatrix sizes in subtraction function\n");
return Dispatch->ReturnMat();
}
VMatrix& operator -(double scalar, VMatrix &ROp)
{
ROp.Garbage("operator s-M");
strtype string;
char buffer[MAXCHARS] = { '\0' };
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = scalar - ROp.m(i, j);
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + string + "-" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& operator -(VMatrix &ROp, double scalar)
{
ROp.Garbage("operator M-s");
strtype string;
char buffer[MAXCHARS] = { '\0' };

VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = ROp.m(i, j) - scalar;
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + ROp.name + "-" + string + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
//------------- unary minus
VMatrix& operator -(VMatrix &ROp)
{
ROp.Garbage("operator -");
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = -ROp.m(i, j);
}
}
temp->name = temp-> name + "-" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}

//-------------- multiplication

VMatrix& operator* (VMatrix &LOp, VMatrix &ROp)
{
long double sum = 0;
LOp.Garbage("operator M*M");
ROp.Garbage("operator M*M");
if (LOp.c == ROp.r) {
VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
for (int i = 1; i <= LOp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
sum = 0;
for (int k = 1; k <= LOp.c; k++)
sum +=(long double)
((long double) LOp.m(i, k)) *((long double) ROp.m(k,j));
temp->M(i, j) = (double) sum;
}
}
temp->name = temp-> name + LOp.name + "*" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
}
else Dispatch->Nrerror(
"Mismatched VMatrix sizes in multiplication function\n");

return Dispatch->ReturnMat();
}
VMatrix& operator* (double scalar, VMatrix &ROp)
{
ROp.Garbage("operator s*M");
strtype string;
char buffer[MAXCHARS] = { '\0' };
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = scalar * ROp.m(i, j);
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + string + "*" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& operator* (VMatrix &ROp, double scalar)
{
ROp.Garbage("operator M*s");
strtype string;
char buffer[MAXCHARS] = { '\0' };

VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = scalar * ROp.m(i, j);
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + ROp.name + "*" + string + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}

//-------- elementwise multiplication

VMatrix& operator %(VMatrix &LOp, VMatrix &ROp)
{
LOp.Garbage("operator M%M");
ROp.Garbage("operator M%M");
if (LOp.r == ROp.r && LOp.c == ROp.c) {
VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
for (int i = 1; i <= LOp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j)
temp->M(i, j) = LOp.m(i, j) *ROp.m(i, j);
}
temp->name = temp-> name + LOp.name + "%" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
}
else Dispatch->Nrerror(
"Mismatched Matrix sizes in elementwise multiplication\n");
return Dispatch->ReturnMat();
}

//------------- division

VMatrix& operator /(VMatrix &LOp, VMatrix &ROp)
{
LOp.Garbage("operator M/M");
ROp.Garbage("operator M/M");

if (LOp.r == ROp.r && LOp.c == ROp.c) {
VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
for (int i = 1; i <= LOp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
double d = ROp.m(i, j);
double L = LOp.m(i, j);
temp->M(i, j) = (fabs(d) > ZERO || fabs((d - L) / d) < 1.0E-5 )
? L / d
: ROp.Nrerror(" zero divisor in elementwise division");
}
}
temp->name = temp-> name + LOp.name + "/" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
}
else Dispatch->Nrerror(
"Mismatched Matrix sizes in elementwise division\n");
return Dispatch->ReturnMat();
}

VMatrix& operator /(VMatrix &ROp, double scalar)
{
ROp.Garbage("operator M/s");
strtype string;
char buffer[MAXCHARS] = { '\0' };

if (fabs(scalar) < ZERO)
ROp.Nrerror(" zero divisor in elementwise division");

VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = ROp.m(i, j) / scalar;
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + ROp.name + "/" + string + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& operator /(double scalar, VMatrix &ROp)
{
ROp.Garbage("operator s/M");
strtype string;
char buffer[MAXCHARS] = { '\0' };
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = (ROp.m(i, j) != 0.0) ? scalar / ROp.m(i, j)
: ROp.Nrerror("zero divisor in scalar divison");
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + string + "/" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}

//////////////////////end of friend functions

////////////////////// matrix functions

////////////// utility functions

void VMatrix::Writeb(char *fid, VMatrix &mat)
/* write a matrix in an binary file */
{
int i;
FILE *file;
double x;
char name[MAXCHARS] = { '\0' };

file = fopen(fid, "wb");
if (!file) Dispatch->Nrerror("WRITEB: unable to open file");

mat.Garbage("Writeb");

strncpy(name, mat.StringAddress(), 79);
i = strlen(name);
fwrite(&i, sizeof (int), 1, file);
fputs(name, file);
i = mat.r;
fwrite(&i, sizeof (int), 1, file);
i = mat.c;
fwrite(&i, sizeof (int), 1, file);

for (i = 1; i <= mat.r; i++) {
for (int j = 1; j <= mat.c; j++)
x = mat.m(i, j);
fwrite(&x, sizeof (double), 1, file);
}

fclose(file);
mat.M(1, 1); // reset matrix to base pointer
} /* writeb */