/[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 1954 - (hide annotations)
Fri Oct 31 03:22:34 2008 UTC (10 years, 10 months ago) by artak
File MIME type: text/plain
File size: 8987 byte(s)
condition is added to check wheter row sum is ZERO and some new coarsening algorithm is also added
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     #define IS_IN_MIS_NOW -2
39     #define IS_IN_MIS -3
40     #define IS_CONNECTED_TO_MIS -4
41 artak 1954 #define ZERO 1.e-10
42 jfenwick 1851
43 artak 1890
44 artak 1931 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
45 jfenwick 1851
46     index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
47 artak 1881 dim_t i,j;
48 artak 1931 /*double threshold=0.05;*/
49 artak 1902 index_t iptr,*index,*where_p,diagptr;
50 artak 1936 bool_t fail;
51 artak 1954 dim_t n=A->numRows;
52     double sum;
53 jfenwick 1851 if (A->pattern->type & PATTERN_FORMAT_SYM) {
54     Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");
55     return;
56     }
57    
58     /* is there any vertex available ?*/
59 artak 1954 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
60 jfenwick 1851
61 artak 1931 #pragma omp parallel for private(i,iptr,index,where_p,diagptr,j) schedule(static)
62     for (i=0;i<n;++i) {
63 jfenwick 1851 if (mis_marker[i]==IS_AVAILABLE) {
64 artak 1902 diagptr=A->pattern->ptr[i];
65     index=&(A->pattern->index[diagptr]);
66     where_p=(index_t*)bsearch(&i,
67     index,
68     A->pattern->ptr[i + 1]-A->pattern->ptr[i],
69     sizeof(index_t),
70     Paso_comparIndex);
71     if (where_p==NULL) {
72     Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
73     } else {
74     diagptr+=(index_t)(where_p-index);
75     }
76 jfenwick 1851 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
77 artak 1902 j=A->pattern->index[iptr]-index_offset;
78 artak 1931 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {
79     mis_marker[j]=IS_CONNECTED_TO_MIS;
80 artak 1954 }
81 jfenwick 1851 }
82 artak 1902 }
83     }
84 artak 1931
85 artak 1954 #pragma omp parallel for private(i,sum,iptr) schedule(static)
86 artak 1931 for (i=0;i<n;++i)
87     if(mis_marker[i]==IS_AVAILABLE)
88 artak 1954 mis_marker[i]=IS_IN_MIS;
89    
90     #pragma omp parallel for private(i,iptr,index,where_p,diagptr) schedule(static)
91 artak 1902 for (i=0;i<n;i++) {
92 artak 1931 if (mis_marker[i]==IS_CONNECTED_TO_MIS) {
93     fail=FALSE;
94 artak 1902 diagptr=A->pattern->ptr[i];
95     index=&(A->pattern->index[diagptr]);
96     where_p=(index_t*)bsearch(&i,
97     index,
98     A->pattern->ptr[i + 1]-A->pattern->ptr[i],
99     sizeof(index_t),
100     Paso_comparIndex);
101     if (where_p==NULL) {
102     Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
103     } else {
104     diagptr+=(index_t)(where_p-index);
105     }
106 jfenwick 1851 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
107 artak 1902 j=A->pattern->index[iptr]-index_offset;
108 artak 1931 if (mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])<-threshold){
109     fail=TRUE;
110 artak 1954 #ifndef _OPENMP
111 artak 1931 break;
112 artak 1954 #endif
113 artak 1883 }
114 jfenwick 1851 }
115 artak 1954 if(!fail) {
116 artak 1931 mis_marker[i]=IS_IN_MIS;
117 artak 1954 }
118 artak 1902 }
119 artak 1954 }
120    
121     #pragma omp parallel for private(i,sum,iptr) schedule(static)
122     for (i=0;i<n;++i)
123     if(mis_marker[i]==IS_IN_MIS)
124     {
125     sum=0.;
126     for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
127     sum+=A->val[iptr];
128     }
129     if(ABS(sum)<ZERO)
130     mis_marker[i]=IS_CONNECTED_TO_MIS;
131     }
132    
133 artak 1902 }
134 jfenwick 1851 /* swap to TRUE/FALSE in mis_marker */
135     #pragma omp parallel for private(i) schedule(static)
136     for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
137     }
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     index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
155 phornby 1917 dim_t i;
156 phornby 1907 index_t iptr;
157 artak 1890 double threshold,min_offdiagonal;
158 artak 1954 dim_t n=A->numRows;
159 artak 1890 if (A->pattern->type & PATTERN_FORMAT_SYM) {
160 artak 1902 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
161 artak 1890 return;
162     }
163    
164     /* is there any vertex available ?*/
165 artak 1902 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
166 artak 1890
167 artak 1902 #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
168 artak 1890 for (i=0;i<n;++i) {
169 artak 1902 min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
170 artak 1890 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
171 artak 1954 if(A->pattern->index[iptr-index_offset] != i){
172 artak 1902 min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
173 artak 1890 }
174     }
175    
176     threshold = theta*min_offdiagonal;
177 artak 1931 mis_marker[i]=IS_IN_MIS;
178 artak 1890 #pragma omp parallel for private(iptr) schedule(static)
179     for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
180 artak 1931 if(A->val[iptr-index_offset] < threshold){
181     mis_marker[i]=IS_CONNECTED_TO_MIS;
182 artak 1936 #ifndef _OPENMP
183     break;
184     #endif
185 artak 1890 }
186     }
187     }
188     }
189     /* swap to TRUE/FALSE in mis_marker */
190     #pragma omp parallel for private(i) schedule(static)
191     for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
192     }
193 artak 1954
194    
195    
196    
197     void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
198     {
199     index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
200     dim_t i,j;
201     index_t iptr;
202     double diag,eps_Aii,val;
203     dim_t n=A->numRows;
204     double* diags=MEMALLOC(n,double);
205    
206     if (A->pattern->type & PATTERN_FORMAT_SYM) {
207     Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
208     return;
209     }
210    
211     if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
212     #pragma omp parallel for private(i,iptr,diag) schedule(static)
213     for (i=0;i<n;++i) {
214     diag = 0;
215     for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
216     if(A->pattern->index[iptr-index_offset] != i){
217     diag+=A->val[iptr-index_offset];
218     }
219     }
220     diags[i]=ABS(diag);
221     }
222    
223     #pragma omp parallel for private(i,iptr,diag) schedule(static)
224     for (i=0;i<n;++i) {
225     eps_Aii = theta*theta*diags[i];
226     mis_marker[i]=IS_CONNECTED_TO_MIS;
227    
228     for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
229     j=A->pattern->index[iptr-index_offset];
230     val=A->val[iptr-index_offset];
231     if(j!= i && val*val>=eps_Aii * diags[j]){
232     mis_marker[i]=IS_IN_MIS;
233     }
234     }
235     }
236     }
237     /* swap to TRUE/FALSE in mis_marker */
238     #pragma omp parallel for private(i) schedule(static)
239     for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
240    
241     MEMFREE(diags);
242     }
243    
244    
245    
246    
247    
248    
249    
250    
251    
252    
253    
254    
255    
256 jfenwick 1851 #undef IS_AVAILABLE
257     #undef IS_IN_MIS_NOW
258     #undef IS_IN_MIS
259 artak 1954 #undef IS_CONNECTED_TO_MIS
260     #undef ZERO

  ViewVC Help
Powered by ViewVC 1.1.26