fmat.h

Go to the documentation of this file.

#ifndef FMAT_H
#define FMAT_H

#include "vec.h"
#include <stdlib.h>
#include <assert.h>

#ifdef TNT_USE_REGIONS
#include "region2d.h"
#endif

template <class T>
class Fortran_matrix 
{
  public:

    typedef         T   value_type;
    typedef         T   element_type;
    typedef         T*  pointer;
    typedef         T*  iterator;
    typedef         T&  reference;
    typedef const   T*  const_iterator;
    typedef const   T&  const_reference;

    Subscript lbound() const { return 1;}
 
        enum AllocMode {INTERNAL_ALLOC, EXTERNAL_ALLOC} amode; 
        // determine if Storage for v_ is allocated when matrix is 
        // created and is destroyed when matrix is destroyed or   
    // v points to externally allocated memory region assumed 
        // to be of the sufficient size (m_*n_)

  protected:
    T* v_;               
    Subscript m_;
    Subscript n_;
    T** col_;           

        int MaxSize; 

 
        // internal helper function to create the array
    // of row pointers

    void initialize(Subscript M, Subscript N)
    {
        // adjust col_[] pointers so that they are 1-offset:
        //   col_[j][i] is really col_[j-1][i-1];
        //
        // v_[] is the internal contiguous array, it is still 0-offset
        //
       if(M == 0 || N == 0)
         {
            m_ =  0;
            n_ =  0;
            v_ =  NULL;
            col_= NULL;
            return ;
         }

                if(amode == INTERNAL_ALLOC) v_ = new T[M*N];
                if(amode == EXTERNAL_ALLOC) assert( (M*N) <= this->MaxSize);

        col_ = new T*[N];

        assert(v_  != NULL);
        assert(col_ != NULL);


        m_ = M;
        n_ = N;  
        T* p = v_ - 1;              
        for (Subscript i=0; i<N; i++)
        {
            col_[i] = p;
            p += M ;
            
        }
        col_ --; 
    }
   
    void copy(const T*  v)
    {
        Subscript N = m_ * n_;
        Subscript i;

#ifdef TNT_UNROLL_LOOPS
        Subscript Nmod4 = N & 4;
        Subscript N4 = N - Nmod4;

        for (i=0; i<N4; i+=4)
        {
            v_[i] = v[i];
            v_[i+1] = v[i+1];
            v_[i+2] = v[i+2];
            v_[i+3] = v[i+3];
        }

        for (i=N4; i< N; i++)
            v_[i] = v[i];
#else

        for (i=0; i< N; i++)
            v_[i] = v[i];
#endif      
    }

    void set(const T& val)
    {
        Subscript N = m_ * n_;
        Subscript i;

#ifdef TNT_UNROLL_LOOPS
        Subscript Nmod4 = N & 4;
        Subscript N4 = N - Nmod4;

        for (i=0; i<N4; i+=4)
        {
            v_[i] = val;
            v_[i+1] = val;
            v_[i+2] = val;
            v_[i+3] = val; 
        }

        for (i=N4; i< N; i++)
            v_[i] = val;
#else

        for (i=0; i< N; i++)
            v_[i] = val;
        
#endif      
    }
    


    void destroy()
    {     
        /* do nothing, if no memory has been previously allocated */
        if (v_ == NULL) return ;

        /* if we are here, then matrix was previously allocated */
                if(amode == INTERNAL_ALLOC)
                {
                        delete [] (v_); // do not delete if externally allocated 
                }
        col_ ++;                // changed back to 0-offset
        delete [] (col_);
    }


  public:

    T* begin() { return v_; }
    const T* begin() const { return v_;}

    T* end() { return v_ + m_*n_; }
    const T* end() const { return v_ + m_*n_; }


    // constructors

    Fortran_matrix() : v_(0), m_(0), n_(0), col_(0),
                               amode(INTERNAL_ALLOC), MaxSize(0) {}

    Fortran_matrix(const Fortran_matrix<T> & A)
    {
                amode=INTERNAL_ALLOC; 
                MaxSize = 0;
        initialize(A.m_, A.n_);
        copy(A.v_);
    }

    Fortran_matrix(Subscript M, Subscript N, const T& value = T(0))
    {
                amode=INTERNAL_ALLOC; MaxSize=0;
        initialize(M,N);
        set(value);
    }

