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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2131 - (show annotations)
Thu Dec 4 06:09:50 2008 UTC (10 years, 9 months ago) by artak
File MIME type: text/plain
File size: 7982 byte(s)
Openmp bug is fixed in Paso_cum_sum and some printouts removed from Solver_AMG
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 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
17 /* Paso: Pattern: Paso_Pattern_coupling
18
19 searches for a maximal independent set MIS in the matrix pattern
20 vertices in the maximal independent set are marked in mis_marker
21 nodes to be considered are marked by -1 on the input in mis_marker
22
23 */
24 /**********************************************************************/
25
26 /* Copyrights by ACcESS Australia 2003,2004,2005 */
27 /* Author: artak@uq.edu.au */
28
29 /**************************************************************/
30
31 #include "PasoUtil.h"
32 #include "Pattern_coupling.h"
33
34
35 /***************************************************************/
36
37 #define IS_AVAILABLE -1
38 #define IS_IN_SET -3
39 #define IS_REMOVED -4
40
41 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
42
43 dim_t i,j;
44 /*double threshold=0.05;*/
45 double sum;
46 index_t iptr,*index,*where_p,*diagptr;
47 bool_t passed=FALSE;
48 dim_t n=A->numRows;
49 diagptr=MEMALLOC(n,index_t);
50
51 if (A->pattern->type & PATTERN_FORMAT_SYM) {
52 Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");
53 return;
54 }
55
56 #pragma omp parallel for private(i) schedule(static)
57 for (i=0;i<n;++i)
58 if(mis_marker[i]==IS_AVAILABLE)
59 mis_marker[i]=IS_IN_SET;
60
61 #pragma omp parallel for private(i,index,where_p) schedule(static)
62 for (i=0;i<n;++i) {
63 diagptr[i]=A->pattern->ptr[i];
64 index=&(A->pattern->index[diagptr[i]]);
65 where_p=(index_t*)bsearch(&i,
66 index,
67 A->pattern->ptr[i + 1]-A->pattern->ptr[i],
68 sizeof(index_t),
69 Paso_comparIndex);
70 if (where_p==NULL) {
71 Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
72 } else {
73 diagptr[i]+=(index_t)(where_p-index);
74 }
75 }
76
77
78
79 /*This loop cannot be parallelized, as order matters here.*/
80 for (i=0;i<n;++i) {
81 if (mis_marker[i]==IS_IN_SET) {
82 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
83 j=A->pattern->index[iptr];
84 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
85 mis_marker[j]=IS_REMOVED;
86 }
87 }
88 }
89 }
90
91
92
93 /*This loop cannot be parallelized, as order matters here.*/
94 for (i=0;i<n;i++) {
95 if (mis_marker[i]==IS_REMOVED) {
96 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
97 j=A->pattern->index[iptr];
98 if (mis_marker[j]==IS_IN_SET) {
99 if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
100 passed=TRUE;
101 }
102 else {
103 passed=FALSE;
104 break;
105 }
106 }
107 }
108 if (passed) {
109 mis_marker[i]=IS_IN_SET;
110 passed=FALSE;
111 }
112 }
113 }
114
115 /* This check is to make sure we dont get some nusty rows which were not removed durring coarsing process.*/
116 /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
117 /*This loop cannot be parallelized, as order matters here.*/
118 for (i=0;i<n;i++) {
119 if (mis_marker[i]==IS_REMOVED) {
120 sum=0;
121 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
122 j=A->pattern->index[iptr];
123 if (mis_marker[j]==IS_REMOVED)
124 sum+=A->val[iptr];
125 }
126 if(ABS(sum)<1.e-12)
127 mis_marker[i]=IS_IN_SET;
128 }
129 }
130
131
132 /* swap to TRUE/FALSE in mis_marker */
133 #pragma omp parallel for private(i) schedule(static)
134 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
135
136 MEMFREE(diagptr);
137 }
138
139 /*
140 *
141 * Return a strength of connection mask using the classical
142 * strength of connection measure by Ruge and Stuben.
143 *
144 * Specifically, an off-diagonal entry A[i.j] is a strong
145 * connection if:
146 *
147 * -A[i,j] >= theta * max( -A[i,k] ) where k != i
148 *
149 * Otherwise, the connection is weak.
150 *
151 */
152 void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
153 {
154 dim_t i;
155 index_t iptr;
156 double threshold,max_offdiagonal;
157 dim_t n=A->numRows;
158 if (A->pattern->type & PATTERN_FORMAT_SYM) {
159 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
160 return;
161 }
162
163 #pragma omp parallel for private(i,iptr,max_offdiagonal,threshold) schedule(static)
164 for (i=0;i<n;++i) {
165 if(mis_marker[i]==IS_AVAILABLE) {
166 /*min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];*/
167 max_offdiagonal = -1.e+15;
168 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
169 if(A->pattern->index[iptr] != i){
170 max_offdiagonal = MAX(max_offdiagonal,-(A->val[iptr]));
171 }
172 }
173
174 threshold = theta*max_offdiagonal;
175 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
176 if(-(A->val[iptr])<threshold && A->pattern->index[iptr]!=i){
177 A->val[iptr]=0.;
178 }
179 }
180 }
181 }
182
183 Paso_Pattern_coup(A,mis_marker,0.05);
184 }
185
186
187
188
189 void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
190 {
191 dim_t i,j;
192 index_t iptr;
193 double diag,eps_Aii,val;
194 dim_t n=A->numRows;
195 double* diags=MEMALLOC(n,double);
196
197 if (A->pattern->type & PATTERN_FORMAT_SYM) {
198 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
199 return;
200 }
201
202
203 #pragma omp parallel for private(i,iptr,diag) schedule(static)
204 for (i=0;i<n;++i) {
205 diag = 0;
206 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
207 if(A->pattern->index[iptr] != i){
208 diag+=A->val[iptr];
209 }
210 }
211 diags[i]=ABS(diag);
212 }
213
214
215 #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
216 for (i=0;i<n;++i) {
217 if (mis_marker[i]==IS_AVAILABLE) {
218 eps_Aii = theta*theta*diags[i];
219 val=0.;
220 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
221 j=A->pattern->index[iptr];
222 val=A->val[iptr];
223 if(j!= i) {
224 if(val*val>=eps_Aii * diags[j]) {
225 A->val[iptr]=val;
226 }
227 else {
228 A->val[iptr]=0.;
229 }
230 }
231 }
232 }
233 }
234
235 Paso_Pattern_coup(A,mis_marker,0.05);
236 MEMFREE(diags);
237 }
238
239
240 #undef IS_AVAILABLE
241 #undef IS_IN_SET
242 #undef IS_REMOVED
243
244 void Paso_Pattern_color1(Paso_SparseMatrix* A, index_t* num_colors, index_t* colorOf) {
245 index_t out=0, *mis_marker=NULL,i;
246 dim_t n=A->numRows;
247 mis_marker=TMPMEMALLOC(n,index_t);
248 if ( !Paso_checkPtr(mis_marker) ) {
249 /* get coloring */
250 #pragma omp parallel for private(i) schedule(static)
251 for (i = 0; i < n; ++i) {
252 colorOf[i]=-1;
253 mis_marker[i]=-1;
254 }
255
256 while (Paso_Util_isAny(n,colorOf,-1) && Paso_noError()) {
257 Paso_Pattern_coup(A,mis_marker,0.05);
258
259 #pragma omp parallel for private(i) schedule(static)
260 for (i = 0; i < n; ++i) {
261 if (mis_marker[i]) colorOf[i]=out;
262 mis_marker[i]=colorOf[i];
263 }
264 ++out;
265 }
266 TMPMEMFREE(mis_marker);
267 *num_colors=out;
268 }
269 }

  ViewVC Help
Powered by ViewVC 1.1.26