Conjugate Gradient

[sueda]

I'm trying to write CG with taco, but I'm having difficulties. Here is my starter code, which doesn't compile yet.

void conjgrad(taco::Tensor<double> &A, taco::Tensor<double> &b, taco::Tensor<double> &x)
{
	taco::Format  dv({taco::Dense});
	int n = b.getDimensions()[0];
	
	taco::Tensor<double> r({n}, dv);
	taco::Tensor<double> p({n}, dv);
	taco::Tensor<double> Ap({n}, dv);

	double rsold, rsnew, alpha;
	taco::IndexVar i, j;
	
	r(i) = b(i) - A(i,j) * x(j); // A can't be const!
	p = r;
	rsold = r(i) * r(i);
	
	int itermax = n;
	double thresh = 1e-10;
	for(int iter = 0; iter < itermax; ++iter) {
		Ap(i) = A(i,j) * p(j);
		alpha = rsold / (p(i) * Ap(i));
		x(i) = x(i) + alpha * p(i);
		r(i) = r(i) - alpha * Ap(i);
		rsnew = r(i) * r(i);
		if(sqrt(rsnew) < thresh) {
			break;
		}
		p(i) = r(i) + (rsnew / rsold) * p(i);
		rsold = rsnew;
	}
}

Here are some comments/issues.

• It complains if I try to declare A and b as const.

• I know that I'm missing evaluate()/assemble()/compile(), but it is difficult to get the correct syntax. If I put evaluate() everywhere, the code would get messy pretty quickly. For instance, assigning the result of a dot product to a double would take multiple lines, which would be cumbersome.

• Since I declare using namespace Eigen in this code, if I declare using namespace taco I get some name collisions.

• I would prefer to say b.getDimension(0) to get the number of rows in b.

• Is there an easy way to declare a Tensor to be of same storage type as another? I want r, p, and Ap to be the same as b without having to create dv.

[stephenchouca]

Thanks for your feedback! 57661f4 and eedb9db should address your first and fourth points respectively. With regards to your last point, the Tensor class does in fact have a getFormat method that returns the format of the tensor, so you can write something like

taco::Tensor<double> r({n}, b.getFormat());