C++ LinearSolver类代码示例

void testBanded()
{
    if (os_) *os_ << "testBanded()\n";

    LinearSolver<> solver;

    ublas::banded_matrix<double> A(2,2,1,1);
    A(0,0) = 1.;
    A(0,1) = 2.;
    A(1,0) = 3.;
    A(1,1) = 4.;

    ublas::vector<double> y(2);
    y(0) = 5.;
    y(1) = 11.;

    ublas::vector<double> x = solver.solve(A, y);

    if (os_) *os_ << "A: " << A << endl;
    if (os_) *os_ << "y: " << y << endl;
    if (os_) *os_ << "x: " << x << endl;

    unit_assert_equal(x(0), 1., 1e-14);
    unit_assert_equal(x(1), 2., 1e-14);
}

void testBandedComplex()
{
    if (os_) *os_ << "testBandedComplex()\n";

    LinearSolver<> solver;

    ublas::banded_matrix< complex<double> > A(2,2,1,1);
    A(0,0) = 1.;
    A(0,1) = 2.;
    A(1,0) = 3.;
    A(1,1) = 4.;

    ublas::vector< complex<double> > y(2);
    y(0) = 5.;
    y(1) = 11.;

    ublas::vector< complex<double> > x = solver.solve(A, y);

    if (os_) *os_ << "A: " << A << endl;
    if (os_) *os_ << "y: " << y << endl;
    if (os_) *os_ << "x: " << x << endl;

    unit_assert(norm(x(0)-1.) < 1e-14);
    unit_assert(norm(x(1)-2.) < 1e-14);
}

void RBF::computeFunctionForData()
{
	switch(acceleration)
	{
		case FastMultipole:
			fmmComputeFunction();
		case None:
		default:
			int n = data->fnc.size();
			printf("Solving linear equations: \n"); fflush(stdout);
			LinearSolver rbfSolver;
			SparseMatrix rbfMatrix(n);
			printf("Constructing matrix ... "); fflush(stdout);
			for(int i=0; i<n; i++)
			{
				for(int j=0; j<n; j++)
				{
					//printf("%d %d ", i,j); fflush(stdout);
					double val = computeKernel(i,j);
					//printf("%lf\n", val); fflush(stdout);
					rbfMatrix.push_back(i,j,val);
				}
			}
			printf("Done\n"); fflush(stdout);
			rbfSolver.setMatrix(&rbfMatrix);
			printf("Running BiCGSTAB Iterations ... "); fflush(stdout);
			rbfSolver.biCGStab(data->fnc, coeff);
			printf("Done\n"); fflush(stdout);
	}
}

void testComplex()
{
    if (os_) *os_ << "testComplex()\n";

    LinearSolver<> solver;

    ublas::matrix< complex<double> > A(2,2);
    A(0,0) = 1;
    A(0,1) = 2;
    A(1,0) = 3;
    A(1,1) = 4;

    ublas::vector< complex<double> > y(2);
    y(0) = 5;
    y(1) = 11;

    ublas::vector< complex<double> > x = solver.solve(A, y);

    if (os_) *os_ << "A: " << A << endl;
    if (os_) *os_ << "y: " << y << endl;
    if (os_) *os_ << "x: " << x << endl;

    unit_assert(x(0) == 1.);
    unit_assert(x(1) == 2.);
}

void testDoubleQR()
{
    if (os_) *os_ << "testDoubleQR()\n";

    LinearSolver<LinearSolverType_QR> solver;

    ublas::matrix<double> A(2,2);
    A(0,0) = 1.;
    A(0,1) = 2.;
    A(1,0) = 3.;
    A(1,1) = 4.;

    ublas::vector<double> y(2);
    y(0) = 5.;
    y(1) = 11.;

    ublas::vector<double> x = solver.solve(A, y);

    if (os_) *os_ << "A: " << A << endl;
    if (os_) *os_ << "y: " << y << endl;
    if (os_) *os_ << "x: " << x << endl;

    if (os_) *os_ << x(0) << " - 1. = " << x(0) - 1. << endl;

    unit_assert_equal(x(0), 1., 1e-14);
    unit_assert_equal(x(1), 2., 1e-14);
}

void MechanicalOperations::m_solveSystem()
{
    LinearSolver* s = ctx->get<LinearSolver>(ctx->getTags(), BaseContext::SearchDown);
    if (!s)
    {
        ctx->serr << "ERROR:  requires a LinearSolver."<<ctx->sendl;
        return;
    }
    s->solveSystem();
}

    boost::numeric::ublas::vector<T> solve(
        const boost::numeric::ublas::matrix<T>& A,
        const boost::numeric::ublas::vector<T>& x)
    {
        LinearSolver<LinearSolverType_QR> solver;

        boost::numeric::ublas::vector<T> y = solver.solve(A, x);

        return y;
    }

void MechanicalOperations::m_setSystemRHVector(core::MultiVecDerivId v)
{
    LinearSolver* s = ctx->get<LinearSolver>(ctx->getTags(), BaseContext::SearchDown);
    if (!s)
    {
        ctx->serr << "ERROR:  requires a LinearSolver."<<ctx->sendl;
        return;
    }
    s->setSystemRHVector(v);
}

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

  // Load the mesh.
  Mesh mesh;
  MeshReaderH2D mloader;
  mloader.load("square_2_elem.mesh", &mesh);

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

  // Create an H1 space with default shapeset.
  H1Space<double> space(&mesh, P_INIT);
  int ndof = Space<double>::get_num_dofs(&space);
  info("ndof = %d", ndof);

  // Initialize the right-hand side.
  CustomRightHandSide rhs_value(K);

  // Initialize the weak formulation.
  CustomWeakForm wf(&rhs_value, BDY_LEFT_RIGHT, K);

  Solution<double> sln; 

  // NON-ADAPTIVE VERSION
  
  // Initialize the linear problem.
  DiscreteProblem<double> dp(&wf, &space);

  // Select matrix solver.
  SparseMatrix<double>* matrix = create_matrix<double>(matrix_solver_type);
  Vector<double>* rhs = create_vector<double>(matrix_solver_type);
  LinearSolver<double>* solver = create_linear_solver<double>(matrix_solver_type, matrix, rhs);

  // Assemble stiffness matrix and rhs.
  dp.assemble(matrix, rhs);
 
  // Solve the linear system of the reference problem. If successful, obtain the solutions.
  if(solver->solve()) Solution<double>::vector_to_solution(solver->get_sln_vector(), &space, &sln);
  else error ("Matrix solver failed.\n");

  // Visualize the solution.
  ScalarView view("Solution", new WinGeom(0, 0, 440, 350));
  view.show(&sln);

  // Calculate error wrt. exact solution.
  CustomExactSolution sln_exact(&mesh, K);
  double err = Global<double>::calc_abs_error(&sln, &sln_exact, HERMES_H1_NORM);
  printf("err = %g, err_squared = %g\n\n", err, err*err);
 
  // Wait for all views to be closed.
  View::wait();
  return 0;
}

void MechanicalOperations::m_setSystemMBKMatrix(double mFact, double bFact, double kFact)
{
    LinearSolver* s = ctx->get<LinearSolver>(ctx->getTags(), BaseContext::SearchDown);
    if (!s)
    {
        ctx->serr << "ERROR:  requires a LinearSolver."<<ctx->sendl;
        return;
    }
    mparams.setMFactor(mFact);
    mparams.setBFactor(bFact);
    mparams.setKFactor(kFact);
    s->setSystemMBKMatrix(&mparams);
}

int LinearRelaxFixed::linearization(const IntervalVector& box, LinearSolver& lp_solver)  {
	int num=0;
	for (int i=0; i<A.nb_rows(); i++) {
		try {
			lp_solver.addConstraint(A[i],LEQ,b[i]);
			num++;
		} catch (LPException&) { }
	}

	return num;
}

void MechanicalOperations::m_print( std::ostream& out )
{
    LinearSolver* s = ctx->get<LinearSolver>(ctx->getTags(), BaseContext::SearchDown);
    if (!s)
    {
        ctx->serr << "ERROR: requires a LinearSolver."<<ctx->sendl;
        return;
    }
    defaulttype::BaseMatrix* m = s->getSystemBaseMatrix();
    if (!m) return;
    //out << *m;
    int ny = m->rowSize();
    int nx = m->colSize();
    out << "[";
    for (int y=0; y<ny; ++y)
    {
        out << "[";
        for (int x=0; x<nx; x++)
            out << ' ' << m->element(x,y);
        out << "]";
    }
    out << "]";
}

