/[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 1975 - (show annotations)
Thu Nov 6 03:07:02 2008 UTC (10 years, 10 months ago) by artak
File MIME type: text/plain
File size: 8853 byte(s)
Some variables are initialized to get rid of compiler warnings
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 index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
44 dim_t i,j;
45 /*double threshold=0.05;*/
46 index_t iptr,*index,*where_p,diagptr;
47 bool_t fail;
48 dim_t n=A->numRows;
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 /* is there any vertex available ?*/
56 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
57
58 #pragma omp parallel for private(i,iptr,index,where_p,diagptr,j) schedule(static)
59 for (i=0;i<n;++i) {
60 if (mis_marker[i]==IS_AVAILABLE) {
61 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 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
74 j=A->pattern->index[iptr]-index_offset;
75 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {
76 mis_marker[j]=IS_REMOVED;
77 }
78 }
79 }
80 }
81
82 #pragma omp parallel for private(i,sum,iptr) schedule(static)
83 for (i=0;i<n;++i)
84 if(mis_marker[i]==IS_AVAILABLE)
85 mis_marker[i]=IS_IN_SET;
86
87 #pragma omp parallel for private(i,iptr,index,where_p,diagptr) schedule(static)
88 for (i=0;i<n;i++) {
89 if (mis_marker[i]==IS_REMOVED) {
90 fail=FALSE;
91 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 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
104 j=A->pattern->index[iptr]-index_offset;
105 if (mis_marker[j]==IS_IN_SET && (A->val[iptr]/A->val[diagptr])<-threshold){
106 fail=TRUE;
107 #ifndef _OPENMP
108 break;
109 #endif
110 }
111 }
112 if(!fail) {
113 mis_marker[i]=IS_IN_SET;
114 }
115 }
116 }
117
118 /* #pragma omp parallel for private(i,sum,iptr) schedule(static)
119 for (i=0;i<n;++i)
120 if(mis_marker[i]==IS_IN_SET)
121 {
122 sum=0.;
123 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
124 sum+=A->val[iptr];
125 }
126 if(ABS(sum)<ZERO)
127 mis_marker[i]=IS_REMOVED;
128 }
129 */
130 }
131 /* swap to TRUE/FALSE in mis_marker */
132 #pragma omp parallel for private(i) schedule(static)
133 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);
134 }
135
136 /*
137 *
138 * Return a strength of connection mask using the classical
139 * strength of connection measure by Ruge and Stuben.
140 *
141 * Specifically, an off-diagonal entry A[i.j] is a strong
142 * connection if:
143 *
144 * -A[i,j] >= theta * max( -A[i,k] ) where k != i
145 *
146 * Otherwise, the connection is weak.
147 *
148 */
149 void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
150 {
151 index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
152 dim_t i;
153 index_t iptr;
154 double threshold,min_offdiagonal;
155 dim_t n=A->numRows;
156 if (A->pattern->type & PATTERN_FORMAT_SYM) {
157 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
158 return;
159 }
160
161 /* is there any vertex available ?*/
162 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
163
164 #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
165 for (i=0;i<n;++i) {
166 min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
167 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
168 if(A->pattern->index[iptr-index_offset] != i){
169 min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
170 }
171 }
172
173 threshold = theta*min_offdiagonal;
174 mis_marker[i]=IS_IN_SET;
175 #pragma omp parallel for private(iptr) schedule(static)
176 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
177 if(A->val[iptr-index_offset] < threshold){
178 mis_marker[i]=IS_REMOVED;
179 #ifndef _OPENMP
180 break;
181 #endif
182 }
183 }
184 }
185 }
186 /* swap to TRUE/FALSE in mis_marker */
187 #pragma omp parallel for private(i) schedule(static)
188 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);
189 }
190
191
192
193
194 void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
195 {
196 index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
197 dim_t i,j;
198 index_t iptr;
199 double diag,eps_Aii,val;
200 dim_t n=A->numRows;
201 double* diags=MEMALLOC(n,double);
202
203 if (A->pattern->type & PATTERN_FORMAT_SYM) {
204 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
205 return;
206 }
207
208 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
209 #pragma omp parallel for private(i,iptr,diag) schedule(static)
210 for (i=0;i<n;++i) {
211 diag = 0;
212 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
213 if(A->pattern->index[iptr-index_offset] != i){
214 diag+=A->val[iptr-index_offset];
215 }
216 }
217 diags[i]=ABS(diag);
218 }
219
220 #pragma omp parallel for private(i,iptr,diag,j,val,eps_Aii) schedule(static)
221 for (i=0;i<n;++i) {
222 eps_Aii = theta*theta*diags[i];
223 mis_marker[i]=IS_REMOVED;
224
225 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
226 j=A->pattern->index[iptr-index_offset];
227 val=A->val[iptr-index_offset];
228 if(j!= i && val*val>=eps_Aii * diags[j]){
229 mis_marker[i]=IS_IN_SET;
230 }
231 }
232 }
233 }
234 /* swap to TRUE/FALSE in mis_marker */
235 #pragma omp parallel for private(i) schedule(static)
236 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET);
237
238 MEMFREE(diags);
239 }
240
241
242 #undef IS_AVAILABLE
243 #undef IS_IN_SET
244 #undef IS_REMOVED

  ViewVC Help
Powered by ViewVC 1.1.26