C++ VectorMatrix类代码示例

VectorMatrixCU32Accessor::VectorMatrixCU32Accessor(VectorMatrix &mat) : mat(mat)
{
    mat.writeLock(1); // 1 = CU32Device!
    data_x = static_cast<CU32Array*>(mat.getArray(1, 0))->ptr();
    data_y = static_cast<CU32Array*>(mat.getArray(1, 1))->ptr();
    data_z = static_cast<CU32Array*>(mat.getArray(1, 2))->ptr();
}

void gradient_cpu(double delta_x, double delta_y, double delta_z, const double *phi, VectorMatrix &field)
{
	const int dim_x = field.dimX();
	const int dim_y = field.dimY();
	const int dim_z = field.dimZ();

	const int phi_sx = 1;
	const int phi_sy = (dim_x+1);
	const int phi_sz = (dim_x+1)*(dim_y+1);

	VectorMatrix::accessor field_acc(field);

	for (int z=0; z<dim_z; ++z)
	for (int y=0; y<dim_y; ++y)
	for (int x=0; x<dim_x; ++x) {
		const int i = phi_sx*x + phi_sy*y + phi_sz*z;

		const double dx = (+ phi[i+1*phi_sx+0*phi_sy+0*phi_sz] - phi[i+0*phi_sx+0*phi_sy+0*phi_sz]
		                   + phi[i+1*phi_sx+1*phi_sy+0*phi_sz] - phi[i+0*phi_sx+1*phi_sy+0*phi_sz]
		                   + phi[i+1*phi_sx+0*phi_sy+1*phi_sz] - phi[i+0*phi_sx+0*phi_sy+1*phi_sz]
		                   + phi[i+1*phi_sx+1*phi_sy+1*phi_sz] - phi[i+0*phi_sx+1*phi_sy+1*phi_sz]) / (4.0 * delta_x);

		const double dy = (+ phi[i+0*phi_sx+1*phi_sy+0*phi_sz] - phi[i+0*phi_sx+0*phi_sy+0*phi_sz]
		                   + phi[i+1*phi_sx+1*phi_sy+0*phi_sz] - phi[i+1*phi_sx+0*phi_sy+0*phi_sz]
		                   + phi[i+0*phi_sx+1*phi_sy+1*phi_sz] - phi[i+0*phi_sx+0*phi_sy+1*phi_sz]
		                   + phi[i+1*phi_sx+1*phi_sy+1*phi_sz] - phi[i+1*phi_sx+0*phi_sy+1*phi_sz]) / (4.0 * delta_y);

		const double dz = (+ phi[i+0*phi_sx+0*phi_sy+1*phi_sz] - phi[i+0*phi_sx+0*phi_sy+0*phi_sz]
		                   + phi[i+1*phi_sx+0*phi_sy+1*phi_sz] - phi[i+1*phi_sx+0*phi_sy+0*phi_sz]
		                   + phi[i+0*phi_sx+1*phi_sy+1*phi_sz] - phi[i+0*phi_sx+1*phi_sy+0*phi_sz]
		                   + phi[i+1*phi_sx+1*phi_sy+1*phi_sz] - phi[i+1*phi_sx+1*phi_sy+0*phi_sz]) / (4.0 * delta_z);

		field_acc.set(x, y, z, Vector3d(dx, dy, dz));
	}
}

VectorVectorConvolution_Simple::VectorVectorConvolution_Simple(const VectorMatrix &lhs, int dim_x, int dim_y, int dim_z)
	: lhs(lhs), dim_x(dim_x), dim_y(dim_y), dim_z(dim_z)
{
	assert(lhs.getShape().getDim(0) == 3);
	exp_x = lhs.getShape().getDim(1);
	exp_y = lhs.getShape().getDim(2);
	exp_z = lhs.getShape().getDim(3);
}

 // must pass in a lvalue
 void copy_to( VectorMatrix<T>& vector_matrix_obj ) {
   vector_matrix_obj.resize( this->nrow_, this->ncol_ );
   const size_t length = this->nrow_ * this->ncol_;
   std :: vector<T> temp_store;
   temp_store.resize(length);
   std :: copy_n( this->store_.begin(), length, temp_store.begin() );
   vector_matrix_obj.set_store() = std :: move( temp_store );
 }

