/[escript]/branches/doubleplusgood/paso/src/Pattern.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/Pattern.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4291 - (show annotations)
Thu Mar 7 08:57:58 2013 UTC (6 years, 1 month ago) by jfenwick
File size: 8091 byte(s)


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17
18 /************************************************************************************/
19
20 /* Paso: Pattern */
21
22 /************************************************************************************/
23
24 /* Author: Lutz Gross, l.gross@uq.edu.au */
25
26 /************************************************************************************/
27
28 #include "Paso.h"
29 #include "Pattern.h"
30
31 /************************************************************************************/
32
33 /* allocates a Pattern */
34
35 Paso_Pattern* Paso_Pattern_alloc(int type, dim_t numOutput, dim_t numInput, index_t* ptr, index_t* index) {
36 Paso_Pattern *out=NULL;
37 index_t index_offset=(type & MATRIX_FORMAT_OFFSET1 ? 1:0);
38 index_t loc_min_index,loc_max_index,min_index=index_offset,max_index=index_offset-1;
39 dim_t i;
40 Esys_resetError();
41
42 if (ptr!=NULL && index != NULL) {
43 #pragma omp parallel private(loc_min_index,loc_max_index,i)
44 {
45 loc_min_index=index_offset;
46 loc_max_index=index_offset-1;
47 if (type & MATRIX_FORMAT_OFFSET1) {
48 #pragma omp for schedule(static)
49 for (i=0;i<numOutput;++i) {
50 if (ptr[i]<ptr[i+1]) {
51 #ifdef USE_QSORTG
52 qsortG(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
53 #else
54 qsort(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
55 #endif
56 loc_min_index=MIN(loc_min_index,index[ptr[i]-1]);
57 loc_max_index=MAX(loc_max_index,index[ptr[i+1]-2]);
58 }
59 }
60 } else {
61 #pragma omp for schedule(static)
62 for (i=0;i<numOutput;++i) {
63 if (ptr[i]<ptr[i+1]) {
64 #ifdef USE_QSORTG
65 qsortG(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
66 #else
67 qsort(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex);
68 #endif
69 loc_min_index=MIN(loc_min_index,index[ptr[i]]);
70 loc_max_index=MAX(loc_max_index,index[ptr[i+1]-1]);
71 }
72 }
73 }
74 #pragma omp critical
75 {
76 min_index=MIN(loc_min_index,min_index);
77 max_index=MAX(loc_max_index,max_index);
78 }
79 }
80 if ( (min_index<index_offset) || (max_index>=numInput+index_offset) ) {
81 Esys_setError(TYPE_ERROR,"Paso_Pattern_alloc: Pattern index out of range.");
82 return NULL;
83 }
84 }
85 out=new Paso_Pattern;
86 if (! Esys_checkPtr(out)) {
87 out->type=type;
88 out->reference_counter=1;
89 out->numOutput=numOutput;
90 out->numInput=numInput;
91 out->ptr=ptr;
92 out->index=index;
93 out->main_iptr = NULL;
94 out->coloring = NULL;
95 out->numColors=-1;
96
97 if (out->ptr == NULL) {
98 out->len=0;
99 } else {
100 out->len=out->ptr[out->numOutput] - index_offset;
101 }
102 }
103 return out;
104 }
105
106 /* returns a reference to in */
107
108 Paso_Pattern* Paso_Pattern_getReference(Paso_Pattern* in) {
109 if (in!=NULL) {
110 ++(in->reference_counter);
111 }
112 return in;
113 }
114
115 /* deallocates a Pattern: */
116
117 void Paso_Pattern_free(Paso_Pattern* in) {
118 if (in!=NULL) {
119 in->reference_counter--;
120 if (in->reference_counter<=0) {
121 delete[] in->ptr;
122 delete[] in->index;
123 delete[] in->main_iptr;
124 delete[] in->coloring;
125 delete in;
126 }
127 }
128 }
129 /* ***********************************************************************************/
130
131 /* some routines which help to get the matrix pattern from elements: */
132
133 /* this routine is used by qsort called in Paso_Pattern_alloc */
134
135 int Paso_comparIndex(const void *index1,const void *index2){
136 index_t Iindex1,Iindex2;
137 Iindex1=*(index_t*)index1;
138 Iindex2=*(index_t*)index2;
139 if (Iindex1<Iindex2) {
140 return -1;
141 } else {
142 if (Iindex1>Iindex2) {
143 return 1;
144 } else {
145 return 0;
146 }
147 }
148 }
149
150 bool_t Paso_Pattern_isEmpty(Paso_Pattern* in) {
151 if (in != NULL) {
152 if ((in->ptr != NULL) && (in->index != NULL)) return FALSE;
153 }
154 return TRUE;
155 }
156
157
158 /* creates a Paso_pattern from a range of indices */
159 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)
160 {
161 const dim_t n=index_list_array->n;
162 Paso_IndexList* index_list = index_list_array->index_list;
163 dim_t *ptr=NULL;
164 register dim_t s,i,itmp;
165 index_t *index=NULL;
166 Paso_Pattern* out=NULL;
167
168 ptr=new index_t[n+1-n0];
169 if (! Esys_checkPtr(ptr) ) {
170 /* get the number of connections per row */
171 #pragma omp parallel for private(i) schedule(static)
172 for(i=n0;i<n;++i) {
173 ptr[i-n0]=Paso_IndexList_count(&index_list[i],range_min,range_max);
174 }
175 /* accumulate ptr */
176 s=0;
177 for(i=n0;i<n;++i) {
178 itmp=ptr[i-n0];
179 ptr[i-n0]=s;
180 s+=itmp;
181 }
182 ptr[n-n0]=s;
183 /* fill index */
184 index=new index_t[ptr[n-n0]];
185 if (! Esys_checkPtr(index)) {
186 #pragma omp parallel for private(i) schedule(static)
187 for(i=n0;i<n;++i) {
188 Paso_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
189 }
190 out=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
191 }
192 }
193 if (! Esys_noError()) {
194 delete[] ptr;
195 delete[] index;
196 Paso_Pattern_free(out);
197 }
198 return out;
199 }
200
201 index_t* Paso_Pattern_borrowMainDiagonalPointer(Paso_Pattern* A)
202 {
203 const dim_t n=A->numOutput;
204 int fail=0;
205 index_t *index,*where_p, i;
206
207 if (A->main_iptr == NULL) {
208 A->main_iptr=new index_t[n];
209 if (! Esys_checkPtr(A->main_iptr) ) {
210 #pragma omp parallel
211 {
212 /* identify the main diagonals */
213 #pragma omp for schedule(static) private(i, index, where_p)
214 for (i = 0; i < n; ++i) {
215 index=&(A->index[A->ptr[i]]);
216 where_p=reinterpret_cast<index_t*>(bsearch(&i,
217 index,
218 (size_t) (A->ptr[i + 1]-A->ptr[i]),
219 sizeof(index_t),
220 Paso_comparIndex));
221
222 if (where_p==NULL) {
223 fail=1;
224 } else {
225 A->main_iptr[i]=A->ptr[i]+(index_t)(where_p-index);
226 }
227 }
228
229 }
230 if (fail > 0) {
231 delete[] A->main_iptr;
232 A->main_iptr=NULL;
233 }
234
235 }
236 }
237 return A->main_iptr;
238 }
239
240
241 dim_t Paso_Pattern_getNumColors(Paso_Pattern* A)
242 {
243 Paso_Pattern_borrowColoringPointer(A); /* make sure numColors is defined */
244 return A->numColors;
245 }
246 index_t* Paso_Pattern_borrowColoringPointer(Paso_Pattern* A)
247 {
248 dim_t n=A->numInput;
249 /* is coloring available ? */
250 if (A->coloring == NULL) {
251
252 A->coloring=new index_t[n];
253 if ( ! Esys_checkPtr(A->coloring)) {
254 Paso_Pattern_color(A,&(A->numColors),A->coloring);
255 if (! Esys_noError()) {
256 delete[] A->coloring;
257 }
258 }
259 }
260 return A->coloring;
261 }
262 dim_t Paso_Pattern_maxDeg(Paso_Pattern* A)
263 {
264 dim_t deg=0, loc_deg=0, i;
265 const dim_t n=A->numInput;
266
267 #pragma omp parallel private(i, loc_deg)
268 {
269 loc_deg=0;
270 #pragma omp for schedule(static)
271 for (i = 0; i < n; ++i) {
272 loc_deg=MAX(loc_deg, A->ptr[i+1]-A->ptr[i]);
273 }
274 #pragma omp critical
275 {
276 deg=MAX(deg, loc_deg);
277 }
278 }
279 return deg;
280 }
281

  ViewVC Help
Powered by ViewVC 1.1.26