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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (hide annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years ago) by jfenwick
File MIME type: text/plain
File size: 12604 byte(s)
Merging dudley and scons updates from branches

1 ksteube 1313
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1313
14 ksteube 1811
15 phornby 1916
16 ksteube 1313 /**************************************************************/
17    
18     /* Paso: Pattern */
19    
20     /**************************************************************/
21    
22 jfenwick 2608 /* Author: Lutz Gross, l.gross@uq.edu.au */
23 ksteube 1313
24     /**************************************************************/
25    
26     #include "Paso.h"
27     #include "Pattern.h"
28    
29     /**************************************************************/
30    
31     /* allocates a Pattern */
32    
33 gross 2551 Paso_Pattern* Paso_Pattern_alloc(int type, dim_t numOutput, dim_t numInput, index_t* ptr, index_t* index) {
34 artak 1976 Paso_Pattern *out=NULL;
35 ksteube 1313 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 artak 1976 dim_t i;
38 jfenwick 3259 Esys_resetError();
39 ksteube 1313
40     if (type & PATTERN_FORMAT_SYM) {
41 jfenwick 3259 Esys_setError(TYPE_ERROR,"Paso_Pattern_alloc: symmetric matrix pattern is not supported yet");
42 ksteube 1313 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 gross 1736 if ( (min_index<index_offset) || (max_index>=numInput+index_offset) ) {
83 jfenwick 3259 Esys_setError(TYPE_ERROR,"Paso_Pattern_alloc: Pattern index out of range.");
84 ksteube 1313 return NULL;
85     }
86     }
87     out=MEMALLOC(1,Paso_Pattern);
88 jfenwick 3259 if (! Esys_checkPtr(out)) {
89 ksteube 1313 out->type=type;
90     out->reference_counter=1;
91     out->numOutput=numOutput;
92 gross 1736 out->numInput=numInput;
93 ksteube 1313 out->ptr=ptr;
94     out->index=index;
95 gross 3005 out->main_iptr = NULL;
96 gross 3094 out->coloring = NULL;
97     out->numColors=-1;
98 gross 2551
99 ksteube 1313 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 gross 3005 MEMFREE(in->main_iptr);
126 gross 3094 MEMFREE(in->coloring);
127 ksteube 1313 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 artak 1881
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 jfenwick 3259 if (! Esys_checkPtr(index_list)) {
172 artak 2107 #pragma omp parallel for private(i) schedule(static)
173 artak 1881 for(i=0;i<A->numOutput;++i) {
174     index_list[i].extension=NULL;
175     index_list[i].n=0;
176     }
177 artak 2107 }
178 artak 1881
179 artak 2107 #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)
180 artak 1881 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 phornby 1923 Paso_IndexList_insertIndex(&(index_list[i]),k);
186 artak 1881 }
187     }
188     }
189    
190 artak 1931 out=Paso_IndexList_createPattern(0, A->numOutput,index_list,0,B->numInput,0);
191 artak 1881
192     /* clean up */
193     if (index_list!=NULL) {
194 artak 2107 #pragma omp parallel for private(i) schedule(static)
195 artak 1881 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 artak 1976 Paso_Pattern *out=NULL;
212     index_t iptrA,iptrB;
213 artak 1890 dim_t i,j,k;
214 artak 1881
215     Paso_IndexList* index_list=NULL;
216    
217     index_list=TMPMEMALLOC(A->numOutput,Paso_IndexList);
218 jfenwick 3259 if (! Esys_checkPtr(index_list)) {
219 artak 2107 #pragma omp parallel for private(i) schedule(static)
220 artak 1881 for(i=0;i<A->numOutput;++i) {
221     index_list[i].extension=NULL;
222     index_list[i].n=0;
223     }
224 artak 2107 }
225    
226     #pragma omp parallel for private(i,iptrA,j,iptrB,k) schedule(static)
227 artak 1902 for(i = 0; i < B->numOutput; i++){
228 artak 1881 iptrA = A->ptr[i],
229     iptrB = B->ptr[i];
230 artak 1890
231 artak 1881 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 phornby 1923 Paso_IndexList_insertIndex(&(index_list[i]),j);
236 artak 1881 iptrA++;
237     } else if (j>k) {
238 phornby 1923 Paso_IndexList_insertIndex(&(index_list[i]),k);
239 artak 1881 iptrB++;
240     } else if (j==k) {
241 phornby 1923 Paso_IndexList_insertIndex(&(index_list[i]),j);
242 artak 1881 iptrB++;
243     iptrA++;
244     }
245     }
246     while(iptrA < A->ptr[i+1]) {
247     j = A->index[iptrA];
248 phornby 1923 Paso_IndexList_insertIndex(&(index_list[i]),j);
249 artak 1881 iptrA++;
250     }
251     while(iptrB < B->ptr[i+1]) {
252     k = B->index[iptrB];
253 phornby 1923 Paso_IndexList_insertIndex(&(index_list[i]),k);
254 artak 1881 iptrB++;
255     }
256     }
257    
258 artak 1931 out=Paso_IndexList_createPattern(0, A->numOutput,index_list,0,A->numInput,0);
259 artak 1881
260    
261     /* clean up */
262     if (index_list!=NULL) {
263 artak 2107 #pragma omp parallel for private(i) schedule(static)
264 artak 1881 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 jfenwick 3259 if (Esys_checkPtr(in->extension)) return;
285 artak 1881 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 jfenwick 3259 if (! Esys_checkPtr(ptr) ) {
351 artak 1881 /* get the number of connections per row */
352 artak 2107 #pragma omp parallel for private(i) schedule(static)
353 artak 1881 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 jfenwick 3259 if (! Esys_checkPtr(index)) {
367 artak 2107 #pragma omp parallel for private(i) schedule(static)
368 artak 1881 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 gross 2551 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
372 artak 1881 }
373     }
374 jfenwick 3259 if (! Esys_noError()) {
375 artak 1881 MEMFREE(ptr);
376     MEMFREE(index);
377     Paso_Pattern_free(out);
378     }
379     return out;
380     }
381 gross 3005
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 jfenwick 3259 if (! Esys_checkPtr(A->main_iptr) ) {
391 gross 3005 #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 gross 3009 if (fail > 0) {
409     MEMFREE(A->main_iptr);
410     A->main_iptr=NULL;
411 gross 3005 }
412 gross 3009
413 gross 3005 }
414     }
415     return A->main_iptr;
416 gross 3009 }
417 gross 3094
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 jfenwick 3259 if ( ! Esys_checkPtr(A->coloring)) {
431 gross 3094 Paso_Pattern_color(A,&(A->numColors),A->coloring);
432 jfenwick 3259 if (! Esys_noError()) {
433 gross 3094 MEMFREE(A->coloring);
434     }
435     }
436     }
437     return A->coloring;
438     }

  ViewVC Help
Powered by ViewVC 1.1.26