finley/src/IndexList.c

/*
 ************************************************************
 *          Copyright 2006 by ACcESS MNRF                   *
 *                                                          *
 *              http://www.access.edu.au                    *
 *       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     */

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

/*  Author: gross@access.edu.au */
/*  Version: $Id$ */

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

#include "IndexList.h"

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

#ifdef PASO_MPI
int Finley_IndexList_localToGlobal(Finley_NodeDistribution *dofDistribution, int localIndex) {
  /*
    get global id of icol
    if icol is internal node (on this CPU): use icol+vtxdist[my_CPU]
    else use indexExternal[icol-numLocal] to get global index of node
    (actually DOF...the NodeDistribution structure should have been called DofDistribution)
  */
  index_t my_CPU=dofDistribution->MPIInfo->rank;
  if (localIndex < dofDistribution->numLocal) {
    localIndex = localIndex + dofDistribution->vtxdist[my_CPU];
  }
  else {
    localIndex = dofDistribution->indexExternal[localIndex-dofDistribution->numLocal];
  }
  return(localIndex);
}
#endif

/**************************************************************/
/* 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_Mesh* mesh, Finley_ElementFile* elements,
                                       bool_t reduce_row_order, index_t* row_Label,
                                       bool_t reduce_col_order, index_t* col_Label) {
  /* 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, num_CPUs = 1, my_CPU = 0;
  dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;
#ifdef PASO_MPI
    num_CPUs = mesh->MPIInfo->size;
    my_CPU = mesh->MPIInfo->rank;
#endif
  /* print_mesh_statistics( mesh, TRUE ); */

  if (elements!=NULL) {
    dim_t NN=elements->ReferenceElement->Type->numNodes;
    index_t id[NN],*row_node,*col_node;
    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;
    }
    if (num_CPUs == 1) {
    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_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
                  for (kc=0;kc<NN_col;kc++) {
                       icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
                       Finley_IndexList_insertIndex(&(index_list[irow]),icol);
                  }
                }
            }
        }
      }
    }
    else {  /* More than one CPU (what's below should also work for one CPU, but let's build confidence in it first) */
#ifdef PASO_MPI
    Finley_NodeDistribution *row_degreeOfFreedomDistribution;
    Finley_NodeDistribution *col_degreeOfFreedomDistribution;
    if (reduce_col_order) {
      col_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;
    }
    else {
      col_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;
    }
    if (reduce_row_order) {
      row_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;
    }
    else {
      row_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;
    }
    /* Not using loop over colors as above */ {
        #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
    for (e=0;e<elements->numElements;e++) {
                for (kr=0;kr<NN_row;kr++) {
                  irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
          if (irow < row_degreeOfFreedomDistribution->numLocal) {
                    for (kc=0;kc<NN_col;kc++) {
               /* Get the local col ID */
                       icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
               /* Convert to global col ID (row ID is saved as local value) */
               icol = Finley_IndexList_localToGlobal(col_degreeOfFreedomDistribution, icol);
                       Finley_IndexList_insertIndex(&(index_list[irow]),icol);
                    }
          }
                }
        }
      }
#endif
    }   /* More than one CPU */
  }
  return;
}

/* 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) {
  if (in==NULL) {
     return 0;
  } else {
     return (in->n)+Finley_IndexList_count(in->extension);
  }
}

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

void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array) {
  dim_t i;
  if (in!=NULL) {
    for (i=0;i<in->n;i++) array[i]=in->index[i];
    Finley_IndexList_toArray(in->extension,&(array[in->n]));
  }
}

/* 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);
  }
}

/*
 * $Log$
 * Revision 1.6  2005/09/15 03:44:22  jgs
 * Merge of development branch dev-02 back to main trunk on 2005-09-15
 *
 * Revision 1.5.2.1  2005/09/07 06:26:18  gross
 * the solver from finley are put into the standalone package paso now
 *
 * Revision 1.5  2005/07/08 04:07:51  jgs
 * Merge of development branch back to main trunk on 2005-07-08
 *
 * Revision 1.4  2004/12/15 07:08:32  jgs
 * *** empty log message ***
 * Revision 1.1.1.1.2.3  2005/06/29 02:34:50  gross
 * some changes towards 64 integers in finley
 *
 * Revision 1.1.1.1.2.2  2004/11/24 01:37:13  gross
 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
 *
 *
 *
 */
