Intrepid
Intrepid_RealSpaceToolsDef_Kokkos.hpp
1 #ifndef INTREPID_REALSPACETOOLSDEF_KOKKOS_HPP
2 #define INTREPID_REALSPACETOOLSDEF_KOKKOS_HPP
3 
4 #ifdef INTREPID_OLD_KOKKOS_CODE
5 #ifdef Kokkos_ENABLE_CXX11
6 #include "Intrepid_RealSpaceTools.hpp"
7 
8 namespace Intrepid{
9 template<class Scalar>
10 template<class Scalar1, class Scalar2,class Layout,class MemorySpace>
11 void RealSpaceTools<Scalar>::inverseTemp(Kokkos::View<Scalar1,Layout,MemorySpace> & inverseMats, const Kokkos::View<Scalar2,Layout,MemorySpace> & inMats){
12 
13  ArrayWrapper<Scalar1,Kokkos::View<Scalar1,Layout,MemorySpace>, Rank<Kokkos::View<Scalar1,Layout,MemorySpace> >::value, false>inverseMatsWrap(inverseMats);
14  ArrayWrapper<Scalar2,Kokkos::View<Scalar2,Layout,MemorySpace>, Rank<Kokkos::View<Scalar2,Layout,MemorySpace> >::value, true>inMatsWrap(inMats);
15 
16  int arrayRank = getrank(inMats);
17 
18  int dim_i0 = 1; // first index dimension (e.g. cell index)
19  int dim_i1 = 1; // second index dimension (e.g. point index)
20  int dim = inMats.dimension(arrayRank-2); // spatial dimension
21 
22  // determine i0 and i1 dimensions
23  switch(arrayRank) {
24  case 4:
25  dim_i0 = inMats.dimension(0);
26  dim_i1 = inMats.dimension(1);
27  switch(dim) {
28  case 3: {
29 
30  Kokkos::parallel_for (dim_i0, [&] (const int i0) {
31 
32  for (int i1=0; i1<dim_i1; i1++) {
33 
34 
35  int i, j, rowID = 0, colID = 0;
36  int rowperm[3]={0,1,2};
37  int colperm[3]={0,1,2}; // Complete pivoting
38  Scalar emax(0);
39 
40  for(i=0; i < 3; ++i){
41  for(j=0; j < 3; ++j){
42  if( std::abs( inMatsWrap(i0,i1,i,j) ) > emax){
43  rowID = i; colID = j; emax = std::abs( inMatsWrap(i0,i1,i,j) );
44  }
45  }
46  }
47 
48  if( rowID ){
49  rowperm[0] = rowID;
50  rowperm[rowID] = 0;
51  }
52  if( colID ){
53  colperm[0] = colID;
54  colperm[colID] = 0;
55  }
56  Scalar B[3][3], S[2][2], Bi[3][3]; // B=rowperm inMat colperm, S=Schur complement(Boo)
57  for(i=0; i < 3; ++i){
58  for(j=0; j < 3; ++j){
59  B[i][j] = inMatsWrap(i0,i1,rowperm[i],colperm[j]);
60  }
61  }
62  B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
63  for(i=0; i < 2; ++i){
64  for(j=0; j < 2; ++j){
65  S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
66  }
67  }
68  Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
69 
70 
71  Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS;
72  Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS;
73 
74  for(j=0; j<2;j++)
75  Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
76  for(i=0; i<2;i++)
77  Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
78 
79  Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
80  Bi[1][1] = Si[0][0];
81  Bi[1][2] = Si[0][1];
82  Bi[2][1] = Si[1][0];
83  Bi[2][2] = Si[1][1];
84  for(i=0; i < 3; ++i){
85  for(j=0; j < 3; ++j){
86  inverseMatsWrap(i0,i1,i,j) = Bi[colperm[i]][rowperm[j]]; // set inverse
87  }
88  }
89  } // for i1
90  }); // for i0
91  break;
92  } // case 3
93 
94  case 2: {
95 
96  Kokkos::parallel_for (dim_i0, [&] (const int i0) {
97 
98 
99  for (int i1=0; i1<dim_i1; i1++) {
100 
101 
102  Scalar determinant = inMatsWrap(i0,i1,0,0)*inMatsWrap(i0,i1,1,1)-inMatsWrap(i0,i1,0,1)*inMatsWrap(i0,i1,1,0);
103 
104  inverseMatsWrap(i0,i1,0,0) = inMatsWrap(i0,i1,1,1) / determinant;
105  inverseMatsWrap(i0,i1,0,1) = - inMatsWrap(i0,i1,0,1) / determinant;
106  //
107  inverseMatsWrap(i0,i1,1,0) = - inMatsWrap(i0,i1,1,0) / determinant;
108  inverseMatsWrap(i0,i1,1,1) = inMatsWrap(i0,i1,0,0) / determinant;
109  } // for i1
110  }); // for i0
111  break;
112  } // case 2
113 
114  case 1: {
115  Kokkos::parallel_for (dim_i0, [&] (const int i0) {
116 
117  for (int i1=0; i1<dim_i1; i1++) {
118  inverseMatsWrap(i0,i1,0,0) = (Scalar)1 / inMatsWrap(i0,i1,0,0);
119  } // for i1
120  }); // for i2
121 
122 
123  break;
124  } // case 1
125 
126  } // switch (dim)
127  break;
128  case 3:
129  dim_i1 = inMats.dimension(0);
130  switch(dim) {
131  case 3: {
132  Kokkos::parallel_for (dim_i1, [&] (const int i1) {
133 
134  int i, j, rowID = 0, colID = 0;
135  int rowperm[3]={0,1,2};
136  int colperm[3]={0,1,2}; // Complete pivoting
137  Scalar emax(0);
138 
139  for(i=0; i < 3; ++i){
140  for(j=0; j < 3; ++j){
141  if( std::abs( inMatsWrap(i1,i,j) ) > emax){
142  rowID = i; colID = j; emax = std::abs( inMatsWrap(i1,i,j) );
143  }
144  }
145  }
146 
147  if( rowID ){
148  rowperm[0] = rowID;
149  rowperm[rowID] = 0;
150  }
151  if( colID ){
152  colperm[0] = colID;
153  colperm[colID] = 0;
154  }
155  Scalar B[3][3], S[2][2], Bi[3][3]; // B=rowperm inMat colperm, S=Schur complement(Boo)
156  for(i=0; i < 3; ++i){
157  for(j=0; j < 3; ++j){
158  B[i][j] = inMatsWrap(i1,rowperm[i],colperm[j]);
159  }
160  }
161  B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
162  for(i=0; i < 2; ++i){
163  for(j=0; j < 2; ++j){
164  S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
165  }
166  }
167  Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
168 
169 
170  Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS;
171  Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS;
172 
173  for(j=0; j<2;j++)
174  Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
175  for(i=0; i<2;i++)
176  Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
177 
178  Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
179  Bi[1][1] = Si[0][0];
180  Bi[1][2] = Si[0][1];
181  Bi[2][1] = Si[1][0];
182  Bi[2][2] = Si[1][1];
183  for(i=0; i < 3; ++i){
184  for(j=0; j < 3; ++j){
185  inverseMatsWrap(i1,i,j) = Bi[colperm[i]][rowperm[j]]; // set inverse
186  }
187  }
188  }); // for i1
189 
190  break;
191  } // case 3
192 
193  case 2: {
194 
195  Kokkos::parallel_for (dim_i1, [&] (const int i1) {
196 
197  Scalar determinant = inMatsWrap(i1,0,0)*inMatsWrap(i1,1,1)-inMatsWrap(i1,0,1)*inMatsWrap(i1,1,0);
198 
199  inverseMatsWrap(i1,0,0) = inMatsWrap(i1,1,1) / determinant;
200  inverseMatsWrap(i1,0,1) = - inMatsWrap(i1,0,1) / determinant;
201  //
202  inverseMatsWrap(i1,1,0) = - inMatsWrap(i1,1,0) / determinant;
203  inverseMatsWrap(i1,1,1) = inMatsWrap(i1,0,0) / determinant;
204  }); // for i1
205 
206  break;
207  } // case 2
208 
209  case 1: {
210 
211  Kokkos::parallel_for (dim_i1, [&] (const int i1) {
212  inverseMatsWrap(i1,0,0) = (Scalar)1 / inMatsWrap(i1,0,0);
213  });
214 
215 
216 
217  break;
218  } // case 1
219 
220  } // switch (dim)
221  break;
222  case 2:
223  switch(dim) {
224  case 3: {
225  int i, j, rowID = 0, colID = 0;
226  int rowperm[3]={0,1,2};
227  int colperm[3]={0,1,2}; // Complete pivoting
228  Scalar emax(0);
229 
230  for(i=0; i < 3; ++i){
231  for(j=0; j < 3; ++j){
232  if( std::abs( inMatsWrap(i,j) ) > emax){
233  rowID = i; colID = j; emax = std::abs( inMatsWrap(i,j) );
234  }
235  }
236  }
237 
238  if( rowID ){
239  rowperm[0] = rowID;
240  rowperm[rowID] = 0;
241  }
242  if( colID ){
243  colperm[0] = colID;
244  colperm[colID] = 0;
245  }
246  Scalar B[3][3], S[2][2], Bi[3][3]; // B=rowperm inMat colperm, S=Schur complement(Boo)
247  for(i=0; i < 3; ++i){
248  for(j=0; j < 3; ++j){
249  B[i][j] = inMatsWrap(rowperm[i],colperm[j]);
250  }
251  }
252  B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
253  for(i=0; i < 2; ++i){
254  for(j=0; j < 2; ++j){
255  S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
256  }
257  }
258  Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
259 
260  Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS;
261  Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS;
262 
263  for(j=0; j<2;j++)
264  Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
265  for(i=0; i<2;i++)
266  Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
267 
268  Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
269  Bi[1][1] = Si[0][0];
270  Bi[1][2] = Si[0][1];
271  Bi[2][1] = Si[1][0];
272  Bi[2][2] = Si[1][1];
273  for(i=0; i < 3; ++i){
274  for(j=0; j < 3; ++j){
275  inverseMatsWrap(i,j) = Bi[colperm[i]][rowperm[j]]; // set inverse
276  }
277  }
278 
279  break;
280  } // case 3
281 
282  case 2: {
283  Scalar determinant = inMatsWrap(0,0)*inMatsWrap(1,1)-inMatsWrap(0,1)*inMatsWrap(1,0);
284 
285  inverseMatsWrap(0,0) = inMatsWrap(1,1) / determinant;
286  inverseMatsWrap(0,1) = - inMatsWrap(0,1) / determinant;
287  //
288  inverseMatsWrap(1,0) = - inMatsWrap(1,0) / determinant;
289  inverseMatsWrap(1,1) = inMatsWrap(0,0) / determinant;
290 
291  break;
292  } // case 2
293 
294  case 1: {
295  inverseMatsWrap(0,0) = (Scalar)1 / inMatsWrap(0,0);
296  break;
297  } // case 1
298 
299  }
300  break;
301  }
302 
303 
304 
305 
306 }
307 
308 }
309 
310 
311 #endif
312 #endif
313 #endif
Intrepid_RealSpaceTools.hpp
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
Generated by   doxygen 1.8.5