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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3303 - (show annotations)
Mon Oct 25 04:33:31 2010 UTC (8 years, 11 months ago) by gross
File MIME type: text/plain
File size: 7866 byte(s)
more clean up work on the AMG
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 /* creates a Paso_pattern from a range of indices */
157 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 {
159 const dim_t n=index_list_array->n;
160 Paso_IndexList* index_list = index_list_array->index_list;
161 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 if (! Esys_checkPtr(ptr) ) {
168 /* get the number of connections per row */
169 #pragma omp parallel for private(i) schedule(static)
170 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 if (! Esys_checkPtr(index)) {
184 #pragma omp parallel for private(i) schedule(static)
185 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 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
189 }
190 }
191 if (! Esys_noError()) {
192 MEMFREE(ptr);
193 MEMFREE(index);
194 Paso_Pattern_free(out);
195 }
196 return out;
197 }
198
199 index_t* Paso_Pattern_borrowMainDiagonalPointer(Paso_Pattern* A)
200 {
201 const dim_t n=A->numOutput;
202 int fail=0;
203 index_t *index,*where_p, i;
204
205 if (A->main_iptr == NULL) {
206 A->main_iptr=MEMALLOC(n,index_t);
207 if (! Esys_checkPtr(A->main_iptr) ) {
208 #pragma omp parallel
209 {
210 /* identify the main diagonals */
211 #pragma omp for schedule(static) private(i, index, where_p)
212 for (i = 0; i < n; ++i) {
213 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
227 }
228 if (fail > 0) {
229 MEMFREE(A->main_iptr);
230 A->main_iptr=NULL;
231 }
232
233 }
234 }
235 return A->main_iptr;
236 }
237
238
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 if ( ! Esys_checkPtr(A->coloring)) {
252 Paso_Pattern_color(A,&(A->numColors),A->coloring);
253 if (! Esys_noError()) {
254 MEMFREE(A->coloring);
255 }
256 }
257 }
258 return A->coloring;
259 }
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