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);
  if(JFNK == 0)
  {
    // 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;
    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;
  }
  else
    jfnk_cg(dp_coarse, space, MATRIX_SOLVER_TOL, MATRIX_SOLVER_MAXITER,
            JFNK_EPSILON, NEWTON_TOL_COARSE, NEWTON_MAX_ITER, matrix_solver);

  // Cleanup.
  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;
  double ftr_errors[MAX_ELEM_NUM];        // This array decides what 
                                          // elements will be refined.

  bool done = false;
  do
  {
    info("---- Adaptivity step %d:", as); 
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  Hermes2D hermes_2D;

  // Load the mesh.
  Mesh mesh;
  H2DReader mloader;
  mloader.load("domain-excentric.mesh", &mesh);
  //mloader.load("domain-concentric.mesh", &mesh);

  // Initial mesh refinements.
  for (int i=0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();
  mesh.refine_towards_boundary("Inner", INIT_BDY_REF_NUM_INNER, false);  // true for anisotropic refinements
  mesh.refine_towards_boundary("Outer", INIT_BDY_REF_NUM_OUTER, false);  // false for isotropic refinements

  // Initialize boundary conditions.
  EssentialBCNonConstX bc_inner_vel_x(std::string("Inner"), VEL, STARTUP_TIME);
  EssentialBCNonConstY bc_inner_vel_y(std::string("Inner"), VEL, STARTUP_TIME);
  DefaultEssentialBCConst bc_outer_vel(std::string("Outer"), 0.0);
  EssentialBCs bcs_vel_x(Hermes::vector<EssentialBoundaryCondition *>(&bc_inner_vel_x, &bc_outer_vel));
  EssentialBCs bcs_vel_y(Hermes::vector<EssentialBoundaryCondition *>(&bc_inner_vel_y, &bc_outer_vel));
  EssentialBCs bcs_pressure;

  // Spaces for velocity components and pressure.
  H1Space xvel_space(&mesh, &bcs_vel_x, P_INIT_VEL);
  H1Space yvel_space(&mesh, &bcs_vel_y, P_INIT_VEL);
#ifdef PRESSURE_IN_L2
  L2Space p_space(&mesh, &bcs_pressure, P_INIT_PRESSURE);
#else
  H1Space p_space(&mesh, &bcs_pressure, P_INIT_PRESSURE);
#endif
  Hermes::vector<Space *> spaces = Hermes::vector<Space *>(&xvel_space, &yvel_space, &p_space);

  // Calculate and report the number of degrees of freedom.
  int ndof = Space::get_num_dofs(spaces);
  info("ndof = %d.", ndof);

  // Define projection norms.
  ProjNormType vel_proj_norm = HERMES_H1_NORM;
#ifdef PRESSURE_IN_L2
  ProjNormType p_proj_norm = HERMES_L2_NORM;
#else
  ProjNormType p_proj_norm = HERMES_H1_NORM;
#endif

  // Solutions for the Newton's iteration and time stepping.
  info("Setting initial conditions.");
  Solution xvel_prev_time, yvel_prev_time, p_prev_time; 
  xvel_prev_time.set_zero(&mesh);
  yvel_prev_time.set_zero(&mesh);
  p_prev_time.set_zero(&mesh);
  Hermes::vector<Solution*> slns = Hermes::vector<Solution*>(&xvel_prev_time, &yvel_prev_time, 
                                                             &p_prev_time);

  // Initialize weak formulation.
  WeakForm* wf = new WeakFormNSNewton(STOKES, RE, TAU, &xvel_prev_time, &yvel_prev_time);

  // Initialize the FE problem.
  DiscreteProblem dp(wf, spaces);

  // 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 views.
  VectorView vview("velocity [m/s]", new WinGeom(0, 0, 600, 500));
  ScalarView pview("pressure [Pa]", new WinGeom(610, 0, 600, 500));
  //vview.set_min_max_range(0, 1.6);
  vview.fix_scale_width(80);
  //pview.set_min_max_range(-0.9, 1.0);
  pview.fix_scale_width(80);
  pview.show_mesh(true);

  // Project the initial condition on the FE space to obtain initial
  // coefficient vector for the Newton's method.
  scalar* coeff_vec = new scalar[Space::get_num_dofs(spaces)];
  // Newton's vector is set to zero (no OG projection needed).
  memset(coeff_vec, 0, ndof * sizeof(double));
  /*
  // This can be used for more complicated initial conditions.
    info("Projecting initial condition to obtain initial vector for the Newton's method.");
    OGProjection::project_global(spaces, slns, coeff_vec, matrix_solver, 
                                 Hermes::vector<ProjNormType>(vel_proj_norm, vel_proj_norm, p_proj_norm));
  */

  // Time-stepping loop:
  char title[100];
  int num_time_steps = T_FINAL / TAU;
  for (int ts = 1; ts <= num_time_steps; ts++)
  {
    current_time += TAU;
    info("---- Time step %d, time = %g:", ts, current_time);

    // Update time-dependent essential BCs.
    info("Updating time-dependent essential BC.");
    Space::update_essential_bc_values(Hermes::vector<Space *>(&xvel_space, &yvel_space, &p_space), current_time);

    // Perform Newton's iteration.
    info("Solving nonlinear problem:");
//.........这里部分代码省略.........

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);

  double elem_errors[MAX_ELEM_NUM];      // This array decides what 
                                         // elements will be refined.
  ElemPtr2 ref_elem_pairs[MAX_ELEM_NUM]; // To store element pairs from the 
                                         // FTR solution. Decides how 
                                         // elements will be hp-refined. 
  for (int i=0; i < MAX_ELEM_NUM; i++) 
  {
    ref_elem_pairs[i][0] = new Element();
    ref_elem_pairs[i][1] = new Element();
  }

  // DOF and CPU convergence graphs.
  SimpleGraph graph_dof_exact, graph_cpu_exact;

  // Adaptivity loop:
  int as = 1;
  bool done = false;
  do
  {
    info("---- Adaptivity step %d:", as);

    // 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;
    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;
//.........这里部分代码省略.........

int main(int argc, char **args) 
{
  // Load the mesh.
  Mesh mesh;
  H3DReader mloader;
  mloader.load("fichera-corner.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 H1 space with default shapeset.
  H1Space 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(bilinear_form<double, double>, bilinear_form<Ord, Ord>, HERMES_SYM, HERMES_ANY);
  wf.add_vector_form(linear_form<double, double>, linear_form<Ord, Ord>, HERMES_ANY);

  // Set exact solution.
  ExactSolution exact(&mesh, fndd);

  // DOF and CPU convergence graphs.
  SimpleGraph graph_dof_est, graph_cpu_est, graph_dof_exact, graph_cpu_exact;

  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

  // 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);

    // 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_H1_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, solutions_for_adapt) * 100;

    // Report results.
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  // Load the mesh.
  Mesh mesh;
  H2DReader mloader;
  mloader.load("GAMM-channel.mesh", &mesh);

  // Perform initial mesh refinements.
  for (int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements(0);
  //mesh.refine_towards_boundary(BDY_SOLID_WALL_BOTTOM, 2);

  // Initialize boundary condition types and spaces with default shapesets.
  L2Space space_rho(&mesh, P_INIT);
  L2Space space_rho_v_x(&mesh, P_INIT);
  L2Space space_rho_v_y(&mesh, P_INIT);
  L2Space space_e(&mesh, P_INIT);
  int ndof = Space::get_num_dofs(Hermes::vector<Space*>(&space_rho, &space_rho_v_x, &space_rho_v_y, &space_e));
  info("ndof: %d", ndof);

  // Initialize solutions, set initial conditions.
  InitialSolutionEulerDensity prev_rho(&mesh, RHO_EXT);
  InitialSolutionEulerDensityVelX prev_rho_v_x(&mesh, RHO_EXT * V1_EXT);
  InitialSolutionEulerDensityVelY prev_rho_v_y(&mesh, RHO_EXT * V2_EXT);
  InitialSolutionEulerDensityEnergy prev_e(&mesh, QuantityCalculator::calc_energy(RHO_EXT, RHO_EXT * V1_EXT, RHO_EXT * V2_EXT, P_EXT, KAPPA));

  // Numerical flux.
  OsherSolomonNumericalFlux num_flux(KAPPA);

  // Initialize weak formulation.
  EulerEquationsWeakFormSemiImplicitMultiComponent wf(&num_flux, KAPPA, RHO_EXT, V1_EXT, V2_EXT, P_EXT, BDY_SOLID_WALL_BOTTOM, BDY_SOLID_WALL_TOP, 
    BDY_INLET, BDY_OUTLET, &prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e, (P_INIT == 0));

  // 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), is_linear);

  // If the FE problem is in fact a FV problem.
  if(P_INIT == 0) 
    dp.set_fvm();  

  // Filters for visualization of Mach number, pressure and entropy.
  MachNumberFilter Mach_number(Hermes::vector<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e), KAPPA);
  PressureFilter pressure(Hermes::vector<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e), KAPPA);
  EntropyFilter entropy(Hermes::vector<MeshFunction*>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e), KAPPA, RHO_EXT, P_EXT);

  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));
  
  ScalarView s1("1", new WinGeom(0, 0, 600, 300));
  ScalarView s2("2", new WinGeom(700, 0, 600, 300));
  ScalarView s3("3", new WinGeom(0, 400, 600, 300));
  ScalarView s4("4", new WinGeom(700, 400, 600, 300));

  // 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);

  // Set up CFL calculation class.
  CFLCalculation CFL(CFL_NUMBER, KAPPA);

  int iteration = 0; double t = 0;
  for(t = 0.0; t < 3.0; t += time_step) {
    info("---- Time step %d, time %3.5f.", iteration++, t);

    // Set the current time step.
    wf.set_time_step(time_step);

    // Assemble the stiffness matrix and rhs.
    info("Assembling the stiffness matrix and right-hand side vector.");
    dp.assemble(matrix, rhs);

    // Solve the matrix problem.
    info("Solving the matrix problem.");
    scalar* solution_vector = NULL;
    if(solver->solve()) {
      solution_vector = solver->get_solution();
      Solution::vector_to_solutions(solution_vector, Hermes::vector<Space *>(&space_rho, &space_rho_v_x, 
      &space_rho_v_y, &space_e), Hermes::vector<Solution *>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));
    }
    else
      error ("Matrix solver failed.\n");

    if(SHOCK_CAPTURING) {
      DiscontinuityDetector discontinuity_detector(Hermes::vector<Space *>(&space_rho, &space_rho_v_x, 
        &space_rho_v_y, &space_e), Hermes::vector<Solution *>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));

      std::set<int> discontinuous_elements = discontinuity_detector.get_discontinuous_element_ids(DISCONTINUITY_DETECTOR_PARAM);

      FluxLimiter flux_limiter(solution_vector, Hermes::vector<Space *>(&space_rho, &space_rho_v_x, 
        &space_rho_v_y, &space_e), Hermes::vector<Solution *>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));

      flux_limiter.limit_according_to_detector(discontinuous_elements);
    }

    CFL.calculate_semi_implicit(Hermes::vector<Solution *>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e), &mesh, time_step);

    // Visualization.
    
//.........这里部分代码省略.........

//.........这里部分代码省略.........
#elif defined WITH_MUMPS
	MatrixSolverType matrix_solver = SOLVER_MUMPS; 
#endif

#ifdef RHS2
	WeakForm wf;
	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>, HERMES_ANY_INT, &fsln);

	// Initialize discrete problem.
	bool is_linear = true;
	DiscreteProblem dp(&wf, &space, is_linear);
#elif defined SYS3
	WeakForm wf(3);
	wf.add_matrix_form(0, 0, biform_1_1<double, scalar>, biform_1_1<Ord, Ord>, HERMES_SYM);
	wf.add_matrix_form(0, 1, biform_1_2<double, scalar>, biform_1_2<Ord, Ord>, HERMES_NONSYM);
	wf.add_vector_form(0, liform_1<double, scalar>, liform_1<Ord, Ord>);

	wf.add_matrix_form(1, 1, biform_2_2<double, scalar>, biform_2_2<Ord, Ord>, HERMES_SYM);
	wf.add_matrix_form(1, 2, biform_2_3<double, scalar>, biform_2_3<Ord, Ord>, HERMES_NONSYM);
	wf.add_vector_form(1, liform_2<double, scalar>, liform_2<Ord, Ord>);

	wf.add_matrix_form(2, 2, biform_3_3<double, scalar>, biform_3_3<Ord, Ord>, HERMES_SYM);

	// Initialize discrete problem.
	bool is_linear = true;
	DiscreteProblem dp(&wf, Hermes::vector<Space *>(&space1, &space2, &space3), is_linear);
#endif
	// Time measurement.
	TimePeriod cpu_time;
	cpu_time.tick();
  
	// 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 stiffness matrix and load vector.
	dp.assemble(matrix, rhs);

	// Solve the linear system. If successful, obtain the solution.
	info("Solving the linear problem.");
	bool solved = solver->solve();

	// Time measurement.
	cpu_time.tick();
	// Print timing information.
	info("Solution and mesh with polynomial orders saved. Total running time: %g s", cpu_time.accumulated());

	// Time measurement.
	TimePeriod sln_time;
	sln_time.tick();

	if (solved) {
#ifdef RHS2
		// Solve the linear system. If successful, obtain the solution.
		info("Solving the linear problem.");
                Solution sln(&mesh1);

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);

  // Initialize solution views (their titles will be2 updated in each time step).
  ScalarView sview_T("", new WinGeom(0, 0, 500, 400));
  sview_T.fix_scale_width(50);
  ScalarView sview_phi("", new WinGeom(0, 500, 500, 400));
  sview_phi.fix_scale_width(50);
  ScalarView sview_T_exact("", new WinGeom(550, 0, 500, 400));
  sview_T_exact.fix_scale_width(50);
  ScalarView sview_phi_exact("", new WinGeom(550, 500, 500, 400));
  sview_phi_exact.fix_scale_width(50);
  char title[100]; // Character array to store the title for an actual view and time step.

  // 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 = true;
    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;
    
    // Show the new time level solution.
//.........这里部分代码省略.........

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, 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;
    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;

    // Test variable.
    int success_test = 1;

    // 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);

//.........这里部分代码省略.........

int main(int argc, char **args)
{
  // Test variable.
  int success_test = 1;

  // Check the number of command-line parameters.
  if (argc < 2) {
    info("Use x, y, z, xy, xz, yz, or xyz as a command-line parameter.");
    error("Not enough command-line parameters.");
  }

  // Determine anisotropy type from the command-line parameter.
  ANISO_TYPE = parse_aniso_type(args[1]);

  // Load the mesh.
  Mesh mesh;
  H3DReader mesh_loader;
  mesh_loader.load("hex-0-1.mesh3d", &mesh);

  // Assign the lowest possible directional polynomial degrees so that the problem's NDOF >= 1.
  assign_poly_degrees();

  // Create an H1 space with default shapeset.
  info("Setting directional polynomial degrees %d, %d, %d.", P_INIT_X, P_INIT_Y, P_INIT_Z);
  H1Space 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(bilinear_form<double, scalar>, bilinear_form<Ord, Ord>, HERMES_SYM, HERMES_ANY);
  wf.add_vector_form(linear_form<double, scalar>, linear_form<Ord, Ord>, HERMES_ANY);

  // Set exact solution.
  ExactSolution exact(&mesh, fndd);

  // Time measurement.
  TimePeriod cpu_time;
  cpu_time.tick();

  // 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");
		  success_test = 0;
	  }

    // 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);

    // 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.
//.........这里部分代码省略.........

 adjMatrix() { matrix=NULL; create_matrix(10); }

int main(int argc, char *argv[]) {
  int rows = 8; 
  int columns = 7; //4;
  char *matrix = create_matrix(rows, columns);
  
  /*
  set_element(matrix, rows, columns, 0, 0, 1);
  set_element(matrix, rows, columns, 0, 1, 0);
  set_element(matrix, rows, columns, 0, 2, 0);
  set_element(matrix, rows, columns, 0, 3, 0);
  
  set_element(matrix, rows, columns, 1, 0, 0);
  set_element(matrix, rows, columns, 1, 1, 0);
  set_element(matrix, rows, columns, 1, 2, 0);
  set_element(matrix, rows, columns, 1, 3, 0);
  
  set_element(matrix, rows, columns, 2, 0, 0);
  set_element(matrix, rows, columns, 2, 1, 0);
  set_element(matrix, rows, columns, 2, 2, 2);
  set_element(matrix, rows, columns, 2, 3, 3);
  */
  
  set_element(matrix, rows, columns, 0, 0, 1);
  set_element(matrix, rows, columns, 0, 1, 0);
  set_element(matrix, rows, columns, 0, 2, 0);
  set_element(matrix, rows, columns, 0, 3, 0);
  set_element(matrix, rows, columns, 0, 4, 0);
  set_element(matrix, rows, columns, 0, 5, 0);
  set_element(matrix, rows, columns, 0, 6, 0);
  
  set_element(matrix, rows, columns, 1, 0, 0);
  set_element(matrix, rows, columns, 1, 1, 0);
  set_element(matrix, rows, columns, 1, 2, 0);
  set_element(matrix, rows, columns, 1, 3, 0);
  set_element(matrix, rows, columns, 1, 4, 0);
  set_element(matrix, rows, columns, 1, 5, 0);
  set_element(matrix, rows, columns, 1, 6, 0);
  
  set_element(matrix, rows, columns, 2, 0, 0);
  set_element(matrix, rows, columns, 2, 1, 0);
  set_element(matrix, rows, columns, 2, 2, 0);
  set_element(matrix, rows, columns, 2, 3, 0);
  set_element(matrix, rows, columns, 2, 4, 0);
  set_element(matrix, rows, columns, 2, 5, 0);
  set_element(matrix, rows, columns, 2, 6, 0);
  
  set_element(matrix, rows, columns, 3, 0, 0);
  set_element(matrix, rows, columns, 3, 1, 0);
  set_element(matrix, rows, columns, 3, 2, 0);
  set_element(matrix, rows, columns, 3, 3, 0);
  set_element(matrix, rows, columns, 3, 4, 0);
  set_element(matrix, rows, columns, 3, 5, 0);
  set_element(matrix, rows, columns, 3, 6, 0);
  
  set_element(matrix, rows, columns, 4, 0, 0);
  set_element(matrix, rows, columns, 4, 1, 0);
  set_element(matrix, rows, columns, 4, 2, 0);
  set_element(matrix, rows, columns, 4, 3, 0);
  set_element(matrix, rows, columns, 4, 4, 0);
  set_element(matrix, rows, columns, 4, 5, 0);
  set_element(matrix, rows, columns, 4, 6, 0);
  
  set_element(matrix, rows, columns, 5, 0, 0);
  set_element(matrix, rows, columns, 5, 1, 0);
  set_element(matrix, rows, columns, 5, 2, 0);
  set_element(matrix, rows, columns, 5, 3, 0);
  set_element(matrix, rows, columns, 5, 4, 0);
  set_element(matrix, rows, columns, 5, 5, 0);
  set_element(matrix, rows, columns, 5, 6, 0);
  
  set_element(matrix, rows, columns, 6, 0, 0);
  set_element(matrix, rows, columns, 6, 1, 0);
  set_element(matrix, rows, columns, 6, 2, 0);
  set_element(matrix, rows, columns, 6, 3, 0);
  set_element(matrix, rows, columns, 6, 4, 0);
  set_element(matrix, rows, columns, 6, 5, 0);
  set_element(matrix, rows, columns, 6, 6, 0);
  
  set_element(matrix, rows, columns, 7, 0, 2);
  set_element(matrix, rows, columns, 7, 1, 0);
  set_element(matrix, rows, columns, 7, 2, 0);
  set_element(matrix, rows, columns, 7, 3, 0);
  set_element(matrix, rows, columns, 7, 4, 0);
  set_element(matrix, rows, columns, 7, 5, 3);
  set_element(matrix, rows, columns, 7, 6, 3);
  
  print_matrix(matrix, rows, columns);
  
  printf("%d\n", search(matrix, rows, columns, 0, 0));
  return 0;
}

int main(int argc, char* argv[])
{
  // Instantiate a class with global functions.
  Hermes2D hermes2d;

  // Load the mesh.
  Mesh mesh, mesh1;
  H2DReader mloader;
  mloader.load("domain.mesh", &mesh);

  // Perform uniform mesh refinement.
  mesh.refine_all_elements();

  // Initialize boundary conditions.
  DefaultEssentialBCConst zero_disp("Bottom", 0.0);
  EssentialBCs bcs(&zero_disp);

  // Create x- and y- displacement space using the default H1 shapeset.
  H1Space u1_space(&mesh, &bcs, P_INIT);
  H1Space u2_space(&mesh, &bcs, P_INIT);
  int ndof = Space::get_num_dofs(Hermes::vector<Space *>(&u1_space, &u2_space));
  info("ndof = %d", ndof);

  // Initialize the weak formulation.
  CustomWeakFormLinearElasticity wf(E, nu, rho*g1, "Top", f0, f1);

  // Initialize the FE problem.
  DiscreteProblem dp(&wf, Hermes::vector<Space *>(&u1_space, &u2_space));

  // 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);

  // Initial coefficient vector for the Newton's method.  
  scalar* coeff_vec = new scalar[ndof];
  memset(coeff_vec, 0, ndof*sizeof(scalar));

  // Perform Newton's iteration.
  bool verbose = true;
  bool jacobian_changed = true;
  if (!hermes2d.solve_newton(coeff_vec, &dp, solver, matrix, rhs, jacobian_changed,
      NEWTON_TOL, NEWTON_MAX_ITER, verbose)) error("Newton's iteration failed.");

  // Translate the resulting coefficient vector into the Solution sln.
  Solution u1_sln, u2_sln;
  Solution::vector_to_solutions(coeff_vec, Hermes::vector<Space *>(&u1_space, &u2_space), 
                                Hermes::vector<Solution *>(&u1_sln, &u2_sln));
  
  // Visualize the solution.
  ScalarView view("Von Mises stress [Pa]", new WinGeom(0, 0, 800, 400));
  double lambda = (E * nu) / ((1 + nu) * (1 - 2*nu));  // First Lame constant.
  double mu = E / (2*(1 + nu));                        // Second Lame constant.
  VonMisesFilter stress(Hermes::vector<MeshFunction *>(&u1_sln, &u2_sln), lambda, mu);
  view.show_mesh(false);
  view.show(&stress, HERMES_EPS_HIGH, H2D_FN_VAL_0, &u1_sln, &u2_sln, 1.5e5);

  // Wait for the view to be closed.
  View::wait();

  // Clean up.
  delete [] coeff_vec;
  delete solver;
  delete matrix;
  delete rhs;

  return 0;
}

int main(void) {
  clock_t begin, end;
  double time_spent;
  begin = clock();

  matrix* a = create_matrix(4, 4);
  value temp_a[16] = { 18, 60, 57, 96,
		       41, 24, 99, 58,
		       14, 30, 97, 66,
		       51, 13, 19, 85 };
  insert_array(temp_a, a);

  matrix* b = create_matrix(4, 4);
  assert(insert_array(temp_a, b));


  //tests check_boundaries
  assert(check_boundaries(1,1,a));
  assert(check_boundaries(4,4,a));
  assert(!check_boundaries(4,5,a));
  assert(!check_boundaries(5,4,a));
  assert(!check_boundaries(0,1,a));
  assert(!check_boundaries(1,0,a));
  assert(!check_boundaries(-1,1,a));
  assert(!check_boundaries(1,-1,a));


  //tests compare_matrices,insert_value and get_value
  assert(compare_matrices(a,b));
  assert(insert_value(10,1,1,b));
  assert(!compare_matrices(a,b));
  assert(get_value(1,1,b)==10);
  assert(insert_value(18,1,1,b));
  assert(compare_matrices(a,b));


  //tests is_matrix
  matrix* c=a;
  assert(compare_matrices(a,c));
  assert(!is_matrix(a,b));
  assert(is_matrix(a,c));


  //tests insert_value by trying to go outside the matrix
  assert(insert_value(1,1,1,c));
  assert(insert_value(2,2,2,c));
  assert(insert_value(3,3,3,c));
  assert(insert_value(4,4,4,c));
  assert(!insert_value(5,5,5,c));
  assert(!insert_value(-1,-1,-1,c));
  assert(!insert_value(-1,-1,1,c));
  assert(!insert_value(-1,1,-1,c));

  //test get_value
  assert(get_value(1,1,c)==1);
  assert(get_value(2,2,c)==2);
  assert(get_value(3,3,c)==3);
  assert(get_value(4,4,c)==4);
  assert(get_value(0,0,c)==0);
  assert(get_value(1,-1,c)==0);
  assert(get_value(-1,1,c)==0);
  assert(get_value(5,5,c)==0);

  //tests insert and get without boundary checks
  insert_value_without_check(4,1,1,c);
  insert_value_without_check(3,2,2,c);
  insert_value_without_check(2,3,3,c);
  insert_value_without_check(1,4,4,c);
  assert(get_value_without_check(1,1,c)==4);
  assert(get_value_without_check(2,2,c)==3);
  assert(get_value_without_check(3,3,c)==2);
  assert(get_value_without_check(4,4,c)==1);

  //tests add_matrices
  value temp_b[16]={
    36,120,114,192,
    82,48,198,116,
    28, 60, 194,132,
    102,26,38,170};
  assert(insert_array(temp_b,a));
  matrix* d = create_matrix(4, 4);
  assert(add_matrices(b,b,d));
  assert(compare_matrices(d,a));

  //tests subtract_matrices
  value temp_c[16]={
    0,0,0,0,
    0,0,0,0,
    0, 0, 0,0,
    0,0,0,0};
  assert(insert_array(temp_c,a));
  assert(subtract_matrices(b,b,d));
  assert(compare_matrices(d,a));

  //tests sum_of_row
  assert(insert_array(temp_a,a));
  assert(sum_of_row(1,a)==231);
  assert(sum_of_row(4,a)==168);
  assert(sum_of_row(0,a)==0);
  assert(sum_of_row(5,a)==0);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  info("Desired number of eigenvalues: %d.", NUMBER_OF_EIGENVALUES);

  // Load the mesh.
  info("Loading and refining mesh...");
  Mesh mesh;
  H3DReader mloader;
  mloader.load("hexahedron.mesh3d", &mesh);

  // Perform initial mesh refinements (optional).
  for (int i = 0; i < INIT_REF_NUM; i++) mesh.refine_all_elements(H3D_H3D_H3D_REFT_HEX_XYZ);

  // Create an H1 space with default shapeset.
  H1Space space(&mesh, bc_types, essential_bc_values, Ord3(P_INIT_X, P_INIT_Y, P_INIT_Z));
  int ndof = Space::get_num_dofs(&space);
  info("ndof: %d.", ndof);

  // Initialize the weak formulation for the left hand side, i.e., H.
  info("Initializing weak form...");
  WeakForm wf_left, wf_right;
  wf_left.add_matrix_form(bilinear_form_left, bilinear_form_left_ord, HERMES_SYM, HERMES_ANY );
  wf_right.add_matrix_form(callback(bilinear_form_right), HERMES_SYM, HERMES_ANY );

  // Initialize matrices and matrix solver.
  SparseMatrix* matrix_left = create_matrix(matrix_solver);
  SparseMatrix* matrix_right = create_matrix(matrix_solver);
  Solver* solver = create_linear_solver(matrix_solver, matrix_left);

  // Assemble the matrices.
  info("Assembling matrices...");
  bool is_linear = true;
  DiscreteProblem dp_left(&wf_left, &space, is_linear);
  dp_left.assemble(matrix_left);
  DiscreteProblem dp_right(&wf_right, &space, is_linear);
  dp_right.assemble(matrix_right);

  // Write matrix_left in MatrixMarket format.
  write_matrix_mm("mat_left.mtx", matrix_left);

  // Write matrix_right in MatrixMarket format.
  write_matrix_mm("mat_right.mtx", matrix_right);

  // Calling Python eigensolver. Solution will be written to "eivecs.dat".
  info("Using eigensolver...");
  char call_cmd[255];
  sprintf(call_cmd, "python solveGenEigenFromMtx.py mat_left.mtx mat_right.mtx %g %d %g %d", 
	  TARGET_VALUE, NUMBER_OF_EIGENVALUES, TOL, MAX_ITER);
  system(call_cmd);

  // Initia