@subsection Linear algebra
The @code{matrix} class can be used with indices to do some simple linear
-algebra (products of vectors and matrices, traces and scalar products):
+algebra (sums and products of vectors and matrices, traces and scalar
+products):
@example
@{
cout << e.simplify_indexed() << endl;
// -> [[ [[2*y+x]], [[4*y+3*x]] ]].i
- e = indexed(A, i, j) * indexed(X, i);
+ e = indexed(A, i, j) * indexed(X, i) + indexed(X, j);
cout << e.simplify_indexed() << endl;
- // -> [[ [[3*y+x,4*y+2*x]] ]].j
+ // -> [[ [[3*y+2*x,5*y+2*x]] ]].j
@}
@end example
-You can of course obtain the same results with the @code{matrix::mul()}
-and @code{matrix::trace()} methods but with indices you don't have to
-worry about transposing matrices.
+You can of course obtain the same results with the @code{matrix::add()},
+@code{matrix::mul()} and @code{matrix::trace()} methods but with indices you
+don't have to worry about transposing matrices.
Matrix indices always start at 0 and their dimension must match the number
of rows/columns of the matrix. Matrices with one row or one column are
vectors and can have one or two indices (it doesn't matter whether it's a
-row or a columnt vector). Other matrices must have two indices.
+row or a column vector). Other matrices must have two indices.
You should be careful when using indices with variance on matrices. GiNaC
doesn't look at the variance and doesn't know that @samp{F~mu~nu} and
return i.hold();
}
+/** Add two indexed expressions. They are guaranteed to be of class indexed
+ * (or a subclass) and their indices are compatible. This function is used
+ * internally by simplify_indexed().
+ *
+ * @param self First indexed expression; it's base object is *this
+ * @param other Second indexed expression
+ * @return sum of self and other
+ * @see ex::simplify_indexed() */
+ex basic::add_indexed(const ex & self, const ex & other) const
+{
+ return self + other;
+}
+
/** Try to contract two indexed expressions that appear in the same product.
* If a contraction exists, the function overwrites one or both of the
* expressions and returns true. Otherwise it returns false. It is
* guaranteed that both expressions are of class indexed (or a subclass)
- * and that at least one dummy index has been found.
+ * and that at least one dummy index has been found. This functions is
+ * used internally by simplify_indexed().
*
* @param self Pointer to first indexed expression; it's base object is *this
* @param other Pointer to second indexed expression
* @param v The complete vector of factors
- * @return true if the contraction was successful, false otherwise */
+ * @return true if the contraction was successful, false otherwise
+ * @see ex::simplify_indexed() */
bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
{
// Do nothing
virtual exvector get_free_indices(void) const;
virtual ex simplify_ncmul(const exvector & v) const;
virtual ex eval_indexed(const basic & i) const;
+ virtual ex add_indexed(const ex & self, const ex & other) const;
virtual bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const;
protected: // non-const functions should be called from class ex only
virtual ex derivative(const symbol & s) const;
// Simplification of sum = sum of simplifications, check consistency of
// free indices in each term
if (is_ex_exactly_of_type(e_expanded, add)) {
- ex sum = _ex0();
+ ex sum = simplify_indexed(e_expanded.op(0), free_indices, sp);
- for (unsigned i=0; i<e_expanded.nops(); i++) {
+ for (unsigned i=1; i<e_expanded.nops(); i++) {
exvector free_indices_of_term;
- sum += simplify_indexed(e_expanded.op(i), free_indices_of_term, sp);
- if (i == 0)
- free_indices = free_indices_of_term;
- else if (!indices_consistent(free_indices, free_indices_of_term))
+ ex term = simplify_indexed(e_expanded.op(i), free_indices_of_term, sp);
+ if (!indices_consistent(free_indices, free_indices_of_term))
throw (std::runtime_error("simplify_indexed: inconsistent indices in sum"));
+ if (is_ex_of_type(sum, indexed) && is_ex_of_type(term, indexed))
+ sum = sum.op(0).bp->add_indexed(sum, term);
+ else
+ sum += term;
}
return sum;
return i.hold();
}
+/** Sum of two indexed matrices. */
+ex matrix::add_indexed(const ex & self, const ex & other) const
+{
+ GINAC_ASSERT(is_ex_of_type(self, indexed));
+ GINAC_ASSERT(is_ex_of_type(other, indexed));
+ GINAC_ASSERT(self.nops() == 2 || self.nops() == 3);
+
+ // Only add two matrices
+ if (is_ex_of_type(other.op(0), matrix)) {
+ GINAC_ASSERT(other.nops() == 2 || other.nops() == 3);
+
+ const matrix &self_matrix = ex_to_matrix(self.op(0));
+ const matrix &other_matrix = ex_to_matrix(other.op(0));
+
+ if (self.nops() == 2 && other.nops() == 2) { // vector + vector
+
+ if (self_matrix.row == other_matrix.row)
+ return indexed(self_matrix.add(other_matrix), self.op(1));
+ else if (self_matrix.row == other_matrix.col)
+ return indexed(self_matrix.add(other_matrix.transpose()), self.op(1));
+
+ } else if (self.nops() == 3 && other.nops() == 3) { // matrix + matrix
+
+ if (self.op(1).is_equal(other.op(1)) && self.op(2).is_equal(other.op(2)))
+ return indexed(self_matrix.add(other_matrix), self.op(1), self.op(2));
+ else if (self.op(1).is_equal(other.op(2)) && self.op(2).is_equal(other.op(1)))
+ return indexed(self_matrix.add(other_matrix.transpose()), self.op(1), self.op(2));
+
+ }
+ }
+
+ // Don't know what to do, return unevaluated sum
+ return self + other;
+}
+
/** Contraction of an indexed matrix with something else. */
bool matrix::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
{
ex evalf(int level=0) const;
ex subs(const lst & ls, const lst & lr) const;
ex eval_indexed(const basic & i) const;
+ ex add_indexed(const ex & self, const ex & other) const;
bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const;
protected:
unsigned return_type(void) const { return return_types::noncommutative; };
git://www.ginac.de / ginac.git / commitdiff