C++ GenericMatrix类代码示例

//-----------------------------------------------------------------------------
void DirichletBC::zero(GenericMatrix& A) const
{
  // Check arguments
  check_arguments(&A, NULL, NULL, 0);

  // A map to hold the mapping from boundary dofs to boundary values
  Map boundary_values;

  // Create local data for application of boundary conditions
  dolfin_assert(_function_space);
  LocalData data(*_function_space);

  // Compute dofs and values
  compute_bc(boundary_values, data, _method);

  // Copy boundary value data to arrays
  std::vector<dolfin::la_index> dofs(boundary_values.size());
  std::size_t i = 0;
  for (auto bv = boundary_values.begin(); bv != boundary_values.end(); ++bv)
    dofs[i++] = bv->first;

  // Modify linear system (A_ii = 1)
  A.zero_local(boundary_values.size(), dofs.data());

  // Finalise changes to A
  A.apply("insert");
}

GenericMatrix IKSolver::pseudoJacobian(std::vector<Joint*> *bones, int type)
{
    GenericMatrix j = jacobian(bones,type);
    //    return j.transpose();
    if(type!=0){
        return j.transpose();
    }
    GenericMatrix j_jt_inverse = (j*j.transpose());
    //j_jt_inverse.debugPrint("pre-inverse");
    j_jt_inverse = j_jt_inverse.inverse();

    if( j_jt_inverse.cols() == 1 ) return j.transpose();

    //j_jt_inverse.debugPrint("jjtinv");

    //bool nan= false;
    // If transpose only

    // Pseudo Inverse
    // Verify if there is inverse
    //    for(int i=0;i<j_jt_inverse.cols();i++){
    //        for(int j=0;j<j_jt_inverse.rows();j++){
    //            if (isnan(j_jt_inverse.get(j,i))){
    //                nan = true;
    //                break;
    //            }
    //        }
    //    }

    //    if(nan){
    //        return j.transpose();
    //    }

    return (j.transpose()*j_jt_inverse);
}

//-----------------------------------------------------------------------------
CoordinateMatrix::CoordinateMatrix(const GenericMatrix& A, bool symmetric,
                                   bool base_one)
  : _symmetric(symmetric), _base_one(base_one)
{
  _size[0] = A.size(0);
  _size[1] = A.size(1);

  // Iterate over local rows
  const std::pair<std::size_t, std::size_t> local_row_range = A.local_range(0);
  if (!_symmetric)
  {
    for (std::size_t i = local_row_range.first; i < local_row_range.second; ++i)
    {
      // Get column and value data for row
      std::vector<std::size_t> columns;
      std::vector<double> values;
      A.getrow(i, columns, values);

      // Insert data at end
      _rows.insert(_rows.end(), columns.size(), i);
      _cols.insert(_cols.end(), columns.begin(), columns.end());
      _vals.insert(_vals.end(), values.begin(), values.end());
    }
    assert(_rows.size() == _cols.size());
  }
  else
  {
    assert(_size[0] == _size[1]);
    for (std::size_t i = local_row_range.first; i < local_row_range.second; ++i)
    {
      // Get column and value data for row
      std::vector<std::size_t> columns;
      std::vector<double> values;
      A.getrow(i, columns, values);

      for (std::size_t j = 0; j < columns.size(); ++j)
      {
        if (columns[j] >= i)
        {
          _rows.push_back(i);
          _cols.push_back(columns[j]);
          _vals.push_back(values[j]);
        }
      }
    }
    assert(_rows.size() == _cols.size());
  }

  // Add 1 for Fortran-style indices
  if (base_one)
  {
    for (std::size_t i = 0; i < _cols.size(); ++i)
    {
      _rows[i]++;
      _cols[i]++;
    }
  }
}

//---------------------------------------------------------------------------
void EigenMatrix::axpy(double a, const GenericMatrix& A,
                       bool same_nonzero_pattern)
{
  // Check for same size
  if (size(0) != A.size(0) or size(1) != A.size(1))
  {
    dolfin_error("EigenMatrix.cpp",
                 "perform axpy operation with Eigen matrix",
                 "Dimensions don't match");
  }

  _matA += (a)*(as_type<const EigenMatrix>(A).mat());
}

GenericMatrix<T> GenericMatrix<T>::negative() const
{
	GenericMatrix foo = *this;
	unsigned dim = foo.rows() * foo.cols();
	for (unsigned el = 0; el < dim; ++el)
	{
		T bar = foo.data()[ el ];
		if (bar > T( 0 ))
			foo.data()[ el ] = T( 0 );
	}

	return foo;
}