KeyFrameGraph::KeyFrameGraph()
: nextEdgeId(0)
{
	typedef g2o::BlockSolver_7_3 BlockSolver;
	typedef g2o::LinearSolverCSparse<BlockSolver::PoseMatrixType> LinearSolver;
	//typedef g2o::LinearSolverPCG<BlockSolver::PoseMatrixType> LinearSolver;
	LinearSolver* solver = new LinearSolver();
	BlockSolver* blockSolver = new BlockSolver(solver);
	g2o::OptimizationAlgorithmLevenberg* algorithm = new g2o::OptimizationAlgorithmLevenberg(blockSolver);
	graph.setAlgorithm(algorithm);

    graph.setVerbose(false); // printOptimizationInfo
	solver->setWriteDebug(true);
	blockSolver->setWriteDebug(true);
	algorithm->setWriteDebug(true);


	totalPoints=0;
	totalEdges=0;
	totalVertices=0;


}

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

      if(iteration == 1) {
        if(CAND_LIST_FLOW == H2D_HP_ANISO)
        {
          prev_rho.set_const((*ref_spaces)[4]->get_mesh(), RHO_EXT);
          prev_rho_v_x.set_const((*ref_spaces)[4]->get_mesh(), RHO_EXT * V1_EXT);
          prev_rho_v_y.set_const((*ref_spaces)[4]->get_mesh(), RHO_EXT * V2_EXT);
          prev_e.set_const((*ref_spaces)[4]->get_mesh(), QuantityCalculator::calc_energy(RHO_EXT, RHO_EXT * V1_EXT, RHO_EXT * V2_EXT, P_EXT, KAPPA));
        }
        prev_c.set_const((*ref_spaces)[4]->get_mesh(), 0.0);
      }

      if(as > 1) {
        delete rsln_rho.get_mesh();
        delete rsln_rho_v_x.get_mesh();
        delete rsln_rho_v_y.get_mesh();
        delete rsln_e.get_mesh();
        delete rsln_c.get_mesh();
      }

      // Report NDOFs.
      info("ndof_coarse: %d, ndof_fine: %d.", 
        Space<double>::get_num_dofs(Hermes::vector<Space<double> *>(&space_rho, &space_rho_v_x, 
        &space_rho_v_y, &space_e, &space_c)), Space<double>::get_num_dofs(*ref_spaces));

      // Very imporant, set the meshes for the flow as the same.
      (*ref_spaces)[1]->get_mesh()->set_seq((*ref_spaces)[0]->get_mesh()->get_seq());
      (*ref_spaces)[2]->get_mesh()->set_seq((*ref_spaces)[0]->get_mesh()->get_seq());
      (*ref_spaces)[3]->get_mesh()->set_seq((*ref_spaces)[0]->get_mesh()->get_seq());

      // Set up the solver, matrix, and rhs according to the solver selection.
      SparseMatrix<double>* matrix = create_matrix<double>(matrix_solver_type);
      Vector<double>* rhs = create_vector<double>(matrix_solver_type);
      LinearSolver<double>* solver = create_linear_solver<double>(matrix_solver_type, matrix, rhs);

      // Initialize the FE problem.
      bool is_linear = true;
      DiscreteProblem<double> dp(wf, *ref_spaces);
      if(SEMI_IMPLICIT)
        static_cast<EulerEquationsWeakFormSemiImplicitCoupled*>(wf)->set_time_step(time_step);
      else
        static_cast<EulerEquationsWeakFormExplicitCoupled*>(wf)->set_time_step(time_step);

      // Assemble 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.");
      if (solver->solve())
        Solution<double>::vector_to_solutions(solver->get_sln_vector(), *ref_spaces, 
        Hermes::vector<Solution<double>*>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e, &rsln_c));
      else
        error ("Matrix solver failed.\n");

      Hermes::vector<Space<double>*> flow_spaces((*ref_spaces)[0], (*ref_spaces)[1], (*ref_spaces)[2], (*ref_spaces)[3]);

      double* flow_solution_vector = new double[Space<double>::get_num_dofs(flow_spaces)];

      OGProjection<double>::project_global(flow_spaces, Hermes::vector<MeshFunction<double> *>(&rsln_rho, &rsln_rho_v_x, &rsln_rho_v_y, &rsln_e), flow_solution_vector);

			FluxLimiter flux_limiter(FluxLimiter::Krivodonova, flow_solution_vector, flow_spaces);

			flux_limiter.limit_according_to_detector();

			flux_limiter.get_limited_solutions(Hermes::vector<Solution<double> *>(&prev_rho, &prev_rho_v_x, &prev_rho_v_y, &prev_e));

