matrix.hpp

ncbi源码

    swap(m_Rows, M.m_Rows);
}


///
/// global operators
///


///
/// addition
///

template <class T, class U>
inline CNcbiMatrix< NCBI_PROMOTE(T,U) >
operator+ (const CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
{
    _ASSERT(m0.GetRows() == m1.GetRows());
    _ASSERT(m0.GetCols() == m1.GetCols());

    CNcbiMatrix<NCBI_PROMOTE(T,U)> res(m0.GetRows(), m0.GetCols());

    typename CNcbiMatrix<NCBI_PROMOTE(T,U)>::iterator res_iter = res.begin();
    typename CNcbiMatrix<NCBI_PROMOTE(T,U)>::iterator res_end  = res.end();
    typename CNcbiMatrix<T>::const_iterator           m0_iter  = m0.begin();
    typename CNcbiMatrix<U>::const_iterator           m1_iter  = m1.begin();

    for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter, ++m1_iter) {
        *res_iter = *m0_iter + *m1_iter;
    }

    return res;
}


template <class T, class U>
inline CNcbiMatrix<T>&
operator+= (CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
{
    _ASSERT(m0.GetRows() == m1.GetRows());
    _ASSERT(m0.GetCols() == m1.GetCols());

    typename CNcbiMatrix<T>::iterator           m0_iter  = m0.begin();
    typename CNcbiMatrix<T>::iterator           m0_end   = m0.end();
    typename CNcbiMatrix<U>::const_iterator     m1_iter  = m1.begin();

    size_t mod = (m0.GetRows() * m0.GetCols()) % 4;
    for ( ;  mod;  --mod, ++m0_iter, ++m1_iter) {
        *m0_iter += *m1_iter;
    }

    for ( ;  m0_iter != m0_end;  ) {
        *m0_iter++ += *m1_iter++;
        *m0_iter++ += *m1_iter++;
        *m0_iter++ += *m1_iter++;
        *m0_iter++ += *m1_iter++;
    }

    return m0;
}


///
/// subtraction
///

template <class T, class U>
inline CNcbiMatrix< NCBI_PROMOTE(T,U) >
operator- (const CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
{
    _ASSERT(m0.GetRows() == m1.GetRows());
    _ASSERT(m0.GetCols() == m1.GetCols());

    CNcbiMatrix< NCBI_PROMOTE(T,U) > res(m0.GetRows(), m0.GetCols());
    typename CNcbiMatrix<NCBI_PROMOTE(T,U)>::iterator res_iter = res.begin();
    typename CNcbiMatrix<NCBI_PROMOTE(T,U)>::iterator res_end  = res.end();
    typename CNcbiMatrix<T>::const_iterator           m0_iter  = m0.begin();
    typename CNcbiMatrix<U>::const_iterator           m1_iter  = m1.begin();

    for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter, ++m1_iter) {
        *res_iter = *m0_iter - *m1_iter;
    }

    return res;
}


template <class T, class U>
inline CNcbiMatrix<T>&
operator-= (CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
{
    _ASSERT(m0.GetRows() == m1.GetRows());
    _ASSERT(m0.GetCols() == m1.GetCols());

    typename CNcbiMatrix<T>::iterator           m0_iter  = m0.begin();
    typename CNcbiMatrix<T>::iterator           m0_end   = m0.end();
    typename CNcbiMatrix<U>::const_iterator     m1_iter  = m1.begin();

    for ( ;  m0_iter != m0_end;  ++m0_iter, ++m1_iter) {
        *m0_iter -= *m1_iter;
    }

    return m0;
}


///
/// multiplication
///
template <class T>
inline CNcbiMatrix< NCBI_PROMOTE(T, size_t) >
operator* (const CNcbiMatrix<T>& m0, size_t val)
{
    CNcbiMatrix< NCBI_PROMOTE(T, size_t) > res(m0.GetRows(), m0.GetCols());

    typename CNcbiMatrix<T>::iterator       res_iter = res.begin();
    typename CNcbiMatrix<T>::iterator       res_end  = res.end();
    typename CNcbiMatrix<T>::const_iterator m0_iter  = m0.begin();

    for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
        *res_iter = *m0_iter * val;
    }

    return res;
}


template <class U>
inline CNcbiMatrix< NCBI_PROMOTE(size_t, U) >
operator* (size_t val, const CNcbiMatrix<U>& m0)
{
    CNcbiMatrix< NCBI_PROMOTE(U, size_t) > res(m0.GetRows(), m0.GetCols());

    typename CNcbiMatrix<U>::iterator       res_iter = res.begin();
    typename CNcbiMatrix<U>::iterator       res_end  = res.end();
    typename CNcbiMatrix<U>::const_iterator m0_iter  = m0.begin();

    for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
        *res_iter = *m0_iter * val;
    }

    return res;
}


template <class T>
inline CNcbiMatrix< NCBI_PROMOTE(T, float) >
operator* (const CNcbiMatrix<T>& m0, float val)
{
    CNcbiMatrix< NCBI_PROMOTE(T, float) > res(m0.GetRows(), m0.GetCols());

    typename CNcbiMatrix<T>::iterator       res_iter = res.begin();
    typename CNcbiMatrix<T>::iterator       res_end  = res.end();
    typename CNcbiMatrix<T>::const_iterator m0_iter  = m0.begin();

    for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
        *res_iter = *m0_iter * val;
    }

    return res;
}


