paso/src/Pattern.c


/*******************************************************
*
* Copyright (c) 2003-2008 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


/**************************************************************/

/* Paso: Pattern */

/**************************************************************/
 
/* Author: gross@access.edu.au */

/**************************************************************/

#include "Paso.h"
#include "Pattern.h"

/**************************************************************/

/* allocates a Pattern  */

Paso_Pattern* Paso_Pattern_alloc(int type, dim_t input_block_size, dim_t output_block_size, dim_t numOutput, dim_t numInput, index_t* ptr, index_t* index) {
  Paso_Pattern*out=NULL;
  index_t index_offset=(type & PATTERN_FORMAT_OFFSET1 ? 1:0);
  index_t loc_min_index,loc_max_index,min_index=index_offset,max_index=index_offset-1;
  dim_t i, sum=0;
  Paso_resetError();

  if (type & PATTERN_FORMAT_SYM) {
    Paso_setError(TYPE_ERROR,"Paso_Pattern_alloc: symmetric matrix pattern is not supported yet");
    return NULL;
  }
  if (ptr!=NULL && index != NULL) {
    #pragma omp parallel private(loc_min_index,loc_max_index,i)
     {
        loc_min_index=index_offset;
        loc_max_index=index_offset-1;
        if (type & PATTERN_FORMAT_OFFSET1) {
           #pragma omp for schedule(static) 
           for (i=0;i<numOutput;++i) {
               if (ptr[i]<ptr[i+1]) {
                 #ifdef USE_QSORTG
                    qsortG(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex); 
                 #else
                    qsort(&(index[ptr[i]-1]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex); 
                 #endif
                 loc_min_index=MIN(loc_min_index,index[ptr[i]-1]);
                 loc_max_index=MAX(loc_max_index,index[ptr[i+1]-2]);
               }
           }
        } else {
           #pragma omp for schedule(static) 
           for (i=0;i<numOutput;++i) {
               if (ptr[i]<ptr[i+1]) {
                 #ifdef USE_QSORTG
                    qsortG(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex); 
                 #else
                    qsort(&(index[ptr[i]]),(size_t)(ptr[i+1]-ptr[i]),sizeof(index_t),Paso_comparIndex); 
                 #endif
                 loc_min_index=MIN(loc_min_index,index[ptr[i]]);
                 loc_max_index=MAX(loc_max_index,index[ptr[i+1]-1]);
               }
           }
        }
        #pragma omp critical
        {
           min_index=MIN(loc_min_index,min_index);
           max_index=MAX(loc_max_index,max_index);
        }
    }
    if ( (min_index<index_offset) || (max_index>=numInput+index_offset) ) {
      Paso_setError(TYPE_ERROR,"Paso_Pattern_alloc: Pattern index out of range.");
      return NULL;
    }
  }
  out=MEMALLOC(1,Paso_Pattern);
  if (! Paso_checkPtr(out)) {
      out->type=type;
      out->reference_counter=1;
      out->numOutput=numOutput;
      out->numInput=numInput;
      out->ptr=ptr;
      out->index=index;
      out->input_block_size=input_block_size;
      out->output_block_size=output_block_size;
      out->block_size=out->input_block_size * out->output_block_size;
      if (out->ptr == NULL) {
          out->len=0;
      } else {
          out->len=out->ptr[out->numOutput] - index_offset;
      }
  }
  #ifdef Paso_TRACE
  printf("Paso_Pattern_alloc: system matrix pattern as been allocated.\n");
  #endif
  return out;
}

/* returns a reference to in */

Paso_Pattern* Paso_Pattern_getReference(Paso_Pattern* in) {
     if (in!=NULL) {
        ++(in->reference_counter);
     }
     return in;
}
  
/* deallocates a Pattern: */

void Paso_Pattern_free(Paso_Pattern* in) {
  if (in!=NULL) {
     in->reference_counter--;
     if (in->reference_counter<=0) {
        MEMFREE(in->ptr);
        MEMFREE(in->index);
        MEMFREE(in);
        #ifdef Paso_TRACE
        printf("Paso_Pattern_free: pattern as been deallocated.\n");
        #endif
     }
   }
}
/* *************************************************************/

/*  some routines which help to get the matrix pattern from elements: */

/*  this routine is used by qsort called in Paso_Pattern_alloc */

int Paso_comparIndex(const void *index1,const void *index2){
   index_t Iindex1,Iindex2;
   Iindex1=*(index_t*)index1;
   Iindex2=*(index_t*)index2;
   if (Iindex1<Iindex2) {
      return -1;
   } else {
      if (Iindex1>Iindex2) {
         return 1;
      } else {
         return 0;
      }
   }
}

bool_t Paso_Pattern_isEmpty(Paso_Pattern* in) {
     if (in != NULL) {
         if ((in->ptr != NULL) && (in->index != NULL)) return FALSE;
     }
     return TRUE;
}


/* computes the pattern coming from matrix-matrix multiplication
*
**/

Paso_Pattern* Paso_Pattern_multiply(int type, Paso_Pattern* A, Paso_Pattern* B) {
  Paso_Pattern*out=NULL;
  index_t index_offset=(type & PATTERN_FORMAT_OFFSET1 ? 1:0);
  index_t iptrA,iptrB;
  dim_t i,j,k;
  Paso_IndexList* index_list=NULL;

  index_list=TMPMEMALLOC(A->numOutput,Paso_IndexList);
  if (! Paso_checkPtr(index_list)) {
  
      #pragma omp parallel private(i)
      {
        #pragma omp for schedule(static)
        for(i=0;i<A->numOutput;++i) {
             index_list[i].extension=NULL;
             index_list[i].n=0;
        }
      }
  }
  
  for(i = 0; i < A->numOutput; i++) {
     for(iptrA = A->ptr[i]; iptrA < A->ptr[i+1]; ++iptrA) {
      j = A->index[iptrA];
      for(iptrB = B->ptr[j]; iptrB < B->ptr[j+1]; ++iptrB) {
        k = B->index[iptrB];
        if(i==k) {
          Finley_IndexList_insertIndex(&(index_list[i]),k);
          break;
        }
     }
    }
  }
    
  out=Paso_IndexList_createPattern(0, A->numOutput,index_list,0,INDEXLIST_LENGTH,0);

  #ifdef Paso_TRACE
  printf("Paso_Pattern_multipy: new pattern has been allocated.\n");
  #endif

 /* clean up */
   if (index_list!=NULL) {
        #pragma omp parallel for private(i) 
        for(i=0;i<A->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
     }
  TMPMEMFREE(index_list);
  
return out;
}


/*
 * Computes the pattern  of C = A binary operation B for CSR matrices A,B
 *
 * Note: we do not check whether A_ij(op)B_ij=0
 *
 */
Paso_Pattern* Paso_Pattern_binop(int type, Paso_Pattern* A, Paso_Pattern* B) {
  Paso_Pattern*out=NULL;
  index_t index_offset=(type & PATTERN_FORMAT_OFFSET1 ? 1:0);
  index_t iptrA,iptrB,*A_row=NULL,*B_row=NULL;
  dim_t i,j,k;

  Paso_IndexList* index_list=NULL;

 index_list=TMPMEMALLOC(A->numOutput,Paso_IndexList);
   if (! Paso_checkPtr(index_list)) {
  
      #pragma omp parallel private(i)
      {
        #pragma omp for schedule(static)
        for(i=0;i<A->numOutput;++i) {
             index_list[i].extension=NULL;
             index_list[i].n=0;
        }
      }
  }
  for(i = 0; i < A->numOutput; i++){
    iptrA = A->ptr[i],
    iptrB = B->ptr[i];
    
    while (iptrA < A->ptr[i+1] && iptrB < B->ptr[i+1]) {
        j = A->index[iptrA];
        k = B->index[iptrB];
        if (j<k) {
           Finley_IndexList_insertIndex(&(index_list[i]),j);
           iptrA++;
        } else if (j>k) {
            Finley_IndexList_insertIndex(&(index_list[i]),k);
            iptrB++;
        } else if (j==k) {
            Finley_IndexList_insertIndex(&(index_list[i]),j);
            iptrB++;
            iptrA++;
        }
    }
    while(iptrA < A->ptr[i+1]) {
        j = A->index[iptrA];
        Finley_IndexList_insertIndex(&(index_list[i]),j);
        iptrA++;
    }
    while(iptrB < B->ptr[i+1]) {
        k = B->index[iptrB];
        Finley_IndexList_insertIndex(&(index_list[i]),k);
        iptrB++;
    }
  }
 
  out=Paso_IndexList_createPattern(0, A->numOutput,index_list,0,INDEXLIST_LENGTH,0);

  #ifdef Paso_TRACE
  printf("Paso_Pattern_binop: new pattern has been allocated.\n");
  #endif

 /* clean up */
   if (index_list!=NULL) {
        #pragma omp parallel for private(i) 
        for(i=0;i<A->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
     }
  TMPMEMFREE(index_list);

  return out;
}

/* inserts row index row into the Paso_IndexList in if it does not exist */

void Paso_IndexList_insertIndex(Paso_IndexList* in, index_t index) {
  dim_t i;
  /* is index in in? */
  for (i=0;i<in->n;i++) {
    if (in->index[i]==index)  return;
  }
  /* index could not be found */
  if (in->n==INDEXLIST_LENGTH) {
     /* if in->index is full check the extension */
     if (in->extension==NULL) {
        in->extension=TMPMEMALLOC(1,Paso_IndexList);
        if (Paso_checkPtr(in->extension)) return;
        in->extension->n=0;
        in->extension->extension=NULL;
     }
     Paso_IndexList_insertIndex(in->extension,index);
  } else {
     /* insert index into in->index*/
     in->index[in->n]=index;
     in->n++;
  }
}

/* counts the number of row indices in the Paso_IndexList in */

dim_t Paso_IndexList_count(Paso_IndexList* in, index_t range_min,index_t range_max) {
  dim_t i;
  dim_t out=0;
  register index_t itmp;
  if (in==NULL) {
     return 0;
  } else {
    for (i=0;i<in->n;i++) {
          itmp=in->index[i];
          if ((itmp>=range_min) && (range_max>itmp)) ++out;
    }
     return out+Paso_IndexList_count(in->extension, range_min,range_max);
  }
}

/* count the number of row indices in the Paso_IndexList in */

void Paso_IndexList_toArray(Paso_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
  dim_t i, ptr;
  register index_t itmp;
  if (in!=NULL) {
    ptr=0;
    for (i=0;i<in->n;i++) {
          itmp=in->index[i];
          if ((itmp>=range_min) && (range_max>itmp)) {
             array[ptr]=itmp+index_offset;
             ptr++;
          }

    }
    Paso_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
  }
}

/* deallocates the Paso_IndexList in by recursive calls */

void Paso_IndexList_free(Paso_IndexList* in) {
  if (in!=NULL) {
    Paso_IndexList_free(in->extension);
    TMPMEMFREE(in);
  }
}

/* creates a Paso_pattern from a range of indices */
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)
{
   dim_t *ptr=NULL;
   register dim_t s,i,itmp;
   index_t *index=NULL;
   Paso_Pattern* out=NULL;

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