return self + other;
}
+/** Multiply an indexed expression with a scalar. This function is used
+ * internally by simplify_indexed().
+ *
+ * @param self Indexed expression; it's base object is *this
+ * @param other Numeric value
+ * @return product of self and other
+ * @see ex::simplify_indexed() */
+ex basic::scalar_mul_indexed(const ex & self, const numeric & 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
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 ex scalar_mul_indexed(const ex & self, const numeric & 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;
}
find_free_and_dummy(un.begin(), un.end(), free_indices, dummy_indices);
+ ex r;
if (something_changed) {
if (non_commutative)
- return ncmul(v);
+ r = ncmul(v);
else
- return mul(v);
+ r = mul(v);
} else
- return e;
+ r = e;
+
+ // Product of indexed object with a scalar?
+ if (is_ex_exactly_of_type(r, mul) && r.nops() == 2
+ && is_ex_exactly_of_type(r.op(1), numeric) && is_ex_of_type(r.op(0), indexed))
+ return r.op(0).op(0).bp->scalar_mul_indexed(r.op(0), ex_to_numeric(r.op(1)));
+ else
+ return r;
}
/** Simplify indexed expression, return list of free indices. */
| exp '!' {$$ = factorial($1);}
| '(' exp ')' {$$ = $2;}
| '[' list_or_empty ']' {$$ = $2;}
- | T_MATRIX_BEGIN matrix T_MATRIX_END {$$ = lst_to_matrix($2);}
+ | T_MATRIX_BEGIN matrix T_MATRIX_END {$$ = lst_to_matrix(ex_to_lst($2));}
;
exprseq : exp {$$ = exprseq($1);}
debugmsg("matrix ctor from unsigned,unsigned,exvector",LOGLEVEL_CONSTRUCT);
}
+/** Construct matrix from (flat) list of elements. If the list has fewer
+ * elements than the matrix, the remaining matrix elements are set to zero.
+ * If the list has more elements than the matrix, the excessive elements are
+ * thrown away. */
+matrix::matrix(unsigned r, unsigned c, const lst & l)
+ : inherited(TINFO_matrix), row(r), col(c)
+{
+ debugmsg("matrix ctor from unsigned,unsigned,lst",LOGLEVEL_CONSTRUCT);
+ m.resize(r*c, _ex0());
+
+ for (unsigned i=0; i<l.nops(); i++) {
+ unsigned x = i % c;
+ unsigned y = i / c;
+ if (y >= r)
+ break; // matrix smaller than list: throw away excessive elements
+ m[y*c+x] = l.op(i);
+ }
+}
+
//////////
// archiving
//////////
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(self.op(0), matrix));
GINAC_ASSERT(is_ex_of_type(other, indexed));
GINAC_ASSERT(self.nops() == 2 || self.nops() == 3);
return self + other;
}
+/** Product of an indexed matrix with a number. */
+ex matrix::scalar_mul_indexed(const ex & self, const numeric & other) const
+{
+ GINAC_ASSERT(is_ex_of_type(self, indexed));
+ GINAC_ASSERT(is_ex_of_type(self.op(0), matrix));
+ GINAC_ASSERT(self.nops() == 2 || self.nops() == 3);
+
+ const matrix &self_matrix = ex_to_matrix(self.op(0));
+
+ if (self.nops() == 2)
+ return indexed(self_matrix.mul(other), self.op(1));
+ else // self.nops() == 3
+ return indexed(self_matrix.mul(other), self.op(1), self.op(2));
+}
+
/** Contraction of an indexed matrix with something else. */
bool matrix::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
{
}
+/** Product of matrix and scalar. */
+matrix matrix::mul(const numeric & other) const
+{
+ exvector prod(row * col);
+
+ for (unsigned r=0; r<row; ++r)
+ for (unsigned c=0; c<col; ++c)
+ prod[r*col+c] = m[r*col+c] * other;
+
+ return matrix(row, col, prod);
+}
+
+
/** operator() to access elements.
*
* @param ro row of element
return k;
}
-/** Convert list of lists to matrix. */
-ex lst_to_matrix(const ex &l)
+ex lst_to_matrix(const lst & l)
{
- if (!is_ex_of_type(l, lst))
- throw(std::invalid_argument("argument to lst_to_matrix() must be a lst"));
-
// Find number of rows and columns
unsigned rows = l.nops(), cols = 0, i, j;
for (i=0; i<rows; i++)
return m;
}
+ex diag_matrix(const lst & l)
+{
+ unsigned dim = l.nops();
+
+ matrix &m = *new matrix(dim, dim);
+ for (unsigned i=0; i<dim; i++)
+ m.set(i, i, l.op(i));
+
+ return m;
+}
+
} // namespace GiNaC
public:
matrix(unsigned r, unsigned c);
matrix(unsigned r, unsigned c, const exvector & m2);
+ matrix(unsigned r, unsigned c, const lst & l);
// functions overriding virtual functions from bases classes
public:
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;
+ ex scalar_mul_indexed(const ex & self, const numeric & other) const;
bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const;
protected:
unsigned return_type(void) const { return return_types::noncommutative; };
matrix add(const matrix & other) const;
matrix sub(const matrix & other) const;
matrix mul(const matrix & other) const;
+ matrix mul(const numeric & other) const;
const ex & operator() (unsigned ro, unsigned co) const;
matrix & set(unsigned ro, unsigned co, ex value);
matrix transpose(void) const;
return static_cast<const matrix &>(*e.bp);
}
-extern ex lst_to_matrix(const ex &l);
+/** Convert list of lists to matrix. */
+extern ex lst_to_matrix(const lst & l);
+
+/** Convert list of diagonal elements to matrix. */
+extern ex diag_matrix(const lst & l);
} // namespace GiNaC