template <class U>
inline CNcbiMatrix< NCBI_PROMOTE(float, U) >
operator* (float val, const CNcbiMatrix<U>& m0)
{
    CNcbiMatrix< NCBI_PROMOTE(U, float) > res(m0.GetRows(), m0.GetCols());

    typename CNcbiMatrix<U>::iterator       res_iter = res.begin();
    typename CNcbiMatrix<U>::iterator       res_end  = res.end();
    typename CNcbiMatrix<U>::const_iterator m0_iter  = m0.begin();

    for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
        *res_iter = *m0_iter * val;
    }

    return res;
}


template <class T>
inline CNcbiMatrix< NCBI_PROMOTE(T, double) >
operator* (const CNcbiMatrix<T>& m0, double val)
{
    CNcbiMatrix< NCBI_PROMOTE(T, double) > res(m0.GetRows(), m0.GetCols());

    typename CNcbiMatrix<T>::iterator       res_iter = res.begin();
    typename CNcbiMatrix<T>::iterator       res_end  = res.end();
    typename CNcbiMatrix<T>::const_iterator m0_iter  = m0.begin();

    for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
        *res_iter = *m0_iter * val;
    }

    return res;
}


template <class U>
inline CNcbiMatrix< NCBI_PROMOTE(double, U) >
operator* (double val, const CNcbiMatrix<U>& m0)
{
    CNcbiMatrix< NCBI_PROMOTE(U, double) > res(m0.GetRows(), m0.GetCols());

    typename CNcbiMatrix<U>::iterator       res_iter = res.begin();
    typename CNcbiMatrix<U>::iterator       res_end  = res.end();
    typename CNcbiMatrix<U>::const_iterator m0_iter  = m0.begin();

    for ( ;  res_iter != res_end;  ++res_iter, ++m0_iter) {
        *res_iter = *m0_iter * val;
    }

    return res;
}


template <class T>
inline CNcbiMatrix<T>&
operator*= (CNcbiMatrix<T>& m0, T val)
{
    typename CNcbiMatrix<T>::iterator m0_iter  = m0.begin();
    typename CNcbiMatrix<T>::iterator m0_end  = m0.end();

    for ( ;  m0_iter != m0_end;  ++m0_iter) {
        *m0_iter *= val;
    }

    return m0;
}


template <class T>
inline CNcbiMatrix<T>&
operator/= (CNcbiMatrix<T>& m0, T val)
{
    val = T(1/val);

    typename CNcbiMatrix<T>::iterator m0_iter  = m0.begin();
    typename CNcbiMatrix<T>::iterator m0_end  = m0.end();

    for ( ;  m0_iter != m0_end;  ++m0_iter) {
        *m0_iter *= val;
    }

    return m0;
}


template <class T, class U>
inline vector< NCBI_PROMOTE(T,U) >
operator* (const CNcbiMatrix<T>& m, const vector<U>& v)
{
    _ASSERT(m.GetCols() == v.size());

    typedef NCBI_PROMOTE(T,U) TPromote;
    vector< TPromote > res(m.GetRows(), TPromote(0));

    for (size_t r = 0;  r < m.GetRows();  ++r) {
        for (size_t i = 0;  i < m.GetCols();  ++i) {
            res[r] += m(r,i) * v[i];
        }
    }

    return res;
}


template <class T, class U>
inline vector< NCBI_PROMOTE(T,U) >
operator* (const vector<T>& v, const CNcbiMatrix<U>& m)
{
    _ASSERT(m.GetRows() == v.size());

    typedef NCBI_PROMOTE(T,U) TPromote;
    vector<TPromote> res(m.GetCols(), TPromote(0));

    for (size_t c = 0;  c < m.GetCols();  ++c) {
        for (size_t i = 0;  i < m.GetRows();  ++i) {
            res[c] += m(i,c) * v[i];
        }
    }

    return res;
}


template <class T, class U>
inline CNcbiMatrix< NCBI_PROMOTE(T,U) >
operator* (const CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
{
    _ASSERT(m0.GetCols() == m1.GetRows());

    typedef NCBI_PROMOTE(T,U) TPromote;
    CNcbiMatrix< TPromote > res(m0.GetRows(), m1.GetCols(), TPromote(0));

#if 1
    for (size_t r = 0; r < m0.GetRows();  ++r) {
        for (size_t c = 0;  c < m1.GetCols();  ++c) {
            for (size_t i = 0;  i < m0.GetCols();  ++i) {
                res(r,c) += m0(r,i) * m1(i,c);
            }
        }
    }
#endif


#if 0
    const vector<T>& lhs = m0.GetData();
    const vector<U>& rhs = m1.GetData();
    size_t col_mod = m0.GetCols() % 4;

    for (size_t r = 0;  r < m0.GetRows();  ++r) {
        size_t r_offs = r * m0.GetCols();

        for (size_t c = 0;  c < m1.GetCols();  ++c) {
            size_t i=0;
            T t0 = 0;
            T t1 = 0;
            T t2 = 0;
            T t3 = 0;

            switch(col_mod) {
            default:
            case 0:
                break;

            case 3:
                t0 += m0.GetData()[r_offs + 2] * m1(2,c);
                ++i;
            case 2:
                t0 += m0.GetData()[r_offs + 1] * m1(1,c);
                ++i;
            case 1:
                t0 += m0.GetData()[r_offs + 0] * m1(0,c);
                ++i;
                break;
            }

            for (;  i < m0.GetCols();  i += 2) {
                t0 += lhs[r_offs + i + 0] * m1(i + 0, c);
                t1 += lhs[r_offs + i + 1] * m1(i + 1, c);
                //t2 += lhs[r_offs + i + 2] * m1(i + 2, c);
                //t3 += lhs[r_offs + i + 3] * m1(i + 3, c);
            }
            res(r,c) = t0 + t1 + t2 + t3;
        }
    }
#endif

    return res;
}


template <class T, class U>
inline CNcbiMatrix<T>&
operator*= (CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
{
    _ASSERT(m0.GetCols() == m1.GetRows());

    m0 = m0 * m1;
    return m0;
}


///
/// comparators
///

template <class T, class U>
inline bool
operator== (const CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
{
    if (m0.GetRows() != m1.GetRows()) {
        return false;
    }
    if (m0.GetCols() != m1.GetCols()) {
        return false;
    }

    for (size_t i = 0;  i < m0.GetData().size();  ++i) {
        if (m0[i] != m1[i]) {
            return false;
        }
    }

    return true;
}


template <class T, class U>
inline bool
operator!= (const CNcbiMatrix<T>& m0, const CNcbiMatrix<U>& m1)
{
    return !(m0 == m1);
}


///
/// stream output
template <class T>
inline CNcbiOstream& operator<<(CNcbiOstream& os, const CNcbiMatrix<T>& M)
{
    for (size_t i = 0;  i < M.GetRows();  ++i) {
        for (size_t j = 0;  j < M.GetCols();  ++j) {
            if (j > 0) {
                os << " ";
            }
            os << M(i, j);
        }
        os << NcbiEndl;
    }
    return os;
}


///
/// stream input
template <class T>
inline CNcbiIstream& operator>>(CNcbiIstream& is, CNcbiMatrix<T>& M)
{
    CNcbiMatrix<T> A;
    string line;
    vector<T> row;
    T entry;
    int linenum = 0;
    while(getline(is, line)) {
        linenum++;
        if (line.empty() || line[0] == '#') {
            continue;
        }
        CNcbiIstrstream iss(line.c_str(), line.size());
        row.clear();
        while(1) {
            iss >> entry;
            if (!iss) {
                break;
            }
            row.push_back(entry);
        }
        if (A.GetCols() == 0) {
            A.Resize(A.GetCols(), row.size());
        }
        if (row.size() != A.GetCols()) {
            NCBI_THROW(CException, eUnknown,
                       "error at line " +
                       NStr::IntToString(linenum) + ": expected " +
                       NStr::IntToString(A.GetCols()) + " columns; got" +
                       NStr::IntToString(row.size()));
        }
        A.Resize(A.GetRows() + 1, A.GetCols());
        for (int i = 0;  i < A.GetCols();  ++i) {
            A(A.GetRows() - 1, i) = row[i];
        }
    }
    M.Swap(A);
    return is;
}



END_NCBI_SCOPE

/*
 * ===========================================================================
 * $Log: matrix.hpp,v $
 * Revision 1000.1  2004/04/16 15:28:41  gouriano
 * PRODUCTION: UPGRADED [CATCHUP_003] Dev-tree R1.6
 *
 * Revision 1.6  2004/03/11 17:31:42  jcherry
 * Sped up adding of rows (greatly speeds reading from streams)
 *
 * Revision 1.5  2004/03/10 14:13:29  dicuccio
 * Corrected include - use ncbistr instead of ncbistre
 *
 * Revision 1.4  2004/02/10 21:17:54  jcherry
 * Made operator>> work for a matrix of any contained type that defines
 * operator>> (not just doubles)
 *
 * Revision 1.3  2004/02/10 14:29:51  jcherry
 * Removed wayward use of protected data in operator>>
 *
 * Revision 1.2  2004/02/10 14:00:40  dicuccio
 * Added comments.  Fixed spacing issues.  Changed extraction operator to be
 * non-friend function
 *
 * Revision 1.1  2004/02/09 17:34:11  dicuccio
 * Moved from gui/math
 *
 * Revision 1.4  2004/02/02 15:55:49  jcherry
 * Added Swap() method and stream reader
 *
 * Revision 1.3  2004/01/29 22:38:03  jcherry
 * Fixed typo
 *
 * Revision 1.2  2004/01/29 21:30:17  jcherry
 * Added stream output operator
 *
 * Revision 1.1  2004/01/23 13:54:15  dicuccio
 * Initial revision
 *
 * ===========================================================================
 */

#endif  // UTIL_MATH___MATRIX__HPP