自学教程：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);}

RealMAST::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);}

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

voidMAST::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());}