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,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 |
} |
119 |
/* swap to TRUE/FALSE in mis_marker */ |
120 |
#pragma omp parallel for private(i) schedule(static) |
121 |
for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET); |
122 |
} |
123 |
|
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 |
dim_t i; |
141 |
index_t iptr; |
142 |
double threshold,min_offdiagonal; |
143 |
dim_t n=A->numRows; |
144 |
if (A->pattern->type & PATTERN_FORMAT_SYM) { |
145 |
Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet"); |
146 |
return; |
147 |
} |
148 |
|
149 |
/* is there any vertex available ?*/ |
150 |
if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) { |
151 |
|
152 |
#pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static) |
153 |
for (i=0;i<n;++i) { |
154 |
min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset]; |
155 |
for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) { |
156 |
if(A->pattern->index[iptr-index_offset] != i){ |
157 |
min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]); |
158 |
} |
159 |
} |
160 |
|
161 |
threshold = theta*min_offdiagonal; |
162 |
mis_marker[i]=IS_IN_SET; |
163 |
#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 |
if(A->val[iptr-index_offset] < threshold){ |
166 |
mis_marker[i]=IS_REMOVED; |
167 |
#ifndef _OPENMP |
168 |
break; |
169 |
#endif |
170 |
} |
171 |
} |
172 |
} |
173 |
} |
174 |
/* swap to TRUE/FALSE in mis_marker */ |
175 |
#pragma omp parallel for private(i) schedule(static) |
176 |
for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET); |
177 |
} |
178 |
|
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 |
#pragma omp parallel for private(i,iptr,diag,j,val,eps_Aii) schedule(static) |
209 |
for (i=0;i<n;++i) { |
210 |
eps_Aii = theta*theta*diags[i]; |
211 |
mis_marker[i]=IS_REMOVED; |
212 |
|
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 |
mis_marker[i]=IS_IN_SET; |
218 |
} |
219 |
} |
220 |
} |
221 |
} |
222 |
/* swap to TRUE/FALSE in mis_marker */ |
223 |
#pragma omp parallel for private(i) schedule(static) |
224 |
for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_SET); |
225 |
|
226 |
MEMFREE(diags); |
227 |
} |
228 |
|
229 |
|
230 |
#undef IS_AVAILABLE |
231 |
#undef IS_IN_SET |
232 |
#undef IS_REMOVED |