// [[Rcpp::export]]
List matrix_generic( GenericMatrix m){
    List output( m.ncol() ) ;
    for( size_t i=0 ; i<4; i++){
        	output[i] = m(i,i) ;
    }
    return output ;
}

T cross(const GenericMatrix<T> &a, const GenericMatrix<T> &b, int dim) {
  casadi_assert_message(a.size1()==b.size1() && a.size2()==b.size2(),"cross(a,b): Inconsistent dimensions. Dimension of a (" << a.dimString() << " ) must equal that of b (" << b.dimString() << ").");
  
  casadi_assert_message(a.size1()==3 || a.size2()==3,"cross(a,b): One of the dimensions of a should have length 3, but got " << a.dimString() << ".");
  casadi_assert_message(dim==-1 || dim==1 || dim==2,"cross(a,b,dim): Dim must be 1, 2 or -1 (automatic).");
  
  
  std::vector<T> ret(3);
  
  
  bool t = a.size1()==3;
  
  if (dim==1) t = true;
  if (dim==2) t = false;
  
  T a1 = t ? a(0,ALL) : a(ALL,0);
  T a2 = t ? a(1,ALL) : a(ALL,1);
  T a3 = t ? a(2,ALL) : a(ALL,2);

  T b1 = t ? b(0,ALL) : b(ALL,0);
  T b2 = t ? b(1,ALL) : b(ALL,1);
  T b3 = t ? b(2,ALL) : b(ALL,2);
    
  ret[0] = a2*b3-a3*b2;
  ret[1] = a3*b1-a1*b3;
  ret[2] = a1*b2-a2*b1;
    
  return t ? vertcat(ret) : horzcat(ret);
  
}

// [[Rcpp::export]]
IntegerVector runit_GenericMatrix_row( GenericMatrix m ){
    GenericMatrix::Row first_row = m.row(0) ;
    IntegerVector out( first_row.size() ) ;
    std::transform(
    	first_row.begin(), first_row.end(),
    	out.begin(),
    	unary_call<SEXP,int>( Function("length" ) ) ) ;
    return out ;
}

// [[Rcpp::export]]
IntegerVector runit_GenericMatrix_column( GenericMatrix m ){
    GenericMatrix::Column col = m.column(0) ;
    IntegerVector out( col.size() ) ;
    std::transform(
    	   col.begin(), col.end(),
    	   out.begin(),
    	   unary_call<SEXP,int>( Function("length" ) )
    ) ;
    return wrap(out) ;
}

GenericMatrix IKSolver::jacobian(std::vector<Joint*>* bones, int type)
{
    //    Joint * effector =  bones->back();
    Joint * effector =  bones->front();
    Joint * joint = effector;
    //int joint_count=0;
    //    if(!(root==NULL)){
    /*while(joint!=NULL){
        //            if(root->getParent()==NULL){
        //                root->DrawObject();
        //            }
        joint_count++;
        joint = joint->parent();
    }*/
    //    }else{
    //        joint_count=1;
    //    }
    //joint = effector;

    GenericMatrix jacobian = GenericMatrix(3,3*bones->size());

    qglviewer::Vec posEffector = effector->globalEffectorPosition();

    for(unsigned int i = 0 ; i < bones->size() ; i++) {
        joint = bones->at(i);
        qglviewer::Vec derivatex,derivatey,derivatez;
        qglviewer::Vec posJoint = joint->globalPosition();
        qglviewer::Vec posrelative = posEffector - posJoint;

        GenericMatrix globalTransform = joint->globalTransformationMatrix().transpose();

        if(type!=2){
            derivatex.setValue(globalTransform.get(0),globalTransform.get(1),globalTransform.get(2));
            derivatey.setValue(globalTransform.get(4),globalTransform.get(5),globalTransform.get(6));
            derivatez.setValue(globalTransform.get(8),globalTransform.get(9),globalTransform.get(10));

            // Cross Product
            derivatex = derivatex^posrelative;
            derivatey = derivatey^posrelative;
            derivatez = derivatez^posrelative;

            double vetx[3] = {derivatex.x,derivatey.x,derivatez.x};
            double vety[3] = {derivatex.y,derivatey.y,derivatez.y};
            double vetz[3] = {derivatex.z,derivatey.z,derivatez.z};


            for(int k=0;k<3;k++){
                jacobian.set( 0 , k+(3*i) , vetx[k] );
                jacobian.set( 1 , k+(3*i) , vety[k] );
                jacobian.set( 2 , k+(3*i) , vetz[k] );
            }
        }else{
            posrelative = ( effector->globalEffectorPosition() - joint->globalPosition() );
            jacobian.set(0, (3*i),      0);
            jacobian.set(0, 1+(3*i),    posrelative.z);
            jacobian.set(0, 2+(3*i),    posrelative.y);
            jacobian.set(1, (3*i),      -posrelative.z);
            jacobian.set(1, 1+(3*i),    0);
            jacobian.set(1, 2+(3*i),    posrelative.x);
            jacobian.set(2, (3*i),      posrelative.y);
            jacobian.set(2, 1+(3*i),    -posrelative.x);
            jacobian.set(2, 2+(3*i),    0);
        }
    }

    //    jacobian.debugPrint("jacobian");
    return jacobian;
}

