C++ OoqpVector类代码示例

void SimpleVector::axdzpy( double alpha, OoqpVector& xvec,
			       OoqpVector& zvec, OoqpVector& select )
{
  assert( n == xvec.length() &&
	  n == zvec.length() );

  SimpleVector & sxvec = dynamic_cast<SimpleVector &>(xvec);
  double * x = sxvec.v;
  SimpleVector & szvec = dynamic_cast<SimpleVector &>(zvec);
  double * z = szvec.v;
  SimpleVector & sselect = dynamic_cast<SimpleVector &>(select);
  double * s   = sselect.v;
  int i;
  if( alpha == 1.0 ) {
    for( i = 0; i < n; i++ ) {
      if( 0.0 != s[i] ) v[i] += x[i] / z[i];
    }
  } else if ( alpha == -1.0 ) {
    for( i = 0; i < n; i++ ) {
      if( 0.0 != s[i] ) v[i] -= x[i] / z[i];
    }
  } else {
    for( i = 0; i < n; i++ ) {
      if( 0.0 != s[i] ) v[i] += alpha * x[i] / z[i];
    }
  }
}

void SimpleVector::axzpy( double alpha, OoqpVector& xvec,
			      OoqpVector& zvec )
{
  assert( n == xvec.length() &&
	  n == zvec.length() );

  SimpleVector & sxvec = dynamic_cast<SimpleVector &>(xvec);
  SimpleVector & szvec = dynamic_cast<SimpleVector &>(zvec);

  double * x = sxvec.v;
  double * z = szvec.v;
  double * lx = x + n;
  double * w = v;

  if( alpha == 1.0 ) {
    while( x < lx ) {
      *w += *x * *z;
      w++; x++; z++;
    }
  } else if ( alpha == -1 ) {
    while( x < lx ) {
      *w -= *x * *z;
      w++; x++; z++;
    }
  } else {
    while( x < lx ) {
      *w += alpha * *x * *z;
      w++; x++; z++;
    }
  }
}

void PetscVector::copyFrom( OoqpVector& v_in )
{
    assert( n == v_in.length() );

    if( v_in.isKindOf( kPetscVector ) ) {
        PetscVector & pet = (PetscVector &) v_in;
        int ierr = VecCopy( pet.pv, pv );
        assert( ierr  == 0);
    } else if ( v_in.isKindOf( kSimpleVector ) ) {
        SimpleVector & sv = (SimpleVector &) v_in;
        this->copyFromArray( sv.elements() );
    } else {
        assert( 0 && "Can't copy from unknown OoqpVector type" );
    }
}

double SimpleVector::dotProductWith( OoqpVector& vec )
{
  assert( n == vec.length() );
  SimpleVector & svec = dynamic_cast<SimpleVector &>(vec);
  double * vvec = svec.v;
  
  double dot1 = 0.0;
  double dot2 = 0.0;

  const int size = 8196;
  int kmax       = n / size;
  int i          = 0;

  int k, imax;
  for( k = 0; k < kmax; k++ ) {
    imax = (k + 1) * 8196;
    for( ; i < imax; i++ ) {
      dot1 += v[i] * vvec[i];
    }
    dot2 += dot1;
    dot1  = 0;
  }
  for( ; i < n; i++ ) {
    dot1 += v[i] * vvec[i];
  }
  
  return dot2 + dot1;
}

void SimpleVector::axpy( double alpha, OoqpVector& vec )
{
  assert( n == vec.length() );
  SimpleVector & sv = dynamic_cast<SimpleVector &>(vec);

  int one = 1;
  daxpy_( &n, &alpha, sv.v, &one, v, &one );
}

void SimpleVector::componentMult( OoqpVector& vec )
{
  assert( n == vec.length() );
  SimpleVector & sv = dynamic_cast<SimpleVector &>(vec);
  double * y = sv.v;
  int i;
  for( i = 0; i < n; i++ ) v[i] *= y[i];
}

