C++ Complex函数代码示例

Complex Complex::operator + (const Complex &c2) const
{
    return Complex(real+c2.real,image+c2.image);
}

Complex Complex::cc() const {
  return Complex(real(), -imag());
}

int ComplexLUDecompose(pcomplex **a, int n, double *vv, int *indx, double *pd)
//	pcomplex		**a;	the matrix whose LU-decomposition is wanted
//	int			n;		order of a
//	double		*vv;	work vector of size n (stores implicit 
//							scaling of each row)
//	int			*indx;	=> row permutation according to partial 
//							pivoting sequence
//	double		*pd;	=> 1 if number of row interchanges was even, 
//							-1 if odd (NULL OK)
{
	int			i, imax, j, k;
	double		big, dum, temp, d;
	pcomplex		sum, cdum;

	d = 1.0;
	imax = 0; // only to shut the compiler up.

	for (i = 0; i < n; i++) {
		big = 0.0;
		for (j = 0; j < n; j++) {
			if ((temp = Cabs(a[i][j])) > big)
				big = temp;
		}
		if (big == 0.0) {
                    printf("singular matrix in routine ComplexLUDecompose\n");
                    return 1;
		}
		vv[i] = 1.0 / big;
	}
	
	for (j = 0; j < n; j++) {
		for (i = 0; i < j; i++) {
			sum = a[i][j];
			for (k = 0; k < i; k++) 
				sum = Csub(sum, Cmul(a[i][k], a[k][j]));
			a[i][j] = sum;
		}
		big = 0.0;
		for (i = j; i < n; i++) {
			sum = a[i][j];
			for (k = 0; k < j; k++)
				sum = Csub(sum, Cmul(a[i][k], a[k][j]));
			a[i][j] = sum;
			dum = vv[i] * Cabs(sum);
			if (dum >= big) {
				big = dum;
				imax = i;
			}
		}
		if (j != imax) {
			for (k = 0; k < n; k++) {
				cdum = a[imax][k];
				a[imax][k] = a[j][k];
				a[j][k] = cdum;
			}	
			d = -d;
			vv[imax] = vv[j];
		}
		indx[j] = imax;
		if (a[j][j].re == 0.0 && a[j][j].im == 0.0)
			a[j][j] = Complex(1.0e-20, 1.0e-20);
		if (j != n - 1){
			cdum = Cdiv(Complex(1.0, 0.0), a[j][j]);
			for (i = j + 1; i < n; i++)
				a[i][j] = Cmul(a[i][j], cdum);
		}
	}

	if (pd != NULL)
		*pd = d;
	return 0;
}

Complex Complex::operator * (const Complex &c)         
{
	double real=re_ * c.re() - im_ * c.im();
	double image=re_ * c.im() + im_ * c.re();
	return Complex(real, image); 
}

// /
const Complex Complex::operator/(const Complex& rhs) {
	// how to divide?
	double real = (this->realpart * rhs.realpart) - (this->imgpart * rhs.imgpart);
	double img = (this->realpart * rhs.imgpart) + (this->imgpart * rhs.realpart);
	return Complex(real, img);
}

void muxImaginaryZeros(T& fromVec,U& toVec)
{
	toVec.resize(fromVec.size());
	for (size_t i=0;  i!=toVec.size(); i++)
		toVec[i] = Complex(fromVec[i],0);
}