double rk_scaled_error_norm_cpu(double h, double eps_abs, double eps_rel, const VectorMatrix &y, const VectorMatrix &y_error)
{
	double norm;
	
	if (false) { // Maximum norm
		double max = 0.0;

		VectorMatrix::const_accessor y_acc(y), y_error_acc(y_error);

		const size_t s = y.size() * 3;
		for (size_t i=0; i<s; ++i) {
			const Vector3d value =       y_acc.get(i);
			const Vector3d error = y_error_acc.get(i);
			{
				const double D0 = eps_rel * std::fabs(value.x) + eps_abs;
				const double r  = std::fabs(error.x) / std::fabs(D0); // scaled error at equation i
				if (r > max) max = r;
			} {
				const double D0 = eps_rel * std::fabs(value.y) + eps_abs;
				const double r  = std::fabs(error.y) / std::fabs(D0); // scaled error at equation i
				if (r > max) max = r;
			} {
				const double D0 = eps_rel * std::fabs(value.z) + eps_abs;
				const double r  = std::fabs(error.z) / std::fabs(D0); // scaled error at equation i
				if (r > max) max = r;
			}
		}
		norm = max;

	} else { // Euclidian norm
		double sum = 0.0;

		VectorMatrix::const_accessor y_acc(y), y_error_acc(y_error);
		const size_t s = y.size();
		for (size_t i=0; i<s; ++i) {
			const Vector3d value =       y_acc.get(i);
			const Vector3d error = y_error_acc.get(i);

			{
				const double D0 = eps_rel * std::fabs(value.x) + eps_abs;
				const double r  = std::fabs(error.x) / std::fabs(D0); // scaled error at equation i
				sum += r*r;
			} {
				const double D0 = eps_rel * std::fabs(value.y) + eps_abs;
				const double r  = std::fabs(error.y) / std::fabs(D0); // scaled error at equation i
				sum += r*r;
			} {
				const double D0 = eps_rel * std::fabs(value.z) + eps_abs;
				const double r  = std::fabs(error.z) / std::fabs(D0); // scaled error at equation i
				sum += r*r;
			}
		}
		norm = std::sqrt(sum / (s*3));
	}

	return norm;
}

Vector3d findExtremum_cpu(VectorMatrix &M, int z_slice, int component)
{
	if (M.getShape().getRank() != 3) {
		throw std::runtime_error("findExtremum: Fixme: Need matrix of rank 3");
	}

	if (component < 0 || component > 2) {
		throw std::runtime_error("findExtremum: Invalid 'component' value, must be 0, 1 or 2.");
	}

	const int dim_x = M.getShape().getDim(0);
	const int dim_y = M.getShape().getDim(1);

	VectorMatrix::const_accessor M_acc(M);

	// Find cell with maximum absolute value
	double max_val = -1.0;
	int max_x = -1, max_y = -1;
	for (int y=1; y<dim_y-1; ++y)
	for (int x=1; x<dim_x-1; ++x) {
		const int val = std::fabs(M_acc.get(x, y, z_slice)[component]);
		if (val > max_val) {
			max_val = val;
			max_x = x;
			max_y = y;
		}
	}
	assert(max_x > 0);
	assert(max_y > 0);
	
	// Refine maximum by fitting to sub-cell precision
	const double xdir_vals[3] = {
		M_acc.get(max_x-1, max_y+0, z_slice)[component],
		M_acc.get(max_x+0, max_y+0, z_slice)[component],
		M_acc.get(max_x+1, max_y+0, z_slice)[component]
	};

	const double ydir_vals[3] = {
		M_acc.get(max_x+0, max_y-1, z_slice)[component],
		M_acc.get(max_x+0, max_y+0, z_slice)[component],
		M_acc.get(max_x+0, max_y+1, z_slice)[component]
	};

	return Vector3d(
		fit(max_x-1, max_x+0, max_x+1, xdir_vals[0], xdir_vals[1], xdir_vals[2]),
		fit(max_y-1, max_y+0, max_y+1, ydir_vals[0], ydir_vals[1], ydir_vals[2]),
		static_cast<double>(z_slice)
	);
}

void VectorMatrix::add(const VectorMatrix &op, double factor)
{
	if (this == &op) {
		scale(1.0 + factor);
	} else if (isUniform() && op.isUniform()) {
		fill(getUniformValue() + op.getUniformValue() * factor);
	} else {
		const int dev = computeStrategy2(op);
		writeLock(dev); op.readLock(dev);
		for (int c=0; c<num_arrays; ++c) {
			matty::getDevice(dev)->add(getArray(dev, c), op.getArray(dev, c), factor);
		}
		writeUnlock(dev); op.readUnlock(dev);
	}
}

void VectorMatrix::assign(const VectorMatrix &op)
{
	if (this == &op) {
		return;
	} else if (op.isUniform()) {
		fill(op.getUniformValue());
	} else {
		const int dev = computeStrategy2(op);
		writeLock(dev); op.readLock(dev);
		for (int c=0; c<num_arrays; ++c) {
			matty::getDevice(dev)->assign(getArray(dev, c), op.getArray(dev, c));
		}
		writeUnlock(dev); op.readUnlock(dev);
	}
}

void rk_combine_result_cpu(
	double h, ButcherTableau &tab,
	const VectorMatrix &k0, const VectorMatrix &k1, const VectorMatrix &k2, const VectorMatrix &k3,
	VectorMatrix &y, VectorMatrix &y_error)
{
	const size_t s = y.size();
	VectorMatrix::accessor y_acc(y), y_error_acc(y_error);

	// tab
	const double  c0 = tab. c[0],  c1 = tab. c[1],  c2 = tab. c[2],  c3 = tab. c[3];
	const double ec0 = tab.ec[0], ec1 = tab.ec[1], ec2 = tab.ec[2], ec3 = tab.ec[3];

	VectorMatrix::const_accessor k0_acc(k0), k1_acc(k1), k2_acc(k2), k3_acc(k3);
	for (size_t i=0; i<s; ++i) {
		const Vector3d y_i = y_acc.get(i);
		y_acc.set(i, 
			y_i + h * (  c0*k0_acc.get(i) 
			           + c1*k1_acc.get(i)
				   + c2*k2_acc.get(i)
				   + c3*k3_acc.get(i))
		);
		y_error_acc.set(i,
			h * (  ec0*k0_acc.get(i) 
			     + ec1*k1_acc.get(i)
			     + ec2*k2_acc.get(i)
			     + ec3*k3_acc.get(i))
		);
	}
}

 void copy_from( const VectorMatrix<T>& vector_matrix_obj ) {
   try {
     const size_t obj_size = vector_matrix_obj.size();
     if( ( obj_size < STACK_DOUBLE_LIMIT ) == false ) {
       throw obj_size;
     }
     this->nrow_ = vector_matrix_obj.nrow();
     this->ncol_ = vector_matrix_obj.ncol();
     std :: vector<T> temp_store = vector_matrix_obj.store();
     std :: copy_n( temp_store.begin(), obj_size, this->store_.begin() );
   }
   catch( size_t n ) {
     std :: cout << " array load exception: " << n << " >= " << STACK_DOUBLE_LIMIT << std :: endl;
     abort();
   }
 }

/*
  @method DistanceList
  @discussion Allocate storage for an n-unit unitList.
  @throws KII_bad_alloc
*/
DistanceList::DistanceList(const FLOAT xlocs[], const FLOAT ylocs[], 
			   const VectorMatrix& radii) 
  : unitList(NULL), numUnits(radii.Size()), numOverlap(0)
{
  // Since this is a symmetric array, there's no need to store
  // information for unit n (all that information is already in the
  // lists for units 1..n-1
  if ((unitList = new ListItem[numUnits-1]) == NULL)
    throw KII_bad_alloc("Failed allocating memory for DistanceList.");

  // Fill the lists. For unit i, we need only store information for
  // connections to units i+1..n, as earlier units already have the
  // information for their connections with i.
  FLOAT tempDist, tempDist2, tempDelta, deltaX, deltaY;
  for (int u1=0; u1<numUnits-1; u1++) {
    unitList[u1].radius = radii[u1];
    for (int u2=u1+1; u2<numUnits; u2++) {
      deltaX = xlocs[u1]-xlocs[u2];
      deltaY = ylocs[u1]-ylocs[u2];
      tempDist2 = deltaX*deltaX + deltaY*deltaY;
      tempDist = sqrt(tempDist2);
      tempDelta = tempDist - (radii[u1]+radii[u2]);
      if (tempDelta < 0.0) {
	unitList[u1].overlapping.push_back(SublistItem(u2, radii[u2], tempDist,
						       tempDist2, tempDelta));
	numOverlap++;
      } else {
	unitList[u1].nonOverlapping.push_back(SublistItem(u2, radii[u2],
							  tempDist, tempDist2, tempDelta));
      }
    }
  }
}

void rk_combine_result(
	double h, ButcherTableau &tab,
	const VectorMatrix &k0, const VectorMatrix &k1, const VectorMatrix &k2,
	const VectorMatrix &k3, const VectorMatrix &k4, const VectorMatrix &k5,
	VectorMatrix &y, VectorMatrix &y_error)
{
	const int s = y.size();
	if (   s != y_error.size()
	    || s != k1.size()
	    || s != k2.size()
	    || s != k3.size()
	    || s != k4.size()
	    || s != k5.size()) throw std::runtime_error("rk_combine_result: Input matrix size mismatch.");
	if (!tab.num_steps == 6) throw std::runtime_error("Need num_steps == 6 in rk_combine_result");

	if (isCudaEnabled()) {
#ifdef HAVE_CUDA
		rk_combine_result_cuda(h, tab, k0, k1, k2, k3, k4, k5, y, y_error, isCuda64Enabled());
#else
		assert(0);
#endif
	} else {
		rk_combine_result_cpu(h, tab, k0, k1, k2, k3, k4, k5, y, y_error);
	}
}

VectorMatrix linearInterpolate(const VectorMatrix &src, Shape dest_dim)
{
	Shape src_dim = src.getShape();

	if (src_dim.getRank() != dest_dim.getRank()) {
		throw std::runtime_error("linearInterpolate: Source and destination matrices need to have the same rank.");
	}

	if (src_dim.getRank() != 3) {
		throw std::runtime_error("linearInterpolate: Fixme: Need to have matrix of rank 3");
	}

	VectorMatrix dest(dest_dim);

	VectorMatrix::      accessor dest_acc(dest);
	VectorMatrix::const_accessor  src_acc(src);

	const bool sing_x = (src_dim.getDim(0) == 1);
	const bool sing_y = (src_dim.getDim(1) == 1);
	const bool sing_z = (src_dim.getDim(2) == 1);

	Vector3d scale(1.0, 1.0, 1.0);
	if (!sing_x) scale.x = double(dest_dim.getDim(0)-1) / double(src_dim.getDim(0)-1);
	if (!sing_y) scale.y = double(dest_dim.getDim(1)-1) / double(src_dim.getDim(1)-1);
	if (!sing_z) scale.z = double(dest_dim.getDim(2)-1) / double(src_dim.getDim(2)-1);

	for (int k=0; k<dest_dim.getDim(2); ++k)
	for (int j=0; j<dest_dim.getDim(1); ++j)
	for (int i=0; i<dest_dim.getDim(0); ++i) {
		// (x,y,z): coordinates of point with indices (i,j,k) in dst matrix
		const double x = i / scale.x;
		const double y = j / scale.y;
		const double z = k / scale.z;

		const double u = x - std::floor(x);
		const double v = y - std::floor(y);
		const double w = z - std::floor(z);

		const int I = std::floor(x);
		const int J = std::floor(y);
		const int K = std::floor(z);

		Vector3d tmp(0.0, 0.0, 0.0);
		if (true                         ) tmp = tmp + (1.0-u) * (1.0-v) * (1.0-w) * src_acc.get(I  , J  , K  );
		if (                      !sing_z) tmp = tmp + (1.0-u) * (1.0-v) *      w  * src_acc.get(I  , J  , K+1);
		if (           !sing_y           ) tmp = tmp + (1.0-u) *      v  * (1.0-w) * src_acc.get(I  , J+1, K  );
		if (           !sing_y && !sing_z) tmp = tmp + (1.0-u) *      v  *      w  * src_acc.get(I  , J+1, K+1);
		if (!sing_x                      ) tmp = tmp +      u  * (1.0-v) * (1.0-w) * src_acc.get(I+1, J  , K  );
		if (!sing_x            && !sing_z) tmp = tmp +      u  * (1.0-v) *      w  * src_acc.get(I+1, J  , K+1);
		if (!sing_x && !sing_y           ) tmp = tmp +      u  *      v  * (1.0-w) * src_acc.get(I+1, J+1, K  );
		if (!sing_x && !sing_y && !sing_z) tmp = tmp +      u  *      v  *      w  * src_acc.get(I+1, J+1, K+1);
		dest_acc.set(i, j, k, tmp);
	}

	return dest;
}

void Transposer_CUDA::copy_pad(const VectorMatrix &M, float *out_x, float *out_y, float *out_z)
{
	// Ifdef HAVE_CUDA_64, we directly support input matrices that
	// are stored with 64 bit precision on the GPU.
#ifdef HAVE_CUDA_64
	const bool M_is_cuda64_bit = M.isCached(2); // 0 = CPU device, 2 = CUDA_64 device
	if (M_is_cuda64_bit) {
		VectorMatrix::const_cu64_accessor M_acc(M);
		// xyz, M -> s1
		cuda_copy_pad_r2r(dim_x, dim_y, dim_z, exp_x, M_acc.ptr_x(), M_acc.ptr_y(), M_acc.ptr_z(), out_x, out_y, out_z);
	} 
	else
#endif
	{
		VectorMatrix::const_cu32_accessor M_acc(M);
		// xyz, M -> s1
		cuda_copy_pad_r2r(dim_x, dim_y, dim_z, exp_x, M_acc.ptr_x(), M_acc.ptr_y(), M_acc.ptr_z(), out_x, out_y, out_z);
	}
}

ConstVectorMatrixAccessor::ConstVectorMatrixAccessor(const VectorMatrix &mat) : mat(mat) 
{
	mat.readLock(0);
	data_x = static_cast<CPUArray*>(mat.getArray(0, 0))->ptr();
	data_y = static_cast<CPUArray*>(mat.getArray(0, 1))->ptr();
	data_z = static_cast<CPUArray*>(mat.getArray(0, 2))->ptr();

	// Precalculate strides
	const int rank = mat.getShape().getRank();
	strides[0] = 1;
	strides[1] = strides[0] * (rank > 0 ? mat.getShape().getDim(0) : 1);
	strides[2] = strides[1] * (rank > 1 ? mat.getShape().getDim(1) : 1);
	strides[3] = strides[2] * (rank > 2 ? mat.getShape().getDim(2) : 1);
}