/[escript]/trunk/paso/src/SystemMatrix_mergeMainAndCouple.c
ViewVC logotype

Contents of /trunk/paso/src/SystemMatrix_mergeMainAndCouple.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3510 - (show annotations)
Fri May 13 06:09:49 2011 UTC (8 years, 4 months ago) by gross
File MIME type: text/plain
File size: 6751 byte(s)
some fixes for the compilation without BOOMERAMG
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /************************************************************/
16 /* Paso: SystemMatrix */
17 /* */
18 /* Merge the MainBlock and CoupleBlock in the matrix */
19 /* Input: SystemMatrix A */
20 /* Output: */
21 /* p_ptr: the pointer to a vector of locations that */
22 /* start a row. */
23 /* p_idx: the pointer to the column indices for each */
24 /* of the rows, and is ordered by rows. */
25 /* p_val: the pointer to the data corresponding */
26 /* directly to the column entries in p_idx. */
27 /************************************************************/
28
29 /* Copyrights by ACcESS Australia 2003 */
30 /* Author: Lin Gao, l.gao@uq.edu.au */
31
32 /************************************************************/
33
34 #include "Paso.h"
35 #include "SystemMatrix.h"
36
37 #ifdef _OPENMP
38 #include <omp.h>
39 #endif
40
41
42 void Paso_SystemMatrix_mergeMainAndCouple(Paso_SystemMatrix* A, index_t** p_ptr, index_t** p_idx, double** p_val){
43 if (A->type & MATRIX_FORMAT_DEFAULT) {
44 Paso_SystemMatrix_mergeMainAndCouple_CSR_OFFSET0(A, p_ptr, p_idx, p_val);
45 } else if (A->type & MATRIX_FORMAT_CSC) {
46 /* CSC part is for PASTIX */
47 if (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) {
48 Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1(A, p_ptr, p_idx, p_val);
49 } else {
50 Esys_setError(SYSTEM_ERROR,"Paso_SystemMatrix_mergeMainAndCouple: CSC with index 0 or block size larger than 1 is not supported. ");
51 }
52 } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
53 Esys_setError(SYSTEM_ERROR,"Paso_SystemMatrix_mergeMainAndCouple: TRILINOS is not supported. ");
54 } else {
55 Esys_setError(SYSTEM_ERROR,"Paso_SystemMatrix_mergeMainAndCouple: CRS is not supported. ");
56 }
57 return;
58 }
59
60
61 void Paso_SystemMatrix_mergeMainAndCouple_CSR_OFFSET0(Paso_SystemMatrix* A, index_t** p_ptr, index_t** p_idx, double** p_val) {
62
63 index_t i, j, i_ub, j_lb, j_ub, row, num_vals, main_num_vals;
64 index_t couple_num_vals, rank, row_offset, ij_ptr=0, idx=0, idx2=0;
65 index_t main_num_rows, couple_num_rows, col_offset;
66 index_t *main_ptr, *main_idx, *couple_ptr, *couple_idx;
67 double *main_val, *couple_val, *rows;
68 Paso_Coupler* coupler;
69
70 if (A->mainBlock->col_block_size!=1 ||
71 A->mainBlock->row_block_size!=1 ||
72 A->col_coupleBlock->col_block_size!=1 ||
73 A->col_coupleBlock->row_block_size!=1) {
74 Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_mergeMainAndCouple_CSR_OFFSET0: requires format with block size 1.");
75 return;
76 }
77
78 if (A->mpi_info->size == 1) {
79 /* initialization */
80 main_num_rows=A->mainBlock->numRows;
81 main_ptr=A->mainBlock->pattern->ptr;
82 main_idx=A->mainBlock->pattern->index;
83 main_val=A->mainBlock->val;
84
85 /* allocate arrays "ptr", "index" and "val" */
86 num_vals = main_ptr[main_num_rows]-1;
87 *p_ptr = MEMALLOC(main_num_rows+1, index_t);
88 *p_idx = MEMALLOC(num_vals, index_t);
89 *p_val = MEMALLOC(num_vals, double);
90
91 #pragma omp for schedule(static) private(i,ij_ptr,j_lb,j_ub)
92 for (i=0; i<main_num_rows; i++) {
93 j_lb = main_ptr[i];
94 j_ub = main_ptr[i+1];
95 (*p_ptr)[i] = j_lb;
96 for (ij_ptr=j_lb; ij_ptr<j_ub; ++ij_ptr) {
97 (*p_idx)[ij_ptr] = main_idx[ij_ptr];
98 (*p_val)[ij_ptr] = main_val[ij_ptr];
99 }
100 }
101 (*p_ptr)[main_num_rows] = main_ptr[main_num_rows];
102 return;
103 }
104
105 main_num_rows=A->mainBlock->numRows;
106 couple_num_rows=A->col_coupleBlock->numRows;
107
108 if (main_num_rows != couple_num_rows) {
109 Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_mergeMainAndCouple_CSR_OFFSET0: number of rows do not match.");
110 return;
111 }
112
113 /* prepare for global coordinates in colCoupleBlock, the results are
114 in coupler->recv_buffer */
115 rows=TMPMEMALLOC(main_num_rows, double);
116 rank = A->mpi_info->rank;
117 row_offset = A->row_distribution->first_component[rank];
118 #pragma omp parallel for private(i) schedule(static)
119 for (i=0; i<main_num_rows; ++i) rows[i]=row_offset+i;
120 coupler= Paso_Coupler_alloc(A->col_coupler->connector, 1);
121 Paso_Coupler_startCollect(coupler, rows);
122
123 /* initalization, including allocate arrays "ptr", "index" and "val" */
124 main_ptr=A->mainBlock->pattern->ptr;
125 main_idx=A->mainBlock->pattern->index;
126 main_val=A->mainBlock->val;
127 couple_ptr=A->col_coupleBlock->pattern->ptr;
128 couple_idx=A->col_coupleBlock->pattern->index;
129 couple_val=A->col_coupleBlock->val;
130 col_offset = A->col_distribution->first_component[rank];
131 main_num_vals = main_ptr[main_num_rows]-main_ptr[0];
132 couple_num_vals = couple_ptr[couple_num_rows]-couple_ptr[0];
133 num_vals = main_num_vals + couple_num_vals;
134 *p_ptr = MEMALLOC(main_num_rows+1, index_t);
135 *p_idx = MEMALLOC(num_vals, index_t);
136 *p_val = MEMALLOC(num_vals, double);
137 (*p_ptr)[0] = 0;
138 i = 0;
139 j = 0;
140
141 Paso_Coupler_finishCollect(coupler);
142
143 /* merge mainBlock and col_coupleBlock */
144 for (row=1; row<=main_num_rows; row++) {
145 i_ub = main_ptr[row];
146 j_ub = couple_ptr[row];
147 while (i < i_ub || j < j_ub) {
148 ij_ptr = i + j;
149 if (j < j_ub) {
150 idx = coupler->recv_buffer[couple_idx[j]];
151 }
152 if (i < i_ub) {
153 idx2 = main_idx[i] + col_offset;
154 }
155 if (j == j_ub || (i < i_ub && idx2 < idx)){
156 (*p_idx)[ij_ptr] = idx2;
157 (*p_val)[ij_ptr] = main_val[i];
158 i++;
159 } else {
160 (*p_idx)[ij_ptr] = idx;
161 (*p_val)[ij_ptr] = couple_val[j];
162 j++;
163 }
164 }
165 (*p_ptr)[row] = ij_ptr+1;
166 }
167
168 TMPMEMFREE(rows);
169 Paso_Coupler_free(coupler);
170 return;
171 }
172
173
174 void Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1(Paso_SystemMatrix* A, index_t** p_ptr, index_t** p_idx, double** p_val) {
175 /*
176 index_t i, j, i_ub, j_ub, col, num_vals, main_num_vals;
177 index_t couple_num_vals, idx, rank, main_offset;
178 index_t main_num_cols=A->mainBlock->pattern->numOutput;
179 index_t couple_num_cols=A->col_coupleBlock->pattern->numOutput;
180 index_t *main_ptr=A->mainBlock->pattern->ptr;
181 index_t *main_idx=A->mainBlock->pattern->index;
182 double *main_val=A->mainBlock->val;
183 index_t *couple_ptr=A->col_coupleBlock->pattern->ptr;
184 index_t *couple_idx=A->col_coupleBlock->pattern->index;
185 double *couple_val=A->col_coupleBlock->val;
186 index_t *couple_global=NULL;
187 Paso_Coupler* coupler=A->col_coupler;
188 */
189
190 return;
191 }
192

  ViewVC Help
Powered by ViewVC 1.1.26