/[escript]/trunk/finley/src/IndexList.cpp
ViewVC logotype

Contents of /trunk/finley/src/IndexList.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4492 - (show annotations)
Tue Jul 2 01:44:11 2013 UTC (5 years, 9 months ago) by caltinay
File size: 6526 byte(s)
Finley changes that were held back while in release mode - moved more stuff
into finley namespace.

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 Finley: Converting an element list into a matrix shape
20
21 *****************************************************************************/
22
23 #include "IndexList.h"
24 #include "ElementFile.h"
25
26 #include "paso/SystemMatrixPattern.h"
27
28 /* Translate from distributed/local array indices to global indices */
29
30 /****************************************************************************/
31 /* inserts the contributions from the element matrices of elements
32 into the row index col. If symmetric is set, only the upper
33 triangle of the matrix is stored.
34 */
35
36 namespace finley {
37
38 void IndexList_insertElements(IndexList* index_list, ElementFile* elements,
39 bool reduce_row_order, const int* row_map,
40 bool reduce_col_order, const int* col_map)
41 {
42 // index_list is an array of linked lists. Each entry is a row (DOF) and
43 // contains the indices to the non-zero columns
44 if (!elements)
45 return;
46
47 const int NN=elements->numNodes;
48 ReferenceElement* refElement =
49 ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, false);
50
51 int NN_row, NN_col, *row_node=NULL, *col_node=NULL, numSub;
52 if (reduce_col_order) {
53 numSub=1;
54 col_node=refElement->Type->linearNodes;
55 NN_col=refElement->LinearBasisFunctions->Type->numShapes * refElement->Type->numSides;
56 } else {
57 numSub=refElement->Type->numSubElements;
58 col_node=refElement->Type->subElementNodes;
59 NN_col=refElement->BasisFunctions->Type->numShapes * refElement->Type->numSides;
60 }
61
62 if (reduce_row_order) {
63 numSub=1;
64 row_node=refElement->Type->linearNodes;
65 NN_row=refElement->LinearBasisFunctions->Type->numShapes * refElement->Type->numSides;
66 } else {
67 numSub=refElement->Type->numSubElements;
68 row_node=refElement->Type->subElementNodes;
69 NN_row=refElement->BasisFunctions->Type->numShapes * refElement->Type->numSides;
70 }
71
72 for (int color=elements->minColor; color<=elements->maxColor; color++) {
73 #pragma omp for
74 for (int e=0; e<elements->numElements; e++) {
75 if (elements->Color[e]==color) {
76 for (int isub=0; isub<numSub; isub++) {
77 for (int kr=0; kr<NN_row; kr++) {
78 const int irow=row_map[elements->Nodes[INDEX2(row_node[INDEX2(kr,isub,NN_row)],e,NN)]];
79 for (int kc=0; kc<NN_col; kc++) {
80 const int icol=col_map[elements->Nodes[INDEX2(col_node[INDEX2(kc,isub,NN_col)],e,NN)]];
81 IndexList_insertIndex(index_list[irow], icol);
82 }
83 }
84 }
85 }
86 }
87 }
88 }
89
90 void IndexList_insertElementsWithRowRangeNoMainDiagonal(IndexList* index_list,
91 int firstRow, int lastRow, ElementFile* elements,
92 int* row_map, int* col_map)
93 {
94 if (!elements)
95 return;
96
97 // this does not resolve macro elements
98 const int NN=elements->numNodes;
99 for (int color=elements->minColor; color<=elements->maxColor; color++) {
100 #pragma omp for
101 for (int e=0; e<elements->numElements; e++) {
102 if (elements->Color[e]==color) {
103 for (int kr=0; kr<NN; kr++) {
104 const int irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
105 if (firstRow<=irow && irow<lastRow) {
106 const int irow_loc=irow-firstRow;
107 for (int kc=0; kc<NN; kc++) {
108 const int icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
109 if (icol != irow)
110 IndexList_insertIndex(index_list[irow_loc], icol);
111 }
112 }
113 }
114 }
115 }
116 }
117 }
118
119 /// inserts a row index into the IndexList in if it does not exist
120 void IndexList_insertIndex(IndexList& in, int index)
121 {
122 if (std::find(in.begin(), in.end(), index) == in.end())
123 in.push_back(index);
124 }
125
126 /// counts the number of row indices in the IndexList in
127 int IndexList_count(const IndexList& in, int range_min, int range_max)
128 {
129 int out=0;
130 for (IndexList::const_iterator it=in.begin(); it!=in.end(); it++) {
131 if (*it >= range_min && range_max > *it)
132 ++out;
133 }
134 return out;
135 }
136
137 /// count the number of row indices in the IndexList in
138 void IndexList_toArray(const IndexList& in, int* array, int range_min, int range_max, int index_offset)
139 {
140 int idx=0;
141 for (IndexList::const_iterator it=in.begin(); it!=in.end(); it++) {
142 if (*it >= range_min && range_max > *it) {
143 array[idx]=(*it)+index_offset;
144 ++idx;
145 }
146 }
147 }
148
149 /// creates a Paso_pattern from a range of indices
150 Paso_Pattern* IndexList_createPattern(int n0, int n,
151 const IndexList* index_list,
152 int range_min, int range_max,
153 int index_offset)
154 {
155 int *ptr=new int[n+1-n0];
156 // get the number of connections per row
157 #pragma omp parallel for
158 for (int i=n0; i<n; ++i) {
159 ptr[i-n0]=IndexList_count(index_list[i], range_min, range_max);
160 }
161 // accumulate ptr
162 int s=0;
163 for (int i=n0; i<n; ++i) {
164 const int itmp=ptr[i-n0];
165 ptr[i-n0]=s;
166 s+=itmp;
167 }
168 ptr[n-n0]=s;
169 // fill index
170 int *index=new int[ptr[n-n0]];
171 #pragma omp parallel for
172 for (int i=n0; i<n; ++i) {
173 IndexList_toArray(index_list[i], &index[ptr[i-n0]], range_min, range_max, index_offset);
174 }
175
176 Paso_Pattern* out=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, n-n0, range_max+index_offset, ptr, index);
177 if (!Finley_noError()) {
178 delete[] ptr;
179 delete[] index;
180 Paso_Pattern_free(out);
181 }
182 return out;
183 }
184
185 } // namespace finley
186

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/lapack2681/finley/src/IndexList.cpp:2682-2741 /branches/pasowrap/finley/src/IndexList.cpp:3661-3674 /branches/py3_attempt2/finley/src/IndexList.cpp:3871-3891 /branches/restext/finley/src/IndexList.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/IndexList.cpp:3669-3791 /branches/stage3.0/finley/src/IndexList.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/IndexList.cpp:3471-3974 /release/3.0/finley/src/IndexList.cpp:2591-2601 /trunk/finley/src/IndexList.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26