本文整理汇总了C++中create_vector函数的典型用法代码示例。如果您正苦于以下问题:C++ create_vector函数的具体用法?C++ create_vector怎么用?C++ create_vector使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了create_vector函数的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: main
int main(int argc, char* argv[])
{
// Time measurement.
TimePeriod cpu_time;
cpu_time.tick();
info("TIME_MAX_ITER = %d", TIME_MAX_ITER);
// Load the mesh file.
Mesh mesh;
H2DReader mloader;
mloader.load("../domain.mesh", &mesh);
// Perform initial mesh refinements.
for (int i=0; i < INIT_GLOB_REF_NUM; i++) mesh.refine_all_elements();
mesh.refine_towards_boundary(1, INIT_BDY_REF_NUM);
// Initialize boundary conditions.
BCTypes bc_types;
bc_types.add_bc_dirichlet(BDY_DIRICHLET);
// Enter Dirichlet boudnary values.
BCValues bc_values;
bc_values.add_zero(BDY_DIRICHLET);
// Create H1 spaces with default shapesets.
H1Space space_T(&mesh, &bc_types, &bc_values, P_INIT);
H1Space space_phi(&mesh, &bc_types, &bc_values, P_INIT);
Hermes::vector<Space*> spaces(&space_T, &space_phi);
// Exact solutions for error evaluation.
ExactSolution T_exact_solution(&mesh, T_exact),
phi_exact_solution(&mesh, phi_exact);
// Solutions in the previous time step.
Solution T_prev_time, phi_prev_time;
Hermes::vector<MeshFunction*> time_iterates(&T_prev_time, &phi_prev_time);
// Solutions in the previous Newton's iteration.
Solution T_prev_newton, phi_prev_newton;
Hermes::vector<Solution*> newton_iterates(&T_prev_newton, &phi_prev_newton);
// Initialize the weak formulation.
WeakForm wf(2);
wf.add_matrix_form(0, 0, jac_TT, jac_TT_ord);
wf.add_matrix_form(0, 1, jac_Tphi, jac_Tphi_ord);
wf.add_vector_form(0, res_T, res_T_ord, HERMES_ANY, &T_prev_time);
wf.add_matrix_form(1, 0, jac_phiT, jac_phiT_ord);
wf.add_matrix_form(1, 1, jac_phiphi, jac_phiphi_ord);
wf.add_vector_form(1, res_phi, res_phi_ord, HERMES_ANY, &phi_prev_time);
// Set initial conditions.
T_prev_time.set_exact(&mesh, T_exact);
phi_prev_time.set_exact(&mesh, phi_exact);
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
solver->set_factorization_scheme(HERMES_REUSE_MATRIX_REORDERING);
// Time stepping.
int t_step = 1;
do {
TIME += TAU;
info("---- Time step %d, t = %g s:", t_step, TIME); t_step++;
info("Projecting to obtain initial vector for the Newton's method.");
scalar* coeff_vec = new scalar[Space::get_num_dofs(spaces)];
OGProjection::project_global(spaces, time_iterates, coeff_vec, matrix_solver);
Solution::vector_to_solutions(coeff_vec, Hermes::vector<Space*>(&space_T, &space_phi),
Hermes::vector<Solution*>(&T_prev_newton, &phi_prev_newton));
// Initialize the FE problem.
bool is_linear = false;
DiscreteProblem dp(&wf, spaces, is_linear);
// Perform Newton's iteration.
info("Newton's iteration...");
bool verbose = false;
if(!solve_newton(coeff_vec, &dp, solver, matrix, rhs, NEWTON_TOL, NEWTON_MAX_ITER, verbose))
error("Newton's iteration failed.");
// Translate the resulting coefficient vector into the Solution sln.
Solution::vector_to_solutions(coeff_vec, spaces, newton_iterates);
delete [] coeff_vec;
// Exact solution for comparison with computational results.
T_exact_solution.update(&mesh, T_exact);
phi_exact_solution.update(&mesh, phi_exact);
// Calculate exact error.
info("Calculating error (exact).");
Hermes::vector<double> exact_errors;
Adapt adaptivity_exact(spaces);
bool solutions_for_adapt = false;
adaptivity_exact.calc_err_exact(Hermes::vector<Solution *>(&T_prev_newton, &phi_prev_newton), Hermes::vector<Solution *>(&T_exact_solution, &phi_exact_solution), &exact_errors, solutions_for_adapt);
double maxerr = std::max(exact_errors[0], exact_errors[1])*100;
//.........这里部分代码省略.........
示例2: fit_ellipsoid
/*******************************************************************************
* int fit_ellipsoid(matrix_t points, vector_t* center, vector_t* lengths)
*
* Fits an ellipsoid to a set of points in 3D space. The principle axes of the
* fitted ellipsoid align with the global coordinate system. Therefore there are
* 6 degrees of freedom defining the ellipsoid: the x,y,z coordinates of the
* centroid and the lengths from the centroid to the surfance in each of the 3
* directions.
*
* matrix_t points is a tall matrix with 3 columns and at least 6 rows. Each row
* must contain the xy&z components of each individual point to be fit. If only
* 6 rows are provided, the resulting ellipsoid will be an exact fit. Otherwise
* the result is a least-squares fit to the overdefined dataset.
*
* vector_t* center is a pointer to a user-created vector which will contain the
* x,y,z position of the centroid of the fit ellipsoid.
*
* vector_t* lengths is a pointer to a user-created vector which will be
* populated with the 3 distances from the surface to the centroid in each of the
* 3 directions.
*******************************************************************************/
int fit_ellipsoid(matrix_t points, vector_t* center, vector_t* lengths){
int i,p;
matrix_t A;
vector_t b;
if(!points.initialized){
printf("ERROR: matrix_t points not initialized\n");
return -1;
}
if(points.cols!=3){
printf("ERROR: matrix_t points must have 3 columns\n");
return -1;
}
p = points.rows;
if(p<6){
printf("ERROR: matrix_t points must have at least 6 rows\n");
return -1;
}
b = create_vector_of_ones(p);
A = create_matrix(p,6);
for(i=0;i<p;i++){
A.data[i][0] = points.data[i][0] * points.data[i][0];
A.data[i][1] = points.data[i][0];
A.data[i][2] = points.data[i][1] * points.data[i][1];
A.data[i][3] = points.data[i][1];
A.data[i][4] = points.data[i][2] * points.data[i][2];
A.data[i][5] = points.data[i][2];
}
vector_t f = lin_system_solve_qr(A,b);
destroy_matrix(&A);
destroy_vector(&b);
// compute center
*center = create_vector(3);
center->data[0] = -f.data[1]/(2*f.data[0]);
center->data[1] = -f.data[3]/(2*f.data[2]);
center->data[2] = -f.data[5]/(2*f.data[4]);
// Solve for lengths
A = create_square_matrix(3);
b = create_vector(3);
// fill in A
A.data[0][0] = (f.data[0] * center->data[0] * center->data[0]) + 1.0;
A.data[0][1] = (f.data[0] * center->data[1] * center->data[1]);
A.data[0][2] = (f.data[0] * center->data[2] * center->data[2]);
A.data[1][0] = (f.data[2] * center->data[0] * center->data[0]);
A.data[1][1] = (f.data[2] * center->data[1] * center->data[1]) + 1.0;
A.data[1][2] = (f.data[2] * center->data[2] * center->data[2]);
A.data[2][0] = (f.data[4] * center->data[0] * center->data[0]);
A.data[2][1] = (f.data[4] * center->data[1] * center->data[1]);
A.data[2][2] = (f.data[4] * center->data[2] * center->data[2]) + 1.0;
// fill in b
b.data[0] = f.data[0];
b.data[1] = f.data[2];
b.data[2] = f.data[4];
// solve for lengths
vector_t scales = lin_system_solve(A, b);
*lengths = create_vector(3);
lengths->data[0] = 1.0/sqrt(scales.data[0]);
lengths->data[1] = 1.0/sqrt(scales.data[1]);
lengths->data[2] = 1.0/sqrt(scales.data[2]);
// cleanup
destroy_vector(&scales);
destroy_matrix(&A);
destroy_vector(&b);
return 0;
}
示例3: main
int main(int argc, char **args)
{
// Time measurement.
TimePeriod cpu_time;
cpu_time.tick();
// Load the mesh.
Mesh mesh;
H3DReader mloader;
mloader.load("lshape_hex.mesh3d", &mesh);
// Perform initial mesh refinement.
for (int i=0; i < INIT_REF_NUM; i++) mesh.refine_all_elements(H3D_H3D_H3D_REFT_HEX_XYZ);
// Create an Hcurl space with default shapeset.
HcurlSpace space(&mesh, bc_types, essential_bc_values, Ord3(P_INIT_X, P_INIT_Y, P_INIT_Z));
// Initialize weak formulation.
WeakForm wf;
wf.add_matrix_form(biform<double, scalar>, biform<Ord, Ord>, HERMES_SYM);
wf.add_matrix_form_surf(biform_surf, biform_surf_ord);
wf.add_vector_form_surf(liform_surf, liform_surf_ord);
// Set exact solution.
ExactSolution exact_sol(&mesh, exact);
// DOF and CPU convergence graphs.
SimpleGraph graph_dof_est, graph_cpu_est, graph_dof_exact, graph_cpu_exact;
// Adaptivity loop.
int as = 1;
bool done = false;
do
{
info("---- Adaptivity step %d:", as);
// Construct globally refined reference mesh and setup reference space.
Space* ref_space = construct_refined_space(&space, 1);
// Initialize discrete problem.
bool is_linear = true;
DiscreteProblem dp(&wf, ref_space, is_linear);
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
// Initialize the preconditioner in the case of SOLVER_AZTECOO.
if (matrix_solver == SOLVER_AZTECOO)
{
((AztecOOSolver*) solver)->set_solver(iterative_method);
((AztecOOSolver*) solver)->set_precond(preconditioner);
// Using default iteration parameters (see solver/aztecoo.h).
}
// Assemble the reference problem.
info("Assembling on reference mesh (ndof: %d).", Space::get_num_dofs(ref_space));
dp.assemble(matrix, rhs);
// Time measurement.
cpu_time.tick();
// Solve the linear system on reference mesh. If successful, obtain the solution.
info("Solving on reference mesh.");
Solution ref_sln(ref_space->get_mesh());
if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), ref_space, &ref_sln);
else error ("Matrix solver failed.\n");
// Time measurement.
cpu_time.tick();
// Project the reference solution on the coarse mesh.
Solution sln(space.get_mesh());
info("Projecting reference solution on coarse mesh.");
OGProjection::project_global(&space, &ref_sln, &sln, matrix_solver, HERMES_HCURL_NORM);
// Time measurement.
cpu_time.tick();
// Output solution and mesh with polynomial orders.
if (solution_output)
{
out_fn_vtk(&sln, "sln", as);
out_orders_vtk(&space, "order", as);
}
// Skip the visualization time.
cpu_time.tick(HERMES_SKIP);
// Calculate element errors and total error estimate.
info("Calculating error estimate and exact error.");
Adapt *adaptivity = new Adapt(&space, HERMES_HCURL_NORM);
bool solutions_for_adapt = true;
double err_est_rel = adaptivity->calc_err_est(&sln, &ref_sln, solutions_for_adapt) * 100;
// Calculate exact error.
solutions_for_adapt = false;
double err_exact_rel = adaptivity->calc_err_exact(&sln, &exact_sol, solutions_for_adapt) * 100;
//.........这里部分代码省略.........
示例4: compute_trajectory
void compute_trajectory(Space *space, DiscreteProblem *dp)
{
info("alpha = (%g, %g, %g, %g), zeta = (%g, %g, %g, %g)",
alpha_ctrl[0], alpha_ctrl[1],
alpha_ctrl[2], alpha_ctrl[3], zeta_ctrl[0],
zeta_ctrl[1], zeta_ctrl[2], zeta_ctrl[3]);
// Newton's loop.
// Fill vector coeff_vec using dof and coeffs arrays in elements.
double *coeff_vec = new double[Space::get_num_dofs(space)];
get_coeff_vector(space, coeff_vec);
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
int it = 1;
while (true)
{
// Obtain the number of degrees of freedom.
int ndof = Space::get_num_dofs(space);
// Assemble the Jacobian matrix and residual vector.
dp->assemble(coeff_vec, matrix, rhs);
// Calculate the l2-norm of residual vector.
double res_l2_norm = get_l2_norm(rhs);
// Info for user.
info("---- Newton iter %d, ndof %d, res. l2 norm %g", it, Space::get_num_dofs(space), res_l2_norm);
// If l2 norm of the residual vector is within tolerance, then quit.
// NOTE: at least one full iteration forced
// here because sometimes the initial
// residual on fine mesh is too small.
if(res_l2_norm < NEWTON_TOL && it > 1) break;
// Multiply the residual vector with -1 since the matrix
// equation reads J(Y^n) \deltaY^{n+1} = -F(Y^n).
for(int i = 0; i < ndof; i++) rhs->set(i, -rhs->get(i));
// Solve the linear system.
if(!solver->solve()) error ("Matrix solver failed.\n");
// Add \deltaY^{n+1} to Y^n.
for (int i = 0; i < ndof; i++) coeff_vec[i] += solver->get_solution()[i];
// If the maximum number of iteration has been reached, then quit.
if (it >= NEWTON_MAX_ITER) error ("Newton method did not converge.");
// Copy coefficients from vector y to elements.
set_coeff_vector(coeff_vec, space);
it++;
}
// Cleanup.
delete matrix;
delete rhs;
delete solver;
delete [] coeff_vec;
}
示例5: QR_decomposition
/*******************************************************************************
* int QR_decomposition(matrix_t A, matrix_t* Q, matrix_t* R)
*
*
*******************************************************************************/
int QR_decomposition(matrix_t A, matrix_t* Q, matrix_t* R){
int i, j, k, s;
int m = A.rows;
int n = A.cols;
vector_t xtemp;
matrix_t Qt, Rt, Qi, F, temp;
if(!A.initialized){
printf("ERROR: matrix not initialized yet\n");
return -1;
}
destroy_matrix(Q);
destroy_matrix(R);
Qt = create_matrix(m,m);
for(i=0;i<m;i++){ // initialize Qt as I
Qt.data[i][i] = 1;
}
Rt = duplicate_matrix(A); // duplicate A to Rt
for(i=0;i<n;i++){ // iterate through columns of A
xtemp = create_vector(m-i); // allocate length, decreases with i
for(j=i;j<m;j++){ // take col of -R from diag down
xtemp.data[j-i] = -Rt.data[j][i];
}
if(Rt.data[i][i] > 0) s = -1; // check the sign
else s = 1;
xtemp.data[0] += s*vector_norm(xtemp); // add norm to 1st element
Qi = create_square_matrix(m); // initialize Qi
F = create_square_matrix(m-i); // initialize shrinking householder_matrix
F = householder_matrix(xtemp); // fill in Househodor
for(j=0;j<i;j++){
Qi.data[j][j] = 1; // fill in partial I matrix
}
for(j=i;j<m;j++){ // fill in remainder (householder_matrix)
for(k=i;k<m;k++){
Qi.data[j][k] = F.data[j-i][k-i];
}
}
// multiply new Qi to old Qtemp
temp = duplicate_matrix(Qt);
destroy_matrix(&Qt);
Qt = multiply_matrices(Qi,temp);
destroy_matrix(&temp);
// same with Rtemp
temp = duplicate_matrix(Rt);
destroy_matrix(&Rt);
Rt = multiply_matrices(Qi,temp);
destroy_matrix(&temp);
// free other allocation used in this step
destroy_matrix(&Qi);
destroy_matrix(&F);
destroy_vector(&xtemp);
}
transpose_matrix(&Qt);
*Q = Qt;
*R = Rt;
return 0;
}
示例6: main
int main()
{
// Time measurement.
TimePeriod cpu_time;
cpu_time.tick();
// Create coarse mesh, set Dirichlet BC, enumerate basis functions.
Space* space = new Space(A, B, NELEM, DIR_BC_LEFT, DIR_BC_RIGHT, P_INIT, NEQ);
// Enumerate basis functions, info for user.
int ndof = Space::get_num_dofs(space);
info("ndof: %d", ndof);
// Initialize the weak formulation.
WeakForm wf;
wf.add_matrix_form(jacobian);
wf.add_vector_form(residual);
// Initialize the FE problem.
bool is_linear = false;
DiscreteProblem *dp_coarse = new DiscreteProblem(&wf, space, is_linear);
// Newton's loop on coarse mesh.
// Fill vector coeff_vec using dof and coeffs arrays in elements.
double *coeff_vec_coarse = new double[Space::get_num_dofs(space)];
get_coeff_vector(space, coeff_vec_coarse);
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix_coarse = create_matrix(matrix_solver);
Vector* rhs_coarse = create_vector(matrix_solver);
Solver* solver_coarse = create_linear_solver(matrix_solver, matrix_coarse, rhs_coarse);
int it = 1;
bool success = false;
while (1) {
// Obtain the number of degrees of freedom.
int ndof_coarse = Space::get_num_dofs(space);
// Assemble the Jacobian matrix and residual vector.
dp_coarse->assemble(coeff_vec_coarse, matrix_coarse, rhs_coarse);
// Calculate the l2-norm of residual vector.
double res_l2_norm = get_l2_norm(rhs_coarse);
// Info for user.
info("---- Newton iter %d, ndof %d, res. l2 norm %g", it, Space::get_num_dofs(space), res_l2_norm);
// If l2 norm of the residual vector is within tolerance, then quit.
// NOTE: at least one full iteration forced
// here because sometimes the initial
// residual on fine mesh is too small.
if(res_l2_norm < NEWTON_TOL_COARSE && it > 1) break;
// Multiply the residual vector with -1 since the matrix
// equation reads J(Y^n) \deltaY^{n+1} = -F(Y^n).
for(int i=0; i < ndof_coarse; i++) rhs_coarse->set(i, -rhs_coarse->get(i));
// Solve the linear system.
if(!solver_coarse->solve())
error ("Matrix solver failed.\n");
// Add \deltaY^{n+1} to Y^n.
for (int i = 0; i < ndof_coarse; i++) coeff_vec_coarse[i] += solver_coarse->get_solution()[i];
// If the maximum number of iteration has been reached, then quit.
if (it >= NEWTON_MAX_ITER) error ("Newton method did not converge.");
// Copy coefficients from vector y to elements.
set_coeff_vector(coeff_vec_coarse, space);
it++;
}
// Cleanup.
delete matrix_coarse;
delete rhs_coarse;
delete solver_coarse;
delete [] coeff_vec_coarse;
delete dp_coarse;
// DOF and CPU convergence graphs.
SimpleGraph graph_dof_est, graph_cpu_est;
SimpleGraph graph_dof_exact, graph_cpu_exact;
// Adaptivity loop:
int as = 1;
bool done = false;
do
{
info("---- Adaptivity step %d:", as);
// Construct globally refined reference mesh and setup reference space.
Space* ref_space = construct_refined_space(space);
// Initialize the FE problem.
bool is_linear = false;
DiscreteProblem* dp = new DiscreteProblem(&wf, ref_space, is_linear);
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
//.........这里部分代码省略.........
示例7: create_beam
int create_beam(int l_size){
beam = (double *)create_vector(l_size);
return 0;
}
示例8: create_noise
int create_noise(int l_size){
noise = (double *)create_vector(l_size);
return 0;
}
示例9: read_beta
int read_beta(){
int number = 2;
int *sizes = create_ivector(number);
ivector_read(&number, proj_size_file, &sizes[0]);
b_lsize = sizes[0]-1;
b_xsize = sizes[1]-1;
b_lvec = create_ivector(b_lsize);
b_xvec = create_vector(b_xsize);
int data_size = (b_lsize+1) * (b_xsize+1);
double **data = create_array(b_lsize+1, b_xsize+1);
char filename[100];
//int xsize_pad = (b_xsize + 7) & ~7;
create_beta(b_lsize, b_xsize);
int n,l,i;
double (*restrict beta)[pmax_prim+1][b_xsize] = (double (*restrict)[pmax_prim+1][b_xsize]) beta_flat;
if(eflag_order_prim!=4)
{
for(n=0;n<pmax_prim+1;n++)
{
char suffix[3] = "";
sprintf(suffix, "%d", n);
filename[0] = '\0';
strcat(filename, proj_data_file);
strcat(filename, "_");
strcat(filename, suffix);
array_read(&data_size, filename, &data[0][0]);
if(n==0)
{
for (i=0; i<b_xsize; i++) b_xvec[i] = data[0][i+1];
for (l=0; l<b_lsize; l++) b_lvec[l] = (int)data[l+1][0];
}
for(l=0; l<b_lsize; l++) beta[l][n][0] = 0.0;
for(i=0; i<b_xsize; i++)
{
for(l=0; l<b_lsize; l++)
{
beta[l][n][i] = data[l+1][i+1];
}
}
}
}
else
{
for(n=0;n<pmax_prim-1;n++)
{
char suffix[3] = "";
sprintf(suffix, "%d", n);
filename[0] = '\0';
strcat(filename, proj_data_file);
strcat(filename, "_");
strcat(filename, suffix);
array_read(&data_size, filename, &data[0][0]);
if(n==0)
{
for (i=0; i<b_xsize; i++) b_xvec[i] = data[0][i+1];
for (l=0; l<b_lsize; l++) b_lvec[l] = (int)data[l+1][0];
}
for(l=0; l<b_lsize; l++) beta[l][n][0] = 0.0;
for(i=0; i<b_xsize; i++)
{
for(l=0; l<b_lsize; l++)
{
beta[l][n][i] = data[l+1][i+1];
}
}
}
filename[0] = '\0';
strcat(filename, proj_data_file);
strcat(filename, "_l1");
array_read(&data_size, filename, &data[0][0]);
for(l=0; l<b_lsize; l++) beta[l][pmax_prim-1][0] = 0.0;
for(i=0; i<b_xsize; i++){
for(l=0; l<b_lsize; l++){
beta[l][pmax_prim-1][i] = data[l+1][i+1];
}
}
filename[0] = '\0';
strcat(filename, proj_data_file);
strcat(filename, "_l2");
//.........这里部分代码省略.........
示例10: create_cl
int create_cl(int l_size){
cl = (double *)create_vector(l_size);
return 0;
}
示例11: main
//------------------------------------------------------------------------------
int main(int argc, char** argv) {
if(argc < 9) {
std::cerr << "usage: " << argv[0]
<< " <platform name> <device type = default | cpu | gpu "
"| acc | all> <device num> <OpenCL source file path>"
" <kernel name> <size> <local size> <vec element width>"
<< std::endl;
exit(EXIT_FAILURE);
}
const int SIZE = atoi(argv[argc - 3]); // number of elements
const int CL_ELEMENT_SIZE = atoi(argv[argc - 1]); // number of per-element
// components
const int CPU_BLOCK_SIZE = 16384; //use block dot product if SIZE divisible
//by this value
const size_t BYTE_SIZE = SIZE * sizeof(real_t);
const int BLOCK_SIZE = atoi(argv[argc - 2]); //local cache for reduction
//equal to local workgroup size
const int REDUCED_SIZE = SIZE / BLOCK_SIZE;
const int REDUCED_BYTE_SIZE = REDUCED_SIZE * sizeof(real_t);
//setup text header that will be prefixed to opencl code
std::ostringstream clheaderStream;
clheaderStream << "#define BLOCK_SIZE " << BLOCK_SIZE << '\n';
clheaderStream << "#define VEC_WIDTH " << CL_ELEMENT_SIZE << '\n';
#ifdef USE_DOUBLE
clheaderStream << "#define DOUBLE\n";
const double EPS = 0.000000001;
#else
const float EPS = 0.00001;
#endif
const bool PROFILE_ENABLE_OPTION = true;
CLEnv clenv = create_clenv(argv[1], argv[2], atoi(argv[3]),
PROFILE_ENABLE_OPTION,
argv[4], argv[5], clheaderStream.str());
cl_int status;
//create input and output matrices
std::vector<real_t> V1 = create_vector(SIZE);
std::vector<real_t> V2 = create_vector(SIZE);
real_t hostDot = std::numeric_limits< real_t >::quiet_NaN();
real_t deviceDot = std::numeric_limits< real_t >::quiet_NaN();
//ALLOCATE DATA AND COPY TO DEVICE
//allocate output buffer on OpenCL device
//the partialReduction array contains a sequence of dot products
//computed on sub-arrays of size BLOCK_SIZE
cl_mem partialReduction = clCreateBuffer(clenv.context,
CL_MEM_WRITE_ONLY,
REDUCED_BYTE_SIZE,
0,
&status);
check_cl_error(status, "clCreateBuffer");
//allocate input buffers on OpenCL devices and copy data
cl_mem devV1 = clCreateBuffer(clenv.context,
CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
BYTE_SIZE,
&V1[0], //<-- copy data from V1
&status);
check_cl_error(status, "clCreateBuffer");
cl_mem devV2 = clCreateBuffer(clenv.context,
CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
BYTE_SIZE,
&V2[0], //<-- copy data from V2
&status);
check_cl_error(status, "clCreateBuffer");
//set kernel parameters
status = clSetKernelArg(clenv.kernel, //kernel
0, //parameter id
sizeof(cl_mem), //size of parameter
&devV1); //pointer to parameter
check_cl_error(status, "clSetKernelArg(V1)");
status = clSetKernelArg(clenv.kernel, //kernel
1, //parameter id
sizeof(cl_mem), //size of parameter
&devV2); //pointer to parameter
check_cl_error(status, "clSetKernelArg(V2)");
status = clSetKernelArg(clenv.kernel, //kernel
2, //parameter id
sizeof(cl_mem), //size of parameter
&partialReduction); //pointer to parameter
check_cl_error(status, "clSetKernelArg(devOut)");
//setup kernel launch configuration
//total number of threads == number of array elements
const size_t globalWorkSize[1] = {SIZE / CL_ELEMENT_SIZE};
//number of per-workgroup local threads
const size_t localWorkSize[1] = {BLOCK_SIZE};
//LAUNCH KERNEL
// make sure all work on the OpenCL device is finished
status = clFinish(clenv.commandQueue);
check_cl_error(status, "clFinish");
cl_event profilingEvent;
timespec kernelStart = {0, 0};
timespec kernelEnd = {0, 0};
clock_gettime(CLOCK_MONOTONIC, &kernelStart);
//launch kernel
status = clEnqueueNDRangeKernel(clenv.commandQueue, //queue
//.........这里部分代码省略.........
示例12: main
int main(int argc, char* argv[])
{
// Load the mesh.
Mesh mesh;
H2DReader mloader;
mloader.load("domain.mesh", &mesh);
// Initial mesh refinements.
for(int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();
// Enter boundary markers.
BCTypes bc_types;
bc_types.add_bc_dirichlet(Hermes::vector<int>(BDY_BOTTOM, BDY_RIGHT, BDY_TOP, BDY_LEFT));
// Enter Dirichlet boundary values.
BCValues bc_values;
bc_values.add_zero(Hermes::vector<int>(BDY_BOTTOM, BDY_RIGHT, BDY_TOP, BDY_LEFT));
// Create an H1 space.
H1Space* phi_space = new H1Space(&mesh, &bc_types, &bc_values, P_INIT);
H1Space* psi_space = new H1Space(&mesh, &bc_types, &bc_values, P_INIT);
int ndof = Space::get_num_dofs(Hermes::vector<Space *>(phi_space, psi_space));
info("ndof = %d.", ndof);
// Initialize previous time level solutions.
Solution phi_prev_time, psi_prev_time;
phi_prev_time.set_exact(&mesh, init_cond_phi);
psi_prev_time.set_exact(&mesh, init_cond_psi);
// Initialize the weak formulation.
WeakForm wf(2);
wf.add_matrix_form(0, 0, callback(biform_euler_0_0));
wf.add_matrix_form(0, 1, callback(biform_euler_0_1));
wf.add_matrix_form(1, 0, callback(biform_euler_1_0));
wf.add_matrix_form(1, 1, callback(biform_euler_1_1));
wf.add_vector_form(0, callback(liform_euler_0), HERMES_ANY, &phi_prev_time);
wf.add_vector_form(1, callback(liform_euler_1), HERMES_ANY, &psi_prev_time);
// Initialize views.
ScalarView view("Psi", new WinGeom(0, 0, 600, 500));
view.fix_scale_width(80);
// Time stepping loop:
int nstep = (int)(T_FINAL/TAU + 0.5);
for(int ts = 1; ts <= nstep; ts++)
{
info("Time step %d:", ts);
info("Solving linear system.");
// Initialize the FE problem.
bool is_linear = true;
DiscreteProblem dp(&wf, Hermes::vector<Space *>(phi_space, psi_space), is_linear);
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
// Assemble the stiffness matrix and right-hand side vector.
info("Assembling the stiffness matrix and right-hand side vector.");
dp.assemble(matrix, rhs);
// Solve the linear system and if successful, obtain the solution.
info("Solving the matrix problem.");
if(solver->solve())
Solution::vector_to_solutions(solver->get_solution(), Hermes::vector<Space *>(phi_space, psi_space), Hermes::vector<Solution *>(&phi_prev_time, &psi_prev_time));
else
error ("Matrix solver failed.\n");
// Show the new time level solution.
char title[100];
sprintf(title, "Time step %d", ts);
view.set_title(title);
view.show(&psi_prev_time);
}
// Wait for all views to be closed.
View::wait();
return 0;
}
示例13: main
int main(int argc, char **args)
{
// Test variable.
int success_test = 1;
for (int i = 0; i < 48; i++) {
for (int j = 0; j < 48; j++) {
info("Config: %d, %d ", i, j);
Mesh mesh;
for (unsigned int k = 0; k < countof(vtcs); k++)
mesh.add_vertex(vtcs[k].x, vtcs[k].y, vtcs[k].z);
unsigned int h1[] = {
hexs[0][i][0] + 1, hexs[0][i][1] + 1, hexs[0][i][2] + 1, hexs[0][i][3] + 1,
hexs[0][i][4] + 1, hexs[0][i][5] + 1, hexs[0][i][6] + 1, hexs[0][i][7] + 1 };
mesh.add_hex(h1);
unsigned int h2[] = {
hexs[1][j][0] + 1, hexs[1][j][1] + 1, hexs[1][j][2] + 1, hexs[1][j][3] + 1,
hexs[1][j][4] + 1, hexs[1][j][5] + 1, hexs[1][j][6] + 1, hexs[1][j][7] + 1 };
mesh.add_hex(h2);
// bc
for (unsigned int k = 0; k < countof(bnd); k++) {
unsigned int facet_idxs[Quad::NUM_VERTICES] = { bnd[k][0] + 1, bnd[k][1] + 1, bnd[k][2] + 1, bnd[k][3] + 1 };
mesh.add_quad_boundary(facet_idxs, bnd[k][4]);
}
mesh.ugh();
// Initialize the space.
H1Space space(&mesh, bc_types, essential_bc_values);
#ifdef XM_YN_ZO
Ord3 ord(4, 4, 4);
#elif defined XM_YN_ZO_2
Ord3 ord(4, 4, 4);
#elif defined X2_Y2_Z2
Ord3 ord(2, 2, 2);
#endif
space.set_uniform_order(ord);
// Initialize the weak formulation.
WeakForm wf;
#ifdef DIRICHLET
wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<Ord, Ord>, HERMES_SYM);
wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>);
#elif defined NEWTON
wf.add_matrix_form(bilinear_form<double, scalar>, bilinear_form<Ord, Ord>, HERMES_SYM);
wf.add_matrix_form_surf(bilinear_form_surf<double, scalar>, bilinear_form_surf<Ord, Ord>);
wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>);
wf.add_vector_form_surf(linear_form_surf<double, scalar>, linear_form_surf<Ord, Ord>);
#endif
// Initialize the FE problem.
bool is_linear = true;
DiscreteProblem dp(&wf, &space, is_linear);
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
// Initialize the preconditioner in the case of SOLVER_AZTECOO.
if (matrix_solver == SOLVER_AZTECOO)
{
((AztecOOSolver*) solver)->set_solver(iterative_method);
((AztecOOSolver*) solver)->set_precond(preconditioner);
// Using default iteration parameters (see solver/aztecoo.h).
}
// Assemble the linear problem.
info("Assembling (ndof: %d).", Space::get_num_dofs(&space));
dp.assemble(matrix, rhs);
// Solve the linear system. If successful, obtain the solution.
info("Solving.");
Solution sln(space.get_mesh());
if(solver->solve()) Solution::vector_to_solution(solver->get_solution(), &space, &sln);
else error ("Matrix solver failed.\n");
ExactSolution ex_sln(&mesh, exact_solution);
// Calculate exact error.
info("Calculating exact error.");
Adapt *adaptivity = new Adapt(&space, HERMES_H1_NORM);
bool solutions_for_adapt = false;
double err_exact = adaptivity->calc_err_exact(&sln, &ex_sln, solutions_for_adapt, HERMES_TOTAL_ERROR_ABS);
if (err_exact > EPS)
{
// Calculated solution is not precise enough.
success_test = 0;
info("failed, error:%g", err_exact);
}
else
info("passed");
// Clean up.
delete matrix;
//.........这里部分代码省略.........
示例14: main
//.........这里部分代码省略.........
wf.add_vector_form_surf(4, callback(linear_form_concentration_inner_edges), H2D_DG_INNER_EDGE,
Hermes::vector<MeshFunction*>(&prev_c, &prev_rho, &prev_rho_v_x, &prev_rho_v_y));
// Initialize the FE problem.
bool is_linear = true;
DiscreteProblem dp(&wf, Hermes::vector<Space*>(&space_rho, &space_rho_v_x, &space_rho_v_y, &space_e, &space_c), is_linear);
// Filters for visualization of pressure and the two components of velocity.
/*
SimpleFilter pressure(calc_pressure_func, Hermes::vector<MeshFunction*>(&sln_rho, &sln_rho_v_x, &sln_rho_v_y, &sln_e));
SimpleFilter u(calc_u_func, Hermes::vector<MeshFunction*>(&sln_rho, &sln_rho_v_x, &sln_rho_v_y, &sln_e));
SimpleFilter w(calc_w_func, Hermes::vector<MeshFunction*>(&sln_rho, &sln_rho_v_x, &sln_rho_v_y, &sln_e));
SimpleFilter Mach_number(calc_Mach_func, Hermes::vector<MeshFunction*>(&sln_rho, &sln_rho_v_x, &sln_rho_v_y, &sln_e));
SimpleFilter entropy_estimate(calc_entropy_estimate_func, Hermes::vector<MeshFunction*>(&sln_rho, &sln_rho_v_x, &sln_rho_v_y, &sln_e));
ScalarView pressure_view("Pressure", new WinGeom(0, 0, 600, 300));
ScalarView Mach_number_view("Mach number", new WinGeom(700, 0, 600, 300));
ScalarView entropy_production_view("Entropy estimate", new WinGeom(0, 400, 600, 300));
VectorView vview("Velocity", new WinGeom(700, 400, 600, 300));
*/
ScalarView s1("w0", new WinGeom(0, 0, 600, 300));
ScalarView s2("w1", new WinGeom(700, 0, 600, 300));
ScalarView s3("w2", new WinGeom(0, 400, 600, 300));
ScalarView s4("w3", new WinGeom(700, 400, 600, 300));
ScalarView s5("Concentration", new WinGeom(350, 200, 600, 300));
// Iteration number.
int iteration = 0;
// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
// Output of the approximate time derivative.
std::ofstream time_der_out("time_der");
for(t = 0.0; t < 3.0; t += TAU) {
info("---- Time step %d, time %3.5f.", iteration++, t);
bool rhs_only = (iteration == 1 ? false : true);
// Assemble stiffness matrix and rhs or just rhs.
if (rhs_only == false) info("Assembling the stiffness matrix and right-hand side vector.");
else info("Assembling the right-hand side vector (only).");
dp.assemble(matrix, rhs, rhs_only);
// Solve the matrix problem.
info("Solving the matrix problem.");
if(solver->solve())
Solution::vector_to_solutions(solver->get_solution(), Hermes::vector<Space *>(&space_rho, &space_rho_v_x,
&space_rho_v_y, &space_e, &space_c), Hermes::vector<Solution *>(&sln_rho, &sln_rho_v_x, &sln_rho_v_y, &sln_e, &sln_c));
else
error ("Matrix solver failed.\n");
// Copy the solutions into the previous time level ones.
prev_rho.copy(&sln_rho);
prev_rho_v_x.copy(&sln_rho_v_x);
prev_rho_v_y.copy(&sln_rho_v_y);
prev_e.copy(&sln_e);
prev_c.copy(&sln_c);
// Visualization.
/*
pressure.reinit();
示例15: create_t_wgt
int create_t_wgt(int l_size){
t_wgt = (double *)create_vector(l_size);
return 0;
}