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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3312 - (hide annotations)
Tue Oct 26 07:54:58 2010 UTC (8 years, 11 months ago) by gross
File MIME type: text/plain
File size: 7863 byte(s)
last step for a clean up version of the AMG
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 gross 3312 index_t index_offset=(type & MATRIX_FORMAT_OFFSET1 ? 1:0);
36 ksteube 1313 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 (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 gross 3312 if (type & MATRIX_FORMAT_OFFSET1) {
46 ksteube 1313 #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 gross 1736 if ( (min_index<index_offset) || (max_index>=numInput+index_offset) ) {
79 jfenwick 3259 Esys_setError(TYPE_ERROR,"Paso_Pattern_alloc: Pattern index out of range.");
80 ksteube 1313 return NULL;
81     }
82     }
83     out=MEMALLOC(1,Paso_Pattern);
84 jfenwick 3259 if (! Esys_checkPtr(out)) {
85 ksteube 1313 out->type=type;
86     out->reference_counter=1;
87     out->numOutput=numOutput;
88 gross 1736 out->numInput=numInput;
89 ksteube 1313 out->ptr=ptr;
90     out->index=index;
91 gross 3005 out->main_iptr = NULL;
92 gross 3094 out->coloring = NULL;
93     out->numColors=-1;
94 gross 2551
95 ksteube 1313 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 gross 3005 MEMFREE(in->main_iptr);
122 gross 3094 MEMFREE(in->coloring);
123 ksteube 1313 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 artak 1881
155    
156     /* creates a Paso_pattern from a range of indices */
157 gross 3283 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)
158 artak 1881 {
159 gross 3283 const dim_t n=index_list_array->n;
160     Paso_IndexList* index_list = index_list_array->index_list;
161 artak 1881 dim_t *ptr=NULL;
162     register dim_t s,i,itmp;
163     index_t *index=NULL;
164     Paso_Pattern* out=NULL;
165    
166     ptr=MEMALLOC(n+1-n0,index_t);
167 jfenwick 3259 if (! Esys_checkPtr(ptr) ) {
168 gross 3283 /* get the number of connections per row */
169 artak 2107 #pragma omp parallel for private(i) schedule(static)
170 artak 1881 for(i=n0;i<n;++i) {
171     ptr[i-n0]=Paso_IndexList_count(&index_list[i],range_min,range_max);
172     }
173     /* accumulate ptr */
174     s=0;
175     for(i=n0;i<n;++i) {
176     itmp=ptr[i-n0];
177     ptr[i-n0]=s;
178     s+=itmp;
179     }
180     ptr[n-n0]=s;
181     /* fill index */
182     index=MEMALLOC(ptr[n-n0],index_t);
183 jfenwick 3259 if (! Esys_checkPtr(index)) {
184 artak 2107 #pragma omp parallel for private(i) schedule(static)
185 artak 1881 for(i=n0;i<n;++i) {
186     Paso_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
187     }
188 gross 3312 out=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
189 artak 1881 }
190     }
191 jfenwick 3259 if (! Esys_noError()) {
192 artak 1881 MEMFREE(ptr);
193     MEMFREE(index);
194     Paso_Pattern_free(out);
195     }
196     return out;
197     }
198 gross 3005
199     index_t* Paso_Pattern_borrowMainDiagonalPointer(Paso_Pattern* A)
200     {
201     const dim_t n=A->numOutput;
202     int fail=0;
203 gross 3283 index_t *index,*where_p, i;
204 gross 3005
205     if (A->main_iptr == NULL) {
206     A->main_iptr=MEMALLOC(n,index_t);
207 jfenwick 3259 if (! Esys_checkPtr(A->main_iptr) ) {
208 gross 3005 #pragma omp parallel
209     {
210     /* identify the main diagonals */
211 gross 3283 #pragma omp for schedule(static) private(i, index, where_p)
212 gross 3005 for (i = 0; i < n; ++i) {
213 gross 3283 index=&(A->index[A->ptr[i]]);
214     where_p=bsearch(&i,
215     index,
216     (size_t) (A->ptr[i + 1]-A->ptr[i]),
217     sizeof(index_t),
218     Paso_comparIndex);
219    
220     if (where_p==NULL) {
221     fail=1;
222     } else {
223     A->main_iptr[i]=A->ptr[i]+(index_t)(where_p-index);
224     }
225     }
226 gross 3005
227     }
228 gross 3009 if (fail > 0) {
229     MEMFREE(A->main_iptr);
230     A->main_iptr=NULL;
231 gross 3005 }
232 gross 3009
233 gross 3005 }
234     }
235     return A->main_iptr;
236 gross 3009 }
237 gross 3283
238 gross 3094
239     dim_t Paso_Pattern_getNumColors(Paso_Pattern* A)
240     {
241     Paso_Pattern_borrowColoringPointer(A); /* make sure numColors is defined */
242     return A->numColors;
243     }
244     index_t* Paso_Pattern_borrowColoringPointer(Paso_Pattern* A)
245     {
246     dim_t n=A->numInput;
247     /* is coloring available ? */
248     if (A->coloring == NULL) {
249    
250     A->coloring=MEMALLOC(n,index_t);
251 jfenwick 3259 if ( ! Esys_checkPtr(A->coloring)) {
252 gross 3094 Paso_Pattern_color(A,&(A->numColors),A->coloring);
253 jfenwick 3259 if (! Esys_noError()) {
254 gross 3094 MEMFREE(A->coloring);
255     }
256     }
257     }
258     return A->coloring;
259 gross 3283 }
260     dim_t Paso_Pattern_maxDeg(Paso_Pattern* A)
261     {
262     dim_t deg=0, loc_deg=0, i;
263     const dim_t n=A->numInput;
264    
265     #pragma omp parallel private(i, loc_deg)
266     {
267     loc_deg=0;
268     #pragma omp for schedule(static)
269     for (i = 0; i < n; ++i) {
270     loc_deg=MAX(loc_deg, A->ptr[i+1]-A->ptr[i]);
271     }
272     #pragma omp critical
273     {
274     deg=MAX(deg, loc_deg);
275     }
276     }
277     return deg;
278     }

  ViewVC Help
Powered by ViewVC 1.1.26