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<MultiVector> fineCoords;
97  ArrayRCP<Scalar> x, y, z;
98  Scalar *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<MultiVector> > (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<SC> Txtemp = Teuchos::arcp<SC>(NumCoords); SC* xtemp = Txtemp.getRawPtr();
154  Teuchos::ArrayRCP<SC> Tytemp = Teuchos::arcp<SC>(NumCoords); SC* ytemp = Tytemp.getRawPtr();
155  Teuchos::ArrayRCP<SC> Tztemp = Teuchos::arcp<SC>(NumCoords); SC* 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  SC xfirst = xtemp[index]; SC 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<MultiVector> > (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, Scalar *xvals, Scalar *yvals, Scalar *zvals, const Teuchos::Comm<int>& comm) const {
257 
258  LO Nnodes, NVertLines, MyNode;
259  LO NumCoords, next; //, subindex, subnext;
260  SC xfirst, yfirst;
261  SC *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<SC> Txtemp = Teuchos::arcp<SC>(NumCoords); xtemp = Txtemp.getRawPtr();
305  Teuchos::ArrayRCP<SC> Tytemp = Teuchos::arcp<SC>(NumCoords); ytemp = Tytemp.getRawPtr();
306  Teuchos::ArrayRCP<SC> Tztemp = Teuchos::arcp<SC>(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>
360  void LineDetectionFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sort_coordinates(LO numCoords, LO* OrigLoc, Scalar* xvals, Scalar* yvals, Scalar* zvals, Scalar* xtemp, Scalar* ytemp, Scalar* ztemp, bool flipXY) const {
361 
362  if( flipXY == false ) { // for line-smoothing
363  for (LO i = 0; i < numCoords; i++) xtemp[i]= xvals[i];
364  } else { // for semi-coarsening
365  for (LO i = 0; i < numCoords; i++) xtemp[i]= yvals[i];
366  }
367  for (LO i = 0; i < numCoords; i++) OrigLoc[i]= i;
368 
369  ML_az_dsort2(xtemp,numCoords,OrigLoc);
370  if( flipXY == false ) { // for line-smoothing
371  for (LO i = 0; i < numCoords; i++) ytemp[i]= yvals[OrigLoc[i]];
372  } else {
373  for (LO i = 0; i < numCoords; i++) ytemp[i]= xvals[OrigLoc[i]];
374  }
375 
376  LO index = 0;
377 
378  while ( index < numCoords ) {
379  SC xfirst = xtemp[index];
380  LO next = index+1;
381  while ( (next != numCoords) && (xtemp[next] == xfirst))
382  next++;
383  ML_az_dsort2(&(ytemp[index]),next-index,&(OrigLoc[index]));
384  for (LO i = index; i < next; i++) ztemp[i]= zvals[OrigLoc[i]];
385  /* One final sort so that the ztemps are in order */
386  LO subindex = index;
387  while (subindex != next) {
388  SC yfirst = ytemp[subindex];
389  LO subnext = subindex+1;
390  while ( (subnext != next) && (ytemp[subnext] == yfirst)) subnext++;
391  ML_az_dsort2(&(ztemp[subindex]),subnext-subindex,&(OrigLoc[subindex]));
392  subindex = subnext;
393  }
394  index = next;
395  }
396 
397  }
398 
399  /* Sort coordinates and additional array accordingly (if provided). This is an expert routine borrowed from ML. Do not change.*/
400  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
401  void LineDetectionFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ML_az_dsort2(Scalar dlist[], LocalOrdinal N, LocalOrdinal list2[]) const {
402  LO l, r, j, i, flag;
403  LO RR2;
404  SC dRR, dK;
405 
406  // note: we use that routine for sorting coordinates only. No complex coordinates are assumed...
407  typedef Teuchos::ScalarTraits<SC> STS;
408 
409  if (N <= 1) return;
410 
411  l = N / 2 + 1;
412  r = N - 1;
413  l = l - 1;
414  dRR = dlist[l - 1];
415  dK = dlist[l - 1];
416 
417  if (list2 != NULL) {
418  RR2 = list2[l - 1];
419  while (r != 0) {
420  j = l;
421  flag = 1;
422 
423  while (flag == 1) {
424  i = j;
425  j = j + j;
426 
427  if (j > r + 1)
428  flag = 0;
429  else {
430  if (j < r + 1)
431  if (STS::real(dlist[j]) > STS::real(dlist[j - 1])) j = j + 1;
432 
433  if (STS::real(dlist[j - 1]) > STS::real(dK)) {
434  dlist[ i - 1] = dlist[ j - 1];
435  list2[i - 1] = list2[j - 1];
436  }
437  else {
438  flag = 0;
439  }
440  }
441  }
442  dlist[ i - 1] = dRR;
443  list2[i - 1] = RR2;
444 
445  if (l == 1) {
446  dRR = dlist [r];
447  RR2 = list2[r];
448  dK = dlist[r];
449  dlist[r ] = dlist[0];
450  list2[r] = list2[0];
451  r = r - 1;
452  }
453  else {
454  l = l - 1;
455  dRR = dlist[ l - 1];
456  RR2 = list2[l - 1];
457  dK = dlist[l - 1];
458  }
459  }
460  dlist[ 0] = dRR;
461  list2[0] = RR2;
462  }
463  else {
464  while (r != 0) {
465  j = l;
466  flag = 1;
467  while (flag == 1) {
468  i = j;
469  j = j + j;
470  if (j > r + 1)
471  flag = 0;
472  else {
473  if (j < r + 1)
474  if (STS::real(dlist[j]) > STS::real(dlist[j - 1])) j = j + 1;
475  if (STS::real(dlist[j - 1]) > STS::real(dK)) {
476  dlist[ i - 1] = dlist[ j - 1];
477  }
478  else {
479  flag = 0;
480  }
481  }
482  }
483  dlist[ i - 1] = dRR;
484  if (l == 1) {
485  dRR = dlist [r];
486  dK = dlist[r];
487  dlist[r ] = dlist[0];
488  r = r - 1;
489  }
490  else {
491  l = l - 1;
492  dRR = dlist[ l - 1];
493  dK = dlist[l - 1];
494  }
495  }
496  dlist[ 0] = dRR;
497  }
498 
499  }
500 } //namespace MueLu
501 
502 #endif // MUELU_LINEDETECTIONFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
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.
User data are always kept. This flag is set automatically when Level::Set("data", data) is used...
LocalOrdinal LO
Namespace for MueLu classes and methods.
void sort_coordinates(LO numCoords, LO *OrigLoc, Scalar *xvals, Scalar *yvals, Scalar *zvals, Scalar *xtemp, Scalar *ytemp, Scalar *ztemp, bool flipXY=false) const
static const NoFactory * get()
Additional warnings.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
LO ML_compute_line_info(LO LayerId[], LO VertLineId[], LO Ndof, LO DofsPerNode, LO MeshNumbering, LO NumNodesPerVertLine, SC *xvals, SC *yvals, SC *zvals, const Teuchos::Comm< int > &comm) const
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
void ML_az_dsort2(SC dlist[], LO N, LO list2[]) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Scalar SC
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()