utils.h

很好用的库

    while (not_at(A_i, (*A_kk).rend())) {
      std::cout << *A_i;
      ++A_i;
      std::cout << ",";
    }
    std::cout << "]";
    std::cout << "," << std::endl;
    ++A_kk;
  } 
  std::cout << "]" << std::endl;
}

///
template <class Matrix>
void print_major(const Matrix& A)
{
  typename Matrix::const_iterator A_kk;
  typename Matrix::OneD::const_iterator A_i;
  A_kk = A.begin();
  std::cout << "[" << std::endl;
  while (not_at(A_kk, A.end())) {
    std::cout << "[";
    A_i = (*A_kk).begin();
    while (not_at(A_i, (*A_kk).end())) {
      std::cout << *A_i;
      ++A_i;
      std::cout << ",";
    }
    std::cout << "]";
    std::cout << "," << std::endl;
    ++A_kk;
  } 
  std::cout << "]" << std::endl;
}

#if 1
template <class MatrixA, class MatrixB>
inline bool
matrix_equal(const MatrixA& A, const MatrixB& B)
{
  typedef typename matrix_traits<MatrixA>::size_type Int;
  if (A.nrows() != Int(B.nrows()) || A.ncols() != Int(B.ncols()))
    return false;

  for (Int i = 0; i < A.nrows(); ++i)
    for (Int j = 0; j < A.ncols(); ++j)
      if (A(i,j) != B(i,j)) {
#if !defined(_MSVCPP_)
		  std::cout << "(" << i << "," << j << ") " 
	     << A(i,j) << " != " << B(i,j) << std::endl;
#endif
		  return false;
      }

  return true;
}

#else

template <class MatrixA, class MatrixB>
inline bool
__matrix_equal(const MatrixA& A, const MatrixB& B, row_tag)
{
  typename MatrixA::const_row_2Diterator A_ii = A.begin_rows();
  while (not_at(A_ii, A.end_rows())) {
    typename MatrixA::RowVector::const_iterator A_j = (*A_ii).begin();
     while (not_at(A_j, (*A_ii).end())) {
       if (abs(*A_j - B(A_ii.index(), A_j.index())) > 0.00001) {
#if !defined(_MSVCPP_)
	 std::cout << "(" << A_ii.index() << "," << A_j.index() << ")="
	      << *A_j << " != " << B(A_ii.index(), A_j.index()) << std::endl;
#endif
	return false;
       }
      ++A_j;
    }
    ++A_ii;
  }
  return true;
}


template <class MatrixA, class MatrixB>
inline bool
__matrix_equal(const MatrixA& A, const MatrixB& B, column_tag)
{
  typename MatrixA::const_column_2Diterator A_jj = A.begin_columns();
  while (A_jj < A.end_columns()) {
    typename MatrixA::ColumnVector::const_iterator A_i = (*A_jj).begin();
     while (A_i != (*A_jj).end()) {
       if (abs(*A_i - B(A_i.index(), A_jj.index())) > 0.00001) {
#if !defined(_MSVCPP_)
		   std::cout << "(" << A_i.index() << "," << A_jj.index() << ")=" 
	      << *A_i << " != " << B(A_i.index(), A_jj.index()) << std::endl;
#endif
		   return false;
       }
      ++A_i;
    }
    ++A_jj;
  }
  return true;
}

template <class MatrixA, class MatrixB>
inline bool
matrix_row_equal(const MatrixA& A, const MatrixB& B) {
  return __matrix_equal(A, B, row_tag());
}

template <class MatrixA, class MatrixB>
inline bool
matrix_column_equal(const MatrixA& A, const MatrixB& B) {
  return __matrix_equal(A, B, column_tag());
}

///
template <class MatrixA, class MatrixB>
inline bool
matrix_equal(const MatrixA& A, const MatrixB& B)
{
  typedef typename MatrixA::orientation orien;
  return __matrix_equal(A, B, orien());
}
#endif

/*********************************************************/
/* Functions to get random numbers to fill the matrices  */
/*********************************************************/
  

inline float
make_rand_element(float)
{
  float r = float(rand());
  return r/float(RAND_MAX)*10.0;
}

inline double
make_rand_element(double)
{
  double r = double(rand());
  return r/double(RAND_MAX)*10.0L;
}

inline complex<float>
make_rand_element(complex<float>)
{
  return complex<float>(make_rand_element(float()), 
			make_rand_element(float()));
}

inline complex<double>
make_rand_element(complex<double>)
{
  return complex<double>(make_rand_element(double()),
			 make_rand_element(double()));
}

#if 0
template <class Orien>
inline void
print_matrix(const matrix<Orien>& mA)
{
#if 0
  const Orien& A = mA.two_d();
  Orien::const_row_2Diterator row_iter = A.begin_rows();
  std::cout << "[";
  while (row_iter < A.end_rows()) {
    Orien::RowVector::const_iterator ri = (*row_iter).begin();
    std::cout << "[";
     while (ri < (*row_iter).end()) {
      std::cout << *ri;
      std::cout << ",";
      ++ri;
    }
    std::cout << "]," << std::endl;;
    ++row_iter;
  }
  std::cout << "]" << std::endl;
#endif

#if 1
  // good one
  const Orien& A = mA.two_d();
  typename Orien::const_row_2Diterator row_iter = A.begin_rows();
  while (row_iter < A.end_rows()) {
    typename Orien::RowVector::const_iterator ri = (*row_iter).begin();
     while (ri < (*row_iter).end()) {
       if (*ri != 0)
	 std::cout << "(" << row_iter.index() << "," << ri.index()
	      << ") = " << *ri << std::endl;;
      ++ri;
    }
    ++row_iter;
  }
#endif

#if 0
  int i, j;
  for (i=0; i < mA.nrows(); ++i)
    for (j=0; j < mA.ncols(); ++j)
       if (mA(i,j) != 0)
	 std::cout << "(" << i << "," << j
	      << ") = " << mA(i,j) << std::endl;;
#endif

#if 0
  int i,j;
  std::cout << "[";
  for (i=0; i < A.nrows(); ++i) {
    std::cout << "[";
    for (j=0; j < A.ncols(); ++j) {
      std::cout << A(i,j);
      if (j < A.ncols() - 1)
	std::cout << ",";
    }
    std::cout << "]";
    if (i < A.nrows() - 1)
      std::cout << "," << std::endl;
  }
  std::cout << "]" << std::endl;
#endif
}

template <class Orien>
inline void
dense_init(matrix<Orien>& A)
{
  int i, j;
  for(i = 0; i < A.nrows(); ++i)
    for (j = 0; j < A.ncols(); ++j)
      A(i,j) = i * A.ncols() + j;
}

#endif

template <class Mat>
inline void
insert_zero_matrix(Mat& A)
{
  typedef typename matrix_traits<Mat>::size_type Int;
  // idea, for sparse matrices, just zero the non zeroes
  // which should remove them
  Int i, j;
  for (i = 0; i < A.nrows(); ++i) {
    Int first = MTL_MAX(0, int(i) - int(A.sub()));
    Int last = MTL_MIN(int(A.ncols()), int(i) + int(A.super()) + 1);
    for (j = 0; j < A.ncols(); ++j)
      if (j >= first && j < last)
	A(i,j) = 0.0;
  }

}

template <class Matrix>
inline void
zero_matrix(Matrix& A){
  typedef typename matrix_traits<Matrix>::value_type T;
  typename Matrix::iterator oneD_iter = A.begin();
  while (oneD_iter < A.end()) {
    typename Matrix::OneD::iterator i = (*oneD_iter).begin();
    while (i < (*oneD_iter).end()) {
      *i = T(0);
      ++i;
    }
    ++oneD_iter;
  }
}
  
#if 0
template <class Orien>
inline void
sparse_init(matrix<Orien>& A, int c)
{
  int i;
  for (i=0; i < A.nrows(); ++i)
    A(i, i) = 4;
  for (i=0; i < A.nrows() - 1; ++i) {
    A(i, i + 1) = -1.0;
    A(i + 1, i) = -1.0;
  }
  for (i=0; i < A.nrows() - 2; ++i) {
    A(i, i + 2) = -1.0;
    A(i + 2, i) = -1.0;
  }
  for (i=0; i < A.nrows() - 3; ++i) {
    A(i, i + 3) = -1.0;
    A(i + 3, i) = -1.0;
  }
  for (i=0; i < A.nrows() - c; ++i) {
    A(i, i + c) = -1.0;
    A(i + c, i) = -1.0;
  }
}


//
//  mat_mat_multiply
//

template <class OA, class OB, class OC>
inline void
calc_correct(const matrix<OA>& mA, const matrix<OB>& mB,
	     matrix<OC>& mC2)
{
  int i, j, k;
  for (k = 0; k < mA.ncols(); ++k)
    for (j = 0; j < mB.ncols(); ++j)
      for (i = 0; i < mA.nrows(); ++i) {
	mC2(i,j) += mA(i,k) * mB(k,j);	 
      }
}

template <class OA>
inline int calc_nonzeroes(const matrix<OA>& mA)
{
  int nonzeroes = 0;

  const OA& A = mA.two_d();

  typename OA::const_row_2Diterator row_iter = A.begin_rows();
  while (row_iter < A.end_rows()) {
    typename OA::RowVector::const_iterator ri = (*row_iter).begin();
    while (ri.not_at((*row_iter).end())) {
      if (*ri != 0.0)
	++nonzeroes;
      ++ri;
    } 
    ++row_iter;
  }

  return nonzeroes;
}
#endif


} /* namespace mtl */

#endif