/*
$Id: utils.cc,v 1.3 1996/10/11 11:30:11 roitzsch Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include <stdarg.h>
#include "utils.h"
static double z_seconds();
MultiCompNode MCNode;
//-------------------------------------------------------------------------
void strToLower(char* s) { while(*s != '\0') { *s = tolower(*s); ++s; } }
void strToUpper(char* s) { while(*s != '\0') { *s = toupper(*s); ++s; } }
//-------------------------------------------------------------------------
Timer:: Timer()
{
t0 = z_seconds();
time(&abst0);
}
float Timer:: cpu(int print)
{
double seconds = z_seconds()-t0;
if (print)
{
cout << " " << Form("%1.2f", seconds) << " cpu sec.";
int minutes = int(seconds)/60;
float restSeconds = seconds - minutes*60;
cout << " (" << minutes << ":" << Form("%1.2f",restSeconds) << ")\n";
}
return seconds;
}
long Timer:: total(int print)
{
time_t seconds;
time(&seconds);
seconds = seconds - abst0;
if (print) cout << " " << seconds << " sec. (" << seconds/60 << ":"
<< seconds%60 << ")\n";
return seconds;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
Parser:: Parser(const char* file, const char* endToken0, const char commentFlag0)
: commentFlag(commentFlag0), fp(0)
{
fp = fopen(file, "r");
if (fp == 0)
{
cout << "\n*** Parser: " << file << " not found\n";
cout.flush(); abort();
}
endToken = new char[strlen(endToken0)+1];
strcpy(endToken, endToken0);
fileName = new char[strlen(file)+1];
strcpy(fileName, file);
}
//-------------------------------------------------------------------------
Parser:: ~Parser() { delete fileName; delete endToken; if (fp) fclose(fp); }
//-------------------------------------------------------------------------
void Parser:: readError(const char* function, char* word)
{
cout << "\n*** Parser." << function << ": read error on word '" << word
<< "' in file " << fileName << "\n";
cout.flush(); abort();
}
//-------------------------------------------------------------------------
int Parser:: skipRestOfLine()
{
int c;
while ((c=getc(fp)) != '\n') if (c==EOF) {status = eof; return status;}
status = ok;
return status;
}
//-------------------------------------------------------------------------
int Parser:: getString(char* word)
{
while (fscanf(fp,"%s",word) != EOF)
{
strToLower(word);
if (word[0] == commentFlag)
{
if (skipRestOfLine() == eof) { status = eof; return status; }
}
else
{
if (!strcmp(word,endToken)) status = end;
else status = ok;
return status;
}
}
status = eof;
return status;
}
//-------------------------------------------------------------------------
int Parser:: findWord(const char* target, int completeMatch)
{
char word[WORDLENGTH], s[WORDLENGTH];
strcpy(s, target);
strToLower(s);
if (completeMatch)
{
while (getString(word) != eof)
if (!strcmp(word,s)) { status = ok; return status; }
}
else
{
int length = strlen(target);
while (getString(word) != eof)
if (!strncmp(word,s,length)) { status = ok; return status; }
}
return status;
}
//-------------------------------------------------------------------------
int Parser:: findChar(const char target)
{
char ch;
int c;
ch = getChar();
c = ch;
while (c != EOF)
{
if (ch == target) { status = ok; return status; }
ch = getChar();
c = ch;
}
return status;
}
//-------------------------------------------------------------------------
int Parser:: getValue(char* c)
{
char word[WORDLENGTH];
if (getString(word) != ok) return status;
if (strlen(word) != 1) readError("getValue(char)", word);
if (sscanf(word,"%c",c) != 1) readError("getValue(char)", word);
status = ok;
return status;
}
int Parser:: getValue(int* i)
{
char word[WORDLENGTH];
if (getString(word) != ok) return status;
if (sscanf(word,"%d",i) != 1) readError("getValue(int)", word);
status = ok;
return status;
}
int Parser:: getValue(float* x)
{
char word[WORDLENGTH];
if (getString(word) != ok) return status;
if (sscanf(word,"%f",x) != 1) readError("getValue(float)", word);
status = ok;
return status;
}
int Parser:: getValue(double* x)
{
char word[WORDLENGTH];
if (getString(word) != ok) return status;
if (sscanf(word,"%lf",x) != 1) readError("getValue(double)", word);
status = ok;
return status;
}
//-------------------------------------------------------------------------
char Parser:: getChar()
{
int c;
if ((c=getc(fp)) == EOF) status = eof;
else status = ok;
return char(c);
}
int Parser:: getLine(char* line)
{
int count = 0;
int c;
while ((c=getc(fp)) != '\n')
{
if (c==EOF) { line[count++] = '\n'; status = eof; return status; }
line[count++] = char(c);
}
line[count++] = '\n';
status = ok;
return status;
}
//-------------------------------------------------------------------------
void Parser:: rewind() { ::rewind(fp); }
long Parser:: ftell () { return ::ftell(fp); }
int Parser:: fseek (long offset, int mode) { return ::fseek(fp,offset,mode); }
//-------------------------------------------------------------------------
BufferedParser:: BufferedParser(const char* file, const char* endToken0,
const char commentFlag0)
: Parser(file, endToken0, commentFlag0), buffer(0), pos(0)
{
if (::fseek(fp,0,endOfFile) != 0)
{
cout << "\n*** BufferedParser: fseek failed on " << fileName << "\n";
fclose(fp);
cout.flush(); abort();
}
bufferLength = ::ftell(fp);
::rewind(fp);
++bufferLength;
buffer = new char[bufferLength];
if (buffer == 0)
{
cout << "\n*** BufferedParser: not enough memory to store "
<< bufferLength << " bytes for " << fileName << "\n";
fclose(fp);
cout.flush(); abort();
}
int size = 1;
fread(buffer, size, bufferLength-1, fp);
buffer[bufferLength-1] = '\0';
pos = buffer;
fclose(fp);
fp = 0;
}
//-------------------------------------------------------------------------
BufferedParser:: ~BufferedParser() { delete buffer; }
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
int BufferedParser:: getValue(char* c)
{
char word[WORDLENGTH];
posBefore = pos;
if (getString(word) != ok) return status;
if (strlen(word) != 1)
{
status=end; pos=posBefore; return status;
}
if (sscanf(word,"%c",c) != 1)
{
status=end; pos=posBefore; return status;
}
status = ok;
return status;
}
int BufferedParser:: getSpecValue(char* c)
{
char word[WORDLENGTH];
posBefore = pos;
if (getString(word) != ok) {
pos=posBefore; return status;
}
if (strlen(word) != 1)
{
status=end; pos=posBefore; return status;
}
if (sscanf(word,"%c",c) != 1)
{
status=end; pos=posBefore; return status;
}
if ( (*c != 'c') && (*c != 's') && (*c != 'p')
&& (*c != 'C') && (*c != 'S') && (*c != 'P') )
{
status=end; pos=posBefore; return status;
}
status = ok;
return status;
}
int BufferedParser:: getValue(int* i)
{
char word[WORDLENGTH];
if (getString(word) != ok) return status;
if (sscanf(word,"%d",i) != 1) readError("getValue(int)", word);
status = ok;
return status;
}
int BufferedParser:: getValue(float* x)
{
char word[WORDLENGTH];
if (getString(word) != ok) return status;
if (sscanf(word,"%f",x) != 1) readError("getValue(float)", word);
status = ok;
return status;
}
int BufferedParser:: getValue(double* x)
{
char word[WORDLENGTH];
if (getString(word) != ok) return status;
if (sscanf(word,"%lf",x) != 1) readError("getValue(double)", word);
status = ok;
return status;
}
//-------------------------------------------------------------------------
int BufferedParser:: getString(char* word)
{
// while (sscanf(pos,"%s",word) == 1) // beware of EOF
// {
// while (isspace(*pos)) ++pos; // skip spaces
// pos += strlen(word);
//
// strToLower(word);
//
// if (word[0] == commentFlag)
// {
// if (skipRestOfLine() == eof) { status = eof; return status; }
// }
// else
// {
// if (!strcmp(word,endToken)) status = end;
// else status = ok;
// return status;
// }
// }
// status = eof;
// return status;
// This code snippet is a performance hack to replace the
// code above and gains a factor 60 speedup in some, if not
// many, situations.
while (1) {
char* w = word;
while (isspace(*pos)) ++pos; // skip leading whitespace
while (*pos && !isspace(*pos)) *w++ = *pos++; // copy string
*w = 0; // terminate copied string
// copied string empty? Must have been eof... (rather eos here ;-)
if (0 == w-word) { status = eof; return status; }
if (word[0] == commentFlag) { // a comment! skip it.
if (skipRestOfLine()==eof) { status = eof; return status; }
} else { // A token! Wow!
strToLower(word);
if (!strcmp(word,endToken)) status = end;
else status = ok;
return status;
}
}
// never get here...
assert(0);
return 0;
}
//-------------------------------------------------------------------------
int BufferedParser:: skipRestOfLine()
{
while (*pos != '\n')
{
if (*pos=='\0') { status = eof; return status; }
else ++pos;
}
++pos;
status = ok;
return status;
}
char BufferedParser:: getChar()
{
char c;
if ((c = *pos) == '\0') status = eof;
else {
++pos;
status = ok;
}
return c;
}
int BufferedParser:: getLine(char* line)
{
int count = 0;
while (*pos != '\n')
{
if (*pos == '\0') { line[++count] = '\n'; status = eof; return status; }
line[++count] = *pos;
++pos;
}
line[++count] = '\n';
++pos;
status = ok;
return status;
}
//-------------------------------------------------------------------------
void BufferedParser:: rewind() { pos = buffer; }
long BufferedParser:: ftell() { return long(pos-buffer); }
int BufferedParser:: fseek(long offset, int mode)
{
if (mode == beginning) pos = buffer + offset;
else if (mode == current) pos = pos + offset;
else if (mode == endOfFile) pos = buffer + bufferLength + offset;
else
{
cout << "\n*** BufferedParser: fseek called on unknown mode "
<< mode <<"\n";
cout.flush(); abort();
}
if (pos < buffer || pos > buffer+bufferLength)
{
cout << "\n*** BufferedParser: fseek requested position outside file "
<< "\n";
cout.flush(); abort();
}
return 1;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
// invert matrix A (which is destroyed in LUDecompose)
void invert(Matrix<Real>& A, Matrix<Real>& AInv)
{
int i,j;
const int dim = A.cHigh();
if (dim == 1)
{
AInv(1,1) = 1.0/A(1,1);
}
else if (dim == 2)
{
AInv(1,1) = A(2,2); AInv(1,2) = -A(1,2);
AInv(2,1) = -A(2,1); AInv(2,2) = A(1,1);
Real det = A(1,1)*A(2,2) - A(1,2)*A(2,1);
for (i=1; i<=dim; ++i)
for (j=1; j<=dim; ++j) AInv(i,j) /= det;
}
else
{
Vector<int> index(dim);
Vector<Real> col(dim);
Decompose(A, index);
for (j=1; j<=dim; ++j)
{
for (i=1; i<=dim; ++i) col[i] = 0.0;
col[j] = 1.0;
FBSubst(A, index, col);
for (i=1; i<=dim; ++i) AInv(i,j) = col[i];
}
}
};
//-------------------------------------------------------------------------
// --LU decomposition with partial pivoting
// (~ Crout's algorithm taken from Numerical Recipes)
void Decompose(Matrix<Real>& A, Vector<int>& index)
{
int i, imax=0, j, k;
Real big, dum, sum, temp;
Real d = 1.0;
const Real tiny = machMin(Real(0.0));
const int dim = index.high();
Vector<Real> vv(dim);
for (i=1; i<=dim; i++)
{
big=0.0;
for (j=1; j<=dim; j++)
{
if (Abs(temp=A(i,j)) > Abs(big)) big = temp;
}
if (big == 0.0) {
cout << "\n*** Decompose: Singular Matrix\n"; cout.flush(); abort();
}
vv[i] = 1.0/big;
}
for (j=1; j<=dim; j++)
{
for (i=1; i<j; i++)
{
sum = A(i,j);
for (k=1; k<i; k++) sum -= A(i,k)*A(k,j);
A(i,j) = sum;
}
big=0.0;
for (i=j; i<=dim; i++)
{
sum = A(i,j);
for (k=1; k<j; k++) sum -= A(i,k)*A(k,j);
A(i,j) = sum;
if (Abs(dum=vv[i]*sum) >= Abs(big))
{
big=dum;
imax=i;
}
}
if (j != imax)
{
for (k=1; k<=dim; k++)
{
dum = A(imax,k);
A(imax,k) = A(j,k);
A(j,k) = dum;
}
d = -d;
vv[imax]=vv[j];
}
index[j] = imax;
if (A(j,j) == 0.0) A(j,j) = tiny;
if (j != dim)
{
dum=1.0/A(j,j);
for (i=j+1; i<=dim; i++) A(i,j) *= dum;
}
}
};
//-------------------------------------------------------------------------
// LU-FB-Substitution
// solve lhs = A**-1 * rhs; A = L*U
// FB Substitution may be done in place, i.e. FBSubst(x,x)
void FBSubst(Matrix<Real>& A, Vector<int>& index, Vector<Real>& x)
{
int i, ii = 0, ip, j;
Real sum;
const int dim = index.high();
for (i=1; i<=dim; i++)
{
ip = index[i];
sum = x[ip];
x[ip] = x[i];
if (ii) for (j=ii; j<=i-1; j++) sum -= A(i,j) * x[j];
else if (Abs(sum)) ii = i;
x[i] = sum;
}
for (i=dim; i>=1; i--)
{
sum = x[i];
for (j=i+1; j<=dim; j++) sum -= A(i,j) * x[j];
x[i] = sum / A(i,i);
}
};
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
// invert matrix A (which is destroyed in LUDecompose)
void invert(Matrix<Complex>& A, Matrix<Complex>& AInv)
{
int i,j;
const int dim = A.cHigh();
if (dim == 1)
{
AInv(1,1) = 1.0/A(1,1);
}
else if (dim == 2)
{
AInv(1,1) = A(2,2); AInv(1,2) = -A(1,2);
AInv(2,1) = -A(2,1); AInv(2,2) = A(1,1);
Complex det = A(1,1)*A(2,2) - A(1,2)*A(2,1);
for (i=1; i<=dim; ++i)
for (j=1; j<=dim; ++j) AInv(i,j) /= det;
}
else
{
Vector<int> index(dim);
Vector<Complex> col(dim);
Decompose(A, index);
for (j=1; j<=dim; ++j)
{
for (i=1; i<=dim; ++i) col[i] = 0.0;
col[j] = 1.0;
FBSubst(A, index, col);
for (i=1; i<=dim; ++i) AInv(i,j) = col[i];
}
}
};
//-------------------------------------------------------------------------
// --LU decomposition with partial pivoting
// (~ Crout's algorithm taken from Numerical Recipes)
void Decompose(Matrix<Complex>& A, Vector<int>& index)
{
int i, imax=0, j, k;
Complex big, dum, sum, temp;
Complex d = 1.0;
const Real tiny = machMin(Real(0.0));
const int dim = index.high();
Vector<Complex> vv(dim);
for (i=1; i<=dim; i++)
{
big=0.0;
for (j=1; j<=dim; j++)
{
if (Abs(temp=A(i,j)) > Abs(big)) big = temp;
}
if (big == 0.0) {
cout << "\n*** Decompose: Singular Matrix\n"; cout.flush(); abort();
}
vv[i] = 1.0/big;
}
for (j=1; j<=dim; j++)
{
for (i=1; i<j; i++)
{
sum = A(i,j);
for (k=1; k<i; k++) sum -= A(i,k)*A(k,j);
A(i,j) = sum;
}
big=0.0;
for (i=j; i<=dim; i++)
{
sum = A(i,j);
for (k=1; k<j; k++) sum -= A(i,k)*A(k,j);
A(i,j) = sum;
if (Abs(dum=vv[i]*sum) >= Abs(big))
{
big=dum;
imax=i;
}
}
if (j != imax)
{
for (k=1; k<=dim; k++)
{
dum = A(imax,k);
A(imax,k) = A(j,k);
A(j,k) = dum;
}
d = -d;
vv[imax]=vv[j];
}
index[j] = imax;
if (A(j,j) == 0.0) A(j,j) = tiny;
if (j != dim)
{
dum=1.0/A(j,j);
for (i=j+1; i<=dim; i++) A(i,j) *= dum;
}
}
};
//-------------------------------------------------------------------------
// LU-FB-Substitution
// solve lhs = A**-1 * rhs; A = L*U
// FB Substitution may be done in place, i.e. FBSubst(x,x)
void FBSubst(Matrix<Complex>& A, Vector<int>& index, Vector<Complex>& x)
{
int i, ii = 0, ip, j;
Complex sum;
const int dim = index.high();
for (i=1; i<=dim; i++)
{
ip = index[i];
sum = x[ip];
x[ip] = x[i];
if (ii) for (j=ii; j<=i-1; j++) sum -= A(i,j) * x[j];
else if (Abs(sum)) ii = i;
x[i] = sum;
}
for (i=dim; i>=1; i--)
{
sum = x[i];
for (j=i+1; j<=dim; j++) sum -= A(i,j) * x[j];
x[i] = sum / A(i,i);
}
};
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
static const int stringSize = 2000;
static char String[stringSize];
char *Form(const char *format ...)
{
va_list va;
va_start(va, format);
if (strlen(format) > 100)
cout <<
"\n--- Warning --- Form: standard output string may be too small\n";
for (int i=0; i<stringSize; ++i) String[i] = '\0';
vsprintf(String, format, va);
va_end(va);
return String;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
#if GETHRVTIME
/*
* It uses the system call gethrtime(3C), which is accurate to
* nanoseconds.
*/
#include <sys/time.h>
static double z_seconds() {
return ( gethrvtime() / 1e9 );
}
#elif MACOS
#include <time.h>
static double z_seconds() {
double t = clock(), cps = CLOCKS_PER_SEC;
return t/cps;
}
#else
#include <sys/times.h>
#include <unistd.h>
#ifndef CLK_TCK
#define CLK_TCK 60
#endif
static double z_seconds()
{
struct tms rusage;
times(&rusage);
return (double)rusage.tms_utime / CLK_TCK;
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1