MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 sort_coordinates(LO numCoords, LO *OrigLoc, Scalar *xvals, Scalar *yvals, Scalar *zvals, Scalar *xtemp, Scalar *ytemp, Scalar *ztemp, bool flipXY=false) const
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
T * getRawPtr() const
User data are always kept. This flag is set automatically when Level::Set(&quot;data&quot;, data) is used...
LocalOrdinal LO
static const NoFactory * get()
void Build(Level &currentLevel) const
Build method.
Additional warnings.
T * getRawPtr() const
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
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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 ML_az_dsort2(SC dlist[], LO N, LO list2[]) const
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()
void DeclareInput(Level &currentLevel) const
Input.