C++ GenericMatrix::cols方法代码示例

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

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

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

int ILSConjugateGradientsLanczos::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;

  // CG-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
  //

  // 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 r(b.size(),0.0); // current residual
  Vector *v_new = new Vector(x.size(),0.0); // new Lanczos vector v_j
  Vector *v_old = new Vector(x.size(),0.0); // Lanczos vector v_{j-1} of the iteration before 
  Vector *v_older = 0; // Lanczos vector v_{j-2} of the iteration before 
  Vector *c_new = new Vector(x.size(),0.0); // current update vector c_j for the solution x
  Vector *c_old = 0; // update vector of iteration before

  // declare some helpers
  double d_new = 0; // current element of diagonal matrix D normally obtained from Cholesky factorization of tridiagonal matrix T, where T consists alpha and beta as below
  double d_old = 0; // corresponding element of the iteration before
  double l_new = 0; // current element of lower unit bidiagonal matrix L normally obtained from Cholesky factorization of tridiagonal matrix T
  double p_new = 0; // current element of vector p, where p is the solution of the modified linear system
  double p_old = 0; // corresponding element of the iteration before
  double alpha = 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
  
  // 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
  d_new=alpha; // d_1 = alpha_1, d_1 is the first element of diagonal matrix D
  p_new = beta/d_new; // p_1 = beta_1 / d_1
  
  *c_new = *v_new; // c_1 = v_1
  Ac = Av; // A*c_1 = A*v_1
  
  // store current solution
  Vector current_x = (p_new*(*c_new)); // first approx. of x: x_1 = p_1 * c_1

  // calculate current residual
  r = b - (p_new*Ac);
  double res = r.scalarProduct(r);
  
  // store minimal residual
  double res_min = res;
  
  // store optimal solution in output variable x
  x = current_x;
  
  double delta_x = fabs(p_new) * c_new->normL2();
  if ( verbose ) {
    cerr << "ILSConjugateGradientsLanczos: iteration 1 / " << maxIterations << endl;
    if ( current_x.size() <= 20 )
      cerr << "ILSConjugateGradientsLanczos: current solution " << current_x << endl;
    cerr << "ILSConjugateGradientsLanczos: delta_x = " << delta_x << endl;
    cerr << "ILSConjugateGradientsLanczos: residual = " << r.scalarProduct(r) << endl;
  }  
  
  // start with second iteration
  uint j = 2;
  while (j <= maxIterations )
  {
    // prepare d and p for next iteration
    d_old = d_new;
    p_old = p_new;
    
    // prepare vectors v_older, v_old, v_new for next iteration
    if ( v_older == 0) v_older = v_old;
    else {
      
      delete v_older;
      v_older = v_old;
    }
    v_old = v_new;
    v_new = new Vector(v_old->size(),0.0);
    
    // prepare vectors c_old, c_new for next iteration
    if ( c_old == 0 ) c_old = c_new;
    else {
      
      delete c_old;
      c_old = c_new;
    }
    c_new = new Vector(c_old->size(),0.0);
    
    //start next iteration:
//.........这里部分代码省略.........