void SimpleVector::componentDiv ( OoqpVector& vec )
{
  assert( n == vec.length() );
  double * pv = v, *lv = v + n;

  SimpleVector & sv = dynamic_cast<SimpleVector &>(vec);
  double * y = sv.v;

  for( ; pv < lv; pv++, y++ ) *pv /= *y;
}

void SimpleVector::addSomeConstants( double c, OoqpVector& select )
{
  SimpleVector & sselect = dynamic_cast<SimpleVector &>(select);
  double * map = sselect.v;

  int i;
  assert( n == select.length() );
  for( i = 0; i < n; i++ ) {
    if( map[i] ) v[i] += c;
  }
}

void SimpleVector::selectNonZeros( OoqpVector& select )
{
  SimpleVector & sselect = dynamic_cast<SimpleVector &>(select);
  double * map = sselect.v;

  assert( n == select.length() );
  int i;
  for( i = 0; i < n; i++ ) {
    if( 0.0 == map[i] ) v[i] = 0.0;
  }
}

void SimpleVector::divideSome( OoqpVector& div, OoqpVector& select )
{
  if( n == 0 ) return;

  SimpleVector & sselect = dynamic_cast<SimpleVector &>(select);
  double * map = sselect.v;

  SimpleVector & sdiv = dynamic_cast<SimpleVector &>(div);
  double * q   = sdiv.v;
  assert( n == div.length() && n == select.length() );

  double * lmap = map + n;
  double * w = v;
  while( map < lmap ) {
    if( 0 != *map ) {
      *w  /= *q;
    }
    map++;
    w++;
    q++;
  }
}

void DenseStorage::fromGetDiagonal( int idiag, OoqpVector& vec )
{ 
  int k;
  int extent = vec.length();
  
  assert( idiag + extent <= n );
  assert( idiag + extent <= m );

  SimpleVector &  sv = (SimpleVector &) vec;

  for ( k = idiag; k < idiag + extent; k++ ) {
    sv[k] = M[k][k];
  }
}

int SimpleVector::somePositive( OoqpVector& select )
{
  SimpleVector & sselect = dynamic_cast<SimpleVector &>(select);
  double * map = sselect.v;

  assert( n == select.length() );

  int i;
  for( i = 0; i < n; i++ ) {
    if( 0.0 != map[i] && v[i] <= 0 ) {
      cout << "Element " << i << " is nonpositive: " << v[i] << endl;
      return 0;
    }
  }
  return 1;
}

int SimpleVector::matchesNonZeroPattern( OoqpVector& select )
{
  SimpleVector & sselect = dynamic_cast<SimpleVector &>(select);
  double * map = sselect.v;

  double * lmap = map + n;
  assert( n == select.length() );

  double *w = v;
  while( map < lmap ) {
    if( *map == 0.0 && *w != 0.0  ) return 0;
    map++;
    w++;
  }

  return 1;
}

void DenseStorage::setToDiagonal( OoqpVector& vec )
{ 
  int i,k;

  int extent = vec.length();

  assert( extent <= n );
  assert( extent <= m );

  SimpleVector &  sv = (SimpleVector &) vec;
  for( i = 0; i < m; i++ ) {
    for( k = 0; k < n; k++ ) {
      M[i][k] = 0.0;
    }
  }

  for ( k = 0; k < extent; k++ ) {
    M[k][k] = sv[k];
  }
}

double SimpleVector::stepbound(OoqpVector & pvec, double maxStep )
{
  assert( n == pvec.length() );

  SimpleVector & spvec = dynamic_cast<SimpleVector &>(pvec);
  double * p = spvec.v;
  double * w = v;
  double bound = maxStep;

  int i;
  for( i = 0; i < n; i++ ) {
    double temp = p[i];
    if( w[i] >= 0 && temp < 0 ) {
      temp = -w[i]/temp;
      if( temp < bound ) {
	bound = temp;
      }
    }
  }
  return bound;
}