/[escript]/trunk/cusplibrary/cusp/gallery/poisson.h
ViewVC logotype

Contents of /trunk/cusplibrary/cusp/gallery/poisson.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5148 - (show annotations)
Mon Sep 15 01:25:23 2014 UTC (5 years, 9 months ago) by caltinay
File MIME type: text/plain
File size: 5357 byte(s)
Merging ripley diagonal storage + CUDA support into trunk.
Options file version has been incremented due to new options
'cuda' and 'nvccflags'.

1 /*
2 * Copyright 2008-2009 NVIDIA Corporation
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17 /*! \file poisson.h
18 * \brief Poisson matrix generators
19 */
20
21 #pragma once
22
23 #include <cusp/detail/config.h>
24
25 #include <cusp/gallery/stencil.h>
26
27 namespace cusp
28 {
29 namespace gallery
30 {
31 /*! \addtogroup gallery Matrix Gallery
32 * \addtogroup poisson Poisson
33 * \ingroup gallery
34 * \{
35 */
36
37 /*! \p poisson5pt: Create a matrix representing a Poisson problem
38 * discretized on an \p m by \p n grid with the standard 5-point
39 * finite-difference stencil.
40 *
41 * \param matrix output
42 * \param m number of grid rows
43 * \param n number of grid columns
44 * \tparam MatrixType matrix container
45 *
46 * \code
47 * #include <cusp/gallery/poisson.h>
48 * #include <cusp/coo_matrix.h>
49 * #include <cusp/print.h>
50 *
51 * int main(void)
52 * {
53 * cusp::coo_matrix<int, float, cusp::device_memory> A;
54 *
55 * // create a matrix for a Poisson problem on a 4x4 grid
56 * cusp::gallery::poisson5pt(A, 4, 4);
57 *
58 * // print matrix
59 * cusp::print(A);
60 *
61 * return 0;
62 * }
63 * \endcode
64 *
65 */
66 template <typename MatrixType>
67 void poisson5pt( MatrixType& matrix, size_t m, size_t n)
68 {
69 CUSP_PROFILE_SCOPED();
70
71 typedef typename MatrixType::index_type IndexType;
72 typedef typename MatrixType::value_type ValueType;
73 typedef thrust::tuple<IndexType,IndexType> StencilIndex;
74 typedef thrust::tuple<StencilIndex,ValueType> StencilPoint;
75
76 cusp::array1d<StencilPoint, cusp::host_memory> stencil;
77 stencil.push_back(StencilPoint(StencilIndex( 0, -1), -1));
78 stencil.push_back(StencilPoint(StencilIndex( -1, 0), -1));
79 stencil.push_back(StencilPoint(StencilIndex( 0, 0), 4));
80 stencil.push_back(StencilPoint(StencilIndex( 1, 0), -1));
81 stencil.push_back(StencilPoint(StencilIndex( 0, 1), -1));
82
83 cusp::gallery::generate_matrix_from_stencil(matrix, stencil, StencilIndex(m,n));
84 }
85
86 template <typename MatrixType>
87 void poisson9pt( MatrixType& matrix, size_t m, size_t n)
88 {
89 CUSP_PROFILE_SCOPED();
90
91 typedef typename MatrixType::index_type IndexType;
92 typedef typename MatrixType::value_type ValueType;
93 typedef thrust::tuple<IndexType,IndexType> StencilIndex;
94 typedef thrust::tuple<StencilIndex,ValueType> StencilPoint;
95
96 cusp::array1d<StencilPoint, cusp::host_memory> stencil;
97 stencil.push_back(StencilPoint(StencilIndex( -1, -1), -1));
98 stencil.push_back(StencilPoint(StencilIndex( 0, -1), -1));
99 stencil.push_back(StencilPoint(StencilIndex( 1, -1), -1));
100 stencil.push_back(StencilPoint(StencilIndex( -1, 0), -1));
101 stencil.push_back(StencilPoint(StencilIndex( 0, 0), 8));
102 stencil.push_back(StencilPoint(StencilIndex( 1, 0), -1));
103 stencil.push_back(StencilPoint(StencilIndex( -1, 1), -1));
104 stencil.push_back(StencilPoint(StencilIndex( 0, 1), -1));
105 stencil.push_back(StencilPoint(StencilIndex( 1, 1), -1));
106
107 cusp::gallery::generate_matrix_from_stencil(matrix, stencil, StencilIndex(m,n));
108 }
109
110 template <typename MatrixType>
111 void poisson7pt( MatrixType& matrix, size_t m, size_t n, size_t k)
112 {
113 CUSP_PROFILE_SCOPED();
114
115 typedef typename MatrixType::index_type IndexType;
116 typedef typename MatrixType::value_type ValueType;
117 typedef thrust::tuple<IndexType,IndexType,IndexType> StencilIndex;
118 typedef thrust::tuple<StencilIndex,ValueType> StencilPoint;
119
120 cusp::array1d<StencilPoint, cusp::host_memory> stencil;
121 stencil.push_back(StencilPoint(StencilIndex( 0, 0, -1), -1));
122 stencil.push_back(StencilPoint(StencilIndex( 0, -1, 0), -1));
123 stencil.push_back(StencilPoint(StencilIndex(-1, 0, 0), -1));
124 stencil.push_back(StencilPoint(StencilIndex( 0, 0, 0), 6));
125 stencil.push_back(StencilPoint(StencilIndex( 1, 0, 0), -1));
126 stencil.push_back(StencilPoint(StencilIndex( 0, 1, 0), -1));
127 stencil.push_back(StencilPoint(StencilIndex( 0, 0, 1), -1));
128
129 cusp::gallery::generate_matrix_from_stencil(matrix, stencil, StencilIndex(m,n,k));
130 }
131
132 template <typename MatrixType>
133 void poisson27pt( MatrixType& matrix, size_t m, size_t n, size_t l)
134 {
135 CUSP_PROFILE_SCOPED();
136
137 typedef typename MatrixType::index_type IndexType;
138 typedef typename MatrixType::value_type ValueType;
139 typedef thrust::tuple<IndexType,IndexType,IndexType> StencilIndex;
140 typedef thrust::tuple<StencilIndex,ValueType> StencilPoint;
141
142 cusp::array1d<StencilPoint, cusp::host_memory> stencil;
143 for( IndexType k = -1; k <= 1; k++ )
144 for( IndexType j = -1; j <= 1; j++ )
145 for( IndexType i = -1; i <= 1; i++ )
146 stencil.push_back(StencilPoint(StencilIndex( i, j, k), (i==0 && j==0 && k==0) ? 26 : -1));
147
148 cusp::gallery::generate_matrix_from_stencil(matrix, stencil, StencilIndex(m,n,l));
149 }
150 /*! \}
151 */
152
153 } // end namespace gallery
154 } // end namespace cusp
155

  ViewVC Help
Powered by ViewVC 1.1.26