int ILSConjugateGradients::solveLin ( const GenericMatrix & gm, const Vector & b, Vector & x )
{
  Timer t;

  if ( timeAnalysis )
    t.start();

  if ( b.size() != gm.rows() ) {
    fthrow(Exception, "Size of vector b (" << b.size() << ") mismatches with the size of the given GenericMatrix (" << gm.rows() << ").");
  }

  if ( x.size() != gm.cols() )
  {
    x.resize(gm.cols());
    x.set(0.0); // bad initial solution, but whatever
  }


  // CG-Method: http://www.netlib.org/templates/templates.pdf
  //

  Vector a;

  // compute r^0 = b - A*x^0
  gm.multiply( a, x ); 
  Vector r = b - a;
  
  if ( timeAnalysis ) {
    t.stop();
    cerr << "r = " << r << endl;
    cerr << "ILSConjugateGradients: TIME " << t.getSum() << " " << r.normL2() << " " << r.normInf() << endl;
    t.start();
  }

  Vector rold ( r.size(), 0.0 );

  // store optimal values 
  double res_min = r.scalarProduct(r);
  Vector current_x = x;
  
  double rhoold = 0.0;
  Vector z ( r.size() );
  Vector p ( z.size() );

  uint i = 1;
  while ( i <= maxIterations )
  {
    // pre-conditioned vector, currently M=I, i.e. no pre-condition
    // otherwise set z = M * r
    if ( jacobiPreconditioner.size() != r.size() )
      z = r;
    else {
      // use simple Jacobi pre-conditioning
      for ( uint jj = 0 ; jj < z.size() ; jj++ )
        z[jj] = r[jj] / jacobiPreconditioner[jj];
    }

    double rho = z.scalarProduct( r );

    if ( verbose ) {
      cerr << "ILSConjugateGradients: iteration " << i << " / " << maxIterations << endl;
      if ( current_x.size() <= 20 )
        cerr << "ILSConjugateGradients: current solution " << current_x << endl;
    }

    if ( i == 1 ) {
      p = z;
    } else {
      double beta;
      if ( useFlexibleVersion ) {
        beta = ( rho - z.scalarProduct(rold) ) / rhoold;
      } else {
        beta = rho / rhoold;
      }
      p = z + beta * p;
    }
    Vector q ( gm.rows() );
    // q = A*p
    gm.multiply ( q, p );

    // sp = p^T A p
    // if A is next to zero this gets nasty, because we divide by sp
    // later on
    double sp = p.scalarProduct(q);
    if ( fabs(sp) < 1e-20 ) {
      // we achieved some kind of convergence, at least this
      // is a termination condition used in the wiki article
      if ( verbose )
        cerr << "ILSConjugateGradients: p^T*q is quite small" << endl;
      break;
    }
    double alpha = rho / sp;
  
    current_x = current_x + alpha * p;

    rold = r;
    r = r - alpha * q;
    
    double res = r.scalarProduct(r);
    double resMax = r.normInf();
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
void DirichletBC::zero_columns(GenericMatrix& A,
                               GenericVector& b,
                               double diag_val) const
{
  // Check arguments
  check_arguments(&A, &b, NULL, 1);

  // A map to hold the mapping from boundary dofs to boundary values
  Map bv_map;
  get_boundary_values(bv_map);

  // Create lookup table of dofs
  //const std::size_t nrows = A.size(0); // should be equal to b.size()
  const std::size_t ncols = A.size(1); // should be equal to max possible dof+1

  std::pair<std::size_t, std::size_t> rows = A.local_range(0);

  std::vector<char> is_bc_dof(ncols);
  std::vector<double> bc_dof_val(ncols);
  for (Map::const_iterator bv = bv_map.begin(); bv != bv_map.end(); ++bv)
  {
    is_bc_dof[bv->first] = 1;
    bc_dof_val[bv->first] = bv->second;
  }

  // Scan through all columns of all rows, setting to zero if
  // is_bc_dof[column]. At the same time, we collect corrections to
  // the RHS

  std::vector<std::size_t> cols;
  std::vector<double> vals;
  std::vector<double> b_vals;
  std::vector<dolfin::la_index> b_rows;

  for (std::size_t row = rows.first; row < rows.second; row++)
  {
    // If diag_val is nonzero, the matrix is a diagonal block
    // (nrows==ncols), and we can set the whole BC row
    if (diag_val != 0.0 && is_bc_dof[row])
    {
      A.getrow(row, cols, vals);
      for (std::size_t j = 0; j < cols.size(); j++)
        vals[j] = (cols[j] == row)*diag_val;
      A.setrow(row, cols, vals);
      A.apply("insert");
      b.setitem(row, bc_dof_val[row]*diag_val);
    }
    else // Otherwise, we scan the row for BC columns
    {
      A.getrow(row, cols, vals);
      bool row_changed = false;
      for (std::size_t j = 0; j < cols.size(); j++)
      {
        const std::size_t col = cols[j];

        // Skip columns that aren't BC, and entries that are zero
        if (!is_bc_dof[col] || vals[j] == 0.0)
          continue;

        // We're going to change the row, so make room for it
        if (!row_changed)
        {
          row_changed = true;
          b_rows.push_back(row);
          b_vals.push_back(0.0);
        }

        b_vals.back() -= bc_dof_val[col]*vals[j];
        vals[j] = 0.0;
      }
      if (row_changed)
      {
        A.setrow(row, cols, vals);
        A.apply("insert");
      }
    }
  }

  b.add_local(&b_vals.front(), b_rows.size(), &b_rows.front());
  b.apply("add");
}

//.........这里部分代码省略.........
      CroutMatrix B1; B1 = B;
      D(2) = determinant(B1) / x - 1;
      C = A * Inverter2(B1) - IdentityMatrix(7);
      Clean(C, 0.000000001); Print(C);
      // Do it again with release
      B.release(); B1 = B;
      D(2) = B.nrows(); D(3) = B.ncols(); D(4) = B.size();
      D(5) = determinant(B1) / x - 1;
      B1.release();
      C = A * Inverter2(B1) - IdentityMatrix(7);
      D(6) = B1.nrows(); D(7) = B1.ncols(); D(8) = B1.size();
      Clean(C, 0.000000001); Print(C);
      // see if we get an implicit invert
      B1 = -A; 
      D(9) = determinant(B1) / x + 1; // odd number of rows - sign will change 
      C = -A * Inverter2(B1) - IdentityMatrix(7);
      Clean(C, 0.000000001); Print(C);
      // check for_return
      B = LU1(A); B1 = LU2(A); CroutMatrix B2 = LU3(A);
      C = A * B.i() - IdentityMatrix(7); Clean(C, 0.000000001); Print(C);
      D(10) = (B == B1 ? 0 : 1) + (B == B2 ? 0 : 1);
      // check lengths
      D(13) = B.size()-49;
      // check release(2)
      B1.release(2);
      B2 = B1; D(15) = B == B2 ? 0 : 1;
      CroutMatrix B3 = B1; D(16) = B == B3 ? 0 : 1;
      D(17) = B1.size();
      // some oddments
      B1 = B; B1 = B1.i(); C = A - B1.i(); Clean(C, 0.000000001); Print(C);
      B1 = B; B1.release(); B1 = B1; B2 = B1;
      D(19) = B == B1 ? 0 : 1; D(20) = B == B2 ? 0 : 1;
      B1.cleanup(); B2 = B1; D(21) = B1.size(); D(22) = B2.size();
      GenericMatrix GM = B; C = A.i() - GM.i(); Clean(C, 0.000000001); Print(C);
      B1 = GM; D(23) = B == B1 ? 0 : 1;
      B1 = A * 0; B2 = B1; D(24) = B2.is_singular() ? 0 : 1;
      // check release again - see if memory moves
      const Real* d = B.const_data();
      const int* i = B.const_data_indx();
      B1 = B;
      const Real* d1 = B1.const_data();
      const int* i1 = B1.const_data_indx();
      B1.release(); B2 = B1;
      const Real* d2 = B2.const_data();
      const int* i2 = B2.const_data_indx();
      D(25) = (d != d1 ? 0 : 1) + (d1 == d2 ? 0 : 1)
         + (i != i1 ? 0 : 1) + (i1 == i2 ? 0 : 1);
  
      Clean(D, 0.000000001); Print(D);
   }

   {
      Tracer et1("Stage 8");
      // Same for BandLUMatrix
      BandMatrix A(7,3,2);
      A.Row(1) <<  3 <<  2 << -1;
      A.Row(2) << -8 <<  7 <<  2 <<  0;
      A.Row(3) <<  2 << -2 <<  3 <<  1 <<  9;
      A.Row(4) << -1 <<  5 <<  2 <<  2 <<  5 << -1;
      A.Row(5)       << -4 <<  1 <<  9 << -8 <<  7 <<  5;
      A.Row(6)             <<  5 << -1 << -2 <<  5 <<  1;
      A.Row(7)                   <<  8 << -1 <<  2 <<  2;
      RowVector D(30); D = 0;
      Real x = determinant(A);
      BandLUMatrix B = A;
      D(1) = determinant(B) / x - 1;

// Base case for all divergence computations. 
// Compute divergence of vector field u.
void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A,
                       const FunctionSpace& DGscalar, 
                       const FunctionSpace& CRvector)
{
  std::shared_ptr<const GenericDofMap>
    CR1_dofmap = CRvector.dofmap(),
    DG0_dofmap = DGscalar.dofmap();

  // Figure out about the local dofs of DG0 
  std::pair<std::size_t, std::size_t>
  first_last_dof = DG0_dofmap->ownership_range();
  std::size_t first_dof = first_last_dof.first;
  std::size_t last_dof = first_last_dof.second;
  std::size_t n_local_dofs = last_dof - first_dof;

  // Get topological dimension so that we know what Facet is
  const Mesh mesh = *DGscalar.mesh();
  std::size_t tdim = mesh.topology().dim(); 
  std::size_t gdim = mesh.geometry().dim();
  
  std::vector<std::size_t> columns;
  std::vector<double> values;
  
  // Fill the values
  for(CellIterator cell(mesh); !cell.end(); ++cell)
  {
    auto dg_dofs = DG0_dofmap->cell_dofs(cell->index());
    // There is only one DG0 dof per cell
    dolfin::la_index cell_dof = dg_dofs[0];
    
    Point cell_mp = cell->midpoint();
    double cell_volume = cell->volume();
    std::size_t local_facet_index = 0;      
    
    auto cr_dofs = CR1_dofmap->cell_dofs(cell->index());
    
    for(FacetIterator facet(*cell); !facet.end(); ++facet)
    {
      double facet_measure=0;
      if(tdim == 2)
        facet_measure = Edge(mesh, facet->index()).length();
      else if(tdim == 3)
        facet_measure = Face(mesh, facet->index()).area();
      // Tdim 1 will not happen because CR is not defined there
      
      Point facet_normal = facet->normal();

      // Flip the normal if it is not outer already
      Point facet_mp = facet->midpoint();
      double sign = (facet_normal.dot(facet_mp - cell_mp) > 0.0) ? 1.0 : -1.0;
      facet_normal *= (sign*facet_measure/cell_volume);
      
      // Dofs of CR on the facet, local order
      std::vector<std::size_t> facet_dofs;
      CR1_dofmap->tabulate_facet_dofs(facet_dofs, local_facet_index);
      
      for (std::size_t j = 0 ; j < facet_dofs.size(); j++)
      {   
        columns.push_back(cr_dofs[facet_dofs[j]]);
        values.push_back(facet_normal[j]);
      }        
      local_facet_index += 1;
    }
    M.setrow(cell_dof, columns, values);
    columns.clear();
    values.clear();
  }
  M.apply("insert");
  //std::shared_ptr<GenericMatrix> Cp = MatMatMult(M, A);
  //return Cp;
}

int ILSSymmLqLanczos::solveLin ( const GenericMatrix & gm, const Vector & b, Vector & x )
{
  if ( b.size() != gm.rows() ) {
    fthrow(Exception, "Size of vector b (" << b.size() << ") mismatches with the size of the given GenericMatrix (" << gm.rows() << ").");
  }

  if ( x.size() != gm.cols() )
  {
    x.resize(gm.cols());
    x.set(0.0); // bad initial solution, but whatever
  }

//   if ( verbose ) cerr << "initial solution: " << x << endl;

  // SYMMLQ-Method based on Lanczos vectors: implementation based on the following:
  //
  // C.C. Paige and M.A. Saunders: "Solution of sparse indefinite systems of linear equations". SIAM Journal on Numerical Analysis, p. 617--629, vol. 12, no. 4, 1975
  // 
  // http://www.netlib.org/templates/templates.pdf
  //
  
  // declare some helpers
  double gamma = 0.0;
  double gamma_bar = 0.0;
  double alpha = 0.0; // alpha_j = v_j^T * A * v_j for new Lanczos vector v_j
  double beta = b.normL2(); // beta_1 = norm(b), in general beta_j = norm(v_j) for new Lanczos vector v_j
  double beta_next = 0.0; // beta_{j+1}
  double c_new = 0.0;
  double c_old = -1.0;
  double s_new = 0.0;
  double s_old = 0.0;
  double z_new = 0.0;
  double z_old = 0.0;
  double z_older = 0.0;
  double delta_new = 0.0;
  double epsilon_next = 0.0;

  // init some helping vectors
  Vector Av(b.size(),0.0); // Av = A * v_j
  Vector Ac(b.size(),0.0); // Ac = A * c_j
  Vector *v_new = new Vector(b.size(),0.0); // new Lanczos vector v_j
  Vector *v_old = 0; // Lanczos vector of the iteration before: v_{j-1}
  Vector *v_next = new Vector(b.size(),0.0); // Lanczos vector of the next iteration: v_{j+1}
  Vector *w_new = new Vector(b.size(),0.0); 
  Vector *w_bar = new Vector(b.size(),0.0); 
  Vector x_L (b.size(),0.0); 
//   Vector x_C (b.size(),0.0); // x_C is a much better approximation than x_L (according to the paper mentioned above)
// NOTE we store x_C in output variable x and only update this solution if the residual decreases (we are able to calculate the residual of x_C without calculating x_C)
  
  // first iteration + initialization, where b will be used as the first Lanczos vector
  *v_new = (1/beta)*b; // init v_1, v_1 = b / norm(b)
  gm.multiply(Av,*v_new); // Av = A * v_1
  alpha = v_new->scalarProduct(Av); // alpha_1 = v_1^T * A * v_1  
  gamma_bar = alpha; // (gamma_bar_1 is equal to alpha_1 in ILSConjugateGradientsLanczos)
  *v_next = Av - (alpha*(*v_new));
  beta_next = v_next->normL2();
  v_next->normalizeL2();
  
  gamma = sqrt( (gamma_bar*gamma_bar) + (beta_next*beta_next) );
  c_new = gamma_bar/gamma;
  s_new = beta_next/gamma;
  
  z_new = beta/gamma;
  
  *w_bar = *v_new;
  
  *w_new = c_new*(*w_bar) + s_new*(*v_next);
  *w_bar = s_new*(*w_bar) - c_new*(*v_next);
  
  x_L = z_new*(*w_new); // first approximation of x
  
  // calculate current residual of x_C
  double res_x_C = (beta*beta)*(s_new*s_new)/(c_new*c_new);
  
  // store minimal residual
  double res_x_C_min = res_x_C;
  
  // store optimal solution x_C in output variable x instead of additional variable x_C
  x = x_L + (z_new/c_new)*(*w_bar); // x_C = x_L + (z_new/c_new)*(*w_bar); 
  
  // calculate delta of x_L
  double delta_x_L = fabs(z_new) * w_new->normL2();
  if ( verbose ) {
    cerr << "ILSSymmLqLanczos: iteration 1 / " << maxIterations << endl;
    if ( x.size() <= 20 )
      cerr << "ILSSymmLqLanczos: current solution x_L: " << x_L << endl;
    cerr << "ILSSymmLqLanczos: delta_x_L = " << delta_x_L << endl;
  }
  
  // start with second iteration
  uint j = 2;
  while (j <= maxIterations )
  {
  
    // prepare next iteration
    if ( v_old == 0 ) v_old = v_new;
    else {
      
      delete v_old;
      v_old = v_new;
//.........这里部分代码省略.........