finley/src/IndexList.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
*
*******************************************************/


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

/* Finley: Converting an element list into a matrix shape     */

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

#include "IndexList.h"

/* Translate from distributed/local array indices to global indices */

/**************************************************************/
/* inserts the contributions from the element matrices of elements
   into the row index col. If symmetric is set, only the upper
   triangle of the matrix is stored. */

void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,
                                       bool_t reduce_row_order, index_t* row_map,
                                       bool_t reduce_col_order, index_t* col_map) {
  /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
  index_t color, *id=NULL;
  dim_t e,kr,kc,NN_row,NN_col,i,icol,irow, NN, *row_node=NULL,*col_node=NULL;
  if (elements!=NULL) {
    NN=elements->numNodes;
    id=TMPMEMALLOC(NN, index_t);
    if (! Finley_checkPtr(id) ) {
       for (i=0;i<NN;i++) id[i]=i;
       if (reduce_col_order) {
          col_node=elements->ReferenceElement->Type->linearNodes;
          NN_col=elements->LinearReferenceElement->Type->numNodes;
       } else {
          col_node=id;
          NN_col=elements->ReferenceElement->Type->numNodes;
       }
       if (reduce_row_order) {
          row_node=elements->ReferenceElement->Type->linearNodes;
          NN_row=elements->LinearReferenceElement->Type->numNodes;
       } else {
          row_node=id;
          NN_row=elements->ReferenceElement->Type->numNodes;
       }
       for (color=elements->minColor;color<=elements->maxColor;color++) {
           #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
           for (e=0;e<elements->numElements;e++) {
               if (elements->Color[e]==color) {
                   for (kr=0;kr<NN_row;kr++) {
                     irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
                     for (kc=0;kc<NN_col;kc++) {
                          icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
                          Finley_IndexList_insertIndex(&(index_list[irow]),icol);
                     }
                   }
               }
           }
       }
       TMPMEMFREE(id);
    }
  }
  return;
}

void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
                                                 Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
{
  index_t color;
  dim_t e,kr,kc,icol,irow, NN;
  if (elements!=NULL) {
    NN=elements->numNodes;
    for (color=elements->minColor;color<=elements->maxColor;color++) {
           #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
           for (e=0;e<elements->numElements;e++) {
               if (elements->Color[e]==color) {
                   for (kr=0;kr<NN;kr++) {
                     irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
                     if ((firstRow<=irow) && (irow < lastRow)) {
                          irow-=firstRow;
                          for (kc=0;kc<NN;kc++) {
                              icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
                              Finley_IndexList_insertIndex(&(index_list[irow]),icol);
                          }
                      }
                  }
               }
           }
    }
  }
}
void Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
                                                              Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
{
  index_t color;
  dim_t e,kr,kc,icol,irow, NN,irow_loc;
  if (elements!=NULL) {
    NN=elements->numNodes;
    for (color=elements->minColor;color<=elements->maxColor;color++) {
           #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
           for (e=0;e<elements->numElements;e++) {
               if (elements->Color[e]==color) {
                   for (kr=0;kr<NN;kr++) {
                     irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
                     if ((firstRow<=irow) && (irow < lastRow)) {
                          irow_loc=irow-firstRow;
                          for (kc=0;kc<NN;kc++) {
                              icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
                              if (icol != irow) Finley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
                          }
                      }
                  }
               }
           }
    }
  }
}

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

void Finley_IndexList_insertIndex(Finley_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,Finley_IndexList);
        if (Finley_checkPtr(in->extension)) return;
        in->extension->n=0;
        in->extension->extension=NULL;
     }
     Finley_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 Finley_IndexList in */

dim_t Finley_IndexList_count(Finley_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+Finley_IndexList_count(in->extension, range_min,range_max);
  }
}

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

void Finley_IndexList_toArray(Finley_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++;
          }

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

/* deallocates the Finley_IndexList in by recursive calls */

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

/* creates a Paso_pattern from a range of indices */
Paso_Pattern* Finley_IndexList_createPattern(dim_t n0, dim_t n,Finley_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 (! Finley_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]=Finley_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 (! Finley_checkPtr(index)) {
              #pragma omp parallel for schedule(static)
              for(i=n0;i<n;++i) {
                  Finley_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 (! Finley_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	/* Finley: Converting an element list into a matrix shape */
18
19	/**************************************************************/
20
21	#include "IndexList.h"
22
23	/* Translate from distributed/local array indices to global indices */
24
25	/**************************************************************/
26	/* inserts the contributions from the element matrices of elements
27	into the row index col. If symmetric is set, only the upper
28	triangle of the matrix is stored. */
29
30	void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,
31	bool_t reduce_row_order, index_t* row_map,
32	bool_t reduce_col_order, index_t* col_map) {
33	/* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
34	index_t color, *id=NULL;
35	dim_t e,kr,kc,NN_row,NN_col,i,icol,irow, NN, row_node=NULL,col_node=NULL;
36	if (elements!=NULL) {
37	NN=elements->numNodes;
38	id=TMPMEMALLOC(NN, index_t);
39	if (! Finley_checkPtr(id) ) {
40	for (i=0;i<NN;i++) id[i]=i;
41	if (reduce_col_order) {
42	col_node=elements->ReferenceElement->Type->linearNodes;
43	NN_col=elements->LinearReferenceElement->Type->numNodes;
44	} else {
45	col_node=id;
46	NN_col=elements->ReferenceElement->Type->numNodes;
47	}
48	if (reduce_row_order) {
49	row_node=elements->ReferenceElement->Type->linearNodes;
50	NN_row=elements->LinearReferenceElement->Type->numNodes;
51	} else {
52	row_node=id;
53	NN_row=elements->ReferenceElement->Type->numNodes;
54	}
55	for (color=elements->minColor;color<=elements->maxColor;color++) {
56	#pragma omp for private(e,irow,kr,kc,icol) schedule(static)
57	for (e=0;e<elements->numElements;e++) {
58	if (elements->Color[e]==color) {
59	for (kr=0;kr<NN_row;kr++) {
60	irow=row_map[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
61	for (kc=0;kc<NN_col;kc++) {
62	icol=col_map[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
63	Finley_IndexList_insertIndex(&(index_list[irow]),icol);
64	}
65	}
66	}
67	}
68	}
69	TMPMEMFREE(id);
70	}
71	}
72	return;
73	}
74
75	void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
76	Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
77	{
78	index_t color;
79	dim_t e,kr,kc,icol,irow, NN;
80	if (elements!=NULL) {
81	NN=elements->numNodes;
82	for (color=elements->minColor;color<=elements->maxColor;color++) {
83	#pragma omp for private(e,irow,kr,kc,icol) schedule(static)
84	for (e=0;e<elements->numElements;e++) {
85	if (elements->Color[e]==color) {
86	for (kr=0;kr<NN;kr++) {
87	irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
88	if ((firstRow<=irow) && (irow < lastRow)) {
89	irow-=firstRow;
90	for (kc=0;kc<NN;kc++) {
91	icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
92	Finley_IndexList_insertIndex(&(index_list[irow]),icol);
93	}
94	}
95	}
96	}
97	}
98	}
99	}
100	}
101	void Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
102	Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
103	{
104	index_t color;
105	dim_t e,kr,kc,icol,irow, NN,irow_loc;
106	if (elements!=NULL) {
107	NN=elements->numNodes;
108	for (color=elements->minColor;color<=elements->maxColor;color++) {
109	#pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
110	for (e=0;e<elements->numElements;e++) {
111	if (elements->Color[e]==color) {
112	for (kr=0;kr<NN;kr++) {
113	irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
114	if ((firstRow<=irow) && (irow < lastRow)) {
115	irow_loc=irow-firstRow;
116	for (kc=0;kc<NN;kc++) {
117	icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
118	if (icol != irow) Finley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
119	}
120	}
121	}
122	}
123	}
124	}
125	}
126	}
127
128	/* inserts row index row into the Finley_IndexList in if it does not exist */
129
130	void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {
131	dim_t i;
132	/* is index in in? */
133	for (i=0;i<in->n;i++) {
134	if (in->index[i]==index) return;
135	}
136	/* index could not be found */
137	if (in->n==INDEXLIST_LENGTH) {
138	/* if in->index is full check the extension */
139	if (in->extension==NULL) {
140	in->extension=TMPMEMALLOC(1,Finley_IndexList);
141	if (Finley_checkPtr(in->extension)) return;
142	in->extension->n=0;
143	in->extension->extension=NULL;
144	}
145	Finley_IndexList_insertIndex(in->extension,index);
146	} else {
147	/* insert index into in->index*/
148	in->index[in->n]=index;
149	in->n++;
150	}
151	}
152
153	/* counts the number of row indices in the Finley_IndexList in */
154
155	dim_t Finley_IndexList_count(Finley_IndexList* in, index_t range_min,index_t range_max) {
156	dim_t i;
157	dim_t out=0;
158	register index_t itmp;
159	if (in==NULL) {
160	return 0;
161	} else {
162	for (i=0;i<in->n;i++) {
163	itmp=in->index[i];
164	if ((itmp>=range_min) && (range_max>itmp)) ++out;
165	}
166	return out+Finley_IndexList_count(in->extension, range_min,range_max);
167	}
168	}
169
170	/* count the number of row indices in the Finley_IndexList in */
171
172	void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
173	dim_t i, ptr;
174	register index_t itmp;
175	if (in!=NULL) {
176	ptr=0;
177	for (i=0;i<in->n;i++) {
178	itmp=in->index[i];
179	if ((itmp>=range_min) && (range_max>itmp)) {
180	array[ptr]=itmp+index_offset;
181	ptr++;
182	}
183
184	}
185	Finley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
186	}
187	}
188
189	/* deallocates the Finley_IndexList in by recursive calls */
190
191	void Finley_IndexList_free(Finley_IndexList* in) {
192	if (in!=NULL) {
193	Finley_IndexList_free(in->extension);
194	TMPMEMFREE(in);
195	}
196	}
197
198	/* creates a Paso_pattern from a range of indices */
199	Paso_Pattern* Finley_IndexList_createPattern(dim_t n0, dim_t n,Finley_IndexList* index_list,index_t range_min,index_t range_max,index_t index_offset)
200	{
201	dim_t *ptr=NULL;
202	register dim_t s,i,itmp;
203	index_t *index=NULL;
204	Paso_Pattern* out=NULL;
205
206	ptr=MEMALLOC(n+1-n0,index_t);
207	if (! Finley_checkPtr(ptr) ) {
208	/* get the number of connections per row */
209	#pragma omp parallel for schedule(static) private(i)
210	for(i=n0;i<n;++i) {
211	ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);
212	}
213	/* accumulate ptr */
214	s=0;
215	for(i=n0;i<n;++i) {
216	itmp=ptr[i-n0];
217	ptr[i-n0]=s;
218	s+=itmp;
219	}
220	ptr[n-n0]=s;
221	/* fill index */
222	index=MEMALLOC(ptr[n-n0],index_t);
223	if (! Finley_checkPtr(index)) {
224	#pragma omp parallel for schedule(static)
225	for(i=n0;i<n;++i) {
226	Finley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
227	}
228	out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,n-n0,range_max+index_offset,ptr,index);
229	}
230	}
231	if (! Finley_noError()) {
232	MEMFREE(ptr);
233	MEMFREE(index);
234	Paso_Pattern_free(out);
235	}
236	return out;
237	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision