ViennaCL - The Vienna Computing Library  1.7.1
Free open-source GPU-accelerated linear algebra and solver library.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iterative-armadillo.cpp
Go to the documentation of this file.
1 /* =========================================================================
2  Copyright (c) 2010-2016, Institute for Microelectronics,
3  Institute for Analysis and Scientific Computing,
4  TU Wien.
5  Portions of this software are copyright by UChicago Argonne, LLC.
6 
7  -----------------
8  ViennaCL - The Vienna Computing Library
9  -----------------
10 
11  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
12 
13  (A list of authors and contributors can be found in the PDF manual)
14 
15  License: MIT (X11), see file LICENSE in the base directory
16 ============================================================================= */
17 
27 // System headers
28 #include <iostream>
29 
30 #ifndef NDEBUG
31  #define NDEBUG
32 #endif
33 
34 // Armadillo headers (disable BLAS and LAPACK to avoid linking issues)
35 #define ARMA_DONT_USE_BLAS
36 #define ARMA_DONT_USE_LAPACK
37 #include <armadillo>
38 
39 // IMPORTANT: Must be set prior to any ViennaCL includes if you want to use ViennaCL algorithms on Armadillo objects
40 #define VIENNACL_WITH_ARMADILLO 1
41 
42 // ViennaCL headers
43 #include "viennacl/linalg/cg.hpp"
47 
48 
49 // Some helper functions for this tutorial:
50 #include "vector-io.hpp"
51 
56 int main(int, char *[])
57 {
58  typedef float ScalarType;
59 
64  std::vector<std::map<unsigned int, ScalarType> > stl_matrix;
65  std::cout << "Reading matrix (this might take some time)..." << std::endl;
66  if (!viennacl::io::read_matrix_market_file(stl_matrix, "../examples/testdata/mat65k.mtx"))
67  {
68  std::cout << "Error reading Matrix file. Make sure you run from the build/-folder." << std::endl;
69  return EXIT_FAILURE;
70  }
71 
72  // Copy over to Armadillo sparse matrix by putting the indices into a matrix and the values into a vector:
73  std::size_t num_nnz = 0;
74  for (std::size_t i=0; i<stl_matrix.size(); ++i)
75  num_nnz += stl_matrix[i].size();
76 
77  arma::Mat<arma::uword> arma_indices(2, num_nnz);
78  arma::Col<ScalarType> arma_values(num_nnz);
79 
80  std::size_t index = 0;
81  for (std::size_t i=0; i<stl_matrix.size(); ++i)
82  {
83  for (std::map<unsigned int, ScalarType>::const_iterator it = stl_matrix[i].begin(); it != stl_matrix[i].end(); ++it)
84  {
85  arma_indices(0, index) = i;
86  arma_indices(1, index) = it->first;
87  arma_values(index) = it->second;
88  ++index;
89  }
90  }
91  std::cout << "Done: reading matrix" << std::endl;
92 
93 
94 
98  arma::SpMat<ScalarType> arma_matrix(arma_indices, arma_values, 65025, 65025);
99  arma::Col<ScalarType> arma_rhs;
100  arma::Col<ScalarType> arma_result;
101  arma::Col<ScalarType> residual;
102 
106  if (!readVectorFromFile("../examples/testdata/rhs65025.txt", arma_rhs))
107  {
108  std::cout << "Error reading RHS file" << std::endl;
109  return EXIT_FAILURE;
110  }
111 
112  if (!readVectorFromFile("../examples/testdata/result65025.txt", arma_result))
113  {
114  std::cout << "Error reading Result file" << std::endl;
115  return EXIT_FAILURE;
116  }
117 
121  std::cout << "----- Running CG -----" << std::endl;
122  arma_result = viennacl::linalg::solve(arma_matrix, arma_rhs, viennacl::linalg::cg_tag());
123 
124  residual = arma_matrix * arma_result - arma_rhs;
125  std::cout << "Relative residual: " << viennacl::linalg::norm_2(residual) / viennacl::linalg::norm_2(arma_rhs) << std::endl;
126 
130  std::cout << "----- Running BiCGStab -----" << std::endl;
131  arma_result = viennacl::linalg::solve(arma_matrix, arma_rhs, viennacl::linalg::bicgstab_tag());
132 
133  residual = arma_matrix * arma_result - arma_rhs;
134  std::cout << "Relative residual: " << viennacl::linalg::norm_2(residual) / viennacl::linalg::norm_2(arma_rhs) << std::endl;
135 
139  std::cout << "----- Running GMRES -----" << std::endl;
140  arma_result = viennacl::linalg::solve(arma_matrix, arma_rhs, viennacl::linalg::gmres_tag());
141 
142  residual = arma_matrix * arma_result - arma_rhs;
143  std::cout << "Relative residual: " << viennacl::linalg::norm_2(residual) / viennacl::linalg::norm_2(arma_rhs) << std::endl;
144 
148  std::cout << std::endl;
149  std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;
150  std::cout << std::endl;
151 }
152 
T norm_2(std::vector< T, A > const &v1)
Definition: norm_2.hpp:96
A reader and writer for the matrix market format is implemented here.
int main()
Definition: bisect.cpp:91
The stabilized bi-conjugate gradient method is implemented here.
A tag for the solver GMRES. Used for supplying solver parameters and for dispatching the solve() func...
Definition: gmres.hpp:49
VectorT solve(MatrixT const &matrix, VectorT const &rhs, bicgstab_tag const &tag, PreconditionerT const &precond)
Definition: bicgstab.hpp:496
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:239
Implementations of the generalized minimum residual method are in this file.
bool readVectorFromFile(const std::string &filename, VectorType &vec)
Definition: vector-io.hpp:104
The conjugate gradient method is implemented here.
float ScalarType
Definition: fft_1d.cpp:42
A tag for the conjugate gradient Used for supplying solver parameters and for dispatching the solve()...
Definition: cg.hpp:48
A tag for the stabilized Bi-conjugate gradient solver. Used for supplying solver parameters and for d...
Definition: bicgstab.hpp:47
long read_matrix_market_file(MatrixT &mat, const char *file, long index_base=1)
Reads a sparse matrix from a file (MatrixMarket format)