/[escript]/trunk/paso/src/Pattern.c
ViewVC logotype

Contents of /trunk/paso/src/Pattern.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3283 - (show annotations)
Mon Oct 18 22:39:28 2010 UTC (9 years ago) by gross
File MIME type: text/plain
File size: 10050 byte(s)
AMG reengineered, BUG is Smoother fixed.


1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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
18 /* Paso: Pattern */
19
20 /**************************************************************/
21
22 /* Author: Lutz Gross, l.gross@uq.edu.au */
23
24 /**************************************************************/
25
26 #include "Paso.h"
27 #include "Pattern.h"
28
29 /**************************************************************/
30
31 /* allocates a Pattern */
32
33 Paso_Pattern* Paso_Pattern_alloc(int type, dim_t numOutput, dim_t numInput, index_t* ptr, index_t* index) {
34 Paso_Pattern *out=NULL;
35 index_t index_offset=(type & PATTERN_FORMAT_OFFSET1 ? 1:0);
36 index_t loc_min_index,loc_max_index,min_index=index_offset,max_index=index_offset-1;
37 dim_t i;
38 Esys_resetError();
39
40 if (ptr!=NULL && index != NULL) {
41 #pragma omp parallel private(loc_min_index,loc_max_index,i)
42 {
43 loc_min_index=index_offset;
44 loc_max_index=index_offset-1;
45 if (type & PATTERN_FORMAT_OFFSET1) {
46 #pragma omp for schedule(static)
47 for (i=0;i<numOutput;++i) {
48 if (ptr[i]<ptr[i+1]) {
49 #ifdef USE_QSORTG
50 qsortG(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
51 #else
52 qsort(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
53 #endif
54 loc_min_index=MIN(loc_min_index,index[ptr[i]-1]);
55 loc_max_index=MAX(loc_max_index,index[ptr[i+1]-2]);
56 }
57 }
58 } else {
59 #pragma omp for schedule(static)
60 for (i=0;i<numOutput;++i) {
61 if (ptr[i]<ptr[i+1]) {
62 #ifdef USE_QSORTG
63 qsortG(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
64 #else
65 qsort(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
66 #endif
67 loc_min_index=MIN(loc_min_index,index[ptr[i]]);
68 loc_max_index=MAX(loc_max_index,index[ptr[i+1]-1]);
69 }
70 }
71 }
72 #pragma omp critical
73 {
74 min_index=MIN(loc_min_index,min_index);
75 max_index=MAX(loc_max_index,max_index);
76 }
77 }
78 if ( (min_index<index_offset) || (max_index>=numInput+index_offset) ) {
79 Esys_setError(TYPE_ERROR,"Paso_Pattern_alloc: Pattern index out of range.");
80 return NULL;
81 }
82 }
83 out=MEMALLOC(1,Paso_Pattern);
84 if (! Esys_checkPtr(out)) {
85 out->type=type;
86 out->reference_counter=1;
87 out->numOutput=numOutput;
88 out->numInput=numInput;
89 out->ptr=ptr;
90 out->index=index;
91 out->main_iptr = NULL;
92 out->coloring = NULL;
93 out->numColors=-1;
94
95 if (out->ptr == NULL) {
96 out->len=0;
97 } else {
98 out->len=out->ptr[out->numOutput] - index_offset;
99 }
100 }
101 return out;
102 }
103
104 /* returns a reference to in */
105
106 Paso_Pattern* Paso_Pattern_getReference(Paso_Pattern* in) {
107 if (in!=NULL) {
108 ++(in->reference_counter);
109 }
110 return in;
111 }
112
113 /* deallocates a Pattern: */
114
115 void Paso_Pattern_free(Paso_Pattern* in) {
116 if (in!=NULL) {
117 in->reference_counter--;
118 if (in->reference_counter<=0) {
119 MEMFREE(in->ptr);
120 MEMFREE(in->index);
121 MEMFREE(in->main_iptr);
122 MEMFREE(in->coloring);
123 MEMFREE(in);
124 }
125 }
126 }
127 /* *************************************************************/
128
129 /* some routines which help to get the matrix pattern from elements: */
130
131 /* this routine is used by qsort called in Paso_Pattern_alloc */
132
133 int Paso_comparIndex(const void *index1,const void *index2){
134 index_t Iindex1,Iindex2;
135 Iindex1=*(index_t*)index1;
136 Iindex2=*(index_t*)index2;
137 if (Iindex1<Iindex2) {
138 return -1;
139 } else {
140 if (Iindex1>Iindex2) {
141 return 1;
142 } else {
143 return 0;
144 }
145 }
146 }
147
148 bool_t Paso_Pattern_isEmpty(Paso_Pattern* in) {
149 if (in != NULL) {
150 if ((in->ptr != NULL) && (in->index != NULL)) return FALSE;
151 }
152 return TRUE;
153 }
154
155
156 /* computes the pattern coming from matrix-matrix multiplication
157 *
158 **/
159
160 Paso_Pattern* Paso_Pattern_multiply(int type, Paso_Pattern* A, Paso_Pattern* B) {
161 Paso_Pattern*out=NULL;
162 index_t iptrA,iptrB;
163 dim_t i,j,k;
164 Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(A->numOutput);
165
166 #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)
167 for(i = 0; i < A->numOutput; i++) {
168 for(iptrA = A->ptr[i]; iptrA < A->ptr[i+1]; ++iptrA) {
169 j = A->index[iptrA];
170 for(iptrB = B->ptr[j]; iptrB < B->ptr[j+1]; ++iptrB) {
171 k = B->index[iptrB];
172 Paso_IndexListArray_insertIndex(index_list,i,k);
173 }
174 }
175 }
176
177 out=Paso_Pattern_fromIndexListArray(0, index_list,0,B->numInput,0);
178
179 /* clean up */
180 Paso_IndexListArray_free(index_list);
181
182 return out;
183 }
184
185
186
187 /*
188 * Computes the pattern of C = A binary operation B for CSR matrices A,B
189 *
190 * Note: we do not check whether A_ij(op)B_ij=0
191 *
192 */
193 Paso_Pattern* Paso_Pattern_binop(int type, Paso_Pattern* A, Paso_Pattern* B) {
194 Paso_Pattern *out=NULL;
195 index_t iptrA,iptrB;
196 dim_t i,j,k;
197
198 Paso_IndexListArray* index_list = Paso_IndexListArray_alloc(A->numOutput);
199
200 #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)
201 for(i = 0; i < B->numOutput; i++){
202 iptrA = A->ptr[i],
203 iptrB = B->ptr[i];
204
205 while (iptrA < A->ptr[i+1] && iptrB < B->ptr[i+1]) {
206 j = A->index[iptrA];
207 k = B->index[iptrB];
208 if (j<k) {
209 Paso_IndexListArray_insertIndex(index_list,i,j);
210 iptrA++;
211 } else if (j>k) {
212 Paso_IndexListArray_insertIndex(index_list,i,k);
213 iptrB++;
214 } else if (j==k) {
215 Paso_IndexListArray_insertIndex(index_list,i,j);
216 iptrB++;
217 iptrA++;
218 }
219 }
220 while(iptrA < A->ptr[i+1]) {
221 j = A->index[iptrA];
222 Paso_IndexListArray_insertIndex(index_list,i,j);
223 iptrA++;
224 }
225 while(iptrB < B->ptr[i+1]) {
226 k = B->index[iptrB];
227 Paso_IndexListArray_insertIndex(index_list,i,k);
228 iptrB++;
229 }
230 }
231
232 out=Paso_Pattern_fromIndexListArray(0, index_list,0,A->numInput,0);
233
234
235 /* clean up */
236 Paso_IndexListArray_free(index_list);
237
238 return out;
239 }
240
241 /* creates a Paso_pattern from a range of indices */
242 Paso_Pattern* Paso_Pattern_fromIndexListArray(dim_t n0, Paso_IndexListArray* index_list_array,index_t range_min,index_t range_max,index_t index_offset)
243 {
244 const dim_t n=index_list_array->n;
245 Paso_IndexList* index_list = index_list_array->index_list;
246 dim_t *ptr=NULL;
247 register dim_t s,i,itmp;
248 index_t *index=NULL;
249 Paso_Pattern* out=NULL;
250
251 ptr=MEMALLOC(n+1-n0,index_t);
252 if (! Esys_checkPtr(ptr) ) {
253 /* get the number of connections per row */
254 #pragma omp parallel for private(i) schedule(static)
255 for(i=n0;i<n;++i) {
256 ptr[i-n0]=Paso_IndexList_count(&index_list[i],range_min,range_max);
257 }
258 /* accumulate ptr */
259 s=0;
260 for(i=n0;i<n;++i) {
261 itmp=ptr[i-n0];
262 ptr[i-n0]=s;
263 s+=itmp;
264 }
265 ptr[n-n0]=s;
266 /* fill index */
267 index=MEMALLOC(ptr[n-n0],index_t);
268 if (! Esys_checkPtr(index)) {
269 #pragma omp parallel for private(i) schedule(static)
270 for(i=n0;i<n;++i) {
271 Paso_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
272 }
273 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
274 }
275 }
276 if (! Esys_noError()) {
277 MEMFREE(ptr);
278 MEMFREE(index);
279 Paso_Pattern_free(out);
280 }
281 return out;
282 }
283
284 index_t* Paso_Pattern_borrowMainDiagonalPointer(Paso_Pattern* A)
285 {
286 const dim_t n=A->numOutput;
287 int fail=0;
288 index_t *index,*where_p, i;
289
290 if (A->main_iptr == NULL) {
291 A->main_iptr=MEMALLOC(n,index_t);
292 if (! Esys_checkPtr(A->main_iptr) ) {
293 #pragma omp parallel
294 {
295 /* identify the main diagonals */
296 #pragma omp for schedule(static) private(i, index, where_p)
297 for (i = 0; i < n; ++i) {
298 index=&(A->index[A->ptr[i]]);
299 where_p=bsearch(&i,
300 index,
301 (size_t) (A->ptr[i + 1]-A->ptr[i]),
302 sizeof(index_t),
303 Paso_comparIndex);
304
305 if (where_p==NULL) {
306 fail=1;
307 } else {
308 A->main_iptr[i]=A->ptr[i]+(index_t)(where_p-index);
309 }
310 }
311
312 }
313 if (fail > 0) {
314 MEMFREE(A->main_iptr);
315 A->main_iptr=NULL;
316 }
317
318 }
319 }
320 return A->main_iptr;
321 }
322
323
324 dim_t Paso_Pattern_getNumColors(Paso_Pattern* A)
325 {
326 Paso_Pattern_borrowColoringPointer(A); /* make sure numColors is defined */
327 return A->numColors;
328 }
329 index_t* Paso_Pattern_borrowColoringPointer(Paso_Pattern* A)
330 {
331 dim_t n=A->numInput;
332 /* is coloring available ? */
333 if (A->coloring == NULL) {
334
335 A->coloring=MEMALLOC(n,index_t);
336 if ( ! Esys_checkPtr(A->coloring)) {
337 Paso_Pattern_color(A,&(A->numColors),A->coloring);
338 if (! Esys_noError()) {
339 MEMFREE(A->coloring);
340 }
341 }
342 }
343 return A->coloring;
344 }
345 dim_t Paso_Pattern_maxDeg(Paso_Pattern* A)
346 {
347 dim_t deg=0, loc_deg=0, i;
348 const dim_t n=A->numInput;
349
350 #pragma omp parallel private(i, loc_deg)
351 {
352 loc_deg=0;
353 #pragma omp for schedule(static)
354 for (i = 0; i < n; ++i) {
355 loc_deg=MAX(loc_deg, A->ptr[i+1]-A->ptr[i]);
356 }
357 #pragma omp critical
358 {
359 deg=MAX(deg, loc_deg);
360 }
361 }
362 return deg;
363 }

  ViewVC Help
Powered by ViewVC 1.1.26