Real
MAST::ComplexAssemblyBase::residual_l2_norm(const libMesh::NumericVector<Real>& real,
                                            const libMesh::NumericVector<Real>& imag) {
    
    START_LOG("complex_solve()", "Residual-L2");

    MAST::NonlinearSystem& nonlin_sys = _system->system();
    
    // iterate over each element, initialize it and get the relevant
    // analysis quantities
    RealVectorX
    sol,
    vec_re;
    RealMatrixX
    mat_re;
    
    ComplexVectorX
    delta_sol,
    vec;
    ComplexMatrixX
    mat;
    
    std::vector<libMesh::dof_id_type> dof_indices;
    const libMesh::DofMap& dof_map = nonlin_sys.get_dof_map();
    
    
    std::unique_ptr<libMesh::NumericVector<Real> >
    residual_re(nonlin_sys.solution->zero_clone().release()),
    residual_im(nonlin_sys.solution->zero_clone().release()),
    localized_base_solution,
    localized_real_solution(build_localized_vector(nonlin_sys, real).release()),
    localized_imag_solution(build_localized_vector(nonlin_sys, imag).release());
    
    
    if (_base_sol)
        localized_base_solution.reset(build_localized_vector(nonlin_sys,
                                                              *_base_sol).release());

    
    // if a solution function is attached, initialize it
    //if (_sol_function)
    //    _sol_function->init( X);
    
    
    libMesh::MeshBase::const_element_iterator       el     =
    nonlin_sys.get_mesh().active_local_elements_begin();
    const libMesh::MeshBase::const_element_iterator end_el =
    nonlin_sys.get_mesh().active_local_elements_end();
    
    MAST::ComplexAssemblyElemOperations&
    ops = dynamic_cast<MAST::ComplexAssemblyElemOperations&>(*_elem_ops);
    
    for ( ; el != end_el; ++el) {
        
        const libMesh::Elem* elem = *el;
        
        dof_map.dof_indices (elem, dof_indices);
        
        ops.init(*elem);
        
        // get the solution
        unsigned int ndofs = (unsigned int)dof_indices.size();
        sol.setZero(ndofs);
        delta_sol.setZero(ndofs);
        vec.setZero(ndofs);
        mat.setZero(ndofs, ndofs);
        
        // set zero velocity for base solution
        ops.set_elem_velocity(sol);
        
        // set the value of the base solution, if provided
        if (_base_sol)
            for (unsigned int i=0; i<dof_indices.size(); i++)
                sol(i) = (*localized_base_solution)(dof_indices[i]);
        ops.set_elem_solution(sol);
        
        // set the value of the small-disturbance solution
        for (unsigned int i=0; i<dof_indices.size(); i++)
            delta_sol(i) = Complex((*localized_real_solution)(dof_indices[i]),
                                   (*localized_imag_solution)(dof_indices[i]));
        
        ops.set_elem_complex_solution(delta_sol);
        
        
//        if (_sol_function)
//            ops.attach_active_solution_function(*_sol_function);
        
        
        // perform the element level calculations
        ops.elem_calculations(false,
                                             vec, mat);
        ops.clear_elem();
        
//        ops.detach_active_solution_function();
        
        // add to the real part of the residual
        vec_re  =  vec.real();
        DenseRealVector v;
        MAST::copy(v, vec_re);
        dof_map.constrain_element_vector(v, dof_indices);
//.........这里部分代码省略.........

void Biquad::setAllpassPole(const Complex& pole) {
  Complex zero = Complex(1, 0) / pole;
  setZeroPolePairs(zero, pole);
}

void
MAST::ComplexAssemblyBase::
residual_and_jacobian_field_split (const libMesh::NumericVector<Real>& X_R,
                                   const libMesh::NumericVector<Real>& X_I,
                                   libMesh::NumericVector<Real>& R_R,
                                   libMesh::NumericVector<Real>& R_I,
                                   libMesh::SparseMatrix<Real>&  J_R,
                                   libMesh::SparseMatrix<Real>&  J_I) {

    libmesh_assert(_system);
    libmesh_assert(_discipline);
    libmesh_assert(_elem_ops);

    MAST::NonlinearSystem& nonlin_sys = _system->system();
    
    R_R.zero();
    R_I.zero();
    J_R.zero();
    J_I.zero();
    
    // iterate over each element, initialize it and get the relevant
    // analysis quantities
    RealVectorX    sol, vec_re;
    RealMatrixX    mat_re;
    ComplexVectorX delta_sol, vec;
    ComplexMatrixX mat;
    
    std::vector<libMesh::dof_id_type> dof_indices;
    const libMesh::DofMap& dof_map = _system->system().get_dof_map();
    
    
    std::unique_ptr<libMesh::NumericVector<Real> >
    localized_base_solution,
    localized_real_solution,
    localized_imag_solution;
    
    
    // localize the base solution, if it was provided
    if (_base_sol)
        localized_base_solution.reset(build_localized_vector(nonlin_sys,
                                                              *_base_sol).release());
    
    
    // localize sol to real vector
    localized_real_solution.reset(build_localized_vector(nonlin_sys,
                                                              X_R).release());
    // localize sol to imag vector
    localized_imag_solution.reset(build_localized_vector(nonlin_sys,
                                                              X_I).release());
    
    
    // if a solution function is attached, initialize it
    //if (_sol_function)
    //    _sol_function->init( X);
    
    
    libMesh::MeshBase::const_element_iterator       el     =
    nonlin_sys.get_mesh().active_local_elements_begin();
    const libMesh::MeshBase::const_element_iterator end_el =
    nonlin_sys.get_mesh().active_local_elements_end();
    
    MAST::ComplexAssemblyElemOperations&
    ops = dynamic_cast<MAST::ComplexAssemblyElemOperations&>(*_elem_ops);

    for ( ; el != end_el; ++el) {
        
        const libMesh::Elem* elem = *el;
        
        dof_map.dof_indices (elem, dof_indices);
        
        ops.init(*elem);
        
        // get the solution
        unsigned int ndofs = (unsigned int)dof_indices.size();
        sol.setZero(ndofs);
        delta_sol.setZero(ndofs);
        vec.setZero(ndofs);
        mat.setZero(ndofs, ndofs);
        
        // first set the velocity to be zero
        ops.set_elem_velocity(sol);
        
        // next, set the base solution, if provided
        if (_base_sol)
            for (unsigned int i=0; i<dof_indices.size(); i++)
                sol(i) = (*localized_base_solution)(dof_indices[i]);
        
        ops.set_elem_solution(sol);
        
        // set the value of the small-disturbance solution
        for (unsigned int i=0; i<dof_indices.size(); i++)
            delta_sol(i) = Complex((*localized_real_solution)(dof_indices[i]),
                                   (*localized_imag_solution)(dof_indices[i]));
        
        ops.set_elem_complex_solution(delta_sol);
        
        
//        if (_sol_function)
//            physics_elem->attach_active_solution_function(*_sol_function);
        
//.........这里部分代码省略.........

void
MAST::ComplexAssemblyBase::
residual_and_jacobian_blocked (const libMesh::NumericVector<Real>& X,
                               libMesh::NumericVector<Real>& R,
                               libMesh::SparseMatrix<Real>&  J,
                               MAST::Parameter* p) {

    libmesh_assert(_system);
    libmesh_assert(_discipline);
    libmesh_assert(_elem_ops);

    START_LOG("residual_and_jacobian()", "ComplexSolve");
    
    MAST::NonlinearSystem& nonlin_sys = _system->system();
    
    R.zero();
    J.zero();
    
    // iterate over each element, initialize it and get the relevant
    // analysis quantities
    RealVectorX
    sol;
    ComplexVectorX
    delta_sol,
    vec;
    ComplexMatrixX
    mat,
    dummy;

    // get the petsc vector and matrix objects
    Mat
    jac_bmat = dynamic_cast<libMesh::PetscMatrix<Real>&>(J).mat();
    
    PetscInt ierr;
    
    std::vector<libMesh::dof_id_type> dof_indices;
    const libMesh::DofMap& dof_map = nonlin_sys.get_dof_map();
    const std::vector<libMesh::dof_id_type>&
    send_list = nonlin_sys.get_dof_map().get_send_list();
    
    
    
    std::unique_ptr<libMesh::NumericVector<Real> >
    localized_base_solution,
    localized_complex_sol(libMesh::NumericVector<Real>::build(nonlin_sys.comm()).release());
    
    // prepare a send list for localization of the complex solution
    std::vector<libMesh::dof_id_type>
    complex_send_list(2*send_list.size());
    
    for (unsigned int i=0; i<send_list.size(); i++) {
        complex_send_list[2*i  ] = 2*send_list[i];
        complex_send_list[2*i+1] = 2*send_list[i]+1;
    }

    localized_complex_sol->init(2*nonlin_sys.n_dofs(),
                                2*nonlin_sys.n_local_dofs(),
                                complex_send_list,
                                false,
                                libMesh::GHOSTED);
    X.localize(*localized_complex_sol, complex_send_list);
    
    // localize the base solution, if it was provided
    if (_base_sol)
        localized_base_solution.reset(build_localized_vector(nonlin_sys,
                                                             *_base_sol).release());
    
    

    // if a solution function is attached, initialize it
    //if (_sol_function)
    //    _sol_function->init( X);
    
    
    libMesh::MeshBase::const_element_iterator       el     =
    nonlin_sys.get_mesh().active_local_elements_begin();
    const libMesh::MeshBase::const_element_iterator end_el =
    nonlin_sys.get_mesh().active_local_elements_end();
    
    MAST::ComplexAssemblyElemOperations&
    ops = dynamic_cast<MAST::ComplexAssemblyElemOperations&>(*_elem_ops);

    for ( ; el != end_el; ++el) {
        
        const libMesh::Elem* elem = *el;
        
        dof_map.dof_indices (elem, dof_indices);
        
        ops.init(*elem);
        
        // get the solution
        unsigned int ndofs = (unsigned int)dof_indices.size();
        sol.setZero(ndofs);
        delta_sol.setZero(ndofs);
        vec.setZero(ndofs);
        mat.setZero(ndofs, ndofs);
        
        // first set the velocity to be zero
        ops.set_elem_velocity(sol);
        
//.........这里部分代码省略.........

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

  try {
    TCLAP::CmdLine cmd("CEVAL Algorithm", ' ', "0.1");
    cmd.add(poly);
    cmd.add(x_min);
    cmd.add(x_max);
    cmd.add(y_min);
    cmd.add(y_max);
    cmd.add(use_rb);
    cmd.add(display);
    cmd.add(no_use_inclusion);
    cmd.add(min_box_size);
    cmd.add(max_box_size);
    cmd.add(random_poly);
    cmd.add(rand_degree);
    cmd.add(rand_seed);

    // Parse the args.
    cmd.parse(argc, argv);

    // Get the value parsed by each arg.
    //string name = nameArg.getValue();
  } catch (TCLAP::ArgException &e) {
    cerr << "Error : " << e.error() << endl;
    cerr << "Processing arg : " << e.argId() << endl;
    return -1;
  }

  Polynomial<PolyType> a;
  if (random_poly.getValue()) {
    benchmark::GenerateRandom(&a, rand_degree.getValue(), 10 ,rand_seed.getValue());
  } else {
    benchmark::GetPoly(poly.getValue().c_str(), &a);
  }

  double max_box_size_d(max_box_size.getValue());
  double min_box_size_d(min_box_size.getValue());

  Complex min, max;
  if (use_rb.getValue()) {
    double min_r, max_r;
    benchmark::GetRootBounds(a, &min_r, &max_r);
    min = Complex(min_r, min_r);
    max = Complex(max_r, max_r);
    if (display.getValue()) {
      cout << "Min : " << min << ",Max : " << max << endl;
    }
  } else {
    double x_min_d(x_min.getValue()), x_max_d(x_max.getValue());
    double y_min_d(y_min.getValue()), y_max_d(y_max.getValue());
    min = Complex(x_min_d, y_min_d);
    max = Complex(x_max_d, y_max_d);
  }

  Box *b = new Box(min, max);
  Box *b_copy = new Box(min, max);

  // start timing code.
  struct timeval start;
  struct timeval end;
  gettimeofday(&start, NULL);
  // end timing code.

  Predicates p(a);
  Algorithm algo(p, b, min_box_size_d, !no_use_inclusion.getValue(), display.getValue());
  algo.Run();

  if (no_use_inclusion.getValue()) {
    algo.AttemptIsolation();
  }

  // start timing code.
  gettimeofday(&end, NULL);
  cout << ",time=" <<
      (end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec);


  if (display.getValue()) {
    cout << "Operated on Bounding box : " << min << "," << max << endl;
    cout << "With polynomial : " << endl;

    //a.dump();
    cout << endl;

    cout << "--------------------------" << endl;
    cout << "Number of roots:" << algo.output()->size() << endl;
    list<const Disk *>::const_iterator it = algo.output()->begin();
    while (it != algo.output()->end()) {
      cout << "m= " << (*it)->centre << ", r= " << (*it)->radius << endl;
      ++it;
    }

    display_funcs::SetDisplayParams(b_copy, algo.reject(), algo.output_boxes(),
        algo.ambiguous());

    startGlutLoop(argc, argv);
  } else {
    cout << ",output=" << algo.output()->size() << endl;
  }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		  olds = news;
		  in_file >> news;

		  /*
		   * No good style, i know...
		   */
		  goto go_and_find_the_next_dataset;
		}

	    }

	  /*
	   * Check the location of the dataset.
	   */
	  if (dataset_location != 1)
	    {
	      libMesh::err << "ERROR: Currently only Data at nodes is supported."
			    << std::endl;
	      libmesh_error();
	    }


	  /*
	   * Now get the foreign node id number and the respective nodal data.
	   */
	  int f_n_id;
	  std::vector<Number> values;

	  while(true)
	    {
	      in_file >> f_n_id;

	      /*
	       * if node_nr = -1 then we have reached the end of the dataset.
	       */
	      if (f_n_id==-1)
		  break;

	      /*
	       * Resize the values vector (usually data in three
	       * principle directions, i.e. NVALDC = 3).
	       */
	      values.resize(NVALDC);

	      /*
	       * Read the meshdata for the respective node.
	       */
	      for (unsigned int data_cnt=0; data_cnt<NVALDC; data_cnt++)
		{
		  /*
		   * Check what data type we are reading.
		   * 2,4: Real
		   * 5,6: Complex
		   * other data types are not supported yet.
		   * As again, these floats may also be written
		   * using a 'D' instead of an 'e'.
		   */
		  if (data_type == 2 || data_type == 4)
		    {
		      std::string buf;
		      in_file >> buf;
		      MeshDataUnvHeader::need_D_to_e(buf);
#ifdef LIBMESH_USE_COMPLEX_NUMBERS
		      values[data_cnt] = Complex(std::atof(buf.c_str()), 0.);
#else
		      values[data_cnt] = std::atof(buf.c_str());
#endif
		    }

		  else if(data_type == 5 || data_type == 6)

		    {
#ifdef LIBMESH_USE_COMPLEX_NUMBERS
		      Real re_val, im_val;

		      std::string buf;
		      in_file >> buf;

		      if (MeshDataUnvHeader::need_D_to_e(buf))
		        {
			  re_val = std::atof(buf.c_str());
			  in_file >> buf;
			  MeshDataUnvHeader::need_D_to_e(buf);
			  im_val = std::atof(buf.c_str());
			}
		      else
		        {
			  re_val = std::atof(buf.c_str());
			  in_file >> im_val;
			}

		      values[data_cnt] = Complex(re_val,im_val);
#else

		      libMesh::err << "ERROR: Complex data only supported" << std::endl
				    << "when libMesh is configured with --enable-complex!"
				    << std::endl;
		      libmesh_error();
#endif
		    }

  void Arnoldi<SCAL>::Calc (int numval, Array<Complex> & lam, int numev, 
                            Array<shared_ptr<BaseVector>> & hevecs, 
                            const BaseMatrix * pre) const
  { 
    static Timer t("arnoldi");    
    static Timer t2("arnoldi - orthogonalize");    
    static Timer t3("arnoldi - compute large vectors");

    RegionTimer reg(t);

    auto hv  = a.CreateVector();
    auto hv2 = a.CreateVector();
    auto hva = a.CreateVector();
    auto hvm = a.CreateVector();
   
    int n = hv.FV<SCAL>().Size();    
    int m = min2 (numval, n);


    Matrix<SCAL> matH(m);
    Array<shared_ptr<BaseVector>> abv(m);
    for (int i = 0; i < m; i++)
      abv[i] = a.CreateVector();

    auto mat_shift = a.CreateMatrix();
    mat_shift->AsVector() = a.AsVector() - shift*b.AsVector();  
    shared_ptr<BaseMatrix> inv;
    if (!pre)
      inv = mat_shift->InverseMatrix (freedofs);
    else
      {
        auto itso = make_shared<GMRESSolver<double>> (*mat_shift, *pre);
        itso->SetPrintRates(1);
        itso->SetMaxSteps(2000);
        inv = itso;
      }

    hv.SetRandom();
    hv.SetParallelStatus (CUMULATED);
    FlatVector<SCAL> fv = hv.FV<SCAL>();
    if (freedofs)
      for (int i = 0; i < hv.Size(); i++)
	if (! (*freedofs)[i] ) fv(i) = 0;

    t2.Start();
    // matV = SCAL(0.0);   why ?
    matH = SCAL(0.0);

    *hv2 = *hv;
    SCAL len = sqrt (S_InnerProduct<SCAL> (*hv, *hv2)); // parallel
    *hv /= len;
    
    for (int i = 0; i < m; i++)
      {
	cout << IM(1) << "\ri = " << i << "/" << m << flush;
	/*
	for (int j = 0; j < n; j++)
	  matV(i,j) = hv.FV<SCAL>()(j);
	*/
	*abv[i] = *hv;

	*hva = b * *hv;
	*hvm = *inv * *hva;

	for (int j = 0; j <= i; j++)
	  {
            /*
            SCAL sum = 0.0;
	    for (int k = 0; k < n; k++)
	      sum += hvm.FV<SCAL>()(k) * matV(j,k);
	    matH(j,i) = sum;
	    for (int k = 0; k < n; k++)
	      hvm.FV<SCAL>()(k) -= sum * matV(j,k);
            */
            /*
            SCAL sum = 0.0;
            FlatVector<SCAL> abvj = abv[j] -> FV<SCAL>();
            FlatVector<SCAL> fv_hvm = hvm.FV<SCAL>();
	    for (int k = 0; k < n; k++)
	      sum += fv_hvm(k) * abvj(k);
	    matH(j,i) = sum;
	    for (int k = 0; k < n; k++)
	      fv_hvm(k) -= sum * abvj(k);
            */

	    matH(j,i) = S_InnerProduct<SCAL> (*hvm, *abv[j]);
	    *hvm -= matH(j,i) * *abv[j];
	  }
		
	*hv = *hvm;
	*hv2 = *hv;
	SCAL len = sqrt (S_InnerProduct<SCAL> (*hv, *hv2));
	if (i<m-1) matH(i+1,i) = len; 
	
	*hv /= len;
      }
      
    t2.Stop();
    t2.AddFlops (double(n)*m*m);
    cout << "n = " << n << ", m = " << m << " n*m*m = " << n*m*m << endl;
//.........这里部分代码省略.........

Complex operator+(double d, const Complex & c1)
{
	return Complex(c1.getReal() + d, c1.getImagionary());
}