f_matrix.h

Go to the documentation of this file.

//-------------------------------------------------------------
// A M Bilgic, 1/97
// last update 97/10/04 AMB
//             98/05/07 AMB
//
// $Id: f_matrix.h,v 1.23.2.61 2022/11/03 17:28:10 garloff Exp $
//-------------------------------------------------------------

#ifndef TBCI_F_MATRIX_H
#define TBCI_F_MATRIX_H  

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

// Avoid -fguiding-decls
#if !defined(NO_GD) && !defined(AUTO_DECL)
# include "tbci/f_matrix_gd.h"
#endif

NAMESPACE_TBCI

#ifdef HAVE_GCC295_FRIEND_BUG
# define _VEC getvec()
# define _ENDVEC getendvec()
# define _DIM size()
# define _ROW rows()
# define _COL columns()
# define _FAC getfac()
#else
# define _VEC vec
# define _ENDVEC endvec
# define _DIM dim
# define _ROW row
# define _COL col
# define _FAC fac
#endif



template <typename T> class CSCMatrix;
template <typename T> class F_Matrix;
template <typename T> class F_TMatrix;
template <typename T> class F_TSMatrix;
template <typename T> class BdMatrix;
template <typename T> class Tensor;

template <typename T> class Matrix;
template <typename T> class TMatrix;

#ifdef PRAGMA_I
# pragma interface "f_matrix.h"
#endif

//-----------------------------------------------------
//-----------------------------------------------------

template <typename T>
class F_TMatrix : public Matrix_Sig<T>
{
 protected:
  typedef T mat_t /*TALIGN(sizeof(T))*/;
  unsigned long dim;
  unsigned int  col, row;
  T ** mat;
  T *vec, *endvec;

  int set_ptrs();

 public:
  typedef T value_type;
  typedef T element_type;
  typedef T aligned_value_type TALIGN(MIN_ALIGN2);

  friend class F_Matrix<T>;
  friend class F_TSMatrix<T>;
  friend class Tensor<T>;
  //friend class CSCMatrix<T>;

#ifdef HAVE_GCC295_FRIEND_BUG
  // Provide access to the internals to workaround gcc-2.95 hostility to friends
  T* const getvec() const  { return vec; }
  T* const getendvec() const  { return endvec; }
#endif

  //protected:
 public:
  explicit F_TMatrix (const unsigned = 0);
  F_TMatrix (const unsigned, const unsigned);
  F_TMatrix (const T&, const unsigned, const unsigned);
  F_TMatrix (const Vector<T>&, const enum rowcolvec = colvec);
  //aliasing
  F_TMatrix (const F_TMatrix<T>&);
  F_TMatrix (F_TSMatrix<T>&);
  //copying
  F_TMatrix (const F_Matrix<T>&);
  // real destructor
  void destroy ();

  //public:
  // empty destructor: F_TMatrix will never be destroyed but explicitly !!!
  ~F_TMatrix ()  {}

  // conversion
  operator TMatrix<T>();
  // allow instantiation (Matrix_Sig)
  /*virtual*/ static const char* mat_info ()  { return ("F_TMatrix"); }

  //public:
  // assignment
  F_TMatrix<T>& operator = (const F_Matrix<T>&);
  F_TMatrix<T>& operator = (const F_TMatrix<T>&);
  F_TMatrix<T>& operator = (F_TSMatrix<T>);
  F_TMatrix<T>& operator = (const T&);
  F_TMatrix<T>& alias (const F_Matrix<T>& m)
     { destroy (); dim = m.dim; col = m.col; row = m.row; 
       mat = m.mat; vec = m.vec; endvec = m.endvec; return *this; }

  F_TMatrix<T>& operator += (F_TMatrix<T>);       // arg will be destroyed
  F_TMatrix<T>& operator += (const F_Matrix<T>&);  // arg wont be   "  "
  F_TMatrix<T>& operator += (const T&);
  F_TMatrix<T>& operator -= (F_TMatrix<T>);
  F_TMatrix<T>& operator -= (const F_Matrix<T>&);
  F_TMatrix<T>& operator -= (const T&);

  F_TSMatrix<T> operator *= (const T&);
  F_TSMatrix<T> operator /= (const T&);
  
  // operations
  F_TMatrix<T>& operator - ();
  F_TMatrix<T>& herm ();  //returns the the hermitian adjungated (Tronsposed and conjugate complex)
  F_TMatrix<T>& conj ();
  F_TMatrix<T>& trans();

  F_TMatrix<T>& operator + (F_TMatrix<T>);              // arg will be destroyed
  F_TMatrix<T>& operator + (F_TSMatrix<T>);             // arg will be destroyed
  F_TMatrix<T>& operator + (const F_Matrix<T>&);        // arg wont be destroyed
  F_TMatrix<T>& operator + (const T&);
  F_TMatrix<T>& operator - (F_TMatrix<T>);
  F_TMatrix<T>& operator - (F_TSMatrix<T>);
  F_TMatrix<T>& operator - (const F_Matrix<T>&);
  F_TMatrix<T>& operator - (const T&);
  F_TSMatrix<T> operator * (const T&);
  F_TSMatrix<T> operator / (const T&);

  F_TMatrix<T> operator * (const F_Matrix<T>&);
  F_TMatrix<T> operator * (F_TMatrix<T>);
  F_TMatrix<T> operator * (F_TSMatrix<T>);
  // TODO: Implement transMult

  TVector<T> operator * (const Vector<T>& v)
    { F_Matrix<T> m (*this); return (m * v); }
  TVector<T> operator * (TVector <T>& tv)
    { F_Matrix<T> m (*this); return (m * tv); }

  // friend operations
#ifndef HAVE_GCC295_FRIEND_BUG
  friend NOINST F_TMatrix<T>  FRIEND_TBCI__ operator + FGD (const T&, F_TMatrix<T>);
  friend NOINST F_TMatrix<T>  FRIEND_TBCI__ operator + FGD (const T&, const F_Matrix<T>&);
  friend NOINST F_TMatrix<T>  FRIEND_TBCI__ operator - FGD (const T&, F_TMatrix<T>);
  friend NOINST F_TMatrix<T>  FRIEND_TBCI__ operator - FGD (const T&, const F_Matrix<T>&);
  friend NOINST F_TSMatrix<T> FRIEND_TBCI__ operator * FGD (const T&, F_TMatrix<T>);
  friend NOINST F_TSMatrix<T> FRIEND_TBCI__ operator * FGD (const T&, const F_Matrix<T>&);
#endif
  
  // Access
  T& operator () (const unsigned int, const unsigned int) const;
  //T operator () (const unsigned, const unsigned) const;
  T* get_fortran_matrix () const  { return vec; }
  const T& get (const unsigned r, const unsigned c) const  { return mat[c][r]; }

  bool operator == (const F_Matrix<T>& m);
  bool operator != (const F_Matrix<T>& m)
    { return !(*this == m); }

  bool operator == (F_TMatrix<T> tm)
    { F_Matrix<T> m (tm); return (*this == m); }
  bool operator != (F_TMatrix<T> tm)
    { return !(*this == tm); }

  bool operator == (F_TSMatrix<T>);
  bool operator != (F_TSMatrix<T> ts)
    { return !(*this == ts); }

  //public:
  unsigned int  columns () const { return col; }
  unsigned int  rows () const { return row; }
  unsigned long size () const { return dim; }
  //TVector<T> operator () (const unsigned int) const;
  TVector<T> get_row (const unsigned int) const;
  TVector<T> get_col (const unsigned int) const;
  void set_row (const Vector<T>&, const unsigned int);
  void set_col (const Vector<T>&, const unsigned int);
  void setval(const T val, const unsigned int r, const unsigned int c)
    {  this->operator() (r, c) = val;  }
  T& setval(const unsigned int r, const unsigned int c)
    {  return this->operator() (r, c);  }

  F_TMatrix<T>& resize (const unsigned int, const unsigned int);
  F_TMatrix<T>& resize (const unsigned int d)  { return this->resize (d, d); }
  F_TMatrix<T>& resize (const T&, const unsigned int, const unsigned int);
  F_TMatrix<T>& fill   (const T& = 0);
  F_TMatrix<T>& clear  ();
  F_TMatrix<T>& setunit(const T& = 1);

  F_TMatrix<T>& swap (F_TMatrix<T>&);
  F_TMatrix<T>  transposed_copy() const;
  F_TMatrix<T>& transpose();

  //friend double fabs FGD (const F_Matrix<T>&);
  //friend double fabs FGD (const F_TMatrix<T>&);
  double fabs () const;

  T trace () const;
  
}; 


template <typename T>
inline void F_TMatrix<T>::destroy ()
{       /* Skip if (dim) ? */
        if (dim) {
                if (vec) TBCIDELETE(T , vec, dim);
                if (mat) TBCIDELETE(T*, mat, col); 
                dim = 0;
        }
}

template <typename T>
F_TMatrix<T>::operator TMatrix<T>()
{
  TMatrix<T> tm(row, col);
  for (unsigned int c = 0; c < col; ++c)
    for (unsigned int r = 0; r < row; ++r)
      tm.setval(this->get(r, c), r, c);
  destroy();
  return tm;
}

template <typename T>
inline int F_TMatrix<T>::set_ptrs()
{
  if (mat && vec) {
    T* ptr = vec;
    for (unsigned int i=0; i<col; i++) {
        mat[i] = ptr;
        ptr += row;
    }
    endvec = ptr;
    return 0;
  } else {
    if (vec) { TBCIDELETE(T, vec, dim); vec = (T*) 0; }
    if (mat) { TBCIDELETE(T*,mat, col); mat = (T**)0; }
    dim = 0; col = 0; row = 0;
    endvec = (T*)0;
    return -1;
  }
}


template <typename T>
inline F_TMatrix<T>::F_TMatrix (const unsigned d)
  : dim((unsigned long)d*d), col(d), row(d),
    mat(d>0?NEW(T*,d):0), vec(d>0?NEW(T,d*d):0)
{
        set_ptrs();
}


template <typename T>
inline F_TMatrix<T>::F_TMatrix (const unsigned r, const unsigned c)
  : dim(((unsigned long)r)*c), col(c), row(r),
    mat(dim>0?NEW(T*,c):0), vec(dim>0?NEW(T,dim):0)
{
        set_ptrs();
}

template <typename T>
inline F_TMatrix<T>::F_TMatrix (const T& val, const unsigned r, const unsigned c)
  : dim(((unsigned long)r)*c), col(c), row(r),
    mat(dim>0?NEW(T*,c):0), vec(dim>0?NEW(T,dim):0)
{
        if (LIKELY(!set_ptrs()))
                TBCIFILL (vec, val, T, dim);
}


template <typename T>
inline F_TMatrix<T>::F_TMatrix (const Vector<T>& v, const enum rowcolvec r)
  : dim(v.dim), col(r==rowvec?1:v.dim), row(r==rowvec?v.dim:1),
    mat(v.dim?NEW(T*,col):0), vec(v.dim>0?NEW(T,v.dim):0) 
{
        if (LIKELY(!set_ptrs()))
                TBCICOPY (vec, v.vec, T, dim);
}

// data is duplicated
template <typename T>
inline F_TMatrix<T>::F_TMatrix (const F_Matrix<T>& m)
  : dim(m.dim), col(m.col), row(m.row),
    mat(m.dim>0?NEW(T*,m.col):0),
    vec(m.dim>0?NEW(T,m.dim):0)
{
        if (LIKELY(!set_ptrs()))
                TBCICOPY (vec, m.vec, T, dim);
}

// copy constructor only copies pointers !
template <typename T>
inline F_TMatrix<T>::F_TMatrix (const F_TMatrix<T>& tm)
  : dim (tm.dim), col (tm.col), row (tm.row), 
    mat (tm.mat), vec (tm.vec), endvec (tm.endvec) 
{}

template <typename T>
inline F_TMatrix<T>::F_TMatrix (F_TSMatrix<T>& ts)
  : dim (ts.dim), col (ts.col), row (ts.row)
{
  ts.eval (); vec = ts.vec; mat = ts.mat; 
  endvec = ts.endvec; ts.mut = false;
}


template <typename T>
inline T& F_TMatrix<T>::operator () (const unsigned int r, const unsigned int c) const
{
  BCHK(r>=row, MatErr, Illegal row index, r, mat[0][0]);
  BCHK(c>=col, MatErr, illegal col index, c, mat[0][0]);
  return mat[c][r];
}

/*
template <typename T>
inline T F_TMatrix<T>::operator () (const unsigned i, const unsigned j) const
{
  BCHK(i<0 || i>=row, MatErr, Illegal row index, i, mat[0][0]);
  BCHK(j<0 || j>=col, MatErr, illegal col index, j, mat[0][0]);
  return mat[j][i];
}
 */


//copy
template <typename T>
inline F_TMatrix<T>& F_TMatrix<T>::operator = (const F_Matrix<T>& a)
{
  BCHK(dim != a.dim, MatErr, Different sizes in assignment, a.dim, *this );
  TBCICOPY (vec, a.vec, T, dim);
  return *this;
}

//alias (do we need it?)
template <typename T>
inline F_TMatrix<T>& F_TMatrix<T>::operator = (const F_TMatrix<T>& a)
{
  BCHK(row != a.row || col != a.col, MatErr, Different sizes in assignment, a.dim, *this );
  if (vec != a.vec) /* Is == possible ? */
  {
        this->destroy();
        dim = a.dim; col = a.col; row = a.row;
        mat = a.mat; vec = a.vec; endvec = a.endvec;
  }
  return *this;
}

template <typename T>
inline F_TMatrix<T>& F_TMatrix<T>::operator = (F_TSMatrix<T> a)
{
  BCHK(col != a.col || row != a.row, MatErr, Different sizes in assignment, a.dim, *this );
  a.eval (this);
  return *this;
}

//fill
template <typename T>
inline F_TMatrix<T>& F_TMatrix<T>::operator = (const T& val)
{
  return this->fill (val);
}


template <typename T>
inline F_TMatrix<T>& F_TMatrix<T>::setunit (const T& fac)
{
  TBCIFILL (vec, 0, T, dim);
  BCHK(col!=row, MatErr, setunit makes only sense on quadratic matrices, col, *this);
  const unsigned iend = MIN(row,col);
  for (REGISTER unsigned i = 0; i < iend; i++) mat[i][i] = fac;
  return *this;
}

template <typename T>
inline F_TMatrix<T>& F_TMatrix<T>::clear ()
{
  if (dim) TBCIFILL (vec, 0, T, dim);
  return *this;
}

#if 0
template <typename T>
inline TVector<T> F_TMatrix<T>::operator () (const unsigned int i) const
{
  TVector<T> v(col);
  BCHK(i<0 || i>= row, MatErr, Illegal row index, i, v);
  TBCICOPY (v.vec, mat[i], T, col);
  return v;
}
#endif

template <typename T>
inline void F_TMatrix<T>::set_col (const Vector<T>& v, const unsigned int c)
{
  BCHKNR(c >= col, MatErr, Illegal col index in set_col, c);
  BCHKNR(v.dim != (unsigned long)row, MatErr, Different sizes in set_row, v.dim);
  TBCICOPY(mat[c], v.vec, T, v.dim);
}

template <typename T>
inline TVector<T> F_TMatrix<T>::get_col (const unsigned int c) const  
{ 
  TVector<T> v(row);
  //BCHK(c >= col, MatErr, Illegal col index in get_col, c, v);
  TBCICOPY (v.vec, mat[c], T, row);
  return v; 
}

template <typename T>
inline void F_TMatrix<T>::set_row (const Vector<T>& v, const unsigned int r)
{
  BCHKNR(r >= row, MatErr, Illegal row index in set_row, r);
  for (REGISTER unsigned int j=0; j<col; j++) 
    mat[j][r] = v.vec[j];
}

template <typename T>
inline TVector<T> F_TMatrix<T>::get_row (const unsigned int r) const
{
  TVector<T> v(col);
  BCHK(r >= row, MatErr, Illegal row index in get_row, r, v);
  for (REGISTER unsigned int j=0; j<col; j++) 
    v.vec[j] = mat[j][r];
  return v;
}


#ifndef LAPACK_INLINE
# define LAPACK_INLINE  
#endif
template <typename T>
LAPACK_INLINE F_TMatrix<T>& F_TMatrix<T>::resize (const unsigned int r, const unsigned int c)
{
  if (r == row && c == col) 
    return *this;
  if (mat) 
    TBCIDELETE(T*,mat,col);
  T* vv = vec;
  unsigned long dd =dim;
  dim = ((unsigned long)r)*c;
  col = c; row = r;
  if (dim) { 
        mat = NEW(T*,c); vec = NEW(T,dim);
  } else { 
        mat = (T**)0; vec = (T*)0; endvec = (T*)0;
  }
  if (!set_ptrs())
        TBCICOPY(vec, vv, T, MIN(dim,dd));
  if (dd && vv)
        TBCIDELETE(T, vv, dd);
  return *this;
}

template <typename T>
F_TMatrix<T>& F_TMatrix<T>::resize (const T& v, const unsigned int r, const unsigned int c)
{
  if (r == row && c == col) {
        TBCIFILL(vec, v, T, dim);
        return *this;
  }
  if (mat) { 
          TBCIDELETE(T*,mat,col);
          TBCIDELETE(T,vec, dim);
  }
  dim = r*c;
  row = r; col = c;
  if (dim) { 
        mat = NEW(T*,c); vec = NEW(T,dim);
  } else { 
        mat = (T**)0; vec = (T*)0; endvec = (T*)0; 
  }
  if (LIKELY(!set_ptrs()))
        TBCIFILL (vec, v, T, dim);
  return *this;
}


template <typename T>
inline F_TMatrix<T>& F_TMatrix<T>::fill (const T& val)
{
  TBCIFILL(vec, val, T, dim);   
  return *this;
}

template <typename T>
inline T F_TMatrix<T>::trace () const
{
   ALIGN3(T tr, vec[0], MIN_ALIGN2);
   BCHK (row != col, MatErr, Trace over non quadratic matrix, 0, tr=0);
   for (REGISTER unsigned i = 1; i < row; i++) 
        tr += mat[i][i];
   return tr;
}

/************* Arithmetics *****************/ 


template <typename T>
F_TMatrix<T>& F_TMatrix<T>::trans()
{
  if (rows() == columns())
  {
    for(unsigned int i=0; i<rows(); i++)
      for(unsigned int j=i+1; j<columns(); j++)
        SWAP(mat[i][j], mat[j][i]);
  } else
  {
    STD__ cout << "Not implemented!" << STD__ endl;
  }
  return (*this);
}


template <typename T>
F_TMatrix<T>& F_TMatrix<T>::herm ()
{
  STD__ cout << "WRONG: F_TMatrix<T>::herm()" << STD__ endl;
  if (rows() == columns())
  {
    for(unsigned int i=0; i<rows(); i++)
    {
//      mat[i][i]=(mat[i][i]).conj();
      for(unsigned int j=i+1; j<columns(); j++)
      {
//      mat[i][j]=(mat[i][j]).conj();
//      mat[j][i]=(mat[j][i]).conj();
        SWAP(mat[i][j], mat[j][i]);
      }
    }
  } else
  {
    STD__ cout << "Not implemented!" << STD__ endl;
  }
  return (*this);
}


template <typename T>
F_TMatrix<T>& F_TMatrix<T>::conj ()
{ 
  STD__ cout << "Wrong: F_TMatrix<T>::conf()" << STD__ endl;
//  for (REGISTER T* ptr=vec; ptr<endvec; ptr++) *ptr = (*ptr).conj();
  return *this;
}

// unary -
template <typename T>
inline F_TMatrix<T>& F_TMatrix<T>::operator - ()
{
  for (REGISTER T* ptr=vec; ptr<endvec; ptr++) *ptr = -(*ptr);
  return *this;
}

template <typename T>
bool F_TMatrix<T>::operator == (const F_Matrix<T>& m)
{
  F_Matrix<T> th (*this);
  if (row != m.row || col != m.col) return false;
  if (vec == m.vec || dim == 0) return true;
  if (TBCICOMP (vec, m.vec, T, dim)) return false;
  /*else*/ return true;
}

template <typename T>
bool F_TMatrix<T>::operator == (F_TSMatrix<T> ts)
{
  F_Matrix<T> th (*this);
  if (row != ts.row || col != ts.col) return (ts.destroy(), false);
  if (dim == 0) { /*ts.destroy();*/ return true; }
  if (ts.fac == (T)1)
    {
      if (vec == ts.vec) return true;
      if (TBCICOMP (vec, ts.vec, T, dim)) { ts.destroy (); return false; }
    }
  else
    {
      if (vec == ts.vec) return false;
      for (REGISTER T *p1 = vec, *p2 = ts.vec; p1 < endvec; p1++, p2++)
        if (*p1 != *p2 * ts.fac) { ts.destroy (); return false; }
    }
  ts.destroy (); return true;
}


#define F_TMFORALL_M(op)                                                \
template <typename T>                                                   \
inline F_TMatrix<T>& F_TMatrix<T>::operator op (const F_Matrix<T>& a)   \
{                                                                       \
  BCHK(row != a.row, MatErr, number of rows differ in op, a.row, *this ); \
  BCHK(col != a.col, MatErr, number of cols differ in op, a.col, *this ); \
  for (REGISTER T *p1 = vec, *p2 = a.vec; p1 < endvec; p1++, p2++) *p1 op *p2; \
  return *this;                                                         \
}

F_TMFORALL_M(+=)
F_TMFORALL_M(-=)


// Changing a F_TMatrix to a F_Matrix destroys it ...
#define F_TMFORALL_TM(op)                                       \
template <typename T>                                           \
inline F_TMatrix<T>& F_TMatrix<T>::operator op (F_TMatrix<T> a)\
{                                                               \
  return this->operator op (F_Matrix<T>(a));                    \
}

F_TMFORALL_TM(+=)
F_TMFORALL_TM(-=)


#define F_TMFORALL_TS(op)                                               \
template <typename T>                                                   \
inline F_TMatrix<T>& F_TMatrix<T>::operator op (F_TSMatrix<T> ts)       \
{                                                                       \
  BCHK(row != ts.row, MatErr, number of rows differ in op, ts.row, (ts.destroy(),*this) ); \
  BCHK(col != ts.col, MatErr, number of cols differ in op, ts.col, (ts.destroy(),*this) ); \
  for (REGISTER T *p1 = vec, *p2 = ts.vec; p1 < endvec; p1++, p2++)     \
     *p1 op##= *p2 * ts.fac;                                            \
  ts.destroy (); return *this;                                          \
}

F_TMFORALL_TS(+)
F_TMFORALL_TS(-)


//Change these to operate on diag only ?
#define F_TMFORALL_T(op)                                        \
template <typename T>                                           \
inline F_TMatrix<T>& F_TMatrix<T>::operator op (const T& a)     \
{                                                               \
  for (REGISTER T* ptr = vec; ptr < endvec; ptr++) *ptr op a;   \
  return *this;                                                 \
}

F_TMFORALL_T(+=)
F_TMFORALL_T(-=)
//F_TMFORALL_T(*=)
//F_TMFORALL_T(/=)


template <typename T>
inline F_TSMatrix<T> F_TMatrix<T>::operator *= (const T& a)
{ return F_TSMatrix<T> (*this, a); }

template <typename T>
inline F_TSMatrix<T> F_TMatrix<T>::operator /= (const T& a)
{ 
  BCHK(a==(T)0, MatErr, F_Matrix division by zero, 0, F_TSMatrix<T> (*this));
  return F_TSMatrix<T> (*this, (T)1/a); 
}


// temporaries can be changed, because we can't refer them
#define F_STDDEF_TMM(op)                                                \
template <typename T>                                                   \
inline F_TMatrix<T>& F_TMatrix<T>::operator op (const F_Matrix<T>& a)   \
{ return this->operator op##= (a); }

F_STDDEF_TMM(+)
F_STDDEF_TMM(-)


#define F_STDDEF_TMTM(op)                                       \
template <typename T>                                           \
inline F_TMatrix<T>& F_TMatrix<T>::operator op (F_TMatrix<T> a) \
{ return this->operator op##= (F_Matrix<T>(a)); }

F_STDDEF_TMTM(+)
F_STDDEF_TMTM(-)


#define F_STDDEF_TMT(op)                                        \
template <typename T>                                           \
inline F_TMatrix<T>& F_TMatrix<T>::operator op (const T& a)     \
{ return this->operator op##= (a); }

F_STDDEF_TMT(+)
F_STDDEF_TMT(-)
//STDDEF_TMT(*)

template <typename T>
inline F_TSMatrix<T> F_TMatrix<T>::operator * (const T& a)
{ return F_TSMatrix<T> (*this, a); }

template <typename T>
inline F_TSMatrix<T> F_TMatrix<T>::operator / (const T& b)
{
  BCHK((b == (T)0), MatErr, F_Matrix division by 0, 0, F_TSMatrix<T> (*this)); 
  return F_TSMatrix<T> (*this, (T)1/b);
}

//friends

//Change to operate on diag only ?
#define F_STDDEF_TTM(op)                                        \
INST(template <typename T> class F_TMatrix friend F_TMatrix<T> operator op (const T&, F_TMatrix<T>);)           \
template <typename T>                                           \
inline F_TMatrix<T> operator op (const T& a, F_TMatrix<T> b)    \
{                                                               \
  for (REGISTER T* ptr=b._VEC; ptr<b._ENDVEC; ptr++)            \
        *ptr = a op *ptr;                                       \
  return b;                                                     \
}

F_STDDEF_TTM(+)
F_STDDEF_TTM(-)
//STDDEF_TTM(*)
//STDDEF_TTM(/)

INST(template <typename T> class F_TMatrix friend F_TSMatrix<T> operator * (const T&, F_TMatrix<T>);)
template <typename T>
inline F_TSMatrix<T> operator * (const T& a, F_TMatrix<T> b)
{ return F_TSMatrix<T> (b, a); } 
// Commutativity is assumed here !


//Change to operate on diag only ?
#define F_STDDEF_TM(op)                                 \
INST(template <typename T> class F_TMatrix friend F_TMatrix<T> operator op (const T&, const F_Matrix<T>&);)                                     \
template <typename T>                                   \
inline F_TMatrix<T> operator op (const T& a, const F_Matrix<T>& b)      \
{                                                       \
  F_TMatrix<T> c (b._ROW, b._COL);                      \
  for (REGISTER T *p1=c._VEC, *p2=b._VEC; p1<c._ENDVEC; p1++, p2++)     \
    *p1 = a op *p2;                                     \
  return c;                                             \
}

F_STDDEF_TM(+)
F_STDDEF_TM(-)
//STDDEF_TM(*)
//STDDEF_TM(/)


INST(template <typename T> class F_TMatrix friend F_TSMatrix<T> operator * (const T&, const F_Matrix<T>&);)
template <typename T>
inline F_TSMatrix<T> operator * (const T& a, const F_Matrix<T>& b)
{ return F_TSMatrix<T> (b, a); }
// Commutativity is assumed here !


template <typename T>
inline double F_TMatrix<T>::fabs () const 
{ F_Matrix<T> m (*this); return m.fabs (); }

//-------------------------------------------------------------
//-------------------------------------------------------------
template <typename T>
class F_TSMatrix : public Matrix_Sig<T> //: public F_TMatrix
{
 protected:
  unsigned long dim;
  unsigned int col, row;
  T** mat;
  T *vec, *endvec; 

  void destroy (); 
  T fac ALIGN(MIN_ALIGN); //factor
  bool mut;  //mutability: does vec point to changeable (non-refered) data ?
  
  void clone (bool = false, F_TMatrix<T>* = 0);
#ifdef HAVE_GCC295_FRIEND_BUG
 public:
  T* const getvec () const     { return vec; }
  T* const getendvec () const  { return endvec; }
  T& getfac ()  { return fac; }
  const T& getfac () const     { return fac; }
#endif
  void detach (F_TMatrix<T>* = 0);

 public:
  typedef T value_type;
  typedef T element_type;
  typedef T aligned_value_type TALIGN(MIN_ALIGN2);
  
  friend class F_TMatrix<T>;
  friend class F_Matrix<T>; 

  F_TSMatrix () : dim(0), col(0), row(0), mat((T**)0), vec((T*)0), endvec((T*)0) {}
  // empty destructor: F_TSMatrix will never be destroyed but explicitly !!!
  ~F_TSMatrix ()  {}
  
  F_TSMatrix (const F_TMatrix<T>& tm, const T& f = (T)1) 
    : dim(tm.dim), col(tm.col), row(tm.row),
      mat(tm.mat), vec(tm.vec), endvec(tm.endvec), fac(f), mut(true)  {}
  F_TSMatrix (const F_Matrix<T>& m, const T& f = (T)1)
    : dim(m.dim), col(m.col), row(m.row),
      mat(m.mat), vec(m.vec), endvec(m.endvec), fac(f), mut(false)  {}
  F_TSMatrix (const F_TSMatrix<T>& ts)
    : dim(ts.dim), col(ts.col), row(ts.row),
      mat(ts.mat), vec(ts.vec), endvec(ts.endvec), fac(ts.fac), mut(ts.mut)  {}

  T operator () (const unsigned r, const unsigned c)
    { T val = mat[c][r] * fac; destroy (); return val; }
  T get (const unsigned r, const unsigned c) const  { return mat[c][r] * fac; }
  
  F_TSMatrix<T>& eval (F_TMatrix<T>* = 0);
  operator TMatrix<T>();

  // allow instantiation (Matrix_Sig)
  /*virtual*/ static const char* mat_info ()  { return ("F_TSMatrix"); }

  // These operations are very cheap now
  F_TSMatrix<T>& operator = (const F_TSMatrix<T>& ts)
     { destroy (); dim = ts.dim; col = ts.col; row = ts.row;
       mat = ts.mat; vec = ts.vec; endvec = ts.endvec; 
       fac = ts.fac; mut = ts.mut; return *this; }
  F_TSMatrix<T>& operator = (const F_TMatrix<T>& tm)
     { destroy (); dim = tm.dim; col = tm.col; row = tm.row;
       mat = tm.mat; vec = tm.vec; endvec = tm.endvec; 
       fac = (T)1; mut = true; return *this; }
   
  F_TSMatrix<T>& operator *= (const T& f)  { fac *= f; return *this; }
  F_TSMatrix<T>& operator /= (const T& f)  { fac /= f; return *this; }
  F_TSMatrix<T>& operator * (const T& f)  { fac *= f; return *this; }
  F_TSMatrix<T>& operator / (const T& f)  { fac /= f; return *this; }

  F_TSMatrix<T>& operator - ()  { fac *= (T)(-1); return *this; } 
  
  // Now for the more complicated
  F_TMatrix<T> operator + (const F_Matrix<T>&);
  F_TMatrix<T> operator + (const F_TMatrix<T>&);
  F_TMatrix<T> operator + (F_TSMatrix<T>);
  F_TMatrix<T> operator + (const T&);
  F_TMatrix<T> operator - (const F_Matrix<T>&);
  F_TMatrix<T> operator - (const F_TMatrix<T>&);
  F_TMatrix<T> operator - (F_TSMatrix<T>);
  F_TMatrix<T> operator - (const T&);

#ifndef HAVE_GCC295_FRIEND_BUG
  friend NOINST F_TSMatrix<T> FRIEND_TBCI__ operator * FGD (const T&, F_TSMatrix<T>);
  friend NOINST F_TMatrix<T> FRIEND_TBCI__ operator + FGD (const T&, F_TSMatrix<T>);
  friend NOINST F_TMatrix<T> FRIEND_TBCI__ operator - FGD (const T&, F_TSMatrix<T>);
#endif
              
  F_TMatrix<T> operator * (const F_Matrix<T>& m);
  F_TMatrix<T> operator * (const CSCMatrix<T>& m);
  TVector<T> operator * (const Vector<T>& v);

  bool operator == (const F_Matrix<T>&);
  bool operator != (const F_Matrix<T>& m) { return !(*this == m); }

  bool operator == (F_TSMatrix<T>);
  bool operator != (F_TSMatrix<T> ts) { return !(*this == ts); }

  bool operator == (F_TMatrix<T> tm) { F_Matrix<T> m(tm); return (*this == m); }
  bool operator != (F_TMatrix<T> tm) { return !(*this == tm); }


  //Friends
  //friend double fabs FGD (F_TSMatrix<T>&);
  double fabs ();
};


template <typename T>
F_TMatrix<T> F_TSMatrix<T>::operator * (const CSCMatrix<T>& m2)
{       
        F_TMatrix<T> m1(*this);
        F_TMatrix<T> erg(0, m1.rows(), m2.columns());
        for (unsigned int i=0; i<m1.rows(); i++)
                for (unsigned int j=0; j<m2.columns(); j++)
                {
                        T val = (T)0;
      // FIXME: Do we multiply by fac twice here?
                        for (unsigned int x=m2.pcol[j]; x<m2.pcol[j+1]; x++) 
        val += m1(i, m2.irow[x] )*m2.comp[x];
                        erg.setval(fac*val, i, j);
                }

        // Really destroy m1!
        mut = true;
        destroy();
        return erg;
}



template <typename T>
inline void F_TSMatrix<T>::destroy ()
{
  if (!mut || !dim) return;
  if (vec) { TBCIDELETE(T,vec,dim); vec = (T*)0; }
  endvec = (T*)0; dim = 0;
  if (mat) TBCIDELETE(T*,mat,col);
}
  
template <typename T>
inline void F_TSMatrix<T>::detach (F_TMatrix<T>* tm)
{
  if ((!tm && mut) || !dim) return;
  if (!tm) //mut = false
     {
        vec = NEW(T ,dim); if (!vec) { dim = 0; mat = (T**)0; endvec = (T*)0; return; }
        mat = NEW(T*,col); if (!mat) { TBCIDELETE(T,vec,dim); dim = 0; vec = (T*)0; endvec = (T*)0; return; }
        for (REGISTER unsigned c = 0; c < row; c++) mat[c] = vec + c*row;
        endvec = vec + dim; 
     }
   else //tm = true, mut = ?
     {
        vec = tm->vec; mat = tm->mat; endvec = tm->endvec;
     }
   mut = true;
}

template <typename T>
void F_TSMatrix<T>::clone (bool evl, F_TMatrix<T>* tm)
{
  if ((mut && !tm)|| !dim) return;
  T* oldv = vec; T** oldm = mat; bool omut = mut;
  detach (tm);
  if (evl && fac != (T)1)
    {
      for (REGISTER T *p1 = vec, *p2 = oldv; p1 < endvec; p1++, p2++)
        *p1 = *p2 * fac;
      fac = (T)1;
    }
  else TBCICOPY (vec, oldv, T, dim);
  if (tm && omut)
     { TBCIDELETE(T,oldv,dim); TBCIDELETE(T*,oldm,col); }
}



//creates mutable evaluated F_TSMatrix (can be easily converted to F_TMatrix)
template <typename T>
inline F_TSMatrix<T>& F_TSMatrix<T>::eval (F_TMatrix<T>* tm)
{
  if (mut && !tm)
    {
      if (fac != (T)1)
        {
          for (REGISTER T* ptr = vec; ptr < endvec; ptr++)
            *ptr *= fac;
          fac = (T)1;
        }
    }
  else
     clone (true, tm);
   return (*this);
}
          

// Arithmetics

#define F_STDDEF_TSM(op)                                                \
template <typename T>                                                   \
inline F_TMatrix<T> F_TSMatrix<T>::operator op (const F_Matrix<T>& m)   \
{                                                                       \
  BCHK(row != m.row || col != m.col, MatErr, Op op on diff size Mats, m.row, F_TMatrix<T> (*this)); \
  F_TSMatrix<T> ts (*this); ts.detach ();                               \
  for (REGISTER T *p1 = ts.vec, *p2 = vec, *p3 = m.vec; p2 < endvec; p1++, p2++, p3++) \
    *p1 = *p2 * fac op *p3;                                             \
  ts.fac = (T)1; return F_TMatrix<T> (ts);                              \
}

F_STDDEF_TSM(+)
F_STDDEF_TSM(-)


#define F_STDDEF_TSTM(op)                                               \
template <typename T>                                                   \
inline F_TMatrix<T> F_TSMatrix<T>::operator op (const F_TMatrix<T>& tm) \
{                                                                       \
  BCHK(row != tm.row || col != tm.col, MatErr, Op op on diff size Mats, tm.row, tm); \
  for (REGISTER T *p1 = tm.vec, *p2 = vec; p2 < endvec; p1++, p2++)     \
    *p1 = *p2 * fac op *p1;                                             \
  destroy (); return tm;                                                \
}

F_STDDEF_TSTM(+)
F_STDDEF_TSTM(-)


#define F_STDDEF_TSTS(op)                                               \
template <typename T>                                                   \
inline F_TMatrix<T> F_TSMatrix<T>::operator op (F_TSMatrix<T> ts)       \
{                                                                       \
  BCHK(row != ts.row || col != ts.col, MatErr, Op op on diff size Mats, ts.row, F_TMatrix<T> (*this)); \
  F_TSMatrix<T> tm;                                                     \
  if (!mut && ts.mut) tm = ts;                                          \
  else { tm = *this; tm.detach (); }                                    \
  for (REGISTER T *p1 = tm.vec, *p2 = vec, *p3 = ts.vec; p2 < endvec; p1++, p2++, p3++) \
    *p1 = *p2 * fac op *p3 * ts.fac;                                    \
  if (!mut && ts.mut) destroy (); else ts.destroy ();                   \
  tm.fac = (T)1; return F_TMatrix<T> (tm);                              \
}

F_STDDEF_TSTS(+)
F_STDDEF_TSTS(-)


#define F_STDDEF_TST(op)                                        \
template <typename T>                                           \
inline F_TMatrix<T> F_TSMatrix<T>::operator op (const T& a)     \
{                                                               \
  F_TSMatrix<T> ts (*this); ts.detach ();                       \
  for (REGISTER T *p1 = ts.vec, *p2 = vec; p2 < endvec; p1++, p2++) \
    *p1 = *p2 * fac op a;                                       \
  ts.fac = (T)1; return F_TMatrix<T> (ts);                      \
}

F_STDDEF_TST(+)
F_STDDEF_TST(-)

// operate on diag only?
#define F_STDDEF_TTS(op)                                        \
INST(template <typename T> class F_TSMatrix friend F_TMatrix<T> operator op (const T&, F_TSMatrix<T>);)                                         \
template <typename T>                                           \
inline F_TMatrix<T> operator op (const T& a, F_TSMatrix<T> ts)  \
{                                                               \
  F_TSMatrix<T> tm (ts); tm.detach ();                          \
  for (REGISTER T *p1 = tm._VEC, *p2 = ts._VEC; p1 < tm._ENDVEC; p1++, p2++) \
    *p1 = a op *p2 * ts._FAC;                                   \
  tm._FAC = (T)1; return F_TMatrix<T> (tm);                     \
}
  
F_STDDEF_TTS(+)
F_STDDEF_TTS(-)
   

template <typename T>
F_TMatrix<T> F_TSMatrix<T>::operator * (const F_Matrix<T>& m)
{
  F_TMatrix<T> ftm (row, m.col);
  BCHK(col != m.row, MatErr, Illegal sizes in F_Matrix multiplication, 0, (destroy(), ftm));
  REGISTER T tmp ALIGN(MIN_ALIGN);
  for (unsigned int r=0; r<row; r++) {
    for (unsigned int c=0; c<m.col; c++) {
      tmp = mat[0][r] * m(0, c);
      for (REGISTER unsigned int l=1; l<col; l++)
        tmp += mat[l][r] * m(l, c);
      ftm.setval(tmp*fac, r, c);
    }
  }
  destroy (); return ftm;
}


template <typename T>
TVector<T> F_TSMatrix<T>::operator * (const Vector<T>& v)
{
  TVector<T> tv (row);
  BCHK (v.size() != col, MatErr, multiply F_Matrix w/ Vector: wrong config, v.size(), (destroy(), tv));
  /* FIXME: We could call do_mat_vec_transmult<>() here */
  for (unsigned int r=0; r<row; r++)
    {
      ALIGN3(REGISTER T el, mat[0][r] * v(0), MIN_ALIGN2);
      for (REGISTER unsigned int c=1; c<col; c++) 
        el += mat[c][r] * v(c);
      tv.set (el * fac, r);
    }
  destroy (); return tv;
}


template <typename T>
bool F_TSMatrix<T>::operator == (const F_Matrix<T>& m)
{
  if (row != m.row || col != m.col) { destroy (); return false; }
  if (dim == 0) { /*destroy();*/ return true; }
  if (fac == (T)1)
    {
      if (vec == m.vec) return true;
      if (TBCICOMP (vec, m.vec, T, dim)) { destroy (); return false; }
    }
  else
    {
      if (vec == m.vec) return false;
      for (REGISTER T *p1 = vec, *p2 = m.vec; p1 < endvec; p1++, p2++)
        if (*p1 * fac != *p2) { destroy (); return false; }
    }
  destroy (); return true;
}


template <typename T>
bool F_TSMatrix<T>::operator == (F_TSMatrix<T> ts)
{
  if (row != ts.row || col != ts.col) { destroy (); ts.destroy (); return false; }
  if (dim == 0) {/*destroy(); ts.destroy();*/ return true; }
  if (fac == ts.fac)
    {
      if (vec == ts.vec) { destroy (); return true; }
      if (TBCICOMP (vec, ts.vec, T, dim)) { destroy (); ts.destroy (); return false; }
    }
  else
    {
      if (vec == ts.vec) { destroy (); return false; }
      for (REGISTER T *p1 = vec, *p2 = ts.vec; p1 < endvec; p1++, p2++)
        if (*p1 * fac != *p2) { destroy (); ts.destroy (); return false; }
    }
  destroy (); ts.destroy (); return true;
}

template <typename T>
F_TSMatrix<T>::operator TMatrix<T>()
{
  TMatrix<T> tm(this->row, this->col);
  for (unsigned int c = 0; c < this->col; ++c)
    for (unsigned int r = 0; r < this->row; ++r)
      tm.setval(this->get(r, c), r, c);
  destroy();
  return tm;
}


template <typename T>
double F_TSMatrix<T>::fabs ()
{
   REGISTER T* ptr = vec;
   double REGISTER rv ALIGN(8) = fabssqr (*ptr++);
   for (; ptr < endvec; ptr++)
     rv += fabssqr (*ptr);
   rv *= fabssqr (fac);
   destroy (); return MATH__ sqrt (rv);
}


INST(template <typename T> class F_TSMatrix friend F_TSMatrix<T> operator * (const T&, F_TSMatrix<T>);)
template <typename T>
inline F_TSMatrix<T> operator * (const T& f, F_TSMatrix<T> ts)
{ ts._FAC = f * ts._FAC; return ts; }

//---------------------------------------------------------------
// Real, user referable, matrix
//---------------------------------------------------------------
template <typename T>
class F_Matrix : public F_TMatrix<T> 
{
  
 public:

  typedef T value_type;
  typedef T element_type;
  typedef T aligned_value_type TALIGN(MIN_ALIGN2);

  friend class F_TSMatrix<T>;

  // constructor, destructor
  explicit F_Matrix (const unsigned d=0) : F_TMatrix<T> (d) {}
  F_Matrix (const unsigned r, const unsigned c) : F_TMatrix<T> (r, c) {}
  F_Matrix (const T& v, const unsigned r, const unsigned c) : F_TMatrix<T> (v, r, c) {}
  //F_Matrix (const Vector<T>& v, const enum rowcolvec r = colvec) 
  //  : F_TMatrix<T> (v, r) {}

  //copy
  F_Matrix (const F_Matrix<T>& m) : F_TMatrix<T> (m) {}
  F_Matrix (const CSCMatrix<T>& m);

  //alias
  F_Matrix (const F_TMatrix<T>& tm) : F_TMatrix<T> (tm) {}
  F_Matrix (F_TSMatrix<T>& ts) : F_TMatrix<T> (ts) {}

  ~F_Matrix () { this->destroy (); }

  // Conversion
  F_Matrix (const Matrix<T>& m);
  operator TMatrix<T>() const;

  // allow instantiation (Matrix_Sig)
  /*virtual*/ static const char* mat_info ()  { return ("F_Matrix"); }
  
  // Access
  T& operator () (const unsigned int, const unsigned int) const;
  //T operator () (const unsigned, const unsigned) const;
  T* get_fortran_matrix() const { return this->vec; }
  const T& get (const unsigned r, const unsigned c) const  
    { return this->mat[c][r]; }

  unsigned int  columns () const  { return this->row; }
  unsigned int  rows    () const  { return this->col; }
  unsigned long size    () const  { return this->dim; }
  //TVector<T> operator () (const unsigned int) const;
  TVector<T> get_row (const unsigned int) const;
  TVector<T> get_col (const unsigned int) const;
  void set_row (const Vector<T>&, const unsigned int);
  void set_col (const Vector<T>&, const unsigned int);
  void setval(const T val, const unsigned int r, const unsigned int c)
    {  this->operator() (r, c) = val;  }

  // basics  
  F_Matrix<T>& operator = (const F_Matrix<T>& m)
    { this->F_TMatrix<T>::operator = (m); return *this; }
  F_Matrix<T>& operator = (const F_TMatrix<T>& tm)
    { this->F_TMatrix<T>::operator = (tm); return *this; }
  F_Matrix<T>& operator = (F_TSMatrix<T> ts)
    { F_TMatrix<T>(this->F_TMatrix<T>::operator = (ts)); return *this; }
  F_Matrix<T>& operator = (const T& v)
    { this->F_TMatrix<T>::operator = (v); return *this; }

  F_Matrix<T>& resize (const unsigned r, const unsigned c)
    { this->F_TMatrix<T>::resize (r, c); return *this; }
  F_Matrix<T>& resize (const unsigned d)  
    { return this->resize (d, d); }
  F_Matrix<T>& resize (const T& v, const unsigned r, const unsigned c)
    { this->F_TMatrix<T>::resize (v, r, c); return *this; }
  F_Matrix<T>& fill (const T& v = 0)
    { this->F_TMatrix<T>::fill (v); return *this; }
  F_Matrix<T>& setunit (const T& f = 1)
    { this->F_TMatrix<T>::setunit (f); return *this; }
  F_Matrix<T>& clear ()
    { this->F_TMatrix<T>::clear (); return *this; }

  // assignment ops (do we really need em ?)
  // Seems we have to make sure, that a F_Matrix and no F_TMatrix
  // is returned; otherwise there's danger to be destroyed
  F_Matrix<T>& operator += (const F_Matrix<T>& a) 
    { this->F_TMatrix<T>::operator += (a); return *this; } 
  F_Matrix<T>& operator -= (const F_Matrix<T>& a)
    { this->F_TMatrix<T>::operator -= (a); return *this; } 
  F_Matrix<T>& operator += (F_TMatrix<T> a) 
    { this->F_TMatrix<T>::operator += (a); return *this; } 
  F_Matrix<T>& operator -= (F_TMatrix<T> a)
    { this->F_TMatrix<T>::operator -= (a); return *this; } 
  F_Matrix<T>& operator += (const T& a)
    { this->F_TMatrix<T>::operator += (a); return *this; } 
  F_Matrix<T>& operator -= (const T& a)
    { this->F_TMatrix<T>::operator -= (a); return *this; } 
  F_Matrix<T>& operator *= (const T& a)
    { this->F_TMatrix<T>::operator *= (a); return *this; } 
  F_Matrix<T>& operator /= (const T& a)
    { this->F_TMatrix<T>::operator /= (a); return *this; } 

  
  // basic operations
  F_TMatrix<T> operator - () const { return -(F_TMatrix<T> (*this)); }
  
  F_TMatrix<T> operator + (const F_Matrix<T>&) const;
  F_TMatrix<T> operator - (const F_Matrix<T>&) const;
  F_TMatrix<T> operator * (const F_Matrix<T>&) const;

  F_TMatrix<T> operator + (F_TMatrix<T>) const;
  F_TMatrix<T> operator - (F_TMatrix<T>) const;
  F_TMatrix<T> operator * (F_TMatrix<T>) const;

  F_TMatrix<T> operator + (F_TSMatrix<T>) const;
  F_TMatrix<T> operator - (F_TSMatrix<T>) const;
  F_TMatrix<T> operator * (F_TSMatrix<T>) const;

  F_TMatrix<T> operator + (const T&) const;
  F_TMatrix<T> operator - (const T&) const;
  F_TSMatrix<T> operator * (const T&) const;
  F_TSMatrix<T> operator / (const T&) const;

  TVector<T> operator * (const Vector<T>&) const;
  TVector<T> operator * (TVector<T>&) const;

  bool operator == (const F_Matrix<T>& m) const;
  bool operator != (const F_Matrix<T>& m) const
    { return !(*this == m); }

  bool operator == (F_TMatrix<T>& tm) const
    { F_Matrix<T> m (tm); return (*this == m); }
  bool operator != (F_TMatrix<T>& tm) const
    { return !(*this == tm); }

  bool operator == (F_TSMatrix<T>&) const;
  bool operator != (F_TSMatrix<T>& ts) const 
    { return !(*this == ts); } 

  // io-streams

  friend STD__ ostream& operator << FGD (STD__ ostream&, const F_Matrix<T>&);
  friend STD__ istream& operator >> FGD (STD__ istream&, F_Matrix<T>&);
  
  // misc
  // returns the the hermitian adjungated (Transposed and conjugate complex)
  F_TMatrix<T> herm () const { return (F_TMatrix<T> (*this)).herm(); }
  F_TMatrix<T> conj () const { return (F_TMatrix<T> (*this)).conj(); }
  F_TMatrix<T> trans() const { return (F_TMatrix<T> (*this)).trans(); }

  double fabs () const;
};


// definitions
  
template <typename T>
inline F_Matrix<T>::F_Matrix (const CSCMatrix<T>& m) 
: F_TMatrix<T>(m.rows(), m.columns())
{
        for (unsigned int r=0; r<m.rows(); r++)
                for (unsigned int c=0; c<m.columns(); c++)
                        operator()(r,c) = m(r,c);
}

template <typename T>
F_Matrix<T>::F_Matrix(const Matrix<T>& m)
: F_TMatrix<T>(m.rows(), m.columns())
{
  for (unsigned int r = 0; r < this->row; ++r)
    for (unsigned int c = 0; c < this->col; ++c)
                        setval(m(r, c), r, c);
}

template <typename T>
F_Matrix<T>::operator TMatrix<T>() const
{
  TMatrix<T> tm(this->row, this->col);
  for (unsigned int c = 0; c < this->col; ++c)
    for (unsigned int r = 0; r < this->row; ++r)
      tm.setval((*this)(r, c), r, c);
  return tm;
}


template <typename T>
STD__ ostream& operator << (STD__ ostream& os, const F_Matrix<T>& m)
{
  for (unsigned int r=0; r<m.row; r++) {
    for (unsigned int c=0; c<m.col; c++) 
      os << m(r,c) << " ";
    os << "\n";
  }
  return os;
}

INST(template <typename T> class F_TMatrix friend STD__ ostream& operator << (STD__ ostream& os, F_TMatrix<T> tm);)
template <typename T>
STD__ ostream& operator << (STD__ ostream& os, F_TMatrix<T> tm)
{
  F_Matrix<T> fm(tm);
        return os << fm;
}

#if 1
INST(template <typename T> class F_TSMatrix friend STD__ ostream& operator << (STD__ ostream& os, F_TSMatrix<T> ts);)
template <typename T>
STD__ ostream& operator << (STD__ ostream& os, F_TSMatrix<T> ts)
{
        return os << F_Matrix<T>(ts);
}
#endif

template <typename T>
STD__ istream& operator >> (STD__ istream& in, F_Matrix<T>& m)
{
  Vector<T> buf(m.row);
  
  for (unsigned int i=0; i<m.col; i++)
    {
      in >> buf;
      m.set_col (buf, i);
    }
  return in;
}


/************** Arithmetics ***************/

//change these to operate on diag only ?
#define F_MSTDDEF_T(op)                                         \
template <typename T>                                           \
inline F_TMatrix<T> F_Matrix<T>::operator op (const T& a) const \
{                                                               \
  F_TMatrix<T> t (this->row, this->col);                        \
  for (REGISTER T *p1=t.vec, *p2=this->vec; p2<this->endvec; ++p1, ++p2) \
    *p1 = *p2 op a;                                             \
  return t;                                                     \
}

F_MSTDDEF_T(+)
F_MSTDDEF_T(-)

template <typename T>
inline F_TSMatrix<T> F_Matrix<T>::operator * (const T& a) const
{ return F_TSMatrix<T> (*this, a); }

template <typename T>
inline F_TSMatrix<T> F_Matrix<T>::operator / (const T& a) const
{ 
   BCHK (a==(T)0, MatErr, Divide Mat by 0, 0, *this);
   return F_TSMatrix<T> (*this, (T)1/a); 
}

#define F_MSTDDEF_M(op)                                         \
template <typename T>                                           \
inline F_TMatrix<T> F_Matrix<T>::operator op (const F_Matrix<T>& a) const \
{                                                               \
  F_TMatrix<T> t (this->row, this->col);                        \
  BCHK(this->dim != a.dim, MatErr, Operator op on diff size matrices, a.dim, t); \
  for (REGISTER T *p1=t.vec, *p2=this->vec, *p3=a.vec; p2<this->endvec; ++p1, ++p2, ++p3) \
    *p1 = *p2 op *p3;                                           \
  return t;                                                     \
}

F_MSTDDEF_M(+)
F_MSTDDEF_M(-)

  
#define F_MSTDDEF_TM(op)                                        \
template <typename T>                                           \
inline F_TMatrix<T> F_Matrix<T>::operator op (F_TMatrix<T> a) const     \
{                                                               \
  BCHK(this->dim != a.dim, MatErr, Applying op on diff. size matrices, a.dim, a); \
  for (REGISTER T *p1=a.vec, *p2=this->vec; p2<this->endvec; ++p1, ++p2) \
    *p1 = *p2 op *p1;                                           \
  return a;                                                     \
}

F_MSTDDEF_TM(+)
F_MSTDDEF_TM(-)

#define F_MSTDDEF_TS(op)                                        \
template <typename T>                                           \
inline F_TMatrix<T> F_Matrix<T>::operator op (F_TSMatrix<T> ts) const \
{                                                               \
  BCHK(this->dim != ts.dim, MatErr, Applying op on diff. size matrices, ts.dim, (ts.destroy(),F_TMatrix<T> (*this)) ); \
  F_TSMatrix<T> tm (ts); tm.detach ();                          \
  for (REGISTER T *p1=tm.vec, *p2=this->vec, *p3=ts.vec; p2<this->endvec; ++p1, ++p2, ++p3) \
    *p1 = *p2 op *p3 * ts.fac;                                  \
  tm.fac = (T)1; return F_TMatrix<T> (tm);                      \
}

F_MSTDDEF_TS(+)
F_MSTDDEF_TS(-)

 
// F_Matrix multiplication
template <typename T>
F_TMatrix<T> F_Matrix<T>::operator * (const F_Matrix<T>& a) const
{
        F_TMatrix<T> ftm(this->row, a.col);
        BCHK(this->col != a.row, MatErr, Illegal sizes in F_Matrix multiplication, 0, ftm);
        for (unsigned int r=0; r<this->row; ++r) {
                for (REGISTER unsigned int c=0; c<a.col; ++c) {
                        ALIGN3(REGISTER T tmp, (*this)(r, 0) * a(0, c), MIN_ALIGN2);
                        for (REGISTER unsigned int l=1; l<this->col; ++l)
                                tmp += (*this)(r, l) * a(l, c);
                        ftm.setval(tmp, r, c);
                }
        }
        return ftm;
}

template <typename T>
F_TMatrix<T> F_Matrix<T>::operator * (F_TSMatrix<T> a) const
{
        F_TMatrix<T> ftm(this->row, a.col);
        BCHK(this->col != a.row, MatErr, Illegal sizes in F_Matrix multiplication, 0, ftm);
        REGISTER T tmp ALIGN(MIN_ALIGN2);
        for (unsigned int r=0; r<this->row; ++r) {
                for (REGISTER unsigned int c=0; c<a.col; ++c) {
                        tmp = (*this)(c, 0) * a.mat[c][0];
                        for (REGISTER unsigned int l=1; l<this->col; ++l) 
                                tmp += (*this)(r, l) * a.mat[c][l];
                        ftm.setval(tmp*a.fac, r, c);
                }
        }
        a.destroy (); return ftm;
}
 

template <typename T>
inline F_TMatrix<T> F_Matrix<T>::operator * (F_TMatrix<T> a) const
{
        F_TMatrix<T> ftm (this->operator * (F_Matrix<T> (a)));
        return ftm;
}

template <typename T>
inline F_TMatrix<T> F_TMatrix<T>::operator * (const F_Matrix<T>& a)
{
        F_Matrix<T> m (*this);
        return (m * a);
}


template <typename T>
inline F_TMatrix<T> F_TMatrix<T>::operator * (F_TMatrix<T> a)
{
        F_Matrix<T> m (*this);
        return (m * (F_Matrix<T> (a)));
}

template <typename T>
inline F_TMatrix<T> F_TMatrix<T>::operator * (F_TSMatrix<T> a)
{
        F_Matrix<T> m (*this);
        return (m * a);
}


template <typename T>
bool F_Matrix<T>::operator == (const F_Matrix<T>& m) const
{
        if (this->row != m.row || this->col != m.col)
                return false;
        if (this->vec == m.vec || this->dim == 0)
                return true;
        if (TBCICOMP (this->vec, m.vec, T, this->dim))
                return false;
        else 
                return true;
}

template <typename T>
bool F_Matrix<T>::operator == (F_TSMatrix<T>& ts) const
{
        if (this->row != ts.row || this->col != ts.col) 
                return (ts.destroy(), false);
        if (this->vec == 0) { 
                /*ts.destroy();*/ return true; 
        }
        if (ts.fac == (T)1) {
                if (this->vec == ts.vec) 
                        return true;
                if (TBCICOMP(this->vec, ts.vec, T, this->dim)) { 
                        ts.destroy (); return false; 
                }
        } else {
                if (this->vec == ts.vec) 
                        return false;
                for (REGISTER T *p1 = this->vec, *p2 = ts.vec; p1 < this->endvec; ++p1, ++p2)
                        if (*p1 != *p2 * ts.fac) { 
                                ts.destroy (); return false; 
                        }
        }
        ts.destroy (); 
        return true;
}


template<typename T>
TVector<T> F_Matrix<T>::operator * (const Vector<T>& v) const
{
        TVector<T> tv (this->row);
        for (unsigned int r=0; r<this->row; ++r) {
                REGISTER T el = 0;
                for (REGISTER unsigned int c=0; c<this->col; ++c) 
                        el += this->operator() (r, c) * v(c);
                tv.set (el, r);
        }
        return tv;
}

template<typename T>
inline TVector<T> F_Matrix<T>::operator * (TVector<T>& c) const
{
        Vector<T> cn (c);
        return (*this * cn);
}

#if 0
template <typename T>
inline TVector<T> F_Matrix<T>::operator () (const unsigned int i) const
{
        TVector<T> v(this->row);
        BCHK(i < 0 || i >= this->col, MatErr, Illegal row index in get_row, i, v);
        for (REGISTER unsigned int j=0; j<this->row; ++j) 
                v.vec[j] = this->mat[j][i];
        return v;
}
#endif

template <typename T>
inline void F_Matrix<T>::set_col (const Vector<T>& v, const unsigned int c)
{
        F_TMatrix<T>::set_col(v, c);
}

template <typename T>
inline TVector<T> F_Matrix<T>::get_col (const unsigned int c) const 
{
        return F_TMatrix<T>::get_col(c); 
}

template <typename T>
inline TVector<T> F_Matrix<T>::get_row (const unsigned int r) const
{
        return F_TMatrix<T>::get_row(r);
}

template <typename T>
inline void F_Matrix<T>::set_row (const Vector<T>& v, const unsigned int r)
{
        F_TMatrix<T>::set_row(v, r);
}

template <typename T>
inline T& F_Matrix<T>::operator () (const unsigned int r, const unsigned int c) const
{
        BCHK(r>=this->row, MatErr, illegal row index, r, this->mat[0][0]);
        BCHK(c>=this->col, MatErr, Illegal col index, c, this->mat[0][0]);
        return this->mat[c][r];
}

/*
template <typename T>
inline T F_Matrix<T>::operator () (const unsigned r, const unsigned c) const
{
        BCHK(r<0 || r>=row, MatErr, Illegal row index, r, mat[0][0]);
        BCHK(c<0 || c>=col, MatErr, illegal col index, c, mat[0][0]);
        return mat[c][r];
}
 */

template <typename T>
double F_Matrix<T>::fabs () const
{
        REGISTER double rv = 0.0;
        for (REGISTER T *ptr = this->vec; ptr < this->endvec; ++ptr)
                rv += fabssqr (*ptr);
        return MATH__ sqrt (rv);
}


template <typename T>
F_TMatrix<T> F_TMatrix<T>::transposed_copy() const
{
        F_TMatrix<T> trm(this->col, this->row);
        for (unsigned r = 0; r < this->row; ++r)
                for (unsigned c = 0; c < this->col; ++c)
                        trm.mat[c][r] = this->mat[r][c];
        // memleak if it's really a F_TMatrix :-(
        //this->destroy();
        //F_Matrix<T> fm(*this);
        return trm;
}


INST(template <typename T> class F_TMatrix friend F_TMatrix<T> transpose(const F_TMatrix<T>& tm);)
template <typename T>
inline F_TMatrix<T> transpose(const F_TMatrix<T>& ftm)
{
        return ftm.transposed_copy();
}

template <typename T>
F_TMatrix<T>& F_TMatrix<T>::transpose()
{
        if (this->col != this->row) {
                // for non quadratic shape, doing inplace swap ops is too hard
                F_Matrix<T> transmat(this->transposed_copy());
                return this->swap(transmat);
        } else {
                // do in place swap ops
                for (unsigned r = 1; r < this->row; ++r)
                        for (unsigned c = 0; c < r; ++c)
                                SWAP(this->mat[r][c], this->mat[c][r]);
                return *this;
        }
}

template <typename T>
inline F_TMatrix<T>& F_TMatrix<T>::swap (F_TMatrix<T>& m)
{
        //BCHK(this->row != m.row || this->col != m.col, MatErr, Swapping Matrices w/ different layout, m.row, *this);
        SWAP(this->dim, m.dim);
        SWAP(this->row, m.row); 
        //SWAP(this->col, m.col);
        unsigned int tcol = this->col; this->col = m.col; m.col = tcol;
        SWAP(this->mat, m.mat);
        SWAP(this->vec, m.vec); SWAP(this->endvec, m.endvec);
        return *this;
}

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

NAMESPACE_END
        
NAMESPACE_CSTD

INST(template <typename T> class TBCI__ F_TMatrix friend double fabs (const TBCI__ F_TMatrix<T>&);)
template <typename T>
inline double fabs (const TBCI__ F_TMatrix<T>& ftm)
{ return ftm.fabs (); }

INST(template <typename T> class TBCI__ F_Matrix friend double fabs (const TBCI__ F_Matrix<T>&);)
template <typename T>
inline double fabs (const TBCI__ F_Matrix<T>& fm)
{ return fm.fabs (); }

INST(template <typename T> class TBCI__ F_TSMatrix friend double fabs (TBCI__ F_TSMatrix<T>&);)
template <typename T>
inline double fabs (/*typename*/ TBCI__ F_TSMatrix<T>& ftsm)
{ return ftsm.fabs (); }

NAMESPACE_CSTD_END

#endif /* TBCI_F_MATRIX_H */