Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
TestAMG.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Stokhos Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <iostream>
43 
44 // CUSP
46 #include <cusp/gallery/poisson.h>
47 #include <cusp/csr_matrix.h>
48 #include <cusp/krylov/blockcg.h>
49 #include <cusp/block_monitor.h>
50 
51 // Utilities
52 #include "Teuchos_CommandLineProcessor.hpp"
53 #include "Teuchos_TimeMonitor.hpp"
54 #include "Teuchos_StandardCatchMacros.hpp"
55 
56 template <typename Orientation,
57  typename IndexType, typename ValueType, typename MemorySpace>
59  const cusp::csr_matrix<IndexType, ValueType, MemorySpace>& A,
60  IndexType nrhs, IndexType max_its, ValueType tol) {
61 
62  // allocate storage for solution (x) and right hand side (b)
63  cusp::array2d<ValueType, MemorySpace, Orientation> x(A.num_rows, nrhs, 0);
64  cusp::array2d<ValueType, MemorySpace, Orientation> b(A.num_rows, nrhs, 1);
65 
66  // set stopping criteria
67  cusp::default_block_monitor<ValueType> monitor(b, max_its, tol, 0.0, false);
68 
69  // setup preconditioner
72 
73  // solve
74  {
75  TEUCHOS_FUNC_TIME_MONITOR("Total Block-CG Solve Time");
76  cusp::krylov::blockcg(A, x, b, monitor, M);
77  }
78 }
79 
80 int main(int argc, char *argv[])
81 {
82  typedef int IndexType;
83  typedef double ValueType;
84  typedef cusp::device_memory MemorySpace;
85  //typedef cusp::row_major Orientation;
86 
87  bool success = true;
88  bool verbose = false;
89  try {
90 
91  // Setup command line options
92  Teuchos::CommandLineProcessor CLP;
93  CLP.setDocString("This test performance of block multiply routines.\n");
94  IndexType n = 32;
95  CLP.setOption("n", &n, "Number of mesh points in the each direction");
96  IndexType nrhs_begin = 32;
97  CLP.setOption("begin", &nrhs_begin,
98  "Staring number of right-hand-sides");
99  IndexType nrhs_end = 512;
100  CLP.setOption("end", &nrhs_end,
101  "Ending number of right-hand-sides");
102  IndexType nrhs_step = 32;
103  CLP.setOption("step", &nrhs_step,
104  "Increment in number of right-hand-sides");
105  IndexType max_its = 100;
106  CLP.setOption("max_iterations", &max_its,
107  "Maximum number of CG iterations");
108  double tol = 1e-6; // has to be double
109  CLP.setOption("tolerance", &tol, "Convergence tolerance");
110  int device_id = 0;
111  CLP.setOption("device", &device_id, "CUDA device ID");
112  CLP.parse( argc, argv );
113 
114  // Set CUDA device
115  cudaSetDevice(device_id);
116  cudaDeviceSetSharedMemConfig(cudaSharedMemBankSizeEightByte);
117 
118  // create 3D Poisson problem
119  cusp::csr_matrix<IndexType, ValueType, MemorySpace> A;
120  cusp::gallery::poisson27pt(A, n, n, n);
121 
122  // Create timers
123  Teuchos::RCP<Teuchos::Time> tm_cg =
124  Teuchos::TimeMonitor::getNewTimer("Total Block-CG Solve Time");
125  Teuchos::RCP<Teuchos::Time> tm_prec =
126  Teuchos::TimeMonitor::getNewTimer("CUSP Block Multilevel Solve");
127  Teuchos::RCP<Teuchos::Time> tm_coarse =
128  Teuchos::TimeMonitor::getNewTimer("CUSP Coarse-grid Solve");
129  Teuchos::RCP<Teuchos::Time> tm_op =
130  Teuchos::TimeMonitor::getNewTimer("CUSP Operator block-apply");
131  Teuchos::RCP<Teuchos::Time> tm_prec_op =
132  Teuchos::TimeMonitor::getNewTimer("CUSP Matrix block-apply");
133 
134  std::cout << "nrhs , num_rows , num_entries , "
135  << "row_cg , row_op , row_prec , row_prec_op , row_coarse , "
136  << "col_cg , col_op , col_prec , col_prec_op , col_coarse"
137  << std::endl;
138 
139  for (IndexType nrhs = nrhs_begin; nrhs <= nrhs_end; nrhs += nrhs_step) {
140 
141  std::cout << nrhs << " , "
142  << A.num_rows << " , " << A.num_entries << " , ";
143 
144  // test row-major storage
145  Teuchos::TimeMonitor::zeroOutTimers();
146  cusp_sa_block_cg<cusp::row_major>(A, nrhs, max_its, tol);
147 
148  std::cout << tm_cg->totalElapsedTime() << " , "
149  << tm_op->totalElapsedTime() << " , "
150  << tm_prec->totalElapsedTime() << " , "
151  << tm_prec_op->totalElapsedTime() << " , "
152  << tm_coarse->totalElapsedTime() << " , ";
153 
154  // test column-major storage
155  Teuchos::TimeMonitor::zeroOutTimers();
156  cusp_sa_block_cg<cusp::column_major>(A, nrhs, max_its, tol);
157 
158  std::cout << tm_cg->totalElapsedTime() << " , "
159  << tm_op->totalElapsedTime() << " , "
160  << tm_prec->totalElapsedTime() << " , "
161  << tm_prec_op->totalElapsedTime() << " , "
162  << tm_coarse->totalElapsedTime() << std::endl;
163 
164  }
165 
166  }
167  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
168 
169  if (success)
170  return 0;
171  return -1;
172 }
void blockcg(LinearOperator &A, Vector &x, Vector &b)
int main(int argc, char *argv[])
Definition: TestAMG.cpp:80
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
void cusp_sa_block_cg(const cusp::csr_matrix< IndexType, ValueType, MemorySpace > &A, IndexType nrhs, IndexType max_its, ValueType tol)
Definition: TestAMG.cpp:58