/[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 2108 - (show annotations)
Fri Nov 28 05:09:23 2008 UTC (10 years, 9 months ago) by gross
File MIME type: text/plain
File size: 7335 byte(s)
some minor changes to PCG and some extra suspicious characters.
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 index_t iptr,*index,*where_p,*diagptr;
46 bool_t passed=FALSE;
47 dim_t n=A->numRows;
48 diagptr=MEMALLOC(n,index_t);
49 /*double sum;*/
50 if (A->pattern->type & PATTERN_FORMAT_SYM) {
51 Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");
52 return;
53 }
54
55 #pragma omp parallel for private(i) schedule(static)
56 for (i=0;i<n;++i)
57 if(mis_marker[i]==IS_AVAILABLE)
58 mis_marker[i]=IS_IN_SET;
59
60 #pragma omp parallel for private(i,index,where_p) schedule(static)
61 for (i=0;i<n;++i) {
62 diagptr[i]=A->pattern->ptr[i];
63 index=&(A->pattern->index[diagptr[i]]);
64 where_p=(index_t*)bsearch(&i,
65 index,
66 A->pattern->ptr[i + 1]-A->pattern->ptr[i],
67 sizeof(index_t),
68 Paso_comparIndex);
69 if (where_p==NULL) {
70 Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
71 } else {
72 diagptr[i]+=(index_t)(where_p-index);
73 }
74 }
75
76
77
78 /*This loop cannot be parallelized, as order matters here.*/
79 for (i=0;i<n;++i) {
80 if (mis_marker[i]==IS_IN_SET) {
81 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
82 j=A->pattern->index[iptr];
83 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
84 mis_marker[j]=IS_REMOVED;
85 }
86 }
87 }
88 }
89
90
91
92 /*This loop cannot be parallelized, as order matters here.*/
93 for (i=0;i<n;i++) {
94 if (mis_marker[i]==IS_REMOVED) {
95 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
96 j=A->pattern->index[iptr];
97 if (mis_marker[j]==IS_IN_SET) {
98 if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
99 passed=TRUE;
100 }
101 else {
102 passed=FALSE;
103 break;
104 }
105 }
106 }
107 if (passed) {
108 mis_marker[i]=IS_IN_SET;
109 passed=FALSE;
110 }
111 }
112 }
113
114
115 /* swap to TRUE/FALSE in mis_marker */
116 #pragma omp parallel for private(i) schedule(static)
117 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
118
119 MEMFREE(diagptr);
120 }
121
122 /*
123 *
124 * Return a strength of connection mask using the classical
125 * strength of connection measure by Ruge and Stuben.
126 *
127 * Specifically, an off-diagonal entry A[i.j] is a strong
128 * connection if:
129 *
130 * -A[i,j] >= theta * max( -A[i,k] ) where k != i
131 *
132 * Otherwise, the connection is weak.
133 *
134 */
135 void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
136 {
137 dim_t i;
138 index_t iptr;
139 double threshold,max_offdiagonal;
140 dim_t n=A->numRows;
141 if (A->pattern->type & PATTERN_FORMAT_SYM) {
142 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
143 return;
144 }
145
146 #pragma omp parallel for private(i,iptr,max_offdiagonal,threshold) schedule(static)
147 for (i=0;i<n;++i) {
148 if(mis_marker[i]==IS_AVAILABLE) {
149 /*min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];*/
150 max_offdiagonal = -1.e+15;
151 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
152 if(A->pattern->index[iptr] != i){
153 max_offdiagonal = MAX(max_offdiagonal,-(A->val[iptr]));
154 }
155 }
156
157 threshold = theta*max_offdiagonal;
158 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
159 if(-(A->val[iptr])<threshold && A->pattern->index[iptr]!=i){
160 A->val[iptr]=0.;
161 }
162 }
163 }
164 }
165
166 Paso_Pattern_coup(A,mis_marker,0.05);
167 }
168
169
170
171
172 void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
173 {
174 dim_t i,j;
175 index_t iptr;
176 double diag,eps_Aii,val;
177 dim_t n=A->numRows;
178 double* diags=MEMALLOC(n,double);
179
180 if (A->pattern->type & PATTERN_FORMAT_SYM) {
181 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
182 return;
183 }
184
185
186 #pragma omp parallel for private(i,iptr,diag) schedule(static)
187 for (i=0;i<n;++i) {
188 diag = 0;
189 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
190 if(A->pattern->index[iptr] != i){
191 diag+=A->val[iptr];
192 }
193 }
194 diags[i]=ABS(diag);
195 }
196
197
198 #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
199 for (i=0;i<n;++i) {
200 if (mis_marker[i]==IS_AVAILABLE) {
201 eps_Aii = theta*theta*diags[i];
202 val=0.;
203 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
204 j=A->pattern->index[iptr];
205 val=A->val[iptr];
206 if(j!= i) {
207 if(val*val>=eps_Aii * diags[j]) {
208 A->val[iptr]=val;
209 }
210 else {
211 A->val[iptr]=0.;
212 }
213 }
214 }
215 }
216 }
217
218 Paso_Pattern_coup(A,mis_marker,0.05);
219 MEMFREE(diags);
220 }
221
222
223 #undef IS_AVAILABLE
224 #undef IS_IN_SET
225 #undef IS_REMOVED
226
227 void Paso_Pattern_color1(Paso_SparseMatrix* A, index_t* num_colors, index_t* colorOf) {
228 index_t out=0, *mis_marker=NULL,i;
229 dim_t n=A->numRows;
230 mis_marker=TMPMEMALLOC(n,index_t);
231 if ( !Paso_checkPtr(mis_marker) ) {
232 /* get coloring */
233 #pragma omp parallel for private(i) schedule(static)
234 for (i = 0; i < n; ++i) {
235 colorOf[i]=-1;
236 mis_marker[i]=-1;
237 }
238
239 while (Paso_Util_isAny(n,colorOf,-1) && Paso_noError()) {
240 Paso_Pattern_coup(A,mis_marker,0.05);
241
242 #pragma omp parallel for private(i) schedule(static)
243 for (i = 0; i < n; ++i) {
244 if (mis_marker[i]) colorOf[i]=out;
245 mis_marker[i]=colorOf[i];
246 }
247 ++out;
248 }
249 TMPMEMFREE(mis_marker);
250 *num_colors=out;
251 }
252 }

  ViewVC Help
Powered by ViewVC 1.1.26