Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88174932
miscm.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Oct 17, 05:43
Size
4 KB
Mime Type
text/x-c++
Expires
Sat, Oct 19, 05:43 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
21715521
Attached To
rLAMMPS lammps
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
Log In to Comment