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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2122 - (hide annotations)
Wed Dec 3 02:52:28 2008 UTC (10 years, 9 months ago) by artak
File MIME type: text/plain
File size: 7984 byte(s)
Codition is added in the corsening process to avoid rows with very small row sum.
1 jfenwick 1851
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 phornby 1914 #include "PasoUtil.h"
32 phornby 1913 #include "Pattern_coupling.h"
33 jfenwick 1851
34    
35     /***************************************************************/
36    
37     #define IS_AVAILABLE -1
38 artak 1975 #define IS_IN_SET -3
39     #define IS_REMOVED -4
40 jfenwick 1851
41 artak 1931 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
42 jfenwick 1851
43 artak 1881 dim_t i,j;
44 artak 1931 /*double threshold=0.05;*/
45 artak 2122 double sum;
46 artak 2107 index_t iptr,*index,*where_p,*diagptr;
47     bool_t passed=FALSE;
48 artak 1954 dim_t n=A->numRows;
49 artak 2107 diagptr=MEMALLOC(n,index_t);
50 artak 2122
51 jfenwick 1851 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 artak 2107 #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 jfenwick 1851
61 artak 2107 #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 artak 1902 }
102 artak 2107 else {
103     passed=FALSE;
104     break;
105 artak 1902 }
106 artak 2107 }
107 artak 1902 }
108 artak 2107 if (passed) {
109     mis_marker[i]=IS_IN_SET;
110     passed=FALSE;
111     }
112 artak 1902 }
113 artak 2107 }
114    
115 artak 2122 /* 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     #pragma omp parallel for private(i,iptr,j,sum) schedule(static)
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 jfenwick 1851 /* swap to TRUE/FALSE in mis_marker */
133     #pragma omp parallel for private(i) schedule(static)
134 artak 2107 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
135    
136     MEMFREE(diagptr);
137 jfenwick 1851 }
138 artak 1890
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 phornby 1917 dim_t i;
155 phornby 1907 index_t iptr;
156 artak 2107 double threshold,max_offdiagonal;
157 artak 1954 dim_t n=A->numRows;
158 artak 1890 if (A->pattern->type & PATTERN_FORMAT_SYM) {
159 artak 1902 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
160 artak 1890 return;
161     }
162    
163 artak 2107 #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 artak 1890 }
172     }
173    
174 artak 2107 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 artak 1890 }
179     }
180 artak 2107 }
181 artak 1890 }
182 artak 2107
183     Paso_Pattern_coup(A,mis_marker,0.05);
184 artak 1890 }
185 artak 1954
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 artak 2107
203 artak 1954 #pragma omp parallel for private(i,iptr,diag) schedule(static)
204     for (i=0;i<n;++i) {
205     diag = 0;
206 artak 2107 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 artak 1954 }
210     }
211     diags[i]=ABS(diag);
212     }
213    
214 artak 2107
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 artak 1954 eps_Aii = theta*theta*diags[i];
219 artak 2107 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 artak 1954 }
231     }
232 artak 2107 }
233     }
234    
235     Paso_Pattern_coup(A,mis_marker,0.05);
236 artak 1954 MEMFREE(diags);
237     }
238    
239    
240 jfenwick 1851 #undef IS_AVAILABLE
241 artak 1975 #undef IS_IN_SET
242     #undef IS_REMOVED
243 artak 2107
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