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());