/[escript]/branches/domexper/paso/src/Pattern.c
ViewVC logotype

Contents of /branches/domexper/paso/src/Pattern.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3234 - (show annotations)
Mon Oct 4 01:46:30 2010 UTC (9 years, 2 months ago) by jfenwick
File MIME type: text/plain
File size: 12604 byte(s)
Some subdirs need to have changes pulled over but all of the unit tests 
except for modellib appear to work

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 (type & PATTERN_FORMAT_SYM) {
41 Esys_setError(TYPE_ERROR,"Paso_Pattern_alloc: symmetric matrix pattern is not supported yet");
42 return NULL;
43 }
44 if (ptr!=NULL && index != NULL) {
45 #pragma omp parallel private(loc_min_index,loc_max_index,i)
46 {
47 loc_min_index=index_offset;
48 loc_max_index=index_offset-1;
49 if (type & PATTERN_FORMAT_OFFSET1) {
50 #pragma omp for schedule(static)
51 for (i=0;i<numOutput;++i) {
52 if (ptr[i]<ptr[i+1]) {
53 #ifdef USE_QSORTG
54 qsortG(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
55 #else
56 qsort(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
57 #endif
58 loc_min_index=MIN(loc_min_index,index[ptr[i]-1]);
59 loc_max_index=MAX(loc_max_index,index[ptr[i+1]-2]);
60 }
61 }
62 } else {
63 #pragma omp for schedule(static)
64 for (i=0;i<numOutput;++i) {
65 if (ptr[i]<ptr[i+1]) {
66 #ifdef USE_QSORTG
67 qsortG(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
68 #else
69 qsort(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
70 #endif
71 loc_min_index=MIN(loc_min_index,index[ptr[i]]);
72 loc_max_index=MAX(loc_max_index,index[ptr[i+1]-1]);
73 }
74 }
75 }
76 #pragma omp critical
77 {
78 min_index=MIN(loc_min_index,min_index);
79 max_index=MAX(loc_max_index,max_index);
80 }
81 }
82 if ( (min_index<index_offset) || (max_index>=numInput+index_offset) ) {
83 Esys_setError(TYPE_ERROR,"Paso_Pattern_alloc: Pattern index out of range.");
84 return NULL;
85 }
86 }
87 out=MEMALLOC(1,Paso_Pattern);
88 if (! Esys_checkPtr(out)) {
89 out->type=type;
90 out->reference_counter=1;
91 out->numOutput=numOutput;
92 out->numInput=numInput;
93 out->ptr=ptr;
94 out->index=index;
95 out->main_iptr = NULL;
96 out->coloring = NULL;
97 out->numColors=-1;
98
99 if (out->ptr == NULL) {
100 out->len=0;
101 } else {
102 out->len=out->ptr[out->numOutput] - index_offset;
103 }
104 }
105 return out;
106 }
107
108 /* returns a reference to in */
109
110 Paso_Pattern* Paso_Pattern_getReference(Paso_Pattern* in) {
111 if (in!=NULL) {
112 ++(in->reference_counter);
113 }
114 return in;
115 }
116
117 /* deallocates a Pattern: */
118
119 void Paso_Pattern_free(Paso_Pattern* in) {
120 if (in!=NULL) {
121 in->reference_counter--;
122 if (in->reference_counter<=0) {
123 MEMFREE(in->ptr);
124 MEMFREE(in->index);
125 MEMFREE(in->main_iptr);
126 MEMFREE(in->coloring);
127 MEMFREE(in);
128 }
129 }
130 }
131 /* *************************************************************/
132
133 /* some routines which help to get the matrix pattern from elements: */
134
135 /* this routine is used by qsort called in Paso_Pattern_alloc */
136
137 int Paso_comparIndex(const void *index1,const void *index2){
138 index_t Iindex1,Iindex2;
139 Iindex1=*(index_t*)index1;
140 Iindex2=*(index_t*)index2;
141 if (Iindex1<Iindex2) {
142 return -1;
143 } else {
144 if (Iindex1>Iindex2) {
145 return 1;
146 } else {
147 return 0;
148 }
149 }
150 }
151
152 bool_t Paso_Pattern_isEmpty(Paso_Pattern* in) {
153 if (in != NULL) {
154 if ((in->ptr != NULL) && (in->index != NULL)) return FALSE;
155 }
156 return TRUE;
157 }
158
159
160 /* computes the pattern coming from matrix-matrix multiplication
161 *
162 **/
163
164 Paso_Pattern* Paso_Pattern_multiply(int type, Paso_Pattern* A, Paso_Pattern* B) {
165 Paso_Pattern*out=NULL;
166 index_t iptrA,iptrB;
167 dim_t i,j,k;
168 Paso_IndexList* index_list=NULL;
169
170 index_list=TMPMEMALLOC(A->numOutput,Paso_IndexList);
171 if (! Esys_checkPtr(index_list)) {
172 #pragma omp parallel for private(i) schedule(static)
173 for(i=0;i<A->numOutput;++i) {
174 index_list[i].extension=NULL;
175 index_list[i].n=0;
176 }
177 }
178
179 #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)
180 for(i = 0; i < A->numOutput; i++) {
181 for(iptrA = A->ptr[i]; iptrA < A->ptr[i+1]; ++iptrA) {
182 j = A->index[iptrA];
183 for(iptrB = B->ptr[j]; iptrB < B->ptr[j+1]; ++iptrB) {
184 k = B->index[iptrB];
185 Paso_IndexList_insertIndex(&(index_list[i]),k);
186 }
187 }
188 }
189
190 out=Paso_IndexList_createPattern(0, A->numOutput,index_list,0,B->numInput,0);
191
192 /* clean up */
193 if (index_list!=NULL) {
194 #pragma omp parallel for private(i) schedule(static)
195 for(i=0;i<A->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
196 }
197 TMPMEMFREE(index_list);
198
199 return out;
200 }
201
202
203
204 /*
205 * Computes the pattern of C = A binary operation B for CSR matrices A,B
206 *
207 * Note: we do not check whether A_ij(op)B_ij=0
208 *
209 */
210 Paso_Pattern* Paso_Pattern_binop(int type, Paso_Pattern* A, Paso_Pattern* B) {
211 Paso_Pattern *out=NULL;
212 index_t iptrA,iptrB;
213 dim_t i,j,k;
214
215 Paso_IndexList* index_list=NULL;
216
217 index_list=TMPMEMALLOC(A->numOutput,Paso_IndexList);
218 if (! Esys_checkPtr(index_list)) {
219 #pragma omp parallel for private(i) schedule(static)
220 for(i=0;i<A->numOutput;++i) {
221 index_list[i].extension=NULL;
222 index_list[i].n=0;
223 }
224 }
225
226 #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)
227 for(i = 0; i < B->numOutput; i++){
228 iptrA = A->ptr[i],
229 iptrB = B->ptr[i];
230
231 while (iptrA < A->ptr[i+1] && iptrB < B->ptr[i+1]) {
232 j = A->index[iptrA];
233 k = B->index[iptrB];
234 if (j<k) {
235 Paso_IndexList_insertIndex(&(index_list[i]),j);
236 iptrA++;
237 } else if (j>k) {
238 Paso_IndexList_insertIndex(&(index_list[i]),k);
239 iptrB++;
240 } else if (j==k) {
241 Paso_IndexList_insertIndex(&(index_list[i]),j);
242 iptrB++;
243 iptrA++;
244 }
245 }
246 while(iptrA < A->ptr[i+1]) {
247 j = A->index[iptrA];
248 Paso_IndexList_insertIndex(&(index_list[i]),j);
249 iptrA++;
250 }
251 while(iptrB < B->ptr[i+1]) {
252 k = B->index[iptrB];
253 Paso_IndexList_insertIndex(&(index_list[i]),k);
254 iptrB++;
255 }
256 }
257
258 out=Paso_IndexList_createPattern(0, A->numOutput,index_list,0,A->numInput,0);
259
260
261 /* clean up */
262 if (index_list!=NULL) {
263 #pragma omp parallel for private(i) schedule(static)
264 for(i=0;i<A->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
265 }
266 TMPMEMFREE(index_list);
267
268 return out;
269 }
270
271 /* inserts row index row into the Paso_IndexList in if it does not exist */
272
273 void Paso_IndexList_insertIndex(Paso_IndexList* in, index_t index) {
274 dim_t i;
275 /* is index in in? */
276 for (i=0;i<in->n;i++) {
277 if (in->index[i]==index) return;
278 }
279 /* index could not be found */
280 if (in->n==INDEXLIST_LENGTH) {
281 /* if in->index is full check the extension */
282 if (in->extension==NULL) {
283 in->extension=TMPMEMALLOC(1,Paso_IndexList);
284 if (Esys_checkPtr(in->extension)) return;
285 in->extension->n=0;
286 in->extension->extension=NULL;
287 }
288 Paso_IndexList_insertIndex(in->extension,index);
289 } else {
290 /* insert index into in->index*/
291 in->index[in->n]=index;
292 in->n++;
293 }
294 }
295
296 /* counts the number of row indices in the Paso_IndexList in */
297
298 dim_t Paso_IndexList_count(Paso_IndexList* in, index_t range_min,index_t range_max) {
299 dim_t i;
300 dim_t out=0;
301 register index_t itmp;
302 if (in==NULL) {
303 return 0;
304 } else {
305 for (i=0;i<in->n;i++) {
306 itmp=in->index[i];
307 if ((itmp>=range_min) && (range_max>itmp)) ++out;
308 }
309 return out+Paso_IndexList_count(in->extension, range_min,range_max);
310 }
311 }
312
313 /* count the number of row indices in the Paso_IndexList in */
314
315 void Paso_IndexList_toArray(Paso_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
316 dim_t i, ptr;
317 register index_t itmp;
318 if (in!=NULL) {
319 ptr=0;
320 for (i=0;i<in->n;i++) {
321 itmp=in->index[i];
322 if ((itmp>=range_min) && (range_max>itmp)) {
323 array[ptr]=itmp+index_offset;
324 ptr++;
325 }
326
327 }
328 Paso_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
329 }
330 }
331
332 /* deallocates the Paso_IndexList in by recursive calls */
333
334 void Paso_IndexList_free(Paso_IndexList* in) {
335 if (in!=NULL) {
336 Paso_IndexList_free(in->extension);
337 TMPMEMFREE(in);
338 }
339 }
340
341 /* creates a Paso_pattern from a range of indices */
342 Paso_Pattern* Paso_IndexList_createPattern(dim_t n0, dim_t n,Paso_IndexList* index_list,index_t range_min,index_t range_max,index_t index_offset)
343 {
344 dim_t *ptr=NULL;
345 register dim_t s,i,itmp;
346 index_t *index=NULL;
347 Paso_Pattern* out=NULL;
348
349 ptr=MEMALLOC(n+1-n0,index_t);
350 if (! Esys_checkPtr(ptr) ) {
351 /* get the number of connections per row */
352 #pragma omp parallel for private(i) schedule(static)
353 for(i=n0;i<n;++i) {
354 ptr[i-n0]=Paso_IndexList_count(&index_list[i],range_min,range_max);
355 }
356 /* accumulate ptr */
357 s=0;
358 for(i=n0;i<n;++i) {
359 itmp=ptr[i-n0];
360 ptr[i-n0]=s;
361 s+=itmp;
362 }
363 ptr[n-n0]=s;
364 /* fill index */
365 index=MEMALLOC(ptr[n-n0],index_t);
366 if (! Esys_checkPtr(index)) {
367 #pragma omp parallel for private(i) schedule(static)
368 for(i=n0;i<n;++i) {
369 Paso_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
370 }
371 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
372 }
373 }
374 if (! Esys_noError()) {
375 MEMFREE(ptr);
376 MEMFREE(index);
377 Paso_Pattern_free(out);
378 }
379 return out;
380 }
381
382 index_t* Paso_Pattern_borrowMainDiagonalPointer(Paso_Pattern* A)
383 {
384 const dim_t n=A->numOutput;
385 int fail=0;
386 index_t i,iptr,iptr_main;
387
388 if (A->main_iptr == NULL) {
389 A->main_iptr=MEMALLOC(n,index_t);
390 if (! Esys_checkPtr(A->main_iptr) ) {
391 #pragma omp parallel
392 {
393 /* identify the main diagonals */
394 #pragma omp for schedule(static) private(i,iptr,iptr_main)
395 for (i = 0; i < n; ++i) {
396 iptr_main=A->ptr[0]-1;
397 for (iptr=A->ptr[i];iptr<A->ptr[i+1]; iptr++) {
398 if (A->index[iptr]==i) {
399 iptr_main=iptr;
400 break;
401 }
402 }
403 A->main_iptr[i]=iptr_main;
404 if (iptr_main==A->ptr[0]-1) fail=1;
405 }
406
407 }
408 if (fail > 0) {
409 MEMFREE(A->main_iptr);
410 A->main_iptr=NULL;
411 }
412
413 }
414 }
415 return A->main_iptr;
416 }
417
418 dim_t Paso_Pattern_getNumColors(Paso_Pattern* A)
419 {
420 Paso_Pattern_borrowColoringPointer(A); /* make sure numColors is defined */
421 return A->numColors;
422 }
423 index_t* Paso_Pattern_borrowColoringPointer(Paso_Pattern* A)
424 {
425 dim_t n=A->numInput;
426 /* is coloring available ? */
427 if (A->coloring == NULL) {
428
429 A->coloring=MEMALLOC(n,index_t);
430 if ( ! Esys_checkPtr(A->coloring)) {
431 Paso_Pattern_color(A,&(A->numColors),A->coloring);
432 if (! Esys_noError()) {
433 MEMFREE(A->coloring);
434 }
435 }
436 }
437 return A->coloring;
438 }

  ViewVC Help
Powered by ViewVC 1.1.26