Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
cusp_sa_blockcg.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 IndexType, typename ValueType, typename MemorySpace,
57 typename Orientation>
58void cusp_sa_block_cg(IndexType nx, IndexType ny, IndexType nz, IndexType nrhs,
59 IndexType max_its, ValueType tol, bool verbose) {
60 // create an empty sparse matrix structure
61 cusp::csr_matrix<IndexType, ValueType, MemorySpace> A;
62
63 // create 2D Poisson problem
64 cusp::gallery::poisson27pt(A, nx, ny, nz);
65 std::cout << "N =" << A.num_rows<< std::endl;
66 std::cout << "nnz of A = " << A.num_entries << std::endl;
67
68 // solve with multiple RHS
69 std::cout << "\nSolving Ax = b with multiple RHS..." << std::endl;
70
71 // allocate storage for solution (x) and right hand side (b)
72 cusp::array2d<ValueType, MemorySpace, Orientation> x(A.num_rows, nrhs, 0);
73 cusp::array2d<ValueType, MemorySpace, Orientation> b(A.num_rows, nrhs, 1);
74
75 std::cout << "numRHS = " << nrhs << std::endl;
76
77 // set stopping criteria
78 cusp::default_block_monitor<ValueType> monitor(b, max_its, tol, 0);
79
80 // setup preconditioner
83 if (verbose) {
84 // print hierarchy information
85 std::cout << "\nPreconditioner statistics" << std::endl;
86 M.print();
87 }
88
89 // solve
90 {
91 TEUCHOS_FUNC_TIME_MONITOR("Total Block-CG Solve Time");
92 cusp::krylov::blockcg(A, x, b, monitor, M);
93 }
94}
95
96// Orientation types
97enum Orient { ROW, COL };
98const int num_orient = 2;
99const Orient orient_values[] = { ROW, COL };
100const char *orient_names[] = { "row", "col" };
101
102int main(int argc, char *argv[])
103{
104 typedef int IndexType;
105 typedef double ValueType;
106 typedef cusp::device_memory MemorySpace;
107 //typedef cusp::row_major Orientation;
108
109 bool success = true;
110 bool verbose = false;
111 try {
112
113 // Setup command line options
114 Teuchos::CommandLineProcessor CLP;
115 CLP.setDocString(
116 "This example runs AMG-preconditioned block-CG with CUSP.\n");
117
118 IndexType nx = 32;
119 CLP.setOption("nx", &nx, "Number of mesh points in the x-direction");
120 IndexType ny = 32;
121 CLP.setOption("ny", &ny, "Number of mesh points in the y-direction");
122 IndexType nz = 32;
123 CLP.setOption("nz", &nz, "Number of mesh points in the z-direction");
124 IndexType nrhs = 32;
125 CLP.setOption("nrhs", &nrhs, "Number of right-hand-sides");
126 Orient orient = ROW;
127 CLP.setOption("orient", &orient, num_orient, orient_values, orient_names,
128 "Orientation of block RHS");
129 IndexType max_its = 100;
130 CLP.setOption("max_iterations", &max_its,
131 "Maximum number of CG iterations");
132 double tol = 1e-6; // has to be double
133 CLP.setOption("tolerance", &tol, "Convergence tolerance");
134 int device_id = 0;
135 CLP.setOption("device", &device_id, "CUDA device ID");
136 CLP.setOption("verbose", "quiet", &verbose, "Verbose output");
137 CLP.parse( argc, argv );
138
139 // Set CUDA device
140 cudaSetDevice(device_id);
141 cudaDeviceSetSharedMemConfig(cudaSharedMemBankSizeEightByte);
142
143 if (orient == ROW)
144 cusp_sa_block_cg<IndexType,ValueType,MemorySpace,cusp::row_major>(
145 nx, ny, nz, nrhs, max_its, tol, verbose);
146 else if (orient == COL)
147 cusp_sa_block_cg<IndexType,ValueType,MemorySpace,cusp::column_major>(
148 nx, ny, nz, nrhs, max_its, tol, verbose);
149
150 Teuchos::TimeMonitor::summarize(std::cout);
151 Teuchos::TimeMonitor::zeroOutTimers();
152 }
153 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
154
155 if (success)
156 return 0;
157 return -1;
158}
int main(int argc, char *argv[])
void cusp_sa_block_cg(IndexType nx, IndexType ny, IndexType nz, IndexType nrhs, IndexType max_its, ValueType tol, bool verbose)
const Orient orient_values[]
const int num_orient
const char * orient_names[]
@ COL
@ ROW
void blockcg(LinearOperator &A, Vector &x, Vector &b)