[Trilinos-Users] matrix/vector recipes?

Baker, Christopher G. bakercg at ornl.gov
Thu Jan 21 11:57:23 MST 2010


Nico, regarding the example of adding two Tpetra vectors, I assume that it was "not quickly realizable" because you couldn't easily find the appropriate class member.
I assume you eventually settled onto
z = x + y
via
z.update(1.0,x,1.0,y,0.0);
and
y = x + y
via
y.update(1.0,x,1.0);

Regarding the difficulty in deciding this, I apologize. Two issues conspire against users here; I shall explain.

The first is that Vector is a subclass of MultiVector, from which it derives most of its functionality. Unfortunately, the current behavior of the Doxygen webpage is to hide the methods that a derived class inherits from its parent. As a result, if you go to Vector's page, you don't see much. You must either realize that MultiVector is the parent and view MultiVector's page, or click on the "List of all members" link on the Vector page. At which point, you will realize that there is no add() method.

The second issue is the balance between providing the functionality that users expect (those operations that are "trivial to describe") and limiting the size of the interface for a class. In the case of Tpetra, there is the further goal, to make Tpetra look as much like Epetra as possible, in order to ease the transition of Epetra users to Tpetra (hence, update() instead of add()). As a result, there is no method MultiVector::add(), and therefore no Vector::add(). One solution to that is to extend the functionality of the classes via non-member arithmetic routines. These are simply wrappers that call the (less numerous) member functions of the Vector class. For example,
z = Tpetra::add(x,y);
would allocate a new Vector, which would hold the addition of x and y.
It has been under discussion for some time to augment Tpetra with a collection of non-member functions. This is similar to the methodology of the Thyra package (in fact, upon renewal of the Thyra/Tpetra adaptors, the Thyra non-member arithmetic interface could be used.) Core development of Tpetra has taken priority, unfortunately. We would be interested in hearing from you, and other Trilinos users, your opinion on this.

Regarding the padded matrix, the approach taken will depend on the desired form of the output. If it is not necessary that the operator ~A maps a contiguous vector to another contiguous vector, then simplest approach is to apply the operator implicitly.
[w z] = [A b; c d] [x y] = [A*x + b*y; c*x + d*y]
via something like
w.scale(b,y); A.multiply(1.0,x,1.0,w);
z = x.dot(z); z += d*y;
However, this requires that you pass around A, b, c, d, x, y, w and z (although you could encapsulate them).
If you must do this via the Operator/Vector interface, it is a matter of creating a new Tpetra::Operator which stores A, b, c and d. In addition, a new class derived from Vector must be created, which contains a Vector pointing to the first Vector with an appropriately-sized Map; and a pointer to the scalar entry at the bottom. The custom Operator would dynamically cast to the custom Vector, extract these, and perform apply() as above.

Note, Thyra already provides this functionality, via its BlockedLinearOpBase and ProductMultiVectorBase. This would be implemented by encapsulating the effect of b*y, c*x and d*y into Operator objects, which would be used to construct a new operator (A is already an Operator):
TildeA = Thyra::block2x2( A, BYOp, CXOp, DYOp );
XYVec = Thyra::defaultProductMultiVector( ..., tuple(x,y) );
WZVec = Thyra::defaultProductMultiVector( ..., tuple(w,z) );

This would require that your code was okay with passing around Thyra::LinearOpBase and Thyra::MultiVectorBase objects, instead of Tpetra::Operator and Tpetra::MultiVector objects. Note, it would be messier to provide the same functionality in Tpetra as the it doesn't have a pure virtual MultiVector. Therefore, I think the only way is the previously mentioned, contiguous memory sub-vector approach.

Lastly, as for a venue for recipes, I don't know of such a thing. It seems that the trilinos-users mailing list is the best place for that, although it has tedious searches and disruptive submissions. Perhaps Trilinos Framework could chime in on this? Time for an official Trilinos wiki?

Chris



On 1/21/10 11:26 AM, "Nico Schlömer" <nico.schloemer at ua.ac.be> wrote:

Hi,

I recently implemented a function that takes a Tpetra::RowMatrix A, two
vectors b,c and a scalar d as input, and spits out the padded matrix ~A

    A b
    c  d

Despite the fact that this might actually be something that "you don't
want to do", it's also something that seems trivial if you come from
Python, MATLAB and the like, but due to Trilinos' sophisticated memory
management is principally requires quite a number of lines of code to be
written (potentially containing bugs).

During everyday programming, I sometimes stumble upon things that are
trivial to describe but not quickly realizable with Trilinos. Take for
example the addition of two Tpetra::Vectors.

Now, is there any place where recipes for simple matrix/vector
operations and other useful snippets are exchanged?

For the particular example above, what I'm doing right now is that I
*copy* over the values of A, b, and c to ~A, which I'd like to avoid.
Any hints here?

Cheers,
Nico

_______________________________________________
Trilinos-Users mailing list
Trilinos-Users at software.sandia.gov
http://software.sandia.gov/mailman/listinfo/trilinos-users






More information about the Trilinos-Users mailing list