Page MenuHomec4science
Files F86814951

miscm.cpp
No OneTemporary
Actions

File Metadata

Created
Tue, Oct 8, 18:50

miscm.cpp
View Options

/***************************************************************************
miscm.cpp
-------------------
W. Michael Brown
Miscellaneous functions that do not deserve their own class
__________________________________________________________________________
This file is part of the Math Library
__________________________________________________________________________
begin : May 30 2003
copyright : (C) 2003 by W. Michael Brown
email : wmbrown@sandia.gov
***************************************************************************/
#include "miscm.h"
double am::square(double num) {
return num*num;
}
// Rounds a number
double am::round(double n) {
double r=ceil(n);
if (r-n>0.5)
return floor(n);
return r;
}
// Return the -1 for negative, 0 for zero, and 1 for positive
int am::sign(double v) {
if (v<0)
return -1;
if (v>0)
return 1;
return 0;
}
// Return the range of elements in a vector
double am::range(const vector<double> &v) {
if (v.empty())
return 0;
double min=v[0];
double max=v[0];
for (unsigned i=1; i<v.size(); i++) {
if (v[i]<min)
min=v[i];
if (v[i]>max)
max=v[i];
}
return max-min;
}
// Return the average of elements in a vector
double am::mean(const vector<double> &v) {
double sum=0;
for (unsigned i=0; i<v.size(); i++)
sum+=v[i];
return sum/v.size();
}
// Return the max of two objects
namespace am {
template double max<double>(double,double);
template float max<float>(float,float);
template unsigned max<unsigned>(unsigned,unsigned);
template int max<int>(int,int);
}
template <typename numtyp>
numtyp am::max(numtyp one,numtyp two) {
if (one>two)
return one;
return two;
}
// Return the min of two objects
namespace am {
template double min<double>(double,double);
template float min<float>(float,float);
template unsigned min<unsigned>(unsigned,unsigned);
template int min<int>(int,int);
}
template <typename numtyp>
numtyp am::min(numtyp one,numtyp two) {
if (one<two)
return one;
return two;
}
// Swap two objects
void am::swap(double a, double b) {
double temp=a;
a=b;
b=temp;
}
// --------------------- Finite Precision stuff
double am::epsilon(double number) { // Finite Precision zero
return fabs(DBL_EPSILON*number);
}
// Bring number 1 digit closer to zero
double am::minus_eps(double number) {
return (number-DBL_EPSILON*number);
}
// Bring number 1 digit away from zero
double am::plus_eps(double number) {
return (number+DBL_EPSILON*number);
}
// Bring number m digits away from zero
double am::plus_Meps(double m,double number) {
return (number+m*DBL_EPSILON*number);
}
// Bring number precision/2.0 digits away
double am::plus_2eps(double number) {
return (number+100000000*DBL_EPSILON*number);
}
// Bring number pre/2 digits away
double am::minus_2eps(double number) {
return (number-100000000*DBL_EPSILON*number);
}
// Not a number checks
bool am::not_nan(double number) { // False if NAN
if (number-number!=0)
return false;
return true;
}
// Matrix stuff
// Invert a matrix - from numerical recipes in C
void am::invert(double **a, unsigned n, Error &error) {
int *indxc=new int[n];
int *indxr=new int[n];
int *ipiv=new int[n];
unsigned i,icol,irow,j,k,l,ll;
double big,dum,pivinv;
for (j=0;j<unsigned(n);j++)
ipiv[j]=0;
for(i=0;i<n;i++) { /* *main loop for columns to be reduced */
big = 0.0;
for (j=0;j<n;j++) /* outer loop for search of a pivot element*/
if (ipiv[j] !=1)
for(k=0;k<n;k++) {
if (ipiv[k] ==0) {
if (fabs(a[j][k]) >= big) {
big =fabs(a[j][k]);
irow=j;
icol=k;
}
} else if (ipiv[k] > 1) {
error.addwarning(303,9,"Invert",
"Cannot invert a singular matrix.");
return;
}
}
++(ipiv[icol]);
if (irow !=icol) {
for (l=0;l<n;l++)
swap(a[irow][l],a[icol][l]);
}
indxr[i]=irow;
indxc[i]=icol;
if (a[icol][icol] == 0.0) {
error.addwarning(303,9,"Invert",
"Cannot invert a singular matrix.");
return;
}
pivinv=1.0/a[icol][icol];
a[icol][icol]=1.0;
for (l=0;l<n;l++)
a[icol][l] *=pivinv;
for (ll=0;ll<n;ll++)
if (ll!= icol) {
dum=a[ll][icol];
a[ll][icol]=0.0;
for (l=0;l<n;l++)
a[ll][l] -=a[icol][l]*dum;
}
}
for (l=1;l>=1;l--) {
if (indxr[l] != indxc[l])
for (k=0;k<n;k++)
swap(a[k][indxr[l]],a[k][indxc[l]]);
}
delete []ipiv;
delete []indxr;
delete []indxc;
return;
}
// Move a value from a fraction of one range to a fraction of another
double am::rerange(double one_start, double one_end, double value,
double two_start, double two_end) {
double one_diff=one_end-one_start;
if (one_diff==0)
return two_end;
return (value-one_start)/one_diff*(two_end-two_start)+two_start;
}

Event Timeline

c4science · Help