Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
56template <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
80int 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 cusp_sa_block_cg(const cusp::csr_matrix< IndexType, ValueType, MemorySpace > &A, IndexType nrhs, IndexType max_its, ValueType tol)
Definition: TestAMG.cpp:58
int main(int argc, char *argv[])
Definition: TestAMG.cpp:80
void blockcg(LinearOperator &A, Vector &x, Vector &b)