MueLu  Version of the Day
MueLu_LineDetectionFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_LINEDETECTIONFACTORY_DEF_HPP
47 #define MUELU_LINEDETECTIONFACTORY_DEF_HPP
48 
49 #include <Xpetra_Matrix.hpp>
50 //#include <Xpetra_MatrixFactory.hpp>
51 
53 
54 #include "MueLu_FactoryManager.hpp"
55 #include "MueLu_Level.hpp"
56 #include "MueLu_MasterList.hpp"
57 #include "MueLu_Monitor.hpp"
58 
59 namespace MueLu {
60 
61  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
63  RCP<ParameterList> validParamList = rcp(new ParameterList());
64 
65 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
66  SET_VALID_ENTRY("linedetection: orientation");
67  SET_VALID_ENTRY("linedetection: num layers");
68 #undef SET_VALID_ENTRY
69 
70  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
71  validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory for coorindates");
72 
73  return validParamList;
74  }
75 
76  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
78  Input(currentLevel, "A");
79 
80  // The factory needs the information about the number of z-layers. While this information is
81  // provided by the user for the finest level, the factory itself is responsible to provide the
82  // corresponding information on the coarser levels. Since a factory cannot be dependent on itself
83  // we use the NoFactory class as generator class, but remove the UserData keep flag, such that
84  // "NumZLayers" is part of the request/release mechanism.
85  // Please note, that this prevents us from having several (independent) CoarsePFactory instances!
86  // TODO: allow factory to dependent on self-generated data for TwoLevelFactories -> introduce ExpertRequest/Release in Level
87  currentLevel.DeclareInput("NumZLayers", NoFactory::get(), this);
88  currentLevel.RemoveKeepFlag("NumZLayers", NoFactory::get(), MueLu::UserData);
89  }
90 
91  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
93  FactoryMonitor m(*this, "Line detection (Ray style)", currentLevel);
94 
95  LO NumZDir = 0;
96  RCP<CoordinateMultiVector> fineCoords;
97  ArrayRCP<coordinate_type> x, y, z;
98  coordinate_type *xptr = NULL, *yptr = NULL, *zptr = NULL;
99 
100  // obtain general variables
101  RCP<Matrix> A = Get< RCP<Matrix> > (currentLevel, "A");
102  LO BlkSize = A->GetFixedBlockSize();
103  RCP<const Map> rowMap = A->getRowMap();
104  LO Ndofs = rowMap->getNodeNumElements();
105  LO Nnodes = Ndofs/BlkSize;
106 
107  // collect information provided by user
108  const ParameterList& pL = GetParameterList();
109  const std::string lineOrientation = pL.get<std::string>("linedetection: orientation");
110 
111  // interpret "line orientation" parameter provided by the user on the finest level
112  if(currentLevel.GetLevelID() == 0) {
113  if(lineOrientation=="vertical")
114  Zorientation_ = VERTICAL;
115  else if (lineOrientation=="horizontal")
116  Zorientation_ = HORIZONTAL;
117  else if (lineOrientation=="coordinates")
118  Zorientation_ = GRID_SUPPLIED;
119  else
120  TEUCHOS_TEST_FOR_EXCEPTION(false, Exceptions::RuntimeError, "LineDetectionFactory: The parameter 'semicoarsen: line orientation' must be either 'vertical', 'horizontal' or 'coordinates'.");
121  }
122 
123  //TEUCHOS_TEST_FOR_EXCEPTION(Zorientation_!=VERTICAL, Exceptions::RuntimeError, "LineDetectionFactory: The 'horizontal' or 'coordinates' have not been tested!!!. Please remove this exception check and carefully test these modes!");
124 
125  // obtain number of z layers (variable over levels)
126  // This information is user-provided on the finest level and transferred to the coarser
127  // levels by the SemiCoarsenPFactor using the internal "NumZLayers" variable.
128  if(currentLevel.GetLevelID() == 0) {
129  if(currentLevel.IsAvailable("NumZLayers", NoFactory::get())) {
130  NumZDir = currentLevel.Get<LO>("NumZLayers", NoFactory::get()); //obtain info
131  GetOStream(Runtime1) << "Number of layers for line detection: " << NumZDir << " (information from Level(0))" << std::endl;
132  } else {
133  // check whether user provides information or it can be reconstructed from coordinates
134  NumZDir = pL.get<LO>("linedetection: num layers");
135  if(NumZDir == -1) {
136  bool CoordsAvail = currentLevel.IsAvailable("Coordinates");
137 
138  if (CoordsAvail == true) {
139  // try to reconstruct the number of layers from coordinates
140  fineCoords = Get< RCP<CoordinateMultiVector> > (currentLevel, "Coordinates");
141  TEUCHOS_TEST_FOR_EXCEPTION(fineCoords->getNumVectors() != 3, Exceptions::RuntimeError, "Three coordinates arrays must be supplied if line detection orientation not given.");
142  x = fineCoords->getDataNonConst(0);
143  y = fineCoords->getDataNonConst(1);
144  z = fineCoords->getDataNonConst(2);
145  xptr = x.getRawPtr();
146  yptr = y.getRawPtr();
147  zptr = z.getRawPtr();
148 
149  LO NumCoords = Ndofs/BlkSize;
150 
151  /* sort coordinates so that we can order things according to lines */
152  Teuchos::ArrayRCP<LO> TOrigLoc= Teuchos::arcp<LO>(NumCoords); LO* OrigLoc= TOrigLoc.getRawPtr();
153  Teuchos::ArrayRCP<coordinate_type> Txtemp = Teuchos::arcp<coordinate_type>(NumCoords); coordinate_type* xtemp = Txtemp.getRawPtr();
154  Teuchos::ArrayRCP<coordinate_type> Tytemp = Teuchos::arcp<coordinate_type>(NumCoords); coordinate_type* ytemp = Tytemp.getRawPtr();
155  Teuchos::ArrayRCP<coordinate_type> Tztemp = Teuchos::arcp<coordinate_type>(NumCoords); coordinate_type* ztemp = Tztemp.getRawPtr();
156 
157  // sort coordinates in {x,y,z}vals (returned in {x,y,z}temp) so that we can order things according to lines
158  // switch x and y coordinates for semi-coarsening...
159  sort_coordinates(NumCoords, OrigLoc, xptr, yptr, zptr, xtemp, ytemp, ztemp, true);
160 
161  /* go through each vertical line and populate blockIndices so all */
162  /* dofs within a PDE within a vertical line correspond to one block.*/
163  LO NumBlocks = 0;
164  LO NumNodesPerVertLine = 0;
165  LO index = 0;
166 
167  while ( index < NumCoords ) {
168  coordinate_type xfirst = xtemp[index]; coordinate_type yfirst = ytemp[index];
169  LO next = index+1;
170  while ( (next != NumCoords) && (xtemp[next] == xfirst) &&
171  (ytemp[next] == yfirst))
172  next++;
173  if (NumBlocks == 0) {
174  NumNodesPerVertLine = next-index;
175  }
176  // the number of vertical lines must be the same on all processors
177  // TAW: Sep 14 2015: or zero as we allow "empty" processors
178  //TEUCHOS_TEST_FOR_EXCEPTION(next-index != NumNodesPerVertLine,Exceptions::RuntimeError, "Error code only works for constant block size now!!!\n");
179  NumBlocks++;
180  index = next;
181  }
182 
183  NumZDir = NumNodesPerVertLine;
184  GetOStream(Runtime1) << "Number of layers for line detection: " << NumZDir << " (information reconstructed from provided node coordinates)" << std::endl;
185  } else {
186  TEUCHOS_TEST_FOR_EXCEPTION(false, Exceptions::RuntimeError, "LineDetectionFactory: BuildP: User has to provide valid number of layers (e.g. using the 'line detection: num layers' parameter).");
187  }
188  } else {
189  GetOStream(Runtime1) << "Number of layers for line detection: " << NumZDir << " (information provided by user through 'line detection: num layers')" << std::endl;
190  }
191  } // end else (user provides information or can be reconstructed) on finest level
192  } else {
193  // coarse level information
194  // TODO get rid of NoFactory here and use SemiCoarsenPFactory as source of NumZLayers instead.
195  if(currentLevel.IsAvailable("NumZLayers", NoFactory::get())) {
196  NumZDir = currentLevel.Get<LO>("NumZLayers", NoFactory::get()); //obtain info
197  GetOStream(Runtime1) << "Number of layers for line detection: " << NumZDir << std::endl;
198  } else {
199  TEUCHOS_TEST_FOR_EXCEPTION(false, Exceptions::RuntimeError, "LineDetectionFactory: BuildP: No NumZLayers variable found. This cannot be.");
200  }
201  }
202 
203  // plausibility check and further variable collection
204  if (Zorientation_ == GRID_SUPPLIED) { // On finest level, fetch user-provided coordinates if available...
205  bool CoordsAvail = currentLevel.IsAvailable("Coordinates");
206 
207  if (CoordsAvail == false) {
208  if (currentLevel.GetLevelID() == 0)
209  throw Exceptions::RuntimeError("Coordinates must be supplied if line detection orientation not given.");
210  else
211  throw Exceptions::RuntimeError("Coordinates not generated by previous invocation of LineDetectionFactory's BuildP() method.");
212  }
213  fineCoords = Get< RCP<CoordinateMultiVector> > (currentLevel, "Coordinates");
214  TEUCHOS_TEST_FOR_EXCEPTION(fineCoords->getNumVectors() != 3, Exceptions::RuntimeError, "Three coordinates arrays must be supplied if line detection orientation not given.");
215  x = fineCoords->getDataNonConst(0);
216  y = fineCoords->getDataNonConst(1);
217  z = fineCoords->getDataNonConst(2);
218  xptr = x.getRawPtr();
219  yptr = y.getRawPtr();
220  zptr = z.getRawPtr();
221  }
222 
223  // perform line detection
224  if (NumZDir > 0) {
225  LO *LayerId, *VertLineId;
226  Teuchos::ArrayRCP<LO> TLayerId = Teuchos::arcp<LO>(Nnodes); LayerId = TLayerId.getRawPtr();
227  Teuchos::ArrayRCP<LO> TVertLineId= Teuchos::arcp<LO>(Nnodes); VertLineId = TVertLineId.getRawPtr();
228 
229  NumZDir = ML_compute_line_info(LayerId, VertLineId, Ndofs, BlkSize,
230  Zorientation_, NumZDir,xptr,yptr,zptr, *(rowMap->getComm()));
231  //it is NumZDir=NCLayers*NVertLines*DofsPerNode;
232 
233  // store output data on current level
234  // The line detection data is used by the SemiCoarsenPFactory and the line smoothers in Ifpack/Ifpack2
235  Set(currentLevel, "CoarseNumZLayers", NumZDir);
236  Set(currentLevel, "LineDetection_Layers", TLayerId);
237  Set(currentLevel, "LineDetection_VertLineIds", TVertLineId);
238  } else {
239  Teuchos::ArrayRCP<LO> TLayerId = Teuchos::arcp<LO>(0);
240  Teuchos::ArrayRCP<LO> TVertLineId = Teuchos::arcp<LO>(0);
241  Teuchos::ArrayRCP<LO> TVertLineIdSmoo= Teuchos::arcp<LO>(0);
242 
243  // store output data on current level
244  // The line detection data is used by the SemiCoarsenPFactory and the line smoothers in Ifpack/Ifpack2
245  Set(currentLevel, "CoarseNumZLayers", NumZDir);
246  Set(currentLevel, "LineDetection_Layers", TLayerId);
247  Set(currentLevel, "LineDetection_VertLineIds", TVertLineId);
248  }
249 
250  // automatically switch to vertical mode on the coarser levels
251  if(Zorientation_ != VERTICAL)
252  Zorientation_ = VERTICAL;
253  }
254 
255  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
256  LocalOrdinal LineDetectionFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ML_compute_line_info(LocalOrdinal LayerId[], LocalOrdinal VertLineId[], LocalOrdinal Ndof, LocalOrdinal DofsPerNode, LocalOrdinal MeshNumbering, LocalOrdinal NumNodesPerVertLine, typename Teuchos::ScalarTraits<Scalar>::coordinateType *xvals, typename Teuchos::ScalarTraits<Scalar>::coordinateType *yvals, typename Teuchos::ScalarTraits<Scalar>::coordinateType *zvals, const Teuchos::Comm<int>& /* comm */) const {
257 
258  LO Nnodes, NVertLines, MyNode;
259  LO NumCoords, next; //, subindex, subnext;
260  coordinate_type xfirst, yfirst;
261  coordinate_type *xtemp, *ytemp, *ztemp;
262  LO *OrigLoc;
263  LO i,j,count;
264  LO RetVal;
265 
266  RetVal = 0;
267  if ((MeshNumbering != VERTICAL) && (MeshNumbering != HORIZONTAL)) {
268  if ( (xvals == NULL) || (yvals == NULL) || (zvals == NULL)) RetVal = -1;
269  }
270  else {
271  if (NumNodesPerVertLine == -1) RetVal = -4;
272  if ( ((Ndof/DofsPerNode)%NumNodesPerVertLine) != 0) RetVal = -3;
273  }
274  if ( (Ndof%DofsPerNode) != 0) RetVal = -2;
275 
276  TEUCHOS_TEST_FOR_EXCEPTION(RetVal == -1, Exceptions::RuntimeError, "Not semicoarsening as no mesh numbering information or coordinates are given\n");
277  TEUCHOS_TEST_FOR_EXCEPTION(RetVal == -4, Exceptions::RuntimeError, "Not semicoarsening as the number of z nodes is not given.\n");
278  TEUCHOS_TEST_FOR_EXCEPTION(RetVal == -3, Exceptions::RuntimeError, "Not semicoarsening as the total number of nodes is not evenly divisible by the number of z direction nodes .\n");
279  TEUCHOS_TEST_FOR_EXCEPTION(RetVal == -2, Exceptions::RuntimeError, "Not semicoarsening as something is off with the number of degrees-of-freedom per node.\n");
280 
281  Nnodes = Ndof/DofsPerNode;
282  for (MyNode = 0; MyNode < Nnodes; MyNode++) VertLineId[MyNode] = -1;
283  for (MyNode = 0; MyNode < Nnodes; MyNode++) LayerId[MyNode] = -1;
284 
285  if (MeshNumbering == VERTICAL) {
286  for (MyNode = 0; MyNode < Nnodes; MyNode++) {
287  LayerId[MyNode]= MyNode%NumNodesPerVertLine;
288  VertLineId[MyNode]= (MyNode- LayerId[MyNode])/NumNodesPerVertLine;
289  }
290  }
291  else if (MeshNumbering == HORIZONTAL) {
292  NVertLines = Nnodes/NumNodesPerVertLine;
293  for (MyNode = 0; MyNode < Nnodes; MyNode++) {
294  VertLineId[MyNode] = MyNode%NVertLines;
295  LayerId[MyNode] = (MyNode- VertLineId[MyNode])/NVertLines;
296  }
297  }
298  else {
299  // coordinates mode: we distinguish between vertical line numbering for semi-coarsening and line smoothing
300  NumCoords = Ndof/DofsPerNode;
301 
302  // reserve temporary memory
303  Teuchos::ArrayRCP<LO> TOrigLoc= Teuchos::arcp<LO>(NumCoords); OrigLoc= TOrigLoc.getRawPtr();
304  Teuchos::ArrayRCP<coordinate_type> Txtemp = Teuchos::arcp<coordinate_type>(NumCoords); xtemp = Txtemp.getRawPtr();
305  Teuchos::ArrayRCP<coordinate_type> Tytemp = Teuchos::arcp<coordinate_type>(NumCoords); ytemp = Tytemp.getRawPtr();
306  Teuchos::ArrayRCP<coordinate_type> Tztemp = Teuchos::arcp<coordinate_type>(NumCoords); ztemp = Tztemp.getRawPtr();
307 
308  // build vertical line info for semi-coarsening
309 
310  // sort coordinates in {x,y,z}vals (returned in {x,y,z}temp) so that we can order things according to lines
311  // switch x and y coordinates for semi-coarsening...
312  sort_coordinates(NumCoords, OrigLoc, xvals, yvals, zvals, xtemp, ytemp, ztemp, /*true*/ true);
313 
314  LO NumBlocks = 0;
315  LO index = 0;
316 
317  while ( index < NumCoords ) {
318  xfirst = xtemp[index]; yfirst = ytemp[index];
319  next = index+1;
320  while ( (next != NumCoords) && (xtemp[next] == xfirst) &&
321  (ytemp[next] == yfirst))
322  next++;
323  if (NumBlocks == 0) {
324  NumNodesPerVertLine = next-index;
325  }
326  // The number of vertical lines must be the same on all processors
327  // TAW: Sep 14, 2015: or zero as we allow for empty processors.
328  //TEUCHOS_TEST_FOR_EXCEPTION(next-index != NumNodesPerVertLine,Exceptions::RuntimeError, "Error code only works for constant block size now!!!\n");
329  count = 0;
330  for (j= index; j < next; j++) {
331  VertLineId[OrigLoc[j]] = NumBlocks;
332  LayerId[OrigLoc[j]] = count++;
333  }
334  NumBlocks++;
335  index = next;
336  }
337  }
338 
339  /* check that everyone was assigned */
340 
341  for (i = 0; i < Nnodes; i++) {
342  if (VertLineId[i] == -1) {
343  GetOStream(Warnings1) << "Warning: did not assign " << i << " to a vertical line?????\n" << std::endl;
344  }
345  if (LayerId[i] == -1) {
346  GetOStream(Warnings1) << "Warning: did not assign " << i << " to a Layer?????\n" << std::endl;
347  }
348  }
349 
350  // TAW: Sep 14 2015: relax plausibility checks as we allow for empty processors
351  //MueLu_maxAll(&comm, NumNodesPerVertLine, i);
352  //if (NumNodesPerVertLine == -1) NumNodesPerVertLine = i;
353  //TEUCHOS_TEST_FOR_EXCEPTION(NumNodesPerVertLine != i,Exceptions::RuntimeError, "Different processors have different z direction line lengths?\n");
354 
355  return NumNodesPerVertLine;
356  }
357 
358  /* Private member function to sort coordinates in arrays. This is an expert routine. Do not use or change.*/
359  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
361  typename Teuchos::ScalarTraits<Scalar>::coordinateType* xvals,
362  typename Teuchos::ScalarTraits<Scalar>::coordinateType* yvals,
363  typename Teuchos::ScalarTraits<Scalar>::coordinateType* zvals,
364  typename Teuchos::ScalarTraits<Scalar>::coordinateType* xtemp,
365  typename Teuchos::ScalarTraits<Scalar>::coordinateType* ytemp,
366  typename Teuchos::ScalarTraits<Scalar>::coordinateType* ztemp,
367  bool flipXY) const {
368 
369  if( flipXY == false ) { // for line-smoothing
370  for (LO i = 0; i < numCoords; i++) xtemp[i]= xvals[i];
371  } else { // for semi-coarsening
372  for (LO i = 0; i < numCoords; i++) xtemp[i]= yvals[i];
373  }
374  for (LO i = 0; i < numCoords; i++) OrigLoc[i]= i;
375 
376  ML_az_dsort2(xtemp,numCoords,OrigLoc);
377  if( flipXY == false ) { // for line-smoothing
378  for (LO i = 0; i < numCoords; i++) ytemp[i]= yvals[OrigLoc[i]];
379  } else {
380  for (LO i = 0; i < numCoords; i++) ytemp[i]= xvals[OrigLoc[i]];
381  }
382 
383  LO index = 0;
384 
385  while ( index < numCoords ) {
386  coordinate_type xfirst = xtemp[index];
387  LO next = index+1;
388  while ( (next != numCoords) && (xtemp[next] == xfirst))
389  next++;
390  ML_az_dsort2(&(ytemp[index]),next-index,&(OrigLoc[index]));
391  for (LO i = index; i < next; i++) ztemp[i]= zvals[OrigLoc[i]];
392  /* One final sort so that the ztemps are in order */
393  LO subindex = index;
394  while (subindex != next) {
395  coordinate_type yfirst = ytemp[subindex];
396  LO subnext = subindex+1;
397  while ( (subnext != next) && (ytemp[subnext] == yfirst)) subnext++;
398  ML_az_dsort2(&(ztemp[subindex]),subnext-subindex,&(OrigLoc[subindex]));
399  subindex = subnext;
400  }
401  index = next;
402  }
403 
404  }
405 
406  /* Sort coordinates and additional array accordingly (if provided). This is an expert routine borrowed from ML. Do not change.*/
407  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
408  void LineDetectionFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ML_az_dsort2(typename Teuchos::ScalarTraits<Scalar>::coordinateType dlist[], LocalOrdinal N, LocalOrdinal list2[]) const {
409  LO l, r, j, i, flag;
410  LO RR2;
411  coordinate_type dRR, dK;
412 
413  // note: we use that routine for sorting coordinates only. No complex coordinates are assumed...
414  typedef Teuchos::ScalarTraits<SC> STS;
415 
416  if (N <= 1) return;
417 
418  l = N / 2 + 1;
419  r = N - 1;
420  l = l - 1;
421  dRR = dlist[l - 1];
422  dK = dlist[l - 1];
423 
424  if (list2 != NULL) {
425  RR2 = list2[l - 1];
426  while (r != 0) {
427  j = l;
428  flag = 1;
429 
430  while (flag == 1) {
431  i = j;
432  j = j + j;
433 
434  if (j > r + 1)
435  flag = 0;
436  else {
437  if (j < r + 1)
438  if (STS::real(dlist[j]) > STS::real(dlist[j - 1])) j = j + 1;
439 
440  if (STS::real(dlist[j - 1]) > STS::real(dK)) {
441  dlist[ i - 1] = dlist[ j - 1];
442  list2[i - 1] = list2[j - 1];
443  }
444  else {
445  flag = 0;
446  }
447  }
448  }
449  dlist[ i - 1] = dRR;
450  list2[i - 1] = RR2;
451 
452  if (l == 1) {
453  dRR = dlist [r];
454  RR2 = list2[r];
455  dK = dlist[r];
456  dlist[r ] = dlist[0];
457  list2[r] = list2[0];
458  r = r - 1;
459  }
460  else {
461  l = l - 1;
462  dRR = dlist[ l - 1];
463  RR2 = list2[l - 1];
464  dK = dlist[l - 1];
465  }
466  }
467  dlist[ 0] = dRR;
468  list2[0] = RR2;
469  }
470  else {
471  while (r != 0) {
472  j = l;
473  flag = 1;
474  while (flag == 1) {
475  i = j;
476  j = j + j;
477  if (j > r + 1)
478  flag = 0;
479  else {
480  if (j < r + 1)
481  if (STS::real(dlist[j]) > STS::real(dlist[j - 1])) j = j + 1;
482  if (STS::real(dlist[j - 1]) > STS::real(dK)) {
483  dlist[ i - 1] = dlist[ j - 1];
484  }
485  else {
486  flag = 0;
487  }
488  }
489  }
490  dlist[ i - 1] = dRR;
491  if (l == 1) {
492  dRR = dlist [r];
493  dK = dlist[r];
494  dlist[r ] = dlist[0];
495  r = r - 1;
496  }
497  else {
498  l = l - 1;
499  dRR = dlist[ l - 1];
500  dK = dlist[l - 1];
501  }
502  }
503  dlist[ 0] = dRR;
504  }
505 
506  }
507 } //namespace MueLu
508 
509 #endif // MUELU_LINEDETECTIONFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
MueLu::DefaultLocalOrdinal LocalOrdinal
void DeclareInput(Level &currentLevel) const
Input.
Timer to be used in factories. Similar to Monitor but with additional timers.
void Build(Level &currentLevel) const
Build method.
void ML_az_dsort2(coordinate_type dlist[], LO N, LO list2[]) const
User data are always kept. This flag is set automatically when Level::Set("data", data) is used...
typename Teuchos::ScalarTraits< SC >::coordinateType coordinate_type
Namespace for MueLu classes and methods.
static const NoFactory * get()
Additional warnings.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
void sort_coordinates(LO numCoords, LO *OrigLoc, coordinate_type *xvals, coordinate_type *yvals, coordinate_type *zvals, coordinate_type *xtemp, coordinate_type *ytemp, coordinate_type *ztemp, bool flipXY=false) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
LO ML_compute_line_info(LO LayerId[], LO VertLineId[], LO Ndof, LO DofsPerNode, LO MeshNumbering, LO NumNodesPerVertLine, coordinate_type *xvals, coordinate_type *yvals, coordinate_type *zvals, const Teuchos::Comm< int > &comm) const
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()