vector_extra.h

Go to the documentation of this file.

/* KG: 1998/08/17 */
/* $Id: vector_extra.h,v 1.9.2.16 2019/05/28 11:13:02 garloff Exp $ */

#ifndef TBCI_VECTOR_EXTRA_H
#define TBCI_VECTOR_EXTRA_H  

#include "tbci/vector.h"
#include "tbci/matrix.h"

// avoid -fguiding-decls (here: unnecessary)
#if !defined(NO_GD) && !defined(AUTO_DECL)
# include "tbci/vector_extra_gd.h"
#endif

NAMESPACE_TBCI

// elementwise multiplication
INST(template <typename T> class TVector  friend TVector<T>   emul (const Vector<T>&, const Vector<T>&);)
template <typename T>
inline TVector<T>   emul (const Vector<T>& a, const Vector<T>& b)
{
  const unsigned long dim = b.size();
  TVector<T> s (dim);
  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, s);
  do_vec_vec_mul<T> (dim, &s.setval(0), &a.getcref(0), &b.getcref(0));
  return s;
}

INST(template <typename T> class TVector friend TVector<T>   emul (TVector<T>, const Vector<T>&);)
template <typename T>
inline TVector<T>   emul (TVector<T> a, const Vector<T>& b)
{
  const unsigned long dim = b.size();
  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, a);
  do_vec_mul_vec<T> (dim, &a.setval(0), &b.getcref(0));
  return a;
}

INST(template <typename T> class TVector friend TVector<T>   emul (const Vector<T>&, TVector<T>);)
template <typename T>
inline TVector<T>   emul (const Vector<T>& a, TVector<T> b)
{
  const unsigned long dim = b.size();
  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, b);
  do_vec_mul_vec<T> (dim, &b.setval(0), &a.getcref(0));
  return b;
}

INST(template <typename T> class TVector friend TVector<T>   emul (TVector<T>, TVector<T>);)
template <typename T>
inline TVector<T>   emul (TVector<T> a, TVector<T> b)
{
  const unsigned long dim = b.size();
  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, a);
  do_vec_mul_vec<T> (dim, &a.setval(0), &b.getcref(0));
  b.destroy (); return a;
}

INST(template <typename T> class TVector  friend TVector<T>  cemul (const Vector<T>&, const Vector<T>&);)
template <typename T>
inline TVector<T>  cemul (const Vector<T>& a, const Vector<T>& b)
{
  const unsigned long dim = b.size();
  TVector<T> s(dim);
  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, s);
  do_vec_vec_cmul<T> (dim, &s.setval(0), &a.getcref(0), &b.getcref(0));
  return s;
}

INST(template <typename T> class TVector  friend TVector<T>  cemul (TVector<T>, const Vector<T>&);)
template <typename T>
inline TVector<T>  cemul (TVector<T> a, const Vector<T>& b)
{
  const unsigned long dim = b.size();
  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, a);
  do_vec_cmul_vec<T> (dim, &a.setval(0), &b.getcref(0));
  return a;
}

INST(template <typename T> class TVector  friend TVector<T>  cemul (const Vector<T>&, TVector<T>);)
template <typename T>
inline TVector<T>  cemul (const Vector<T>& a, TVector<T> b)
{
  const unsigned long dim = b.size();
  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, b);
  do_vec_cmul_vec_inv<T> (dim, &b.setval(0), &a.getcref(0));
  return b;
}

INST(template <typename T> class TVector  friend TVector<T>  cemul (TVector<T>, TVector<T>);)
template <typename T>
inline TVector<T>  cemul (TVector<T> a, TVector<T> b)
{
  const unsigned long dim = b.size();
  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, a);
  do_vec_cmul_vec<T> (dim, &a.setval(0), &b.getcref(0));
  b.destroy (); return a;
}


// elementwise division

INST(template <typename T> class TVector friend TVector<T>   ediv (const Vector<T>&, const Vector<T>&);)
template <typename T>
inline TVector<T>   ediv (const Vector<T>& a, const Vector<T>& b)
{
  TVector<T> s (b.size());
  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), s);
  do_vec_vec_div<T> (b.size(), &s.setval(0), &a.getcref(0), &b.getcref(0));
  return s;
}

INST(template <typename T> class TVector friend TVector<T>   ediv (TVector<T>, const Vector<T>&);)
template <typename T>
inline TVector<T>   ediv (TVector<T> a, const Vector<T>& b)
{
  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), a);
  do_vec_div_vec<T> (b.size(), &a.setval(0), &b.getcref(0));
  return a;
}

INST(template <typename T> class TVector friend TVector<T>   ediv (const Vector<T>&, TVector<T>);)
template <typename T>
inline TVector<T>   ediv (const Vector<T>& a, TVector<T> b)
{
  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), b);
  do_vec_div_vec_inv<T> (b.size(), &b.setval(0), &a.getcref(0));
  return b;
}

INST(template <typename T> class TVector friend TVector<T>   ediv (TVector<T>, TVector<T>);)
template <typename T>
inline TVector<T>   ediv (TVector<T> a, TVector<T> b)
{
  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), a);
  do_vec_div_vec<T> (b.size(), &a.setval(0), &b.getcref(0));
  b.destroy (); return a;
}

INST(template <typename T> class TVector friend TVector<T>  cediv (const Vector<T>&, const Vector<T>&);)
template <typename T>
inline TVector<T>  cediv (const Vector<T>& a, const Vector<T>& b)
{
  TVector<T> s (b.size());
  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), s);
  do_vec_vec_cdiv<T> (b.size(), &s.setval(0), &a.getcref(0), &b.getcref(0));
  return s;
}

INST(template <typename T> class TVector friend TVector<T>  cediv (TVector<T>, const Vector<T>&);)
template <typename T>
inline TVector<T>  cediv (TVector<T> a, const Vector<T>& b)
{
  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), a);
  do_vec_cdiv_vec<T> (b.size(), &a.setval(0), &b.getcref(0));
  return a;
}

INST(template <typename T> class TVector friend TVector<T>  cediv (const Vector<T>&, TVector<T>);)
template <typename T>
inline TVector<T>  cediv (const Vector<T>& a, TVector<T> b)
{
  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), b);
  do_vec_cdiv_vec_inv<T> (a.size(), &b.setval(0), &a.getcref(0));
  return b;
}

INST(template <typename T> class TVector friend TVector<T>  cediv (TVector<T>, TVector<T>);)
template <typename T>
inline TVector<T>  cediv (TVector<T> a, TVector<T> b)
{
  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), a);
  do_vec_cdiv_vec<T> (b.size(), &a.setval(0), &b.getcref(0));
  b.destroy (); return a;
}


// Now for the combinations with TSVector<T>
// Not yet optimized with the unrolled kernels from vec_kern_unr_pref.h

INST(template <typename T> class TSVector friend TVector<T>   emul (TVector<T>, const TSVector<T>&);)
template <typename T>
TVector<T>   emul (TVector<T> tv, const TSVector<T>& ts)
{
   BCHK(tv.size() != ts.size(), VecErr, Wrong dim in emul, ts.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.setval(i) *= ts.get(i);
   ts.destroy (); return tv;
}

INST(template <typename T> class TSVector friend TVector<T>   emul (const TSVector<T>&, TVector<T>);)
template <typename T>
TVector<T>   emul (const TSVector<T>& ts, TVector<T> tv)
{
   BCHK(tv.size() != ts.size(), VecErr, Wrong dim in emul, tv.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (ts.get(i) * tv.get(i), i);
   ts.destroy (); return tv;
}


INST(template <typename T> class TSVector friend TVector<T>   emul (const Vector<T>&, const TSVector<T>&);)
template <typename T>
TVector<T>   emul (const Vector<T>& v, const TSVector<T>& ts)
{
   TVector<T> tv;
   if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
   else tv.resize (ts.size());
   BCHK(v.size() != ts.size(), VecErr, Wrong dim in emul, ts.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (v.get(i) * ts.get(i), i);
   return tv;
}

INST(template <typename T> class TSVector friend TVector<T>   emul (const TSVector<T>&, const Vector<T>&);)
template <typename T>
TVector<T>   emul (const TSVector<T>& ts, const Vector<T>& v)
{
   TVector<T> tv;
   if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
   else tv.resize (ts.size());
   BCHK(v.size() != ts.size(), VecErr, Wrong dim in emul, v.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (ts.get(i) * v.get(i), i);
   return tv;
}

INST(template <typename T> class TSVector friend TVector<T>   emul (const TSVector<T>&, const TSVector<T>&);)
template <typename T>
TVector<T>   emul (const TSVector<T>& ts1, const TSVector<T>& ts2)
{
   TVector<T> tv;
   if (ts1.mut) { tv.setptr(ts1.vecptr()); tv.setsize(ts1.size()); }
   else if (ts2.mut) { tv.setptr(ts2.vecptr()); tv.setsize(ts2.size()); }
        else tv.resize (ts1.size());
   BCHK(ts1.size() != ts2.size(), VecErr, Wrong dim in emul, ts2.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (ts1.get(i), i);
   if (ts1.mut && ts2.mut) ts2.destroy ();
   return tv;
}


INST(template <typename T> class TSVector friend TVector<T>  cemul (TVector<T>, const TSVector<T>&);)
template <typename T>
TVector<T>  cemul (TVector<T> tv, const TSVector<T>& ts)
{
   BCHK(tv.size() != ts.size(), VecErr, Wrong dim in emul, ts.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (CPLX__ conj(tv.get(i)) * ts.get(i), i);
   ts.destroy (); return tv;
}

INST(template <typename T> class TSVector friend TVector<T>  cemul (const TSVector<T>&, TVector<T>);)
template <typename T>
TVector<T>  cemul (const TSVector<T>& ts, TVector<T> tv)
{
   BCHK(tv.size() != ts.size(), VecErr, Wrong dim in emul, tv.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (CPLX__ conj (ts.get(i)) * tv.get(i), i);
   ts.destroy (); return tv;
}

INST(template <typename T> class TSVector friend TVector<T>  cemul (const Vector<T>&, const TSVector<T>&);)
template <typename T>
TVector<T>  cemul (const Vector<T>& v, const TSVector<T>& ts)
{
   TVector<T> tv;
   if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
   else tv.resize (ts.size());
   BCHK(v.size() != ts.size(), VecErr, Wrong dim in emul, ts.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (CPLX__ conj(v.get(i)) * ts.get(i), i);
   return tv;
}


INST(template <typename T> class TSVector friend TVector<T>  cemul (const TSVector<T>&, const Vector<T>&);)
template <typename T>
TVector<T>  cemul (const TSVector<T>& ts, const Vector<T>& v)
{
   TVector<T> tv;
   if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
   else tv.resize (ts.size());
   BCHK(v.size() != ts.size(), VecErr, Wrong dim in emul, v.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (CPLX__ conj (ts.get(i)) * v.get(i), i);
   return tv;
}

INST(template <typename T> class TSVector friend TVector<T>  cemul (const TSVector<T>&, const TSVector<T>&);)
template <typename T>
TVector<T>  cemul (const TSVector<T>& ts1, const TSVector<T>& ts2)
{
   TVector<T> tv;
   if (ts1.mut) { tv.setptr(ts1.vecptr()); tv.setsize(ts1.size()); }
   else if (ts2.mut) { tv.setptr(ts2.vecptr()); tv.setsize(ts2.size()); }
        else tv.resize (ts1.size());
   BCHK(ts1.size() != ts2.size(), VecErr, Wrong dim in emul, ts2.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (CPLX__ conj (ts1.get(i)), i);
   if (ts1.mut && ts2.mut) ts2.destroy ();
   return tv;
}


INST(template <typename T> class TSVector friend  TVector<T>   ediv (TVector<T>, const TSVector<T>&);)
template <typename T>
TVector<T>   ediv (TVector<T> tv, const TSVector<T>& ts)
{
   BCHK(tv.size() != ts.size(), VecErr, Wrong dim in ediv, ts.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.setval(i) /= ts.get(i);
   ts.destroy (); return tv;
}

INST(template <typename T> class TSVector friend  TVector<T>   ediv (const Vector<T>&, const TSVector<T>&);)
template <typename T>
TVector<T>   ediv (const Vector<T>& v, const TSVector<T>& ts)
{
   TVector<T> tv;
   if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
   else tv.resize (ts.size());
   BCHK(v.size() != ts.size(), VecErr, Wrong dim in ediv, ts.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (v.get(i) / ts.get(i), i);
   return tv;
}

INST(template <typename T> class TSVector friend  TVector<T>   ediv (const TSVector<T>&, const TSVector<T>&);)
template <typename T>
TVector<T>   ediv (const TSVector<T>& ts1, const TSVector<T>& ts2)
{
   TVector<T> tv;
   if (ts1.mut) { tv.setptr(ts1.vecptr()); tv.setsize(ts1.size()); }
   else if (ts2.mut) { tv.setptr(ts2.vecptr()); tv.setsize(ts2.size()); }
        else tv.resize (ts1.size());
   BCHK(ts1.size() != ts2.size(), VecErr, Wrong dim in ediv, ts2.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (ts1.get(i), i);
   if (ts1.mut && ts2.mut) ts2.destroy ();
   return tv;
}

INST(template <typename T> class TSVector friend  TVector<T>   ediv (const TSVector<T>&, const Vector<T>&);)
template <typename T>
TVector<T>   ediv (const TSVector<T>& ts, const Vector<T>& v)
{
   TVector<T> tv;
   if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
   else tv.resize (ts.size());
   BCHK(v.size() != ts.size(), VecErr, Wrong dim in ediv, v.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (ts.get(i) / v.get(i), i);
   return tv;
}

INST(template <typename T> class TSVector friend  TVector<T>   ediv (const TSVector<T>&, TVector<T>);)
template <typename T>
TVector<T>   ediv (const TSVector<T>& ts, TVector<T> tv)
{
   BCHK(tv.size() != ts.size(), VecErr, Wrong dim in ediv, tv.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (ts.get(i) / tv.get(i), i);
   ts.destroy (); return tv;
}


INST(template <typename T> class TSVector friend TVector<T>   cediv (TVector<T>, const TSVector<T>&);)
template <typename T>
TVector<T>   cediv (TVector<T> tv, const TSVector<T>& ts)
{
   BCHK(tv.size() != ts.size(), VecErr, Wrong dim in ediv, ts.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (CPLX__ conj(tv.get(i)) / ts.get(i), i);
   ts.destroy (); return tv;
}

INST(template <typename T> class TSVector friend TVector<T>   cediv (const TSVector<T>&, TVector<T>);)
template <typename T>
TVector<T>  cediv (const TSVector<T>& ts, TVector<T> tv)
{
   BCHK(tv.size() != ts.size(), VecErr, Wrong dim in ediv, tv.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (CPLX__ conj (ts.get(i)) / tv.get(i), i);
   ts.destroy (); return tv;
}

INST(template <typename T> class TSVector friend TVector<T>   cediv (const Vector<T>&, const TSVector<T>&);)
template <typename T>
TVector<T>  cediv (const Vector<T>& v, const TSVector<T>& ts)
{
   TVector<T> tv;
   if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
   else tv.resize (ts.size());
   BCHK(v.size() != ts.size(), VecErr, Wrong dim in ediv, ts.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (CPLX__ conj(v.get(i)) / ts.get(i), i);
   return tv;
}


INST(template <typename T> class TSVector friend TVector<T>   cediv (const TSVector<T>&, const Vector<T>&);)
template <typename T>
TVector<T>  cediv (const TSVector<T>& ts, const Vector<T>& v)
{
   TVector<T> tv;
   if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
   else tv.resize (ts.size());
   BCHK(v.size() != ts.size(), VecErr, Wrong dim in ediv, v.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (CPLX__ conj (ts.get(i)) / v.get(i), i);
   return tv;
}

INST(template <typename T> class TSVector friend TVector<T>   cediv (const TSVector<T>&, const TSVector<T>&);)
template <typename T>
TVector<T>  cediv (const TSVector<T>& ts1, const TSVector<T>& ts2)
{
   TVector<T> tv;
   if (ts1.mut) { tv.setptr(ts1.vecptr()); tv.setsize(ts1.size()); }
   else if (ts2.mut) { tv.setptr(ts2.vecptr()); tv.setsize(ts2.size()); }
        else tv.resize (ts1.size());
   BCHK(ts1.size() != ts2.size(), VecErr, Wrong dim in ediv, ts2.size(), tv); 
   for (register unsigned long i = 0; i < tv.size(); i++)
     tv.set (CPLX__ conj (ts1.get(i)), i);
   if (ts1.mut && ts2.mut) ts2.destroy ();
   return tv;
}



INST(template <typename T> class TVector friend TVector<T> ediv (const T&, TVector<T>);)
template <typename T>
inline TVector<T> ediv (const T& a, TVector<T> b)
{
        STD_SMP_TEMPLATE2C(val_div_vec, b.size(), &b.setval(0), a);
        return b;
}
        
INST(template <typename T> class TVector friend TVector<T> ediv (const T&, const Vector<T>&);)
template <typename T>
inline TVector<T> ediv (const T& a, const Vector<T>& b)
{ 
        const unsigned d = b.size();
        TVector<T> res (d);
        STD_SMP_TEMPLATE3VC(val_vec_div, d, &res.setval(0), &b.getcref(0), a);
        return res;
}
        
INST(template <typename T> class TSVector friend TVector<T> ediv (const T&, const TSVector<T>&);)
template <typename T>
inline TVector<T> ediv (const T& a, const TSVector<T>& ts)
{
        TSVector<T> tm (ts); tm.detach ();
        STD_SMP_TEMPLATE4C(val_svc_div, tm.size(), tm.vecptr(), &ts.getcref(0), a, ts.getfac());
        tm.getfacref() = (T)1; return TVector<T> (tm);
}

/* Summing up a vector is seldomly relevant.
 * TODO: Figure a way to use the algorithms below for other reductions, 
 *  such as dot products, fabs, Mat-Vec and Mat-Mat products.
 * Also create SIMD and parallel versions of it.
 */

/* http://en.wikipedia.org/wiki/Kahan_summation_algorithm */
INST(template <typename T> class Vector friend T kahan_sum(const Vector<T>&);)
template <typename T> T kahan_sum(const Vector<T>& v)
{
        T sum = T(0.0); T comp = T(0.0);
        const unsigned long sz = v.size();
        for (register unsigned long i = 0; i < sz; ++i) {
                T y = v(i) - comp;
                T t = sum + y;
                comp = (t - sum) - y;
                sum = t;
        }
        return sum /* - comp */;
}

/* Improved algo */
INST(template <typename T> class Vector friend T exact_sum(const Vector<T>&);)
template <typename T> T exact_sum(const Vector<T>& v)
{
        T sum = T(0.0); T comp = T(0.0);
        const unsigned long sz = v.size();
        for (register unsigned long i = 0; i < sz; ++i) {
                T y = v(i);
                T t = sum + y;
                comp += (t - sum) - y;
                sum = t;
        }
        return sum - comp;
}

// elementwise multiplication / division (with and without taking complex conj.)
// for matrices

# define _VEC getvec()
# define _ENDVEC getendvec()
# define _DIM size()
# define _ROW rows()
# define _COL columns()
# define _FAC getfac()

INST(template <typename T> class TMatrix friend TMatrix<T> emul (const Matrix<T>&, const Matrix<T>&);)
template <typename T>
TMatrix<T> emul (const Matrix<T>& a, const Matrix<T>& b)
{
        const unsigned int rw = a.rows(), cl = a.columns();
        BCHK(rw != b.rows(), MatErr, rows differ in emul, rw, TMatrix<T>(0,0));
        BCHK(cl != b.columns(), MatErr, columns differ in emul, cl, TMatrix<T>(0,0));
        TMatrix<T> res (rw, cl);
        do_vec_vec_mul<T> (a._DIM, res._VEC, a._VEC, b._VEC);
        return res;
}

INST(template <typename T> class TMatrix friend TMatrix<T> cemul (const Matrix<T>&, const Matrix<T>&);)
template <typename T>
TMatrix<T> cemul (const Matrix<T>& a, const Matrix<T>& b)
{
        const unsigned int rw = a .rows(), cl = a.columns();
        BCHK(rw != b.rows(), MatErr, rows differ in cemul, rw, TMatrix<T>(0,0));
        BCHK(cl != b.columns(), MatErr, columns differ in cemul, cl, TMatrix<T>(0,0));
        TMatrix<T> res (rw, cl);
        do_vec_vec_cmul<T> (a._DIM, res._VEC, a._VEC, b._VEC);
        return res;
}

INST(template <typename T> class TMatrix friend TMatrix<T> ediv (const Matrix<T>&, const Matrix<T>&);)
template <typename T>
TMatrix<T> ediv (const Matrix<T>& a, const Matrix<T>& b)
{
        const unsigned int rw = a .rows(), cl = a.columns();
        BCHK(rw != b.rows(), MatErr, rows differ in ediv, rw, TMatrix<T>(0,0));
        BCHK(cl != b.columns(), MatErr, columns differ in ediv, cl, TMatrix<T>(0,0));
        TMatrix<T> res (rw, cl);
        do_vec_vec_div<T> (a._DIM, res._VEC, a._VEC, b._VEC);
        return res;
}

INST(template <typename T> class TMatrix friend TMatrix<T> cediv (const Matrix<T>&, const Matrix<T>&);)
template <typename T>
TMatrix<T> cediv (const Matrix<T>& a, const Matrix<T>& b)
{
        const unsigned int rw = a .rows(), cl = a.columns();
        BCHK(rw != b.rows(), MatErr, rows differ in cediv, rw, TMatrix<T>(0,0));
        BCHK(cl != b.columns(), MatErr, columns differ in cediv, cl, TMatrix<T>(0,0));
        TMatrix<T> res (rw, cl);
        do_vec_vec_cdiv<T> (a._DIM, res._VEC, a._VEC, b._VEC);
        return res;
}


INST(template <typename T> class TMatrix friend TMatrix<T> ediv (const T&, TMatrix<T>);)
template <typename T>
inline TMatrix<T> ediv (const T& a, TMatrix<T> b)
{
        STD_SMP_TEMPLATE2C(val_div_vec, b._DIM, b._VEC, a);
        return b;
}
        
INST(template <typename T> class TMatrix friend TMatrix<T> ediv (const T&, const Matrix<T>&);)
template <typename T>
inline TMatrix<T> ediv (const T& a, const Matrix<T>& b)
{ 
        const unsigned r = b.rows(), c = b.columns();
        TMatrix<T> res (r, c);
        STD_SMP_TEMPLATE3VC(val_vec_div, b._DIM, res._VEC, b._VEC, a);
        return res;
}
        
INST(template <typename T> class TSMatrix friend TMatrix<T> ediv (const T&, const TSMatrix<T>&);)
template <typename T>
inline TMatrix<T> ediv (const T& a, const TSMatrix<T>& ts)
{
        TSMatrix<T> tm (ts); tm.detach ();
        STD_SMP_TEMPLATE4C(val_svc_div, tm._DIM, tm._VEC, ts._VEC, a, ts._FAC);
        tm._FAC = (T)1; return TMatrix<T> (tm);
}

#undef _DIM
#undef _VEC
#undef _ENDVEC
#undef _ROW
#undef _COL
#undef _FAC

NAMESPACE_END

#endif /* TBCI_VECTOR_EXTRA_H */