- simplify_indexed() simplifies c*(M).i.j -> (M').i.j (c: numeric, M: matrix)

  with M'=c*M
- new member function matrix::mul(const numeric &other) to multiply matrix
  with scalar
- lst_to_matrix() takes a "lst &" instead of an "ex &"
- added constructor of matrix from flat list
- added function diag_matrix() to construct matrix from list of diagonal
  elements

diff --git a/ginac/basic.cpp b/ginac/basic.cpp
index 6415f33fde5a79456ebb97ba0e81f36139110c97..1c6fe6fa0d1bf5553cafb38e81cc8f1c991de894 100644 (file)
--- a/ginac/basic.cpp
+++ b/ginac/basic.cpp
@@ -310,6 +310,18 @@ ex basic::add_indexed(const ex & self, const ex & other) const
        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

diff --git a/ginac/basic.h b/ginac/basic.h
index 4277d949e83e5bca42a71b60131c4a0dada6a552..2e1aefdff859567a4da9605f91c04e42dcb17203 100644 (file)
--- a/ginac/basic.h
+++ b/ginac/basic.h
@@ -131,6 +131,7 @@ public: // only const functions please (may break reference counting)
        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;

diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp
index e23709d5f506647280eca2f54fb528f807a5dc8f..86fe6115d9b8863819c44e3a393d9c99cd0b01d4 100644 (file)
--- a/ginac/indexed.cpp
+++ b/ginac/indexed.cpp
@@ -679,13 +679,21 @@ try_again:
        }
        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. */

diff --git a/ginac/input_parser.yy b/ginac/input_parser.yy
index d1e876690967227869fa4bc18007351a443f877c..9c4ffd971ac66b89902631818862e949368f7900 100644 (file)
--- a/ginac/input_parser.yy
+++ b/ginac/input_parser.yy
@@ -115,7 +115,7 @@ exp : T_NUMBER              {$$ = $1;}
        | 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);}

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index f34534afdb1ab3ee39b567bd11c7ed7408010855..650799e55f6353db6807e1d42a8d7904cc55291b 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -95,6 +95,25 @@ matrix::matrix(unsigned r, unsigned c, const exvector & m2)
        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
 //////////
@@ -385,6 +404,7 @@ ex matrix::eval_indexed(const basic & i) const
 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);
 
@@ -416,6 +436,21 @@ ex matrix::add_indexed(const ex & self, const ex & other) const
        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
 {
@@ -581,6 +616,19 @@ matrix matrix::mul(const matrix & other) 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
@@ -1392,12 +1440,8 @@ int matrix::pivot(unsigned ro, unsigned co, bool symbolic)
        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++)
@@ -1415,4 +1459,15 @@ ex lst_to_matrix(const ex &l)
        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

diff --git a/ginac/matrix.h b/ginac/matrix.h
index 53afeb2e6e05c8ab4448ec4f746c3bd39775db82..1041c6433c8001f64f7c79d7533e929c6f870abc 100644 (file)
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -38,6 +38,7 @@ class matrix : public basic
 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:
@@ -53,6 +54,7 @@ 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; };
@@ -68,6 +70,7 @@ public:
        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;
@@ -137,7 +140,11 @@ inline const matrix &ex_to_matrix(const ex &e)
        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