/[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 1980 - (hide annotations)
Thu Nov 6 04:45:08 2008 UTC (10 years, 10 months ago) by artak
File MIME type: text/plain
File size: 8324 byte(s)
sum variable is removed from #pragma
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     index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
44 artak 1881 dim_t i,j;
45 artak 1931 /*double threshold=0.05;*/
46 artak 1902 index_t iptr,*index,*where_p,diagptr;
47 artak 1936 bool_t fail;
48 artak 1954 dim_t n=A->numRows;
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     /* is there any vertex available ?*/
56 artak 1954 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
57 jfenwick 1851
58 artak 1931 #pragma omp parallel for private(i,iptr,index,where_p,diagptr,j) schedule(static)
59     for (i=0;i<n;++i) {
60 jfenwick 1851 if (mis_marker[i]==IS_AVAILABLE) {
61 artak 1902 diagptr=A->pattern->ptr[i];
62     index=&(A->pattern->index[diagptr]);
63     where_p=(index_t*)bsearch(&i,
64     index,
65     A->pattern->ptr[i + 1]-A->pattern->ptr[i],
66     sizeof(index_t),
67     Paso_comparIndex);
68     if (where_p==NULL) {
69     Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
70     } else {
71     diagptr+=(index_t)(where_p-index);
72     }
73 jfenwick 1851 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
74 artak 1902 j=A->pattern->index[iptr]-index_offset;
75 artak 1931 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {
76 artak 1975 mis_marker[j]=IS_REMOVED;
77 artak 1954 }
78 jfenwick 1851 }
79 artak 1902 }
80     }
81 artak 1931
82 artak 1980 #pragma omp parallel for private(i,iptr) schedule(static)
83 artak 1931 for (i=0;i<n;++i)
84     if(mis_marker[i]==IS_AVAILABLE)
85 artak 1975 mis_marker[i]=IS_IN_SET;
86 artak 1954
87     #pragma omp parallel for private(i,iptr,index,where_p,diagptr) schedule(static)
88 artak 1902 for (i=0;i<n;i++) {
89 artak 1975 if (mis_marker[i]==IS_REMOVED) {
90 artak 1931 fail=FALSE;
91 artak 1902 diagptr=A->pattern->ptr[i];
92     index=&(A->pattern->index[diagptr]);
93     where_p=(index_t*)bsearch(&i,
94     index,
95     A->pattern->ptr[i + 1]-A->pattern->ptr[i],
96     sizeof(index_t),
97     Paso_comparIndex);
98     if (where_p==NULL) {
99     Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
100     } else {
101     diagptr+=(index_t)(where_p-index);
102     }
103 jfenwick 1851 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
104 artak 1902 j=A->pattern->index[iptr]-index_offset;
105 artak 1975 if (mis_marker[j]==IS_IN_SET && (A->val[iptr]/A->val[diagptr])<-threshold){
106 artak 1931 fail=TRUE;
107 artak 1954 #ifndef _OPENMP
108 artak 1931 break;
109 artak 1954 #endif
110 artak 1883 }
111 jfenwick 1851 }
112 artak 1954 if(!fail) {
113 artak 1975 mis_marker[i]=IS_IN_SET;
114 artak 1954 }
115 artak 1902 }
116 artak 1954 }
117    
118 artak 1902 }
119 jfenwick 1851 /* swap to TRUE/FALSE in mis_marker */
120     #pragma omp parallel for private(i) schedule(static)
121 artak 1975 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);
122 jfenwick 1851 }
123 artak 1890
124     /*
125     *
126     * Return a strength of connection mask using the classical
127     * strength of connection measure by Ruge and Stuben.
128     *
129     * Specifically, an off-diagonal entry A[i.j] is a strong
130     * connection if:
131     *
132     * -A[i,j] >= theta * max( -A[i,k] ) where k != i
133     *
134     * Otherwise, the connection is weak.
135     *
136     */
137     void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
138     {
139     index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
140 phornby 1917 dim_t i;
141 phornby 1907 index_t iptr;
142 artak 1890 double threshold,min_offdiagonal;
143 artak 1954 dim_t n=A->numRows;
144 artak 1890 if (A->pattern->type & PATTERN_FORMAT_SYM) {
145 artak 1902 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
146 artak 1890 return;
147     }
148    
149     /* is there any vertex available ?*/
150 artak 1902 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
151 artak 1890
152 artak 1902 #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
153 artak 1890 for (i=0;i<n;++i) {
154 artak 1902 min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
155 artak 1890 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
156 artak 1954 if(A->pattern->index[iptr-index_offset] != i){
157 artak 1902 min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
158 artak 1890 }
159     }
160    
161     threshold = theta*min_offdiagonal;
162 artak 1975 mis_marker[i]=IS_IN_SET;
163 artak 1890 #pragma omp parallel for private(iptr) schedule(static)
164     for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
165 artak 1931 if(A->val[iptr-index_offset] < threshold){
166 artak 1975 mis_marker[i]=IS_REMOVED;
167 artak 1936 #ifndef _OPENMP
168     break;
169     #endif
170 artak 1890 }
171     }
172     }
173     }
174     /* swap to TRUE/FALSE in mis_marker */
175     #pragma omp parallel for private(i) schedule(static)
176 artak 1975 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);
177 artak 1890 }
178 artak 1954
179    
180    
181    
182     void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
183     {
184     index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
185     dim_t i,j;
186     index_t iptr;
187     double diag,eps_Aii,val;
188     dim_t n=A->numRows;
189     double* diags=MEMALLOC(n,double);
190    
191     if (A->pattern->type & PATTERN_FORMAT_SYM) {
192     Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
193     return;
194     }
195    
196     if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
197     #pragma omp parallel for private(i,iptr,diag) schedule(static)
198     for (i=0;i<n;++i) {
199     diag = 0;
200     for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
201     if(A->pattern->index[iptr-index_offset] != i){
202     diag+=A->val[iptr-index_offset];
203     }
204     }
205     diags[i]=ABS(diag);
206     }
207    
208 artak 1975 #pragma omp parallel for private(i,iptr,diag,j,val,eps_Aii) schedule(static)
209 artak 1954 for (i=0;i<n;++i) {
210     eps_Aii = theta*theta*diags[i];
211 artak 1975 mis_marker[i]=IS_REMOVED;
212 artak 1954
213     for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
214     j=A->pattern->index[iptr-index_offset];
215     val=A->val[iptr-index_offset];
216     if(j!= i && val*val>=eps_Aii * diags[j]){
217 artak 1975 mis_marker[i]=IS_IN_SET;
218 artak 1954 }
219     }
220     }
221     }
222     /* swap to TRUE/FALSE in mis_marker */
223     #pragma omp parallel for private(i) schedule(static)
224 artak 1975 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);
225 artak 1954
226     MEMFREE(diags);
227     }
228    
229    
230 jfenwick 1851 #undef IS_AVAILABLE
231 artak 1975 #undef IS_IN_SET
232     #undef IS_REMOVED

  ViewVC Help
Powered by ViewVC 1.1.26