//.........这里部分代码省略.........
  ParameterXMLFileReader reader("nox.xml");
  ParameterList solverParams = reader.getParameters();

  Out::root() << "finding periodic solution" << endl;
  NOXSolver solver(solverParams);
  prob.solve(solver);

  /* unfold the solution onto a non-periodic mesh */
      
  Expr uP = unfoldPeriodicDiscreteFunction(u0, "u_p");
  Out::root() << "uP=" << uP << endl;
      
  Mesh unfoldedMesh = DiscreteFunction::discFunc(uP)->mesh();
  DiscreteSpace unfDiscSpace = DiscreteFunction::discFunc(uP)->discreteSpace();

  FieldWriter writer = new MatlabWriter("Floquet.dat");
  writer.addMesh(unfoldedMesh);
  writer.addField("u_p[0]", new ExprFieldWrapper(uP[0]));
  writer.addField("u_p[1]", new ExprFieldWrapper(uP[1]));

  Array<Expr> a(2);
  a[0] = new Sundance::Parameter(0.0, "a1");
  a[1] = new Sundance::Parameter(0.0, "a2");


  Expr bc = EssentialBC(left, v1*(u1-uP[0]-a[0]) + v2*(u2-uP[1]-a[1]), quad);

  NonlinearProblem unfProb(unfoldedMesh, eqn, bc, 
    List(v1,v2), List(u1,u2), uP, vecType);

  unfProb.setEvalPoint(uP);

  LinearOperator<double> J = unfProb.allocateJacobian();
  Vector<double> b = J.domain().createMember();

  LinearSolver<double> linSolver
    = LinearSolverBuilder::createSolver("amesos.xml");
        
  SerialDenseMatrix<int, double> F(a.size(), a.size());

  for (int i=0; i<a.size(); i++)
  {
    Out::root() << "doing perturbed orbit #" << i << endl;
    for (int j=0; j<a.size(); j++) 
    {
      if (i==j) a[j].setParameterValue(1.0);
      else a[j].setParameterValue(0.0);
    }
        
    unfProb.computeJacobianAndFunction(J, b);
    Vector<double> w = b.copy();
    linSolver.solve(J, b, w);
    Expr w_i = new DiscreteFunction(unfDiscSpace, w);

    for (int j=0; j<a.size(); j++)
    {
      Out::root() << "postprocessing" << i << endl;

      writer.addField("w[" + Teuchos::toString(i)
        + ", " + Teuchos::toString(j) + "]", new ExprFieldWrapper(w_i[j]));
      Expr g = Integral(right, w_i[j], quad);
      F(j,i) = evaluateIntegral(unfoldedMesh, g);
    }
  }

  writer.write();

  Out::root() << "Floquet matrix = " << endl
              << F << endl;
        

  Out::root() << "doing eigenvalue analysis" << endl;
  Array<double> ew_r(a.size());
  Array<double> ew_i(a.size());
  int lWork = 6*a.size();
  Array<double> work(lWork);
  int info = 0;
  LAPACK<int, double> lapack;
  lapack.GEEV('N','N', a.size(), F.values(),
    a.size(), &(ew_r[0]), &(ew_i[0]), 0, 1, 0, 1, &(work[0]), lWork,
    &info);

  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,
    std::runtime_error,
    "LAPACK GEEV returned error code =" << info);
      
  Array<double> ew(a.size());
  for (int i=0; i<a.size(); i++)
  {
    ew[i] = sqrt(ew_r[i]*ew_r[i]+ew_i[i]*ew_i[i]);
    Out::root() << setw(5) << i 
                << setw(16) << ew_r[i] 
                << setw(16) << ew_i[i] 
                << setw(16) << ew[i]
                << endl;
  }

  double err = ::fabs(ew[0] - 0.123);
  return SundanceGlobal::checkTest(err, 0.001);
}