   Fortran_matrix(Subscript M, Subscript N,  T* v, 
                const int New_amode=0)
    {
                Set_amode(New_amode);
                if(amode == EXTERNAL_ALLOC)
                {
                        MaxSize=M*N;
                        v_= v;
                }
        initialize(M,N);
        if(amode == INTERNAL_ALLOC) copy(v);
    }


//    Fortran_matrix(Subscript M, Subscript N, const T* v, 
//              const int New_amode=0)
//    {
//              Set_amode(New_amode);
//              if(amode == EXTERNAL_ALLOC)
//              {
//                      MaxSize=M*N;
//                      v_= v;
//              }
//       initialize(M,N);
//        if(amode == INTERNAL_ALLOC) copy(v);
//    }


    Fortran_matrix(Subscript M, Subscript N, char *s)
    {
                amode=INTERNAL_ALLOC;
        initialize(M,N);
        istrstream ins(s);

        Subscript i, j;

        for (i=1; i<=M; i++)
            for (j=1; j<=N; j++)
                ins >> (*this)(i,j);
    }

    // destructor
    ~Fortran_matrix()
    {
        destroy();
    }


    // assignments
    //
    Fortran_matrix<T>& operator=(const Fortran_matrix<T> &A)
    {
        if (v_ == A.v_)
            return *this;

        if (m_ == A.m_  && n_ == A.n_)      // no need to re-alloc
            copy(A.v_);

        else
        {
            destroy();
            initialize(A.m_, A.n_);
            copy(A.v_);
        }

        return *this;
    }
        
    Fortran_matrix<T>& operator=(const T& scalar)
    { 
        set(scalar); 
        return *this;
    }


    Subscript dim(Subscript d) const 
    {
#ifdef TNT_BOUNDS_CHECK
       assert( d >= 1);
        assert( d <= 2);
#endif
        return (d==1) ? m_ : ((d==2) ? n_ : 0); 
    }

    Subscript num_rows() const { return m_; }
    Subscript num_cols() const { return n_; }

    void newsize(Subscript M, Subscript N)
    {
        if (num_rows() == M && num_cols() == N)
            return;

        destroy();
        initialize(M,N);

        return;
    }

   void Set_amode(int mode) 
        { 
                assert(mode == 0 || mode == 1);
                if(mode == 0) amode=INTERNAL_ALLOC; 
                if(mode == 1) amode=EXTERNAL_ALLOC;  
        }

        void SetMaxSize(const int NewMaxSize)
        {
                MaxSize=NewMaxSize;
        }


    inline reference operator()(Subscript i, Subscript j)
    { 
#ifdef TNT_BOUNDS_CHECK
        assert(1<=i);
        assert(i <= m_) ;
        assert(1<=j);
        assert(j <= n_);
#endif
        return col_[j][i]; 
    }

    inline const_reference operator() (Subscript i, Subscript j) const
    {
#ifdef TNT_BOUNDS_CHECK
        assert(1<=i);
        assert(i <= m_) ;
        assert(1<=j);
        assert(j <= n_);
#endif
        return col_[j][i]; 
    }

    inline reference operator[](Subscript i)
    { 
#ifdef TNT_BOUNDS_CHECK
        assert(0 <= i);
        assert(i <  m_*n_) ;
#endif
        return v_[i]; 
    }

    inline const_reference operator[] (Subscript i) const
    {
#ifdef TNT_BOUNDS_CHECK
        assert(0 <= i);
        assert(i <  m_*n_) ;
#endif
        return v_[i]; 
    }

        inline value_type GetVal(Subscript i, Subscript j) const    
    {
#ifdef TNT_BOUNDS_CHECK
        assert(1<=i);
        assert(i <= m_) ;
        assert(1<=j);
        assert(j <= n_);
#endif
        return col_[j][i]; 
        }

        inline void SetVal(Subscript i, Subscript j, const value_type new_val)
        {
#ifdef TNT_BOUNDS_CHECK
        assert(1<=i);
        assert(i <= m_) ;
        assert(1<=j);
        assert(j <= n_);
#endif
        col_[j][i] = new_val; 
        }

        inline value_type GetVal_lidx(Subscript i) const    
    {
#ifdef TNT_BOUNDS_CHECK
        assert(1<=i);
        assert(i <= m_*n_) ;
#endif
        return v_[i]; 
        }

        inline void SetVal_lidx(Subscript i, const value_type new_val)
        {
#ifdef TNT_BOUNDS_CHECK
        assert(1<=i);
        assert(i <= m_*n_) ;
#endif
        v_[i] = new_val; 
        }

 
    //#ifdef __GNUC__
    //template <class T>   
    //#endif  
#if defined(_MSC_VER) || defined(__xlC__)
friend istream& operator>>(istream &s, Fortran_matrix<T> &A);  
#else

friend istream& operator>><>(istream &s, Fortran_matrix<T> &A); 
#endif
        
#ifdef TNT_USE_REGIONS

    typedef Region2D<Fortran_matrix<T> > Region;
    typedef const_Region2D<Fortran_matrix<T>,T > const_Region;

    Region operator()(const Index1D &I, const Index1D &J)
    {
        return Region(*this, I,J);
    }

    const_Region operator()(const Index1D &I, const Index1D &J) const
    {
        return const_Region(*this, I,J);
    }

#endif


        int Print_format(ostream &sout, const char* format) const
        {
                
                char buf[120]; // set max 20 characters per number
                
                sout << m_ << " " << n_ << endl;
                
                for (Subscript i=1; i<= m_; i++)
                {
                        for (Subscript j=1; j<= n_; j++)
                        {
                                sprintf(buf, format, (*this)(i,j) );
                                sout << buf;
                        }
                        sout << endl;
                }
                
                return 0;
        }

    int GetSymmPart(pointer dsymm)
                // get Symmetrical part of the square matrix in low-diagonal form:
        {
                assert(dsymm != NULL);
                assert(m_ == n_);
                
                for (Subscript i=1; i<= m_; i++)
                {
                        for (Subscript j=1; j<= i; j++)
                        {
                                (*dsymm)= ((*this)(i,j) + (*this)(j,i))/2;
                                dsymm++;
                        }
                }
                return 1;               
        }

    int GetASymmPart(pointer dasymm)
    // get AntiSymmetrical part of the square matrix in lower triangular form
        // (diagonal inlcuded)
        {
                assert(dasymm != NULL);
                assert(m_ == n_);
                
                for (Subscript i=1; i<= m_; i++)
                {
                        for (Subscript j=1; j<= i; j++)
                        {
                                (*dasymm)= ((*this)(i,j) - (*this)(j,i))/2;
                                dasymm++;
                        }
                }
                return 1;               
        }

    int AddSymmPart(pointer dsymm)
        // add Symmetrical matrix in low triangular form 
        // to the current square matrix 
        {
                assert(dsymm != NULL);
                assert(m_ == n_);
                
                for (Subscript i=1; i<= m_; i++)
                {
                        for (Subscript j=1; j<= i; j++)
                        {
                                (*this)(i,j)+= (*dsymm);
                                if(i != j) (*this)(j,i)+= (*dsymm);
                                dsymm++;
                        }
                }
                return 1;               
        }

    int AddASymmPart(pointer dasymm)
        // add AntiSymmetrical matrix in low triangular form (diagonal included)
        // to the current square matrix 
        {
                assert(dasymm != NULL);
                assert(m_ == n_);
                
                for (Subscript i=1; i<= m_; i++)
                {
                        for (Subscript j=1; j<= i; j++)
                        {
                                (*this)(i,j)-= (*dasymm);
                                (*this)(j,i)+= (*dasymm);
                                dasymm++;
                        }
                }
                return 1;               
        }

        int Symmetrize()
        // A= 1/2(A+A^T)
        {
                assert(m_ == n_);
                if(m_ == 0)
                        return 1;
                for(Subscript i=1; i <= m_; i++)
                {
                        for (Subscript j=1; j< i; j++)
                        {
                                (*this)(i,j)= 0.5*((*this)(i,j)+ (*this)(j,i));
                        }

                }
                return 1;
        }

        int ASymmetrize()
        // A= 1/2(A-A^T)
        {
                assert(m_ == n_);
                if(m_ == 0)
                        return 1;
                for(Subscript i=1; i <= m_; i++)
                {
                        for (Subscript j=1; j< i; j++)
                        {
                                (*this)(i,j)= 0.5*((*this)(i,j) - (*this)(j,i));
                                (*this)(j,i)=-(*this)(i,j);
                        }

                }
                return 1;
        }

};


/* ***************************  I/O  ********************************/

template <class T>
ostream& operator<<(ostream &s, const Fortran_matrix<T> &A)
{
    Subscript M=A.num_rows();
    Subscript N=A.num_cols();

    s << M << " " << N << endl;

    for (Subscript i=1; i<=M; i++)
    {
        for (Subscript j=1; j<=N; j++)
        {
            s << A(i,j) << " ";
        }
        s << endl;
    }


    return s;
}



template <class T>
istream& operator>>(istream &s, Fortran_matrix<T> &A)
{

    Subscript M, N;

    s >> M >> N;

    if ( !(M == A.m_ && N == A.n_) )
    {
        A.destroy();
        A.initialize(M,N);
    }


    for (Subscript i=1; i<=M; i++)
        for (Subscript j=1; j<=N; j++)
        {
            s >>  A(i,j);
        }

    return s;
}

//*******************[ basic matrix algorithms ]***************************


template <class T>
Fortran_matrix<T> operator+(const Fortran_matrix<T> &A, 
    const Fortran_matrix<T> &B)
{
    Subscript M = A.num_rows();
    Subscript N = A.num_cols();

    assert(M==B.num_rows());
    assert(N==B.num_cols());

    Fortran_matrix<T> tmp(M,N);
    Subscript i,j;

    for (i=1; i<=M; i++)
        for (j=1; j<=N; j++)
            tmp(i,j) = A(i,j) + B(i,j);

    return tmp;
}

template <class T>
Fortran_matrix<T> operator-(const Fortran_matrix<T> &A, 
    const Fortran_matrix<T> &B)
{
    Subscript M = A.num_rows();
    Subscript N = A.num_cols();

    assert(M==B.num_rows());
    assert(N==B.num_cols());

    Fortran_matrix<T> tmp(M,N);
    Subscript i,j;

    for (i=1; i<=M; i++)
        for (j=1; j<=N; j++)
            tmp(i,j) = A(i,j) - B(i,j);

    return tmp;
}

// element-wise multiplication  (use matmult() below for matrix
// multiplication in the linear algebra sense.)
//
//
template <class T>
Fortran_matrix<T> mult_element(const Fortran_matrix<T> &A, 
    const Fortran_matrix<T> &B)
{
    Subscript M = A.num_rows();
    Subscript N = A.num_cols();

    assert(M==B.num_rows());
    assert(N==B.num_cols());

    Fortran_matrix<T> tmp(M,N);
    Subscript i,j;

    for (i=1; i<=M; i++)
        for (j=1; j<=N; j++)
            tmp(i,j) = A(i,j) * B(i,j);

    return tmp;
}


template <class T>
Fortran_matrix<T> transpose(const Fortran_matrix<T> &A)
{
    Subscript M = A.num_rows();
    Subscript N = A.num_cols();

    Fortran_matrix<T> S(N,M);
    Subscript i, j;

    for (i=1; i<=M; i++)
        for (j=1; j<=N; j++)
            S(j,i) = A(i,j);

    return S;
}


    
template <class T>
inline Fortran_matrix<T> matmult(const Fortran_matrix<T>  &A, 
    const Fortran_matrix<T> &B)
{

#ifdef TNT_BOUNDS_CHECK
    assert(A.num_cols() == B.num_rows());
#endif

    Subscript M = A.num_rows();
    Subscript N = A.num_cols();
    Subscript K = B.num_cols();

    Fortran_matrix<T> tmp(M,K);
    T sum;

    for (Subscript i=1; i<=M; i++)
    for (Subscript k=1; k<=K; k++)
    {
        sum = 0;
        for (Subscript j=1; j<=N; j++)
            sum = sum +  A(i,j) * B(j,k);

        tmp(i,k) = sum; 
    }

    return tmp;
}

template <class T>
inline Fortran_matrix<T> operator*(const Fortran_matrix<T> &A, 
    const Fortran_matrix<T> &B)
{
    return matmult(A,B);
}

template <class T>
inline int mat_diff(Fortran_matrix<T>& C, Fortran_matrix<T>  &A, 
    Fortran_matrix<T> &B)
// C= A-B
// C may coincide with A and B
{

    assert(A.num_cols() == B.num_cols());
        assert(A.num_rows() == B.num_rows());

    Subscript M = A.num_rows();
    Subscript N = A.num_cols();

    C.newsize(M,N);         // adjust shape of C, if necessary

    for(Subscript i=1; i <= M; i++)
        {
                for(Subscript j=1; j <= N; j++)
                {
                        C(i,j)=A(i,j)-B(i,j);
                }
        }
        return 0;
}

template <class T>
inline int mat_add(Fortran_matrix<T>& C, Fortran_matrix<T>  &A, 
    Fortran_matrix<T> &B)
// C= A+B
// C may coincide with A and B
{

    assert(A.num_cols() == B.num_cols());
        assert(A.num_rows() == B.num_rows());

    Subscript M = A.num_rows();
    Subscript N = A.num_cols();

    C.newsize(M,N);         // adjust shape of C, if necessary

    for(Subscript i=1; i <= M; i++)
        {
                for(Subscript j=1; j <= N; j++)
                {
                        C(i,j)=A(i,j)+B(i,j);
                }
        }
        return 0;
}



template <class T>
inline int mat_add_unit(Fortran_matrix<T>& C, Fortran_matrix<T>  &A, 
    const T& factor)
// C= A + factor*I , where I is a unit matrix
// C may coincide with A 
{
        Subscript N=A.num_rows();
        assert(A.num_rows() == A.num_cols());
    C.newsize(N,N);    // adjust shape of C, if necessary
        C=A;
        for(Subscript i=1; i<=N; i++)
        {
                C(i,i)=C(i,i)+factor;
        }
        return 0;
}

template <class T>
T SProd(const Fortran_matrix<T>  &A, 
    const Fortran_matrix<T> &B)
// Scalar Product (sum of element pair product of two matricies
{

    assert(A.num_cols() == B.num_cols());
        assert(A.num_rows() == B.num_rows());

    Subscript M = A.num_rows();
    Subscript N = A.num_cols();
        
        T ss;
    
    ss=0;
    for(Subscript i=1; i <= M; i++)
        {
                for(Subscript j=1; j <= N; j++)
                {
                        ss+=A(i,j)*B(i,j);
                }
        }
        return ss;
}


template <class T>
inline int mat_scale(Fortran_matrix<T>& C, Fortran_matrix<T>  &A, 
    const T& factor)
// C= factor*A  , where I is a unit matrix
// C may coincide with A 
{
    Subscript M = A.num_rows();
    Subscript N = A.num_cols();

    Subscript MN = M*N; 

    C.newsize(M,N);    // adjust shape of C, if necessary

    const T* a = A.begin();
    T* t = C.begin();
    T* tend = C.end();

    for (; t < tend; t++, a++)
        *t = *a * factor;

        return 0;
}

template <class T>
inline int mat_copy_scale(Fortran_matrix<T>& C, const Fortran_matrix<T>  &A, 
    const T& factor)
// C= factor*A  , where I is a unit matrix
// C may coincide with A 
{
    Subscript M = A.num_rows();
    Subscript N = A.num_cols();

    Subscript MN = M*N; 

    C.newsize(M,N);    // adjust shape of C, if necessary

    const T* a = A.begin();
    T* t = C.begin();
    T* tend = C.end();

    for (; t < tend; t++, a++)
        *t = *a * factor;

        return 0;
}

template <class T>
inline int matmult(Fortran_matrix<T>& C, const Fortran_matrix<T>  &A, 
    const Fortran_matrix<T> &B)
{
    assert(A.num_cols() == B.num_rows());

    Subscript M = A.num_rows();
    Subscript N = A.num_cols();
    Subscript K = B.num_cols();

    C.newsize(M,K);         // adjust shape of C, if necessary

    T sum; 

    const T* row_i;
    const T* col_k;

    for (Subscript i=1; i<=M; i++)
    {
        for (Subscript k=1; k<=K; k++)
        {
            row_i = &A(i,1);
            col_k = &B(1,k);
            sum = 0;
            for (Subscript j=1; j<=N; j++)
            {
                sum +=  *row_i * *col_k;
                row_i += M;
                col_k ++;
            }
            C(i,k) = sum; 
        }
    }
    return 0;
}


template <class T>
NumVector<T> matmult(const Fortran_matrix<T>  &A, const NumVector<T> &x)
{

#ifdef TNT_BOUNDS_CHECK
    assert(A.num_cols() == x.dim());
#endif

    Subscript M = A.num_rows();
    Subscript N = A.num_cols();

    NumVector<T> tmp(M);
    T sum;

    for (Subscript i=1; i<=M; i++)
    {
        sum = 0;
        for (Subscript j=1; j<=N; j++)
            sum = sum +  A(i,j) * x(j);

        tmp(i) = sum; 
    }

    return tmp;
}

template <class T>
inline NumVector<T> operator*(const Fortran_matrix<T>  &A, const NumVector<T> &x)
{
    return matmult(A,x);
}

template <class T>
inline Fortran_matrix<T> operator*(const Fortran_matrix<T>  &A, const T &x)
{
    Subscript M = A.num_rows();
    Subscript N = A.num_cols();

    Subscript MN = M*N; 

    Fortran_matrix<T> res(M,N);
    const T* a = A.begin();
    T* t = res.begin();
    T* tend = res.end();

    for (t=res.begin(); t < tend; t++, a++)
        *t = *a * x;

    return res;
} 

// my functions

template <class T>
void raw_copy(T* to, const T* from, const int nelem)
// copy elements from the source to destination
{
        assert(nelem > 0);
        assert(to != NULL); 
        assert(from != NULL);

        const T* ref=from;
        T* dest=to;
        int i=nelem;
        for(; i != 0; dest++, ref++, i--)
        {
                *dest=*ref;
        }
}

template <class T>
inline int matmult_T1(Fortran_matrix<T>& C, const Fortran_matrix<T>  &A, 
    const Fortran_matrix<T> &B)
// C = A^T * B 
{

    assert(A.num_rows() == B.num_rows());

    Subscript M = A.num_cols();
    Subscript N = A.num_rows();
    Subscript K = B.num_cols();

    C.newsize(M,K);         // adjust shape of C, if necessary


    T sum; 

    const T* col_i;
    const T* col_k;

    for (Subscript i=1; i<=M; i++)
    {
        for (Subscript k=1; k<=K; k++)
        {
            col_i = &A(1,i);
            col_k = &B(1,k);
            sum = 0;
            for (Subscript j=1; j<=N; j++)
            {
                sum +=  *col_i * *col_k;
                col_i ++;
                col_k ++;
            }
        
            C(i,k) = sum; 
        }

    }

    return 0;
}

template <class T>
inline int matmult_T2(Fortran_matrix<T>& C, const Fortran_matrix<T>  &A, 
    const Fortran_matrix<T> &B)
// C = A * B^T
{

    assert(A.num_cols() == B.num_cols());

    Subscript M = A.num_rows();
    Subscript N = A.num_cols();
    Subscript K = B.num_rows();

    C.newsize(M,K);         // adjust shape of C, if necessary


    T sum; 

    const T* row_i;
    const T* row_k;

    for (Subscript i=1; i<=M; i++)
    {
        for (Subscript k=1; k<=K; k++)
        {
            row_i = &A(i,1);
            row_k = &B(k,1);
            sum = 0;
            for (Subscript j=1; j<=N; j++)
            {
                sum +=  *row_i * *row_k;
                row_i += M;
                row_k += K;
            }
        
            C(i,k) = sum; 
        }

    }

    return 0;
}

template<class NumMatrix, class NumVector>
int mat_symm_from_lin(NumMatrix & smat, const NumVector & slin, Subscript n)
// convert Symmetrical Matrix from linear storage to Square form
// ( indexation of the matrix and the vector are assumed to be 1-based )
{
        smat.newsize(n,n);
        for(Subscript i=1; i <= n ; i++)
        {
                for(Subscript j=1; j <= n ; j++)
                {
              smat(i,j)=((i >= j) ? slin(i*(i-1)/2 + j) : slin(j*(j-1)/2 + i));
                }
        }
        return 1;
}

template<class NumMatrix, class NumVector>
int mat_asymm_from_lin(NumMatrix & smat, NumVector & slin, Subscript n)
// convert AntiSymmetrical Matrix from linear storage to Square form
// ( indexation of the matrix and the vector are assumed to be 1-based )
{
        smat.newsize(n,n);
        for(Subscript i=1; i <= n ; i++)
        {
                for(Subscript j=1; j <= n ; j++)
                {
              smat(i,j)=((i >= j) ? slin(i*(i-1)/2 + j) : -slin(j*(j-1)/2 + i));
                }
        }
        return 1;
}

#endif
// FMAT_H

---