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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4803 - (show annotations)
Wed Mar 26 06:52:28 2014 UTC (5 years, 9 months ago) by caltinay
File size: 6630 byte(s)
Removed obsolete wrappers for malloc and friends.
Paso_Pattern -> paso::Pattern

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

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