1	/*
2	************************************************************
3	* Copyright 2006 by ACcESS MNRF *
4	* *
5	* http://www.access.edu.au *
6	* Primary Business: Queensland, Australia *
7	* Licensed under the Open Software License version 3.0 *
8	* http://www.opensource.org/licenses/osl-3.0.php *
9	* *
10	************************************************************
11	*/
12
13	/**************************************************************/
14
15	/* Finley: Converting an element list into a matrix shape */
16
17	/**************************************************************/
18
19	/* Author: gross@access.edu.au */
20	/* Version: $Id$ */
21
22	/**************************************************************/
23
24	#include "IndexList.h"
25
26	/* Translate from distributed/local array indices to global indices */
27
28	#ifdef PASO_MPI
29	int Finley_IndexList_localToGlobal(Finley_NodeDistribution *dofDistribution, int localIndex) {
30	/*
31	get global id of icol
32	if icol is internal node (on this CPU): use icol+vtxdist[my_CPU]
33	else use indexExternal[icol-numLocal] to get global index of node
34	(actually DOF...the NodeDistribution structure should have been called DofDistribution)
35	*/
36	index_t my_CPU=dofDistribution->MPIInfo->rank;
37	if (localIndex < dofDistribution->numLocal) {
38	localIndex = localIndex + dofDistribution->vtxdist[my_CPU];
39	}
40	else {
41	localIndex = dofDistribution->indexExternal[localIndex-dofDistribution->numLocal];
42	}
43	return(localIndex);
44	}
45	#endif
46
47	/**************************************************************/
48	/* inserts the contributions from the element matrices of elements
49	into the row index col. If symmetric is set, only the upper
50	triangle of the matrix is stored. */
51
52	void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_Mesh* mesh, Finley_ElementFile* elements,
53	bool_t reduce_row_order, index_t* row_Label,
54	bool_t reduce_col_order, index_t* col_Label) {
55	/* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
56	index_t color, num_CPUs = 1, my_CPU = 0;
57	dim_t e,kr,kc,NN_row,NN_col,i,icol,irow;
58	#ifdef PASO_MPI
59	num_CPUs = mesh->MPIInfo->size;
60	my_CPU = mesh->MPIInfo->rank;
61	#endif
62	/* print_mesh_statistics( mesh, TRUE ); */
63
64	if (elements!=NULL) {
65	dim_t NN=elements->ReferenceElement->Type->numNodes;
66	index_t id[NN],row_node,col_node;
67	for (i=0;i<NN;i++) id[i]=i;
68	if (reduce_col_order) {
69	col_node=elements->ReferenceElement->Type->linearNodes;
70	NN_col=elements->LinearReferenceElement->Type->numNodes;
71	} else {
72	col_node=id;
73	NN_col=elements->ReferenceElement->Type->numNodes;
74	}
75	if (reduce_row_order) {
76	row_node=elements->ReferenceElement->Type->linearNodes;
77	NN_row=elements->LinearReferenceElement->Type->numNodes;
78	} else {
79	row_node=id;
80	NN_row=elements->ReferenceElement->Type->numNodes;
81	}
82	if (num_CPUs == 1) {
83	for (color=elements->minColor;color<=elements->maxColor;color++) {
84	#pragma omp for private(e,irow,kr,kc,icol) schedule(static)
85	for (e=0;e<elements->numElements;e++) {
86	if (elements->Color[e]==color) {
87	for (kr=0;kr<NN_row;kr++) {
88	irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
89	for (kc=0;kc<NN_col;kc++) {
90	icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
91	Finley_IndexList_insertIndex(&(index_list[irow]),icol);
92	}
93	}
94	}
95	}
96	}
97	}
98	else { /* More than one CPU (what's below should also work for one CPU, but let's build confidence in it first) */
99	#ifdef PASO_MPI
100	Finley_NodeDistribution *row_degreeOfFreedomDistribution;
101	Finley_NodeDistribution *col_degreeOfFreedomDistribution;
102	if (reduce_col_order) {
103	col_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;
104	}
105	else {
106	col_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;
107	}
108	if (reduce_row_order) {
109	row_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;
110	}
111	else {
112	row_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;
113	}
114	/* Not using loop over colors as above */ {
115	#pragma omp for private(e,irow,kr,kc,icol) schedule(static)
116	for (e=0;e<elements->numElements;e++) {
117	for (kr=0;kr<NN_row;kr++) {
118	irow=row_Label[elements->Nodes[INDEX2(row_node[kr],e,NN)]];
119	if (irow < row_degreeOfFreedomDistribution->numLocal) {
120	for (kc=0;kc<NN_col;kc++) {
121	/* Get the local col ID */
122	icol=col_Label[elements->Nodes[INDEX2(col_node[kc],e,NN)]];
123	/* Convert to global col ID (row ID is saved as local value) */
124	icol = Finley_IndexList_localToGlobal(col_degreeOfFreedomDistribution, icol);
125	Finley_IndexList_insertIndex(&(index_list[irow]),icol);
126	}
127	}
128	}
129	}
130	}
131	#endif
132	} /* More than one CPU */
133	}
134	return;
135	}
136
137	/* inserts row index row into the Finley_IndexList in if it does not exist */
138
139	void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {
140	dim_t i;
141	/* is index in in? */
142	for (i=0;i<in->n;i++) {
143	if (in->index[i]==index) return;
144	}
145	/* index could not be found */
146	if (in->n==INDEXLIST_LENGTH) {
147	/* if in->index is full check the extension */
148	if (in->extension==NULL) {
149	in->extension=TMPMEMALLOC(1,Finley_IndexList);
150	if (Finley_checkPtr(in->extension)) return;
151	in->extension->n=0;
152	in->extension->extension=NULL;
153	}
154	Finley_IndexList_insertIndex(in->extension,index);
155	} else {
156	/* insert index into in->index*/
157	in->index[in->n]=index;
158	in->n++;
159	}
160	}
161
162	/* counts the number of row indices in the Finley_IndexList in */
163
164	dim_t Finley_IndexList_count(Finley_IndexList* in) {
165	if (in==NULL) {
166	return 0;
167	} else {
168	return (in->n)+Finley_IndexList_count(in->extension);
169	}
170	}
171
172	/* count the number of row indices in the Finley_IndexList in */
173
174	void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array) {
175	dim_t i;
176	if (in!=NULL) {
177	for (i=0;i<in->n;i++) array[i]=in->index[i];
178	Finley_IndexList_toArray(in->extension,&(array[in->n]));
179	}
180	}
181
182	/* deallocates the Finley_IndexList in by recursive calls */
183
184	void Finley_IndexList_free(Finley_IndexList* in) {
185	if (in!=NULL) {
186	Finley_IndexList_free(in->extension);
187	TMPMEMFREE(in);
188	}
189	}
190
191	/*
192	* $Log$
193	* Revision 1.6 2005/09/15 03:44:22 jgs
194	* Merge of development branch dev-02 back to main trunk on 2005-09-15
195	*
196	* Revision 1.5.2.1 2005/09/07 06:26:18 gross
197	* the solver from finley are put into the standalone package paso now
198	*
199	* Revision 1.5 2005/07/08 04:07:51 jgs
200	* Merge of development branch back to main trunk on 2005-07-08
201	*
202	* Revision 1.4 2004/12/15 07:08:32 jgs
203	* * empty log message *
204	* Revision 1.1.1.1.2.3 2005/06/29 02:34:50 gross
205	* some changes towards 64 integers in finley
206	*
207	* Revision 1.1.1.1.2.2 2004/11/24 01:37:13 gross
208	* some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
209	*
210	*
211	*
212	*/
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision