/[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 2107 - (hide annotations)
Fri Nov 28 04:39:07 2008 UTC (10 years, 9 months ago) by artak
File MIME type: text/plain
File size: 7337 byte(s)
Current version of AMG uses ILU for relaxation, however it is not stable when schur matrix becames more denser. For example it is not stable for -e paramenter is more 400 in 2D for convection problem 
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 2107 index_t iptr,*index,*where_p,*diagptr;
46     bool_t passed=FALSE;
47 artak 1954 dim_t n=A->numRows;
48 artak 2107 diagptr=MEMALLOC(n,index_t);
49 artak 1975 /*double sum;*/
50 jfenwick 1851 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 artak 2107 #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 jfenwick 1851
60 artak 2107 #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 artak 1902 }
101 artak 2107 else {
102     passed=FALSE;
103     break;
104 artak 1902 }
105 artak 2107 }
106 artak 1902 }
107 artak 2107 if (passed) {
108     mis_marker[i]=IS_IN_SET;
109     passed=FALSE;
110     }
111 artak 1902 }
112 artak 2107 }
113    
114    
115 jfenwick 1851 /* swap to TRUE/FALSE in mis_marker */
116     #pragma omp parallel for private(i) schedule(static)
117 artak 2107 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
118    
119     MEMFREE(diagptr);
120 jfenwick 1851 }
121 artak 1890
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 phornby 1917 dim_t i;
138 phornby 1907 index_t iptr;
139 artak 2107 double threshold,max_offdiagonal;
140 artak 1954 dim_t n=A->numRows;
141 artak 1890 if (A->pattern->type & PATTERN_FORMAT_SYM) {
142 artak 1902 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
143 artak 1890 return;
144     }
145    
146 artak 2107 #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 artak 1890 }
155     }
156    
157 artak 2107 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 artak 1890 }
162     }
163 artak 2107 }
164 artak 1890 }
165 artak 2107
166     Paso_Pattern_coup(A,mis_marker,0.05);
167 artak 1890 }
168 artak 1954
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 artak 2107
186 artak 1954 #pragma omp parallel for private(i,iptr,diag) schedule(static)
187     for (i=0;i<n;++i) {
188     diag = 0;
189 artak 2107 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 artak 1954 }
193     }
194     diags[i]=ABS(diag);
195     }
196    
197 artak 2107
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 artak 1954 eps_Aii = theta*theta*diags[i];
202 artak 2107 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 artak 1954 }
214     }
215 artak 2107 }
216     }
217    
218     Paso_Pattern_coup(A,mis_marker,0.05);
219 artak 1954 MEMFREE(diags);
220     }
221    
222    
223 jfenwick 1851 #undef IS_AVAILABLE
224 artak 1975 #undef IS_IN_SET
225     #undef IS_REMOVED
226 artak 2107
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     }
253    
254    

  ViewVC Help
Powered by ViewVC 1.1.26