Skip to content

GitLab

Commit 54f75cbe authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Fixes for hessian calculation 

+unit tests
parent fd6a6a92
Hide whitespace changes
Inline Side-by-side
Source/Solver/sparse.cpp
View file @ 54f75cbe
......@@ -403,15 +403,19 @@ namespace INMOST
assert(HR.isSorted());
output.Resize(HL.Size()+HR.Size()+JL.Size()*JR.Size());
INMOST_DATA_ENUM_TYPE i = 0, j = 0, k1 = 0, l1 = 0, k2 = 0, l2 = 0, q = 0, kk1 = 0, ll2 = 0, r;
while(kk1 < JL.Size() && JL.GetIndex(kk1) < JR.GetIndex(l1) ) kk1++;
k1 = kk1;
while(ll2 < JR.Size() && JL.GetIndex(k2) > JR.GetIndex(ll2) ) ll2++;
l2 = ll2;
entry candidate[4] = {stub_entry,stub_entry,stub_entry,stub_entry};
if( i < HL.Size() )
candidate[0] = make_entry(HL.GetIndex(i),b*HL.GetValue(i));
if( j < HR.Size() )
candidate[1] = make_entry(HR.GetIndex(j),c*HR.GetValue(j));
if( k1 < JL.Size() && l1 < JR.Size() )
candidate[2] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),0.5*a*JL.GetValue(k1)*JR.GetValue(l1));
candidate[2] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),(JL.GetIndex(k1)==JR.GetIndex(l1)?0.5:1)*a*JL.GetValue(k1)*JR.GetValue(l1));
if( k2 < JL.Size() && l2 < JR.Size() )
candidate[3] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),0.5*a*JL.GetValue(k2)*JR.GetValue(l2));
candidate[3] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),(JL.GetIndex(k2)==JR.GetIndex(l2)?0.5:1)*a*JL.GetValue(k2)*JR.GetValue(l2));
do
{
//pick smallest
......@@ -454,7 +458,7 @@ namespace INMOST
}
}
if( k1 < JL.Size() && l1 < JR.Size() )
candidate[2] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),(k1==l1?0.5:1)*a*JL.GetValue(k1)*JR.GetValue(l1));
candidate[2] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),(JL.GetIndex(k1)==JR.GetIndex(l1)?0.5:1)*a*JL.GetValue(k1)*JR.GetValue(l1));
else
candidate[2] = stub_entry;
}
......@@ -470,7 +474,7 @@ namespace INMOST
}
}
if( k2 < JL.Size() && l2 < JR.Size() )
candidate[3] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),(k2==l2?0.5:1)*a*JL.GetValue(k2)*JR.GetValue(l2));
candidate[3] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),(JL.GetIndex(k2)==JR.GetIndex(l2)?0.5:1)*a*JL.GetValue(k2)*JR.GetValue(l2));
else
candidate[3] = stub_entry;
}
......@@ -489,13 +493,17 @@ namespace INMOST
assert(H.isSorted());
output.Resize(H.Size()+JL.Size()*JR.Size());
INMOST_DATA_ENUM_TYPE i = 0, k1 = 0, l1 = 0, q = 0, k2 = 0, l2 = 0, r, ll2 = 0, kk1 = 0;
while(kk1 < JL.Size() && JL.GetIndex(kk1) < JR.GetIndex(l1) ) kk1++;
k1 = kk1;
while(ll2 < JR.Size() && JL.GetIndex(k2) > JR.GetIndex(ll2) ) ll2++;
l2 = ll2;
entry candidate[3] = {stub_entry,stub_entry,stub_entry};
if( i < H.Size() )
candidate[0] = make_entry(H.GetIndex(i),b*H.GetValue(i));
if( k1 < JL.Size() && l1 < JR.Size() )
candidate[1] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),0.5*a*JL.GetValue(k1)*JR.GetValue(l1));
candidate[1] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),(JL.GetIndex(k1)==JR.GetIndex(l1)?0.5:1)*a*JL.GetValue(k1)*JR.GetValue(l1));
if( k2 < JL.Size() && l2 < JR.Size() )
candidate[2] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),0.5*a*JL.GetValue(k2)*JR.GetValue(l2));
candidate[2] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),(JL.GetIndex(k2)==JR.GetIndex(l2)?0.5:1)*a*JL.GetValue(k2)*JR.GetValue(l2));
do
{
//pick smallest
......@@ -533,7 +541,7 @@ namespace INMOST
}
}
if( k1 < JL.Size() && l1 < JR.Size() )
candidate[1] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),(k1==l1?0.5:1)*a*JL.GetValue(k1)*JR.GetValue(l1));
candidate[1] = make_entry(make_index(JL.GetIndex(k1),JR.GetIndex(l1)),(JL.GetIndex(k1)==JR.GetIndex(l1)?0.5:1)*a*JL.GetValue(k1)*JR.GetValue(l1));
else
candidate[1] = stub_entry;
}
......@@ -549,7 +557,7 @@ namespace INMOST
}
}
if( k2 < JL.Size() && l2 < JR.Size() )
candidate[2] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),(k2==l2?0.5:1)*a*JL.GetValue(k2)*JR.GetValue(l2));
candidate[2] = make_entry(make_index(JL.GetIndex(k2),JR.GetIndex(l2)),(JL.GetIndex(k2)==JR.GetIndex(l2)?0.5:1)*a*JL.GetValue(k2)*JR.GetValue(l2));
else
candidate[2] = stub_entry;
}
......
Tests/CMakeLists.txt
View file @ 54f75cbe
if(USE_AUTODIFF)
add_subdirectory(autodiff_test000)
add_subdirectory(autodiff_test001)
add_subdirectory(autodiff_test002)
endif(USE_AUTODIFF)
if(USE_SOLVER)
......
Tests/autodiff_test000/main.cpp
View file @ 54f75cbe
......@@ -13,6 +13,7 @@ int main(int argc,char ** argv)
double _dx, _dy, _dxdx, _dydy, _dxdy;
unknown x(_x,0), y(_y,1);
hessian_variable f;
variable f2;
if( test == 0 ) //check derivative and hessian of sin(x*x+y*y)
......@@ -34,6 +35,7 @@ int main(int argc,char ** argv)
_dydy = 2*cos(_x*_x+_y*_y)-4*sin(_x*_x+_y*_y)*_y*_y;
_dxdy = -8*sin(_x*_x+_y*_y)*_x*_y;
f = sin(x*x+y*y);
f2 = sin(x*x + y*y);
}
else if( test == 1 )
{
......@@ -55,6 +57,7 @@ int main(int argc,char ** argv)
_dydy = 2*cos(_x*_x+_y*_y+_x*_y)-sin(_x*_x+_y*_y+_x*_y)*(_x+2*_y)*(_x+2*_y);
_dxdy = 2*cos(_x*_x+_y*_y+_x*_y)-2*sin(_x*_x+_y*_y+_x*_y)*(2*_x+_y)*(_x+2*_y);
f = sin(x*x+y*y+x*y);
f2 = sin(x*x + y*y + x*y);
}
else if( test == 2 )
{
......@@ -76,6 +79,7 @@ int main(int argc,char ** argv)
_dydy = -2*sin(_x*_x+_y*_y+_x*_y)-cos(_x*_x+_y*_y+_x*_y)*(_x+2*_y)*(_x+2*_y);
_dxdy = -2*sin(_x*_x+_y*_y+_x*_y)-2*cos(_x*_x+_y*_y+_x*_y)*(2*_x+_y)*(_x+2*_y);
f = cos(x*x+y*y+x*y);
f2 = cos(x*x + y*y + x*y);
}
else if( test == 3 )
{
......@@ -93,6 +97,7 @@ int main(int argc,char ** argv)
_dxdy = 2*cos(_x);
_dydy = 0;
f = sin(x)*y;
f2 = sin(x)*y;
}
else if( test == 4 )
{
......@@ -110,6 +115,7 @@ int main(int argc,char ** argv)
_dxdy = -2*sin(_x);
_dydy = 0;
f = cos(x)*y;
f2 = cos(x)*y;
}
else if( test == 5 )
{
......@@ -127,6 +133,7 @@ int main(int argc,char ** argv)
_dxdy = -6*_x*_y/(4*pow(_x*_x+_x*_y+_y*_y,1.5));
_dydy = 3*_x*_x/(4*pow(_x*_x+_x*_y+_y*_y,1.5));
f = sqrt(x*x+y*y+x*y);
f2 = sqrt(x*x + y*y + x*y);
}
else if( test == 6 )
{
......@@ -136,6 +143,7 @@ int main(int argc,char ** argv)
_dxdy = 8;
_dydy = 6;
f = 2*x*x+3*y*y+4*x*y;
f2 = 2 * x*x + 3 * y*y + 4 * x*y;
}
else if( test == 7 )
{
......@@ -145,6 +153,7 @@ int main(int argc,char ** argv)
_dxdy = -32*(_x-0.5)*(_y-0.5)*sin(2*((_x-0.5)*(_x-0.5)+(_y-0.5)*(_y-0.5)));
_dydy = 4*cos(2*((_x - 0.5)*(_x - 0.5) + (_y - 0.5)*(_y - 0.5))) - 16*(_y - 0.5)*(_y - 0.5)*sin(2.0*((_x - 0.5)*(_x - 0.5) + (_y - 0.5)*(_y - 0.5)));
f = sin(2*((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)));
f2 = sin(2 * ((x - 0.5)*(x - 0.5) + (y - 0.5)*(y - 0.5)));
}
else if( test == 8 )
{
......@@ -154,6 +163,7 @@ int main(int argc,char ** argv)
_dxdy = -(-32*(_x-0.5)*(_y-0.5)*sin(2*((_x-0.5)*(_x-0.5)+(_y-0.5)*(_y-0.5))));
_dydy = -(4*cos(2*((_x - 0.5)*(_x - 0.5) + (_y - 0.5)*(_y - 0.5))) - 16*(_y - 0.5)*(_y - 0.5)*sin(2.0*((_x - 0.5)*(_x - 0.5) + (_y - 0.5)*(_y - 0.5))));
f = 1.0-sin(2*((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)));
f2 = 1.0 - sin(2 * ((x - 0.5)*(x - 0.5) + (y - 0.5)*(y - 0.5)));
}
else if( test == 9 )
{
......@@ -163,6 +173,7 @@ int main(int argc,char ** argv)
_dxdy = 32*(_x-0.5)*(_y-0.5)*cos(2*((_x-0.5)*(_x-0.5)+(_y-0.5)*(_y-0.5)));
_dydy = 4*sin(2*((_x - 0.5)*(_x - 0.5) + (_y - 0.5)*(_y - 0.5))) + 16*(_y - 0.5)*(_y - 0.5)*cos(2.0*((_x - 0.5)*(_x - 0.5) + (_y - 0.5)*(_y - 0.5)));
f = 1.0-cos(2*((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)));
f2 = 1.0 - cos(2 * ((x - 0.5)*(x - 0.5) + (y - 0.5)*(y - 0.5)));
}
else if( test == 10 )
{
......@@ -172,6 +183,7 @@ int main(int argc,char ** argv)
_dxdy = 8*(_y-0.5)*cos(2*((_x-0.5)*(_x-0.5)+(_y-0.5)*(_y-0.5)))-4*(_x-0.5)*_x*sin(2*((_x-0.5)*(_x-0.5)+(_y-0.5)*(_y-0.5)));
_dydy = 4*_x*cos(2*((_x-0.5)*(_x-0.5)+(_y-0.5)*(_y-0.5))) - 16*_x*(_y-0.5)*(_y-0.5)*sin(2*((_x-0.5)*(_x-0.5)+(_y-0.5)*(_y-0.5)));
f = sin(2*((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)))*x;
f2 = sin(2 * ((x - 0.5)*(x - 0.5) + (y - 0.5)*(y - 0.5)))*x;
}
else if( test == 11 )
{
......@@ -189,6 +201,7 @@ int main(int argc,char ** argv)
_dxdy = 2*cos(_x*_x+_y) - 4*_x*_x*sin(_x*_x+_y);
_dydy = -_x*sin(_x*_x+_y);
f = sin(x*x+y)*x;
f2 = sin(x*x + y)*x;
}
else if( test == 12 )
{
......@@ -206,6 +219,7 @@ int main(int argc,char ** argv)
_dxdy = 4*_y*(cos(_x*_x+_y*_y) - 2*_x*_x*sin(_x*_x+_y*_y));
_dydy = 2*_x*(cos(_x*_x+_y*_y)-2*_y*_y*sin(_x*_x+_y*_y));
f = sin(x*x+y*y)*x;
f2 = sin(x*x + y*y)*x;
}
double dx = f.GetRow()[0];
double dy = f.GetRow()[1];
......@@ -214,6 +228,7 @@ int main(int argc,char ** argv)
double dydy = f.GetHessianRow()[Sparse::HessianRow::make_index(1,1)];
bool error = false;
std::cout << "For GetHessian:" << std::endl;
std::cout << std::setw(10) << "derivative " << std::setw(10) << "original " << std::setw(10) << "computed" << std::endl;
std::cout << std::setw(10) << "dx " << std::setw(10) << _dx << std::setw(10) << dx << std::endl;
std::cout << std::setw(10) << "dy " << std::setw(10) << _dy << std::setw(10) << dy << std::endl;
......@@ -226,7 +241,15 @@ int main(int argc,char ** argv)
if( std::abs(dxdy-_dxdy) > 1.0e-9 ) error = true, std::cout << "Error in dxdy: " << std::abs(dxdy-_dxdy) << " original " << _dxdy << " computed " << dxdy << std::endl;
if( std::abs(dydy-_dydy) > 1.0e-9 ) error = true, std::cout << "Error in dydy: " << std::abs(dydy-_dydy) << " original " << _dydy << " computed " << dydy << std::endl;
if( error ) return -1;
dx = f2.GetRow()[0];
dy = f2.GetRow()[1];
std::cout << "For GetJacobian:" << std::endl;
std::cout << std::setw(10) << "derivative " << std::setw(10) << "original " << std::setw(10) << "computed" << std::endl;
std::cout << std::setw(10) << "dx " << std::setw(10) << _dx << std::setw(10) << dx << std::endl;
std::cout << std::setw(10) << "dy " << std::setw(10) << _dy << std::setw(10) << dy << std::endl;
if( error ) return -1;
return 0;
}
Tests/autodiff_test001/main.cpp
View file @ 54f75cbe
......@@ -14,6 +14,7 @@ int main(int argc,char ** argv)
double dx, dy, dz, dxdx, dydy, dzdz, dxdy, dxdz, dydz;
unknown vx(x,0), vy(y,1), vz(z,2);
hessian_variable f;
variable f2;
if( test == 0 ) //check derivative and hessian of sin(x*x+y*y)
......@@ -33,6 +34,7 @@ int main(int argc,char ** argv)
f = sin(8*((vx-0.5)*(vx-0.5)+(vy-0.5)*(vy-0.5)+(vz-0.5)*(vz-0.5)))*(2*vx-1);
f2 = sin(8 * ((vx - 0.5)*(vx - 0.5) + (vy - 0.5)*(vy - 0.5) + (vz - 0.5)*(vz - 0.5)))*(2 * vx - 1);
}
//mixed derivative computed twice: dxdy and dydx
dxdy *= 2;
......@@ -51,6 +53,7 @@ int main(int argc,char ** argv)
bool error = false;
std::cout << "For GetHessian:" << std::endl;
std::cout << std::setw(10) << "derivative " << std::setw(10) << "original " << std::setw(10) << "computed" << std::endl;
std::cout << std::setw(10) << "dx " << std::setw(10) << dx << std::setw(10) << vdx << std::endl;
std::cout << std::setw(10) << "dy " << std::setw(10) << dy << std::setw(10) << vdy << std::endl;
......@@ -70,7 +73,22 @@ int main(int argc,char ** argv)
if( std::abs(dydy-vdydy) > 1.0e-9 ) error = true, std::cout << "Error in dydy: " << std::abs(dydy-vdydy) << " original " << dydy << " computed " << vdydy << std::endl;
if( std::abs(dydz-vdydz) > 1.0e-9 ) error = true, std::cout << "Error in dydz: " << std::abs(dydz-vdydz) << " original " << dydz << " computed " << vdydz << std::endl;
if( std::abs(dzdz-vdzdz) > 1.0e-9 ) error = true, std::cout << "Error in dzdz: " << std::abs(dzdz-vdzdz) << " original " << dzdz << " computed " << vdzdz << std::endl;
if( error ) return -1;
vdx = f2.GetRow()[0];
vdy = f2.GetRow()[1];
vdz = f2.GetRow()[2];
std::cout << "For GetJacobian:" << std::endl;
std::cout << std::setw(10) << "derivative " << std::setw(10) << "original " << std::setw(10) << "computed" << std::endl;
std::cout << std::setw(10) << "dx " << std::setw(10) << dx << std::setw(10) << vdx << std::endl;
std::cout << std::setw(10) << "dy " << std::setw(10) << dy << std::setw(10) << vdy << std::endl;
std::cout << std::setw(10) << "dz " << std::setw(10) << dz << std::setw(10) << vdz << std::endl;
if (std::abs(dx - vdx) > 1.0e-9) error = true, std::cout << "Error in dx: " << std::abs(dx - vdx) << " original " << dx << " computed " << vdx << std::endl;
if (std::abs(dy - vdy) > 1.0e-9) error = true, std::cout << "Error in dy: " << std::abs(dy - vdy) << " original " << dy << " computed " << vdy << std::endl;
if (std::abs(dz - vdz) > 1.0e-9) error = true, std::cout << "Error in dz: " << std::abs(dz - vdz) << " original " << dz << " computed " << vdz << std::endl;
if (error) return -1;
return 0;
}
Tests/autodiff_test002/CMakeLists.txt 0 → 100644
View file @ 54f75cbe
project(autodiff_test002)
set(SOURCE main.cpp)
add_executable(autodiff_test002 ${SOURCE})
target_link_libraries(autodiff_test002 inmost)
if(USE_MPI)
message("linking autodiff_test002 with MPI")
target_link_libraries(autodiff_test002 ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(autodiff_test002 PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
add_test(NAME autodiff_test002_hessian_0 COMMAND $<TARGET_FILE:autodiff_test002> 0)
add_test(NAME autodiff_test002_hessian_1 COMMAND $<TARGET_FILE:autodiff_test002> 1)
add_test(NAME autodiff_test002_hessian_2 COMMAND $<TARGET_FILE:autodiff_test002> 2)
add_test(NAME autodiff_test002_hessian_3 COMMAND $<TARGET_FILE:autodiff_test002> 3)
add_test(NAME autodiff_test002_hessian_4 COMMAND $<TARGET_FILE:autodiff_test002> 4)
add_test(NAME autodiff_test002_hessian_5 COMMAND $<TARGET_FILE:autodiff_test002> 5)
add_test(NAME autodiff_test002_hessian_6 COMMAND $<TARGET_FILE:autodiff_test002> 6)
add_test(NAME autodiff_test002_hessian_7 COMMAND $<TARGET_FILE:autodiff_test002> 7)
Tests/autodiff_test002/main.cpp 0 → 100644
View file @ 54f75cbe
#include <cstdio>
#include <cmath>
#include "inmost.h"
using namespace INMOST;
int main(int argc,char ** argv)
{
int test = 0;
if (argc > 1) test = atoi(argv[1]);
const double pi = 3.1415926535897932384626433832795;
double x = 0.51;
double y = 0.25;
double z = 0.15;
double t = 0.1;
double dx, dy, dz, dt, dxdx, dydy, dzdz, dtdt, dxdy, dxdz, dydz, dxdt, dydt, dzdt;
unknown vx(x,0), vy(y,1), vz(z,2), vt(t,3);
hessian_variable f;
variable f2;
if( test == 0 )
{
dx = exp(-t)*(sin(2*pi*x*y)*sin(2*pi*z)*sin(2*pi*x*z) + 2 * pi*x*y*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x*sin(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 3 * x*x * y*y * z);
dy = exp(-t)*(2 * pi*x *x * cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * x *x*x * y*z);
dz = exp(-t)*(2 * pi*x*sin(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x *x * sin(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + x *x*x * y*y);
dt = -exp(-t)*(x*sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + x *x*x * y *y * z);
dxdx = exp(-t)*(-4 * pi *pi * x*sin(2 * pi*x*y)*z *z * sin(2 * pi*z)*sin(2 * pi*x*z) - 4 * pi *pi * x*y *y * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*
y*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*sin(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 8 * pi *pi * x*y*cos(2 * pi*x*y)*z*
sin(2 * pi*z)*cos(2 * pi*x*z) + 6 * x*y *y * z);
dxdy = exp(-t)*(-4 * pi *pi * x *x * y*sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*x*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x *x *
cos(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 6 * x *x * y*z);
dxdz = exp(-t)*(-4 * pi *pi * x *x * sin(2 * pi*x*y)*z*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*sin(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x*y*
cos(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*x*sin(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + 4 * pi *pi * x *x * y*cos(2 * pi*x*y)*sin(2 * pi*z)*
cos(2 * pi*x*z) + 4 * pi *pi * x*sin(2 * pi*x*y)*z*cos(2 * pi*z)*cos(2 * pi*x*z) + 3 * x *x * y *y);
dydy = exp(-t)*(2 * x *x*x * z - 4 * pi *pi * x *x*x * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z));
dydz = exp(-t)*(4 * pi *pi * x *x * cos(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x *x*x * cos(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + 2 * x *x*x * y);
dzdz = exp(-t)*(-4 * pi *pi * x *x*x * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) - 4 * pi *pi * x*sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 8 * pi *pi * x *x *
sin(2 * pi*x*y)*cos(2 * pi*z)*cos(2 * pi*x*z));
dxdt = -exp(-t)*(sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x*y*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x*
sin(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 3 * x *x * y *y * z);
dydt = -exp(-t)*(2 * pi*x *x * cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * x *x*x * y*z);
dzdt = -exp(-t)*(2 * pi*x*sin(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x *x * sin(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + x *x*x * y *y);
dtdt = exp(-t)*(x*sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + x *x*x * y *y * z);
f = (vx*vx*vx*vy*vy*vz + vx*sin(2 * pi*vx*vz)*sin(2 * pi*vx*vy)*sin(2 * pi*vz))*exp(-vt);
f2 = (vx*vx*vx*vy*vy*vz + vx*sin(2 * pi*vx*vz)*sin(2 * pi*vx*vy)*sin(2 * pi*vz))*exp(-vt);
}
else if (test == 1)
{
dx = t*(sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x*y*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x*sin(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 3 * x*x * y*y * z);
dy = t*(2 * pi*x *x * cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * x *x*x * y*z);
dz = t*(2 * pi*x*sin(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x *x * sin(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + x *x*x * y*y);
dt = (x*sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + x *x*x * y *y * z);
dxdx = t*(-4 * pi *pi * x*sin(2 * pi*x*y)*z *z * sin(2 * pi*z)*sin(2 * pi*x*z) - 4 * pi *pi * x*y *y * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*
y*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*sin(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 8 * pi *pi * x*y*cos(2 * pi*x*y)*z*
sin(2 * pi*z)*cos(2 * pi*x*z) + 6 * x*y *y * z);
dxdy = t*(-4 * pi *pi * x *x * y*sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*x*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x *x *
cos(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 6 * x *x * y*z);
dxdz = t*(-4 * pi *pi * x *x * sin(2 * pi*x*y)*z*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*sin(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x*y*
cos(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*x*sin(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + 4 * pi *pi * x *x * y*cos(2 * pi*x*y)*sin(2 * pi*z)*
cos(2 * pi*x*z) + 4 * pi *pi * x*sin(2 * pi*x*y)*z*cos(2 * pi*z)*cos(2 * pi*x*z) + 3 * x *x * y *y);
dydy = t*(2 * x *x*x * z - 4 * pi *pi * x *x*x * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z));
dydz = t*(4 * pi *pi * x *x * cos(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x *x*x * cos(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + 2 * x *x*x * y);
dzdz = t*(-4 * pi *pi * x *x*x * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) - 4 * pi *pi * x*sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 8 * pi *pi * x *x *
sin(2 * pi*x*y)*cos(2 * pi*z)*cos(2 * pi*x*z));
dxdt = (sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x*y*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x*
sin(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 3 * x *x * y *y * z);
dydt = (2 * pi*x *x * cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * x *x*x * y*z);
dzdt = (2 * pi*x*sin(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x *x * sin(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + x *x*x * y *y);
dtdt = 0;
f = (vx*vx*vx*vy*vy*vz + vx*sin(2 * pi*vx*vz)*sin(2 * pi*vx*vy)*sin(2 * pi*vz))*vt;
f2 = (vx*vx*vx*vy*vy*vz + vx*sin(2 * pi*vx*vz)*sin(2 * pi*vx*vy)*sin(2 * pi*vz))*vt;
}
else if (test == 2)
{
dx = (sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x*y*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x*sin(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 3 * x*x * y*y * z);
dy = (2 * pi*x *x * cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * x *x*x * y*z);
dz = (2 * pi*x*sin(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x *x * sin(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + x *x*x * y*y);
dt = 0;
dxdx = (-4 * pi *pi * x*sin(2 * pi*x*y)*z *z * sin(2 * pi*z)*sin(2 * pi*x*z) - 4 * pi *pi * x*y *y * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*
y*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*sin(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 8 * pi *pi * x*y*cos(2 * pi*x*y)*z*
sin(2 * pi*z)*cos(2 * pi*x*z) + 6 * x*y *y * z);
dxdy = (-4 * pi *pi * x *x * y*sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*x*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x *x *
cos(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 6 * x *x * y*z);
dxdz = (-4 * pi *pi * x *x * sin(2 * pi*x*y)*z*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*sin(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x*y*
cos(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*x*sin(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + 4 * pi *pi * x *x * y*cos(2 * pi*x*y)*sin(2 * pi*z)*
cos(2 * pi*x*z) + 4 * pi *pi * x*sin(2 * pi*x*y)*z*cos(2 * pi*z)*cos(2 * pi*x*z) + 3 * x *x * y *y);
dydy = (2 * x *x*x * z - 4 * pi *pi * x *x*x * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z));
dydz = (4 * pi *pi * x *x * cos(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x *x*x * cos(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + 2 * x *x*x * y);
dzdz = (-4 * pi *pi * x *x*x * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) - 4 * pi *pi * x*sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 8 * pi *pi * x *x *
sin(2 * pi*x*y)*cos(2 * pi*z)*cos(2 * pi*x*z));
dxdt = 0;
dydt = 0;
dzdt = 0;
dtdt = 0;
f = (vx*vx*vx*vy*vy*vz + vx*sin(2 * pi*vx*vz)*sin(2 * pi*vx*vy)*sin(2 * pi*vz));
f2 = (vx*vx*vx*vy*vy*vz + vx*sin(2 * pi*vx*vz)*sin(2 * pi*vx*vy)*sin(2 * pi*vz));
}
else if (test == 3)
{
dx = (sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x*y*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x*sin(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) );
dy = (2 * pi*x *x * cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z));
dz = (2 * pi*x*sin(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*x *x * sin(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z));
dt = 0;
dxdx = (-4 * pi *pi * x*sin(2 * pi*x*y)*z *z * sin(2 * pi*z)*sin(2 * pi*x*z) - 4 * pi *pi * x*y *y * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*
y*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*sin(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z) + 8 * pi *pi * x*y*cos(2 * pi*x*y)*z*
sin(2 * pi*z)*cos(2 * pi*x*z));
dxdy = (-4 * pi *pi * x *x * y*sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*x*cos(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x *x *
cos(2 * pi*x*y)*z*sin(2 * pi*z)*cos(2 * pi*x*z));
dxdz = (-4 * pi *pi * x *x * sin(2 * pi*x*y)*z*sin(2 * pi*z)*sin(2 * pi*x*z) + 2 * pi*sin(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x*y*
cos(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi*x*sin(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z) + 4 * pi *pi * x *x * y*cos(2 * pi*x*y)*sin(2 * pi*z)*
cos(2 * pi*x*z) + 4 * pi *pi * x*sin(2 * pi*x*y)*z*cos(2 * pi*z)*cos(2 * pi*x*z));
dydy = ( - 4 * pi *pi * x *x*x * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z));
dydz = (4 * pi *pi * x *x * cos(2 * pi*x*y)*cos(2 * pi*z)*sin(2 * pi*x*z) + 4 * pi *pi * x *x*x * cos(2 * pi*x*y)*sin(2 * pi*z)*cos(2 * pi*x*z));
dzdz = (-4 * pi *pi * x *x*x * sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) - 4 * pi *pi * x*sin(2 * pi*x*y)*sin(2 * pi*z)*sin(2 * pi*x*z) + 8 * pi *pi * x *x *
sin(2 * pi*x*y)*cos(2 * pi*z)*cos(2 * pi*x*z));
dxdt = 0;
dydt = 0;
dzdt = 0;
dtdt = 0;
f = (vx*sin(2 * pi*vx*vz)*sin(2 * pi*vx*vy)*sin(2 * pi*vz));
f2 = (vx*sin(2 * pi*vx*vz)*sin(2 * pi*vx*vy)*sin(2 * pi*vz));
}
else if (test == 4)
{
dx = sin(2 * pi*z);
dy = 0;
dz = 2 * pi*x*cos(2 * pi*z);
dt = 0;
dxdx = 0;
dxdy = 0;
dxdz = 2 * pi*cos(2 * pi*z);
dydy = 0;
dydz = 0;
dzdz = -4 * pi *pi * x*sin(2 * pi*z);
dxdt = 0;
dydt = 0;
dzdt = 0;
dtdt = 0;
f = (vx*sin(2 * pi*vz));
f2 = (vx*sin(2 * pi*vz));
}
else if (test == 5)
{
dx = sin(2 * pi*x*z) + 2 * pi*x*z*cos(2 * pi*x*z);
dy = 0;
dz = 2 * pi*x *x * cos(2 * pi*x*z);
dt = 0;
dxdx = 4 * pi*z*cos(2 * pi*x*z) - 4 * pi *pi * x*z *z * sin(2 * pi*x*z);
dxdy = 0;
dxdz = 4 * pi*x*cos(2 * pi*x*z) - 4 * pi *pi * x *x * z*sin(2 * pi*x*z);
dydy = 0;
dydz = 0;
dzdz = -4 * pi *pi * x *x*x * sin(2 * pi*x*z);
dxdt = 0;
dydt = 0;
dzdt = 0;
dtdt = 0;
f = (vx*sin(2 * pi*vx*vz));
f2 = (vx*sin(2 * pi*vx*vz));
}
else if (test == 6)
{
dx = sin(2 * pi*x*y)*sin(2 * pi*z) + 2 * pi*x*y*cos(2 * pi*x*y)*sin(2 * pi*z);
dy = 2 * pi*x *x * cos(2 * pi*x*y)*sin(2 * pi*z);
dz = 2 * pi*x*sin(2 * pi*x*y)*cos(2 * pi*z);
dt = 0;
dxdx = 4 * pi*y*cos(2 * pi*x*y)*sin(2 * pi*z) - 4 * pi *pi * x*y *y * sin(2 * pi*x*y)*sin(2 * pi*z);
dxdy = 4 * pi*x*cos(2 * pi*x*y)*sin(2 * pi*z) - 4 * pi *pi * x *x * y*sin(2 * pi*x*y)*sin(2 * pi*z);
dxdz = 2 * pi*sin(2 * pi*x*y)*cos(2 * pi*z) + 4 * pi *pi * x*y*cos(2 * pi*x*y)*cos(2 * pi*z);
dydy = -4 * pi *pi * x *x*x * sin(2 * pi*x*y)*sin(2 * pi*z);
dydz = 4 * pi *pi * x *x * cos(2 * pi*x*y)*cos(2 * pi*z);
dzdz = -4 * pi *pi * x*sin(2 * pi*x*y)*sin(2 * pi*z);
dxdt = 0;
dydt = 0;
dzdt = 0;
dtdt = 0;
f = (vx*sin(2 * pi*vx*vy)*sin(2 * pi*vz));
f2 = (vx*sin(2 * pi*vx*vy)*sin(2 * pi*vz));
}
else if( test == 7 )
{
dx = 16*exp(-t)*(x-0.5)*cos(8*((z-0.5)*(z-0.5)+(y-0.5)*(y-0.5)+(x-0.5)*(x-0.5)));
dy = 16*exp(-t)*(y-0.5)*cos(8*((z-0.5)*(z-0.5)+(y-0.5)*(y-0.5)+(x-0.5)*(x-0.5)));
dz = 16*exp(-t)*cos(8*((z-0.5)*(z-0.5)+(y-0.5)*(y-0.5)+(x-0.5)*(x-0.5)))*(z-0.5);