/[escript]/branches/trilinos_from_5897/paso/src/SystemMatrix_mergeMainAndCouple.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/paso/src/SystemMatrix_mergeMainAndCouple.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5933 - (show annotations)
Wed Feb 17 23:53:30 2016 UTC (21 months, 4 weeks ago) by caltinay
File size: 10404 byte(s)
sync with trunk.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17
18 /****************************************************************************
19 * Paso: SystemMatrix
20 *
21 * Merge the MainBlock and CoupleBlock in the matrix
22 * Input: SystemMatrix A
23 * Output:
24 * p_ptr: the pointer to a vector of locations that start a row.
25 * p_idx: the pointer to the column indices for each of the rows,
26 * ordered by rows.
27 * p_val: the pointer to the data corresponding directly to the
28 * column entries in p_idx.
29 ****************************************************************************/
30
31 /* Copyrights by ACcESS Australia 2003 */
32 /* Author: Lin Gao, l.gao@uq.edu.au */
33
34 /****************************************************************************/
35
36 #include "Paso.h"
37 #include "SystemMatrix.h"
38
39 namespace paso {
40
41 void SystemMatrix::mergeMainAndCouple(index_t** p_ptr, index_t** p_idx, double** p_val) const
42 {
43 if (type & MATRIX_FORMAT_DEFAULT) {
44 mergeMainAndCouple_CSR_OFFSET0(p_ptr, p_idx, p_val);
45 } else if (type & MATRIX_FORMAT_CSC) {
46 /* CSC part is for PASTIX */
47 if (type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) {
48 mergeMainAndCouple_CSC_OFFSET1(p_ptr, p_idx, p_val);
49 } else {
50 Esys_setError(SYSTEM_ERROR, "SystemMatrix::mergeMainAndCouple: CSC with index 0 or block size larger than 1 is not supported.");
51 }
52 } else {
53 Esys_setError(SYSTEM_ERROR,"SystemMatrix::mergeMainAndCouple: CRS is not supported.");
54 }
55 }
56
57 void SystemMatrix::mergeMainAndCouple_CSR_OFFSET0(index_t** p_ptr, index_t** p_idx, double** p_val) const
58 {
59 if (mainBlock->col_block_size != 1 || mainBlock->row_block_size != 1 ||
60 col_coupleBlock->col_block_size != 1 ||
61 col_coupleBlock->row_block_size != 1) {
62 mergeMainAndCouple_CSR_OFFSET0_Block(p_ptr, p_idx, p_val);
63 return;
64 }
65
66 const index_t main_num_rows = mainBlock->numRows;
67 index_t* main_ptr = mainBlock->pattern->ptr;
68 index_t* main_idx = mainBlock->pattern->index;
69 double* main_val = mainBlock->val;
70
71 if (mpi_info->size == 1) {
72 // allocate arrays "ptr", "index" and "val"
73 const index_t num_vals = main_ptr[main_num_rows]-1;
74 *p_ptr = new index_t[main_num_rows+1];
75 *p_idx = new index_t[num_vals];
76 *p_val = new double[num_vals];
77
78 #pragma omp parallel for
79 for (index_t i=0; i<main_num_rows; i++) {
80 const index_t j_lb = main_ptr[i];
81 const index_t j_ub = main_ptr[i+1];
82 (*p_ptr)[i] = j_lb;
83 for (index_t ij_ptr=j_lb; ij_ptr<j_ub; ++ij_ptr) {
84 (*p_idx)[ij_ptr] = main_idx[ij_ptr];
85 (*p_val)[ij_ptr] = main_val[ij_ptr];
86 }
87 }
88 (*p_ptr)[main_num_rows] = main_ptr[main_num_rows];
89 return;
90 }
91
92 const index_t couple_num_rows = col_coupleBlock->numRows;
93
94 if (main_num_rows != couple_num_rows) {
95 Esys_setError(TYPE_ERROR, "SystemMatrix::mergeMainAndCouple_CSR_OFFSET0: number of rows do not match.");
96 return;
97 }
98
99 double* rows = NULL;
100 const index_t rank = mpi_info->rank;
101 Coupler_ptr coupler;
102 if (global_id == NULL) {
103 // prepare for global coordinates in colCoupleBlock, the results are
104 // in coupler->recv_buffer
105 rows = new double[main_num_rows];
106 const index_t row_offset = row_distribution->first_component[rank];
107 #pragma omp parallel for
108 for (index_t i=0; i<main_num_rows; ++i)
109 rows[i] = row_offset+i;
110 coupler.reset(new Coupler(col_coupler->connector, 1));
111 coupler->startCollect(rows);
112 }
113
114 // initialisation, including allocate arrays "ptr", "index" and "val"
115 index_t* couple_ptr = col_coupleBlock->pattern->ptr;
116 index_t* couple_idx = col_coupleBlock->pattern->index;
117 double* couple_val = col_coupleBlock->val;
118 const index_t col_offset = col_distribution->first_component[rank];
119 const index_t main_num_vals = main_ptr[main_num_rows]-main_ptr[0];
120 const index_t couple_num_vals = couple_ptr[couple_num_rows]-couple_ptr[0];
121 const index_t num_vals = main_num_vals + couple_num_vals;
122 *p_ptr = new index_t[main_num_rows+1];
123 *p_idx = new index_t[num_vals];
124 *p_val = new double[num_vals];
125 (*p_ptr)[0] = 0;
126
127 if (global_id == NULL) {
128 coupler->finishCollect();
129 delete[] rows;
130 const index_t num_cols = col_coupleBlock->numCols;
131 global_id = new index_t[num_cols];
132 #pragma omp parallel for
133 for (index_t i=0; i<num_cols; ++i)
134 global_id[i] = coupler->recv_buffer[i];
135 coupler.reset();
136 }
137
138 index_t i = 0;
139 index_t j = 0;
140 index_t idx = 0;
141 index_t idx2 = 0;
142 index_t ij_ptr = 0;
143
144 // merge mainBlock and col_coupleBlock
145 for (index_t row=1; row<=main_num_rows; row++) {
146 const index_t i_ub = main_ptr[row];
147 const index_t j_ub = couple_ptr[row];
148 while (i < i_ub || j < j_ub) {
149 ij_ptr = i + j;
150 if (j < j_ub) {
151 idx = global_id[couple_idx[j]];
152 }
153 if (i < i_ub) {
154 idx2 = main_idx[i] + col_offset;
155 }
156 if (j == j_ub || (i < i_ub && idx2 < idx)){
157 (*p_idx)[ij_ptr] = idx2;
158 (*p_val)[ij_ptr] = main_val[i];
159 i++;
160 } else {
161 (*p_idx)[ij_ptr] = idx;
162 (*p_val)[ij_ptr] = couple_val[j];
163 j++;
164 }
165 }
166 (*p_ptr)[row] = ij_ptr+1;
167 }
168 }
169
170
171 void SystemMatrix::mergeMainAndCouple_CSR_OFFSET0_Block(index_t** p_ptr, index_t** p_idx, double** p_val) const
172 {
173 const index_t main_num_rows = mainBlock->numRows;
174 index_t* main_ptr = mainBlock->pattern->ptr;
175 index_t* main_idx = mainBlock->pattern->index;
176 double* main_val = mainBlock->val;
177
178 if (mpi_info->size == 1) {
179 // allocate arrays "ptr", "index" and "val"
180 const index_t num_vals = main_ptr[main_num_rows]-1;
181 *p_ptr = new index_t[main_num_rows+1];
182 *p_idx = new index_t[num_vals];
183 *p_val = new double[num_vals * block_size];
184
185 #pragma omp parallel for
186 for (index_t i=0; i<main_num_rows; i++) {
187 const index_t j_lb = main_ptr[i];
188 const index_t j_ub = main_ptr[i+1];
189 (*p_ptr)[i] = j_lb;
190 for (index_t ij_ptr=j_lb; ij_ptr<j_ub; ij_ptr++) {
191 (*p_idx)[ij_ptr] = main_idx[ij_ptr];
192 for (index_t ib=0; ib<block_size; ib++) {
193 (*p_val)[ij_ptr*block_size+ib] = main_val[ij_ptr*block_size+ib];
194 }
195 }
196 }
197 (*p_ptr)[main_num_rows] = main_ptr[main_num_rows];
198 return;
199 }
200
201 const index_t couple_num_rows = col_coupleBlock->numRows;
202
203 if (main_num_rows != couple_num_rows) {
204 Esys_setError(TYPE_ERROR, "SystemMatrix_mergeMainAndCouple_CSR_OFFSET0: number of rows do not match.");
205 return;
206 }
207
208 double* rows = NULL;
209 const index_t rank = mpi_info->rank;
210 Coupler_ptr coupler;
211 if (global_id == NULL) {
212 // prepare for global coordinates in colCoupleBlock, the results are
213 // in coupler->recv_buffer
214 rows = new double[main_num_rows];
215 const index_t row_offset = row_distribution->first_component[rank];
216 #pragma omp parallel for
217 for (index_t i=0; i<main_num_rows; ++i)
218 rows[i]=row_offset+i;
219 coupler.reset(new Coupler(col_coupler->connector, 1));
220 coupler->startCollect(rows);
221 }
222
223 // initialisation, including allocate arrays "ptr", "index" and "val"
224 index_t* couple_ptr = col_coupleBlock->pattern->ptr;
225 index_t* couple_idx = col_coupleBlock->pattern->index;
226 double* couple_val = col_coupleBlock->val;
227 const index_t col_offset = col_distribution->first_component[rank];
228 const index_t main_num_vals = main_ptr[main_num_rows]-main_ptr[0];
229 const index_t couple_num_vals = couple_ptr[couple_num_rows]-couple_ptr[0];
230 const index_t num_vals = main_num_vals + couple_num_vals;
231 *p_ptr = new index_t[main_num_rows+1];
232 *p_idx = new index_t[num_vals];
233 *p_val = new double[num_vals*block_size];
234 (*p_ptr)[0] = 0;
235
236 if (global_id == NULL) {
237 coupler->finishCollect();
238 delete[] rows;
239 const index_t num_cols = col_coupleBlock->numCols;
240 global_id = new index_t[num_cols];
241 #pragma omp parallel for
242 for (index_t i=0; i<num_cols; ++i)
243 global_id[i] = coupler->recv_buffer[i];
244 coupler.reset();
245 }
246
247 index_t i = 0;
248 index_t j = 0;
249 index_t idx = 0;
250 index_t idx2 = 0;
251 index_t ij_ptr = 0;
252
253 // merge mainBlock and col_coupleBlock
254 for (index_t row=1; row<=main_num_rows; row++) {
255 const index_t i_ub = main_ptr[row];
256 const index_t j_ub = couple_ptr[row];
257 while (i < i_ub || j < j_ub) {
258 ij_ptr = i + j;
259 if (j < j_ub) {
260 idx = global_id[couple_idx[j]];
261 }
262 if (i < i_ub) {
263 idx2 = main_idx[i] + col_offset;
264 }
265 if (j == j_ub || (i < i_ub && idx2 < idx)){
266 (*p_idx)[ij_ptr] = idx2;
267 for (index_t ib=0; ib<block_size; ib++)
268 (*p_val)[ij_ptr*block_size+ib] = main_val[i*block_size+ib];
269 i++;
270 } else {
271 (*p_idx)[ij_ptr] = idx;
272 for (index_t ib=0; ib<block_size; ib++)
273 (*p_val)[ij_ptr*block_size+ib] = couple_val[j*block_size+ib];
274 j++;
275 }
276 }
277 (*p_ptr)[row] = ij_ptr+1;
278 }
279 }
280
281 void SystemMatrix::mergeMainAndCouple_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val) const
282 {
283 Esys_setError(TYPE_ERROR, "SystemMatrix_mergeMainAndCouple_CSC_OFFSET1: not implemented.");
284 }
285
286 } // namespace paso
287

Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/paso/src/SystemMatrix_mergeMainAndCouple.cpp:5567-5588 /branches/amg_from_3530/paso/src/SystemMatrix_mergeMainAndCouple.cpp:3531-3826 /branches/lapack2681/paso/src/SystemMatrix_mergeMainAndCouple.cpp:2682-2741 /branches/pasowrap/paso/src/SystemMatrix_mergeMainAndCouple.cpp:3661-3674 /branches/py3_attempt2/paso/src/SystemMatrix_mergeMainAndCouple.cpp:3871-3891 /branches/restext/paso/src/SystemMatrix_mergeMainAndCouple.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/SystemMatrix_mergeMainAndCouple.cpp:3669-3791 /branches/stage3.0/paso/src/SystemMatrix_mergeMainAndCouple.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/SystemMatrix_mergeMainAndCouple.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/SystemMatrix_mergeMainAndCouple.cpp:3517-3974 /release/3.0/paso/src/SystemMatrix_mergeMainAndCouple.cpp:2591-2601 /release/4.0/paso/src/SystemMatrix_mergeMainAndCouple.cpp:5380-5406 /trunk/paso/src/SystemMatrix_mergeMainAndCouple.cpp:4257-4344,5898-5932 /trunk/ripley/test/python/paso/src/SystemMatrix_mergeMainAndCouple.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26