FlashGraph-ng
A new frontier in large-scale graph analysis and data mining
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Pages
FM_MultiVector.h
1 #ifndef __FM_MULTIVECTOR_H__
2 #define __FM_MULTIVECTOR_H__
3 
4 /*
5  * Copyright 2014 Open Connectome Project (http://openconnecto.me)
6  * Written by Da Zheng (zhengda1936@gmail.com)
7  *
8  * This file is part of FlashMatrix.
9  *
10  * Licensed under the Apache License, Version 2.0 (the "License");
11  * you may not use this file except in compliance with the License.
12  * You may obtain a copy of the License at
13  *
14  * http://www.apache.org/licenses/LICENSE-2.0
15  *
16  * Unless required by applicable law or agreed to in writing, software
17  * distributed under the License is distributed on an "AS IS" BASIS,
18  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
19  * See the License for the specific language governing permissions and
20  * limitations under the License.
21  */
22 
23 #include <stdlib.h>
24 #include <time.h>
25 
26 #include <boost/format.hpp>
27 
28 #include "log.h"
29 
30 #ifndef TRILINOS_VERSION
31 #define TRILINOS_VERSION 12
32 #endif
33 
34 #include "AnasaziMultiVec.hpp"
35 
36 #include "NUMA_vector.h"
37 #include "matrix_config.h"
38 #include "generic_type.h"
39 
40 #ifdef FM_VERIFY
41 #include "Epetra_MultiVector.h"
42 #include "Epetra_SerialComm.h"
43 #include "Epetra_Map.h"
44 #include "AnasaziEpetraAdapter.hpp"
45 #endif
46 #include "block_dense_matrix.h"
47 
48 #include "mem_matrix_store.h"
49 #include "dense_matrix.h"
50 #include "dotp_matrix_store.h"
51 #include "matrix_stats.h"
52 
53 static int MV_id;
54 
55 namespace fm
56 {
57 
58 namespace eigen
59 {
60 
61 static std::string get_curr_time_str()
62 {
63  time_t curr = time(NULL);
64  std::string ret = ctime(&curr);
65  ret[ret.length() - 1] = 0;
66  return ret;
67 }
68 
69 template<class ScalarType>
70 class FM_MultiVector: public Anasazi::MultiVec<ScalarType>
71 {
72  size_t subspace_size;
73  std::string solver;
74  bool in_mem;
75  std::string name;
76  block_multi_vector::ptr mat;
77 #ifdef FM_VERIFY
78  std::shared_ptr<Anasazi::EpetraMultiVec> ep_mat;
79 #endif
80 
81  bool is_subspace(size_t num_cols) const {
82  // This works for KrylovSchur. TODO I need to test other eigensolvers.
83  return solver == "KrylovSchur" && num_cols >= subspace_size;
84  }
85 
86  FM_MultiVector(const std::string &extra, bool in_mem,
87  const std::string &solver, size_t subspace_size) {
88  char name_buf[128];
89  snprintf(name_buf, sizeof(name_buf), "MV-%d", MV_id++);
90  this->name = std::string(name_buf) + " " + extra;
91  this->in_mem = in_mem;
92  this->solver = solver;
93  this->subspace_size = subspace_size;
94  }
95 public:
96  FM_MultiVector(size_t num_rows, size_t num_cols, size_t block_size,
97  size_t subspace_size, bool in_mem, const std::string &solver) {
98  this->subspace_size = subspace_size;
99  this->solver = solver;
100  this->in_mem = in_mem;
101  // We don't materialize the column matrix.
102  mat = block_multi_vector::create(num_rows, num_cols, block_size,
103  fm::get_scalar_type<ScalarType>(), in_mem, is_subspace(num_cols));
104 
105 #ifdef FM_VERIFY
106  Epetra_SerialComm Comm;
107  Teuchos::RCP<Epetra_Map> Map = Teuchos::rcp (
108  new Epetra_Map((int) num_rows, 0, Comm));
109  ep_mat = std::shared_ptr<Anasazi::EpetraMultiVec>(
110  new Anasazi::EpetraMultiVec(*Map, num_cols));
111 #endif
112 
113  char name_buf[128];
114  snprintf(name_buf, sizeof(name_buf), "MV-%d", MV_id++);
115  this->name = name_buf;
116  }
117 
118  void verify() const {
119 #ifdef FM_VERIFY
120  int len = ep_mat->GlobalLength();
121  assert((size_t) len == mat->get_num_rows());
122  assert(mat->get_num_cols() == (size_t) ep_mat->NumVectors());
123  size_t ep_col_idx = 0;
124  for (size_t i = 0; i < mat->get_num_blocks(); i++) {
125  dense_matrix::ptr block = mat->get_block(i);
126  detail::local_matrix_store::const_ptr portion
127  = block->get_data().get_portion(0);
128  assert(block->get_num_rows() == portion->get_num_rows());
129  assert(block->get_num_cols() == portion->get_num_cols());
130  for (size_t j = 0; j < portion->get_num_cols(); j++) {
131  for (size_t k = 0; k < portion->get_num_rows(); k++) {
132  double v1 = (*ep_mat)[ep_col_idx][k];
133  double v2 = portion->get<ScalarType>(k, j);
134  assert(abs(v1 - v2) == 0);
135  }
136  ep_col_idx++;
137  }
138  }
139 #endif
140  }
141 
142 #ifdef FM_VERIFY
143  Anasazi::EpetraMultiVec &get_ep_mv() {
144  return *ep_mat;
145  }
146 
147  const Anasazi::EpetraMultiVec &get_ep_mv() const {
148  return *ep_mat;
149  }
150 #endif
151 
152  block_multi_vector::ptr get_data() {
153  return mat;
154  }
155 
156  block_multi_vector::ptr get_data() const {
157  return mat;
158  }
159 
160  size_t get_block_size() const {
161  return mat->get_block_size();
162  }
163 
164  void Random () {
165  mat->init_rand<ScalarType>(-1, 1);
166  sync_fm2ep();
167  verify();
168  }
169 
170  void sync_ep2fm() {
171 #ifdef FM_VERIFY
172  assert(mat->get_num_cols() == (size_t) ep_mat->NumVectors());
173  // This is a temporary solution.
174  for (size_t i = 0; i < mat->get_num_blocks(); i++) {
175  fm::detail::mem_matrix_store::ptr store
176  = fm::detail::mem_matrix_store::create(mat->get_num_rows(),
177  mat->get_block_size(), fm::matrix_layout_t::L_COL,
178  mat->get_type(), mat->get_num_nodes());
179  fm::detail::mem_col_matrix_store::ptr col_store
180  = std::static_pointer_cast<fm::detail::mem_col_matrix_store>(store);
181  for (size_t j = 0; j < mat->get_block_size(); j++) {
182  size_t col_idx = i * mat->get_block_size() + j;
183  memcpy(col_store->get_col(j), (*ep_mat)[col_idx],
184  mat->get_num_rows() * mat->get_entry_size());
185  }
186  mat->set_block(i, fm::dense_matrix::create(col_store));
187  }
188 #endif
189  }
190 
191  void sync_fm2ep() {
192 #ifdef FM_VERIFY
193  BOOST_LOG_TRIVIAL(info) << boost::format("There are %1% blocks")
194  % mat->get_num_blocks();
195  size_t ep_col_idx = 0;
196  for (size_t i = 0; i < mat->get_num_blocks(); i++) {
197  dense_matrix::ptr block = mat->get_block(i);
198  if (block == NULL)
199  continue;
200 
201  detail::local_matrix_store::const_ptr portion
202  = block->get_data().get_portion(0);
203  assert(block->get_num_rows() == portion->get_num_rows());
204  assert(block->get_num_cols() == portion->get_num_cols());
205  for (size_t j = 0; j < portion->get_num_cols(); j++) {
206  for (size_t k = 0; k < portion->get_num_rows(); k++)
207  (*ep_mat)[ep_col_idx][k] = portion->get<ScalarType>(k, j);
208  ep_col_idx++;
209  }
210  }
211 #endif
212  }
213 
214  std::string get_name() const {
215  return name;
216  }
217 
219 
222  virtual Anasazi::MultiVec<ScalarType> * Clone(const int numvecs) const {
223  FM_MultiVector<ScalarType> *ret;
224  if (numvecs % mat->get_block_size() == 0)
225  ret = new FM_MultiVector<ScalarType>(mat->get_num_rows(), numvecs,
226  mat->get_block_size(), subspace_size, in_mem, solver);
227  else
228  ret = new FM_MultiVector<ScalarType>(
229  mat->get_num_rows(), numvecs, numvecs, subspace_size,
230  in_mem, solver);
231  BOOST_LOG_TRIVIAL(info) << boost::format("create new %1% (#cols: %2%)")
232  % ret->get_name() % numvecs;
233  return ret;
234  }
235 
238  virtual Anasazi::MultiVec<ScalarType> * CloneCopy () const {
239  num_col_writes_concept += mat->get_num_cols();
240  num_col_reads_concept += mat->get_num_cols();
241 
242  BOOST_LOG_TRIVIAL(info) << boost::format("deep copy %1% (#cols: %2%)")
243  % name % mat->get_num_cols();
244  std::string extra = std::string("(deep copy from ") + get_name() + ")";
245  FM_MultiVector<ScalarType> *ret = new FM_MultiVector<ScalarType>(extra,
246  in_mem, solver, subspace_size);
247  ret->mat = this->mat->clone();
248 #ifdef FM_VERIFY
249  ret->ep_mat = std::shared_ptr<Anasazi::EpetraMultiVec>(
250  dynamic_cast<Anasazi::EpetraMultiVec *>(ep_mat->CloneCopy()));
251 #endif
252  verify();
253  ret->verify();
254  return ret;
255  }
256 
263  virtual Anasazi::MultiVec<ScalarType> * CloneCopy (
264  const std::vector<int>& index) const {
265  num_col_writes_concept += index.size();
266  num_col_reads_concept += index.size();
267 
268  BOOST_LOG_TRIVIAL(info) << boost::format("deep copy sub %1% (#cols: %2%)")
269  % name % index.size();
270  std::string extra = std::string("(deep copy from sub ") + get_name() + ")";
271  FM_MultiVector<ScalarType> *ret = new FM_MultiVector<ScalarType>(extra,
272  in_mem, solver, subspace_size);
273  ret->mat = mat->get_cols(index);
274 #ifdef FM_VERIFY
275  ret->ep_mat = std::shared_ptr<Anasazi::EpetraMultiVec>(
276  dynamic_cast<Anasazi::EpetraMultiVec *>(ep_mat->CloneCopy(index)));
277 #endif
278  verify();
279  ret->verify();
280  return ret;
281  }
282 
289  virtual Anasazi::MultiVec<ScalarType> * CloneViewNonConst (
290  const std::vector<int>& index) {
291  std::string extra = std::string("(") + get_name() + "[" + vec2str(index) + "])";
292  FM_MultiVector<ScalarType> *ret = new FM_MultiVector<ScalarType>(extra,
293  in_mem, solver, subspace_size);
294  BOOST_LOG_TRIVIAL(info) << boost::format("view %1% (#cols: %2%)")
295  % ret->name % index.size();
296  ret->mat = mat->get_cols_mirror(index);
297 #ifdef FM_VERIFY
298  ret->ep_mat = std::shared_ptr<Anasazi::EpetraMultiVec>(
299  dynamic_cast<Anasazi::EpetraMultiVec *>(ep_mat->CloneViewNonConst(index)));
300 #endif
301 // verify();
302 // ret->verify();
303  return ret;
304  }
305 
306  static std::string vec2str(const std::vector<int> &index) {
307  std::vector<int> copy = index;
308  std::sort(copy.begin(), copy.end());
309  typedef std::pair<int, int> range_t;
310  std::vector<range_t> ranges;
311  ranges.push_back(range_t(copy.front(), copy.front()));
312  for (size_t i = 1; i < copy.size(); i++) {
313  if (ranges.back().second + 1 == copy[i])
314  ranges.back().second++;
315  else
316  ranges.push_back(range_t(copy[i], copy[i]));
317  }
318  std::string ret = "";
319  for (size_t i = 0; i < ranges.size(); i++) {
320  if (ranges[i].first == ranges[i].second)
321  ret += itoa(ranges[i].first);
322  else
323  ret += std::string(std::string(itoa(ranges[i].first)) + ":") + itoa(ranges[i].second);
324  }
325  return ret;
326  }
327 
334  virtual const Anasazi::MultiVec<ScalarType> * CloneView (
335  const std::vector<int>& index) const {
336  std::string extra = std::string("(const ") + get_name() + "[" + vec2str(index) + "])";
337  FM_MultiVector<ScalarType> *ret = new FM_MultiVector<ScalarType>(extra,
338  in_mem, solver, subspace_size);
339  if (index.size() == 1 && mat->get_block_size() > 1
340  && (solver == "KrylovSchur" || solver == "Davidson")) {
341  fm::detail::matrix_stats_t orig_stats = detail::matrix_stats;
342  size_t block_idx = index[0] / mat->get_block_size();
343  dense_matrix::ptr block = mat->get_block(block_idx);
344  const dotp_matrix_store *store
345  = dynamic_cast<const dotp_matrix_store *>(block->get_raw_store().get());
346  if (store == NULL) {
347  BOOST_LOG_TRIVIAL(info) << boost::format(
348  "clone view: materialize %1% on the fly")
349  % block->get_data().get_name();
350  dotp_matrix_store::ptr dotp = dotp_matrix_store::create(
351  mat->get_block(block_idx)->get_raw_store());
352  mat->set_block(block_idx, fm::dense_matrix::create(dotp));
353  }
354  fm::detail::matrix_stats.print_diff(orig_stats);
355  }
356  BOOST_LOG_TRIVIAL(info) << boost::format("const view %1% (#cols: %2%)")
357  % ret->name % index.size();
358  ret->mat = mat->get_cols(index);
359 #ifdef FM_VERIFY
360  ret->ep_mat = std::shared_ptr<Anasazi::EpetraMultiVec>(
361  dynamic_cast<Anasazi::EpetraMultiVec *>(ep_mat->CloneViewNonConst(index)));
362 #endif
363 // verify();
364  ret->verify();
365  return ret;
366  }
367 
369 
371  ANASAZI_DEPRECATED virtual int GetVecLength () const {
372  return mat->get_num_rows();
373  }
374 
377  virtual ptrdiff_t GetGlobalLength () const {
378  return mat->get_num_rows();
379  }
380 
382  virtual int GetNumberVecs () const {
383  return mat->get_num_cols();
384  }
385 
387 
389  virtual void MvTimesMatAddMv (ScalarType alpha,
390  const Anasazi::MultiVec<ScalarType>& A,
391  const Teuchos::SerialDenseMatrix<int,ScalarType>& B, ScalarType beta) {
392  struct timeval start, end;
393  gettimeofday(&start, NULL);
394  sync_fm2ep();
395  const FM_MultiVector &fm_A = dynamic_cast<const FM_MultiVector &>(A);
396  BOOST_LOG_TRIVIAL(info) << get_curr_time_str() << ":" << boost::format(
397  "this(%1%) = %2% * A(%3%) * B(%4%x%5%) + %6% * this")
398  % name % alpha % fm_A.name % B.numRows() % B.numCols() % beta;
399  if (alpha != 0)
400  num_col_reads_concept += fm_A.mat->get_num_cols();
401  if (beta != 0)
402  num_col_reads_concept += mat->get_num_cols();
403  num_col_writes_concept += mat->get_num_cols();
404  num_multiply_concept
405  += fm_A.mat->get_num_rows() * fm_A.mat->get_num_cols() * B.numCols();
406 
407  fm::detail::mem_col_matrix_store::ptr Bstore
408  = fm::detail::mem_col_matrix_store::create(
409  B.numRows(), B.numCols(), fm::get_scalar_type<ScalarType>());
410  for (int i = 0; i < B.numRows(); i++) {
411  for (int j = 0; j < B.numCols(); j++)
412  Bstore->set<ScalarType>(i, j, B(i, j));
413  }
414  fm::scalar_variable_impl<ScalarType> alpha_var(alpha);
415  fm::scalar_variable_impl<ScalarType> beta_var(beta);
416  this->mat->assign(*this->mat->gemm(*fm_A.mat, Bstore, alpha_var, beta_var));
417 #ifdef FM_VERIFY
418  this->ep_mat->MvTimesMatAddMv(alpha, *fm_A.ep_mat, B, beta);
419 #endif
420  fm_A.verify();
421  this->verify();
422  gettimeofday(&end, NULL);
423  BOOST_LOG_TRIVIAL(info) << "gemm takes " << time_diff(start, end)
424  << " seconds";
425  // In KrylovSchur, if the input MV isn't in the subspace, we can drop
426  // the dense matrix in MV to save memory.
427  if (solver == "KrylovSchur" && fm_A.mat->get_num_blocks() == 1
428  && fm_A.mat->get_block(0) != NULL
429  && fm_A.mat->get_block(0)->is_in_mem() && beta == 0) {
430  BOOST_LOG_TRIVIAL(info) << "Drop the matrix to save space";
431  fm_A.mat->set_block(0, fm::dense_matrix::const_ptr());
432  }
433  }
434 
436  virtual void MvAddMv ( ScalarType alpha, const Anasazi::MultiVec<ScalarType>& A,
437  ScalarType beta, const Anasazi::MultiVec<ScalarType>& B ) {
438  sync_fm2ep();
439  const FM_MultiVector &fm_A = dynamic_cast<const FM_MultiVector &>(A);
440  const FM_MultiVector &fm_B = dynamic_cast<const FM_MultiVector &>(B);
441  BOOST_LOG_TRIVIAL(info) << get_curr_time_str() << ":" << boost::format(
442  "this(%1%) = %2% * A(%3%) + %4% * B(%5%)")
443  % name % alpha % fm_A.name % beta % fm_B.name;
444  num_col_reads_concept += fm_A.mat->get_num_cols() + fm_B.mat->get_num_cols();
445  num_col_writes_concept += mat->get_num_cols();
446 
447  if (alpha == 1 && beta == 0)
448  this->mat->assign(*fm_A.mat);
449  else if (alpha == 0 && beta == 1)
450  this->mat->assign(*fm_B.mat);
451  else {
452  block_multi_vector::ptr aA = fm_A.mat->multiply_scalar(alpha);
453  block_multi_vector::ptr bB = fm_B.mat->multiply_scalar(beta);
454  this->mat->assign(*aA->add(*bB));
455  }
456 #ifdef FM_VERIFY
457  this->ep_mat->MvAddMv(alpha, *fm_A.ep_mat, beta, *fm_B.ep_mat);
458 #endif
459  fm_A.verify();
460  fm_B.verify();
461  this->verify();
462  }
463 
465  virtual void MvScale ( ScalarType alpha ) {
466  num_col_writes_concept += mat->get_num_cols();
467  num_col_reads_concept += mat->get_num_cols();
468  num_multiply_concept += mat->get_num_rows() * mat->get_num_cols();
469 
470  sync_fm2ep();
471  BOOST_LOG_TRIVIAL(info) << get_curr_time_str() << ":" << boost::format(
472  "this(%1%) *= %2%") % name % alpha;
473  mat->assign(*mat->multiply_scalar<ScalarType>(alpha));
474 #ifdef FM_VERIFY
475  ep_mat->MvScale(alpha);
476 #endif
477  this->verify();
478  }
479 
481  virtual void MvScale ( const std::vector<ScalarType>& alpha ) {
482  num_col_writes_concept += mat->get_num_cols();
483  num_col_reads_concept += mat->get_num_cols();
484  num_multiply_concept += mat->get_num_rows() * mat->get_num_cols();
485 
486  sync_fm2ep();
487  BOOST_LOG_TRIVIAL(info) << get_curr_time_str() << ":" << boost::format(
488  "this(%s) *= vec") % name;
489  mat->assign(*mat->scale_cols<ScalarType>(alpha));
490 #ifdef FM_VERIFY
491  ep_mat->MvScale(alpha);
492 #endif
493  this->verify();
494  }
495 
499  virtual void MvTransMv ( ScalarType alpha, const Anasazi::MultiVec<ScalarType>& A,
500  Teuchos::SerialDenseMatrix<int,ScalarType>& B
501 #ifdef HAVE_ANASAZI_EXPERIMENTAL
502  , ConjType conj = Anasazi::CONJ
503 #endif
504  ) const {
505  struct timeval start, end;
506  gettimeofday(&start, NULL);
507 
508  const FM_MultiVector &fm_A = dynamic_cast<const FM_MultiVector &>(A);
509  BOOST_LOG_TRIVIAL(info) << get_curr_time_str() << ":" << boost::format(
510  "B(%1%x%2%) = %3% * A(%4%)^T * this(%5%)")
511  % B.numRows() % B.numCols() % alpha % fm_A.name % name;
512  num_col_reads_concept += fm_A.mat->get_num_cols() + mat->get_num_cols();
513  num_multiply_concept
514  += fm_A.mat->get_num_cols() * mat->get_num_rows() * mat->get_num_cols();
515 
516  assert((size_t) B.numRows() == fm_A.mat->get_num_cols());
517  assert((size_t) B.numCols() == this->mat->get_num_cols());
518  assert(fm_A.mat->get_num_rows() == this->mat->get_num_rows());
519  this->verify();
520 
521  long double lalpha = alpha;
522  // This is a special case: we compute the dot product of a col vector
523  // with itself.
524  if (fm_A.mat->get_num_cols() == 1 && mat->get_num_cols() == 1) {
525  const sub_dotp_matrix_store *col1
526  = dynamic_cast<const sub_dotp_matrix_store *>(
527  fm_A.mat->get_block(0)->get_raw_store().get());
528  const sub_dotp_matrix_store *col2
529  = dynamic_cast<const sub_dotp_matrix_store *>(
530  mat->get_block(0)->get_raw_store().get());
531  if (col1 && col2
532  && col1->get_orig_store() == col2->get_orig_store()) {
533  B(0, 0) = col1->get_col_dot(0) * lalpha;
534  return;
535  }
536  }
537 
538  fm::dense_matrix::ptr res = mat->MvTransMv(*fm_A.mat);
539  const_cast<FM_MultiVector *>(this)->sync_fm2ep();
540 #if 0
541  ep_mat->MvTransMv(alpha, *fm_A.ep_mat, B);
542  for (int i = 0; i < B.numRows(); i++)
543  for (int j = 0; j < B.numCols(); j++) {
544  printf("%g\n", B(i, j) - (double) (res->get<ScalarType>(i, j) * lalpha));
545  assert(B(i, j) == (double) (res->get<ScalarType>(i, j) * lalpha));
546  }
547 #endif
548  const fm::detail::mem_matrix_store &mem_res
549  = dynamic_cast<const fm::detail::mem_matrix_store &>(res->get_data());
550  for (int i = 0; i < B.numRows(); i++) {
551  for (int j = 0; j < B.numCols(); j++) {
552  B(i, j) = mem_res.get<ScalarType>(i, j) * lalpha;
553  }
554  }
555  gettimeofday(&end, NULL);
556  BOOST_LOG_TRIVIAL(info) << "MvTransMv takes " << time_diff(start, end)
557  << " seconds";
558  }
559 
565  virtual void MvDot ( const Anasazi::MultiVec<ScalarType>& A,
566  std::vector<ScalarType> & b
567 #ifdef HAVE_ANASAZI_EXPERIMENTAL
568  , ConjType conj = Anasazi::CONJ
569 #endif
570  ) const {
571  struct timeval start, end;
572  gettimeofday(&start, NULL);
573 
574  const FM_MultiVector &fm_A = dynamic_cast<const FM_MultiVector &>(A);
575  BOOST_LOG_TRIVIAL(info) << get_curr_time_str() << ":"
576  << boost::format("MvDot(A(%1%), this(%2%))") % fm_A.name % name;
577  b = mat->MvDot(*fm_A.mat);
578  gettimeofday(&end, NULL);
579  BOOST_LOG_TRIVIAL(info) << "MvDot takes " << time_diff(start, end)
580  << " seconds";
581  }
582 
584 
589  virtual void MvNorm (
590  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> & normvec) const {
591  struct timeval start, end;
592  gettimeofday(&start, NULL);
593 
594  num_col_reads_concept += mat->get_num_cols();
595  num_multiply_concept += mat->get_num_rows() * mat->get_num_cols();
596 
597  BOOST_LOG_TRIVIAL(info) << get_curr_time_str() << ":" << boost::format(
598  "norm(%1%)(#cols: %2%)") % name % normvec.size();
599  verify();
600  fm::detail::matrix_stats_t orig_stats = fm::detail::matrix_stats;
601  for (size_t i = 0; i < mat->get_num_blocks(); i++) {
602  BOOST_LOG_TRIVIAL(info) << boost::format(
603  "norm: materialize %1% on the fly")
604  % mat->get_block(i)->get_data().get_name();
605  dotp_matrix_store::ptr dotp
606  = dotp_matrix_store::create(mat->get_block(i)->get_raw_store());
607  mat->set_block(i, fm::dense_matrix::create(dotp));
608  std::vector<ScalarType> col_dots = dotp->get_col_dot_prods();
609  // TODO do I need long double here?
610  for (size_t j = 0; j < col_dots.size(); j++)
611  normvec[i * mat->get_block_size() + j] = std::sqrt(col_dots[j]);
612  }
613  fm::detail::matrix_stats.print_diff(orig_stats);
614  const_cast<FM_MultiVector *>(this)->sync_fm2ep();
615 #ifdef FM_VERIFY
616  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> normvec1(
617  normvec.size());
618  ep_mat->MvNorm(normvec1);
619  for (size_t i = 0; i < normvec.size(); i++) {
620  BOOST_LOG_TRIVIAL(info) << boost::format("%1%: %2%")
621  % i % (normvec[i] - normvec1[i]);
622  assert(normvec[i] == normvec1[i]);
623  }
624 #endif
625  gettimeofday(&end, NULL);
626  BOOST_LOG_TRIVIAL(info) << "MvNorm takes " << time_diff(start, end)
627  << " seconds";
628  }
629 
631 
636  virtual void SetBlock (const Anasazi::MultiVec<ScalarType>& A,
637  const std::vector<int>& index) {
638  struct timeval start, end;
639  gettimeofday(&start, NULL);
640 
641  num_col_reads_concept += index.size();
642  num_col_writes_concept += index.size();
643 
644  assert((size_t) A.GetNumberVecs() == index.size());
645  const FM_MultiVector &fm_A = dynamic_cast<const FM_MultiVector &>(A);
646  BOOST_LOG_TRIVIAL(info) << get_curr_time_str() << ":" << boost::format(
647  "this(%1%)[%2%] = A(%3%)") % name % vec2str(index) % fm_A.name;
648  sync_fm2ep();
649  this->mat->set_block(*fm_A.mat, index);
650 #ifdef FM_VERIFY
651  this->ep_mat->SetBlock(*fm_A.ep_mat, index);
652 #endif
653  fm_A.verify();
654  this->verify();
655  gettimeofday(&end, NULL);
656  BOOST_LOG_TRIVIAL(info) << "SetBlock takes " << time_diff(start, end)
657  << " seconds";
658  }
659 
661  virtual void MvRandom () {
662  struct timeval start, end;
663  gettimeofday(&start, NULL);
664 
665  num_col_writes_concept += mat->get_num_cols();
666 
667  BOOST_LOG_TRIVIAL(info) << get_curr_time_str() << ":" << boost::format(
668  "this(%1%) = random") % name;
669  mat->init_rand<ScalarType>(-1, 1);
670 // ep_mat->MvRandom();
671  sync_fm2ep();
672  this->verify();
673  gettimeofday(&end, NULL);
674  BOOST_LOG_TRIVIAL(info) << "MvRandom takes " << time_diff(start, end)
675  << " seconds";
676  }
677 
679  virtual void MvInit ( ScalarType alpha ) {
680  // TODO
681  assert(0);
682  }
683 
685 
687  virtual void MvPrint ( std::ostream& os ) const {
688  // TODO
689  assert(0);
690  }
691 
692 #ifdef HAVE_ANASAZI_TSQR
693 
717  virtual void factorExplicit (Anasazi::MultiVec<ScalarType>& Q,
718  Teuchos::SerialDenseMatrix<int, ScalarType>& R,
719  const bool forceNonnegativeDiagonal=false) {
720  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "The Anasazi::MultiVec<"
721  << Teuchos::TypeNameTraits<ScalarType>::name() << "> subclass which you "
722  "are using does not implement the TSQR-related method factorExplicit().");
723  }
724 
759  virtual int revealRank (Teuchos::SerialDenseMatrix<int, ScalarType>& R,
760  const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& tol) {
761  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "The Anasazi::MultiVec<"
762  << Teuchos::TypeNameTraits<ScalarType>::name() << "> subclass which you "
763  "are using does not implement the TSQR-related method revealRank().");
764  }
765 
766 #endif // HAVE_ANASAZI_TSQR
767 };
768 
769 }
770 
771 }
772 
773 namespace Anasazi
774 {
775 
776 // \brief Specialization of MultiVecTraits for Belos::MultiVec.
777 //
778 // Anasazi interfaces to every multivector implementation through a
779 // specialization of MultiVecTraits. Thus, we provide a
780 // specialization of MultiVecTraits for the MultiVec run-time
781 // polymorphic interface above.
782 //
783 // \tparam ScalarType The type of entries in the multivector; the
784 // template parameter of MultiVec.
785 template<class ScalarType>
786 class MultiVecTraits<ScalarType, fm::eigen::FM_MultiVector<ScalarType> > {
787 public:
789 
790 
791  // \brief Create a new empty \c MultiVec containing \c numvecs columns.
792  // \return Reference-counted pointer to the new \c MultiVec.
793  static Teuchos::RCP<fm::eigen::FM_MultiVector<ScalarType> > Clone (
794  const fm::eigen::FM_MultiVector<ScalarType>& mv, const int numvecs) {
795  return Teuchos::rcp ((fm::eigen::FM_MultiVector<ScalarType> *)
796  const_cast<fm::eigen::FM_MultiVector<ScalarType>&> (mv).Clone(numvecs));
797  }
798 
799  /*
800  * \brief Creates a new \c Anasazi::MultiVec and copies contents of \c mv
801  * into the new vector (deep copy).
802  \return Reference-counted pointer to the new \c Anasazi::MultiVec.
803  */
804  static Teuchos::RCP<fm::eigen::FM_MultiVector<ScalarType> > CloneCopy(
805  const fm::eigen::FM_MultiVector<ScalarType>& mv) {
806  return Teuchos::rcp((fm::eigen::FM_MultiVector<ScalarType> *)
807  const_cast<fm::eigen::FM_MultiVector<ScalarType>&>(mv).CloneCopy());
808  }
809 
810  /* \brief Creates a new \c Anasazi::MultiVec and copies the selected
811  * contents of \c mv into the new vector (deep copy).
812 
813  The copied vectors from \c mv are indicated by the \c index.size() indices in \c index.
814  \return Reference-counted pointer to the new \c Anasazi::MultiVec.
815  */
816  static Teuchos::RCP<fm::eigen::FM_MultiVector<ScalarType> > CloneCopy(
817  const fm::eigen::FM_MultiVector<ScalarType>& mv,
818  const std::vector<int>& index) {
819  return Teuchos::rcp((fm::eigen::FM_MultiVector<ScalarType> *)
820  const_cast<fm::eigen::FM_MultiVector<ScalarType>&>(
821  mv).CloneCopy(index));
822  }
823 
824  /* \brief Creates a new \c Anasazi::MultiVec that shares the selected
825  * contents of \c mv (shallow copy).
826 
827  The index of the \c numvecs vectors shallow copied from \c mv are indicated
828  by the indices given in \c index.
829  \return Reference-counted pointer to the new \c Anasazi::MultiVec.
830  */
831  static Teuchos::RCP<fm::eigen::FM_MultiVector<ScalarType> > CloneViewNonConst(
832  fm::eigen::FM_MultiVector<ScalarType>& mv,
833  const std::vector<int>& index) {
834  return Teuchos::rcp((fm::eigen::FM_MultiVector<ScalarType> *)
835  mv.CloneViewNonConst(index));
836  }
837 
838  /* \brief Creates a new const \c Anasazi::MultiVec that shares
839  * the selected contents of \c mv (shallow copy).
840 
841  The index of the \c numvecs vectors shallow copied from \c mv are
842  indicated by the indices given in \c index.
843  \return Reference-counted pointer to the new const \c Anasazi::MultiVec.
844  */
845  static Teuchos::RCP<const fm::eigen::FM_MultiVector<ScalarType> > CloneView(
846  const fm::eigen::FM_MultiVector<ScalarType>& mv,
847  const std::vector<int>& index) {
848  return Teuchos::rcp((fm::eigen::FM_MultiVector<ScalarType> *)
849  const_cast<fm::eigen::FM_MultiVector<ScalarType>&>(
850  mv).CloneView(index));
851  }
852 
854 
856 
857 
858  // Obtain the vector length of \c mv.
859  ANASAZI_DEPRECATED static int GetVecLength(
860  const fm::eigen::FM_MultiVector<ScalarType>& mv) {
861  return mv.GetGlobalLength();
862  }
863 
864  // Obtain the number of vectors in \c mv
865  static int GetNumberVecs(const fm::eigen::FM_MultiVector<ScalarType>& mv) {
866  return mv.GetNumberVecs();
867  }
868 
870 
871  // @name Update methods
873 
874  /* \brief Update \c mv with \f$ \alpha AB + \beta mv \f$.
875  */
876  static void MvTimesMatAddMv( ScalarType alpha,
877  const fm::eigen::FM_MultiVector<ScalarType>& A,
878  const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
879  ScalarType beta, fm::eigen::FM_MultiVector<ScalarType>& mv )
880  { mv.MvTimesMatAddMv(alpha, A, B, beta); }
881 
882  /* \brief Replace \c mv with \f$\alpha A + \beta B\f$.
883  */
884  static void MvAddMv( ScalarType alpha,
885  const fm::eigen::FM_MultiVector<ScalarType>& A,
886  ScalarType beta, const fm::eigen::FM_MultiVector<ScalarType>& B,
887  fm::eigen::FM_MultiVector<ScalarType>& mv)
888  { mv.MvAddMv(alpha, A, beta, B); }
889 
890  /* \brief Compute a dense matrix \c B through the matrix-matrix multiply \f$ \alpha A^Tmv \f$.
891  */
892  static void MvTransMv( ScalarType alpha,
893  const fm::eigen::FM_MultiVector<ScalarType>& A,
894  const fm::eigen::FM_MultiVector<ScalarType>& mv,
895  Teuchos::SerialDenseMatrix<int, ScalarType>& B
896 #ifdef HAVE_ANASAZI_EXPERIMENTAL
897  , ConjType conj = Anasazi::CONJ
898 #endif
899  )
900  { mv.MvTransMv(alpha, A, B
901 #ifdef HAVE_ANASAZI_EXPERIMENTAL
902  , conj
903 #endif
904  ); }
905 
906  /* \brief Compute a vector \c b where the components are the individual dot-products of the \c i-th columns of \c A and \c mv, i.e.\f$b[i] = A[i]^H mv[i]\f$.
907  */
908  static void MvDot( const fm::eigen::FM_MultiVector<ScalarType>& mv,
909  const fm::eigen::FM_MultiVector<ScalarType>& A,
910  std::vector<ScalarType> & b
911 #ifdef HAVE_ANASAZI_EXPERIMENTAL
912  , ConjType conj = Anasazi::CONJ
913 #endif
914  )
915  { mv.MvDot( A, b
916 #ifdef HAVE_ANASAZI_EXPERIMENTAL
917  , conj
918 #endif
919  ); }
920 
921  // Scale each element of the vectors in \c *this with \c alpha.
922  static void MvScale ( fm::eigen::FM_MultiVector<ScalarType>& mv, ScalarType alpha )
923  { mv.MvScale( alpha ); }
924 
926  static void MvScale ( fm::eigen::FM_MultiVector<ScalarType>& mv,
927  const std::vector<ScalarType>& alpha )
928  { mv.MvScale( alpha ); }
929 
931  // @name Norm method
933 
934  /* \brief Compute the 2-norm of each individual vector of \c mv.
935  Upon return, \c normvec[i] holds the value of \f$||mv_i||_2\f$, the \c i-th column of \c mv.
936  */
937  static void MvNorm( const fm::eigen::FM_MultiVector<ScalarType>& mv,
938  std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> & normvec)
939  { mv.MvNorm(normvec); }
940 
942  // @name Initialization methods
944  /* \brief Copy the vectors in \c A to a set of vectors in \c mv indicated by the indices given in \c index.
945 
946  The \c numvecs vectors in \c A are copied to a subset of vectors in \c mv indicated by the indices given in \c index,
947  i.e.<tt> mv[index[i]] = A[i]</tt>.
948  */
949  static void SetBlock( const fm::eigen::FM_MultiVector<ScalarType>& A,
950  const std::vector<int>& index,
951  fm::eigen::FM_MultiVector<ScalarType>& mv )
952  { mv.SetBlock(A, index); }
953 
954  /* \brief Replace the vectors in \c mv with random vectors.
955  */
956  static void MvRandom( fm::eigen::FM_MultiVector<ScalarType>& mv )
957  { mv.MvRandom(); }
958 
959  /* \brief Replace each element of the vectors in \c mv with \c alpha.
960  */
961  static void MvInit( fm::eigen::FM_MultiVector<ScalarType>& mv,
962  ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero() )
963  { mv.MvInit(alpha); }
964 
966  // @name Print method
968 
969  // Print the \c mv multi-vector to the \c os output stream.
970  static void MvPrint( const fm::eigen::FM_MultiVector<ScalarType>& mv,
971  std::ostream& os )
972  { mv.MvPrint(os); }
973 
975 
976 #ifdef HAVE_ANASAZI_TSQR
977  // \typedef tsqr_adaptor_type
978  // \brief TSQR adapter for MultiVec.
979  //
980  // Our TSQR adapter for MultiVec calls MultiVec's virtual
981  // methods. If you want to use TSQR with your MultiVec subclass,
982  // you must implement these methods yourself, as the default
983  // implementations throw std::logic_error.
984  typedef details::MultiVecTsqrAdapter<ScalarType> tsqr_adaptor_type;
985 #endif // HAVE_ANASAZI_TSQR
986 
987 #if TRILINOS_VERSION >= 12
988  // Obtain the vector length of \c mv.
989  // \note This method supersedes GetVecLength, which will be deprecated.
990  static ptrdiff_t GetGlobalLength( const fm::eigen::FM_MultiVector<ScalarType>& mv ) {
991  return mv.GetGlobalLength();
992  }
993 #endif
994 };
995 
996 #if TRILINOS_VERSION < 12
997 
998 // \brief An extension of the MultiVecTraits class that adds a new vector length method.
999 // \ingroup anasazi_opvec_interfaces
1000 //
1001 // This traits class provides an additional method to the multivector
1002 // operations for finding the number of rows that is 64-bit compatible.
1003 // The method in this traits class will replace the GetVecLength()
1004 // method, which will be deprecated, and removed in the next major
1005 // Trilinos release. At this time, this traits class will call the
1006 // GetVecLength() method by default for any traits implementation that
1007 // does not specialize this template. However, for 64-bit support this
1008 // template will need to be specialized.
1009 //
1010 // \note You do <i>not</i> need to write a specialization of
1011 // MultiVecTraitsExt if you are using Epetra, Tpetra, or Thyra
1012 // multivectors. Anasazi already provides specializations for
1013 // these types. Just relax and enjoy using the solvers!
1014 template<class ScalarType>
1015 class MultiVecTraitsExt<ScalarType, fm::eigen::FM_MultiVector<ScalarType> > {
1016 public:
1017  // @name New attribute methods
1019 
1020  // Obtain the vector length of \c mv.
1021  // \note This method supersedes GetVecLength, which will be deprecated.
1022  static ptrdiff_t GetGlobalLength( const fm::eigen::FM_MultiVector<ScalarType>& mv ) {
1023  return mv.GetGlobalLength();
1024  }
1025 
1027 };
1028 
1029 #endif
1030 
1031 }
1032 
1033 #endif