From 4c2cac1858a026df776925357fdb4c6992df497a Mon Sep 17 00:00:00 2001
From: Ansgar Burchardt <Ansgar.Burchardt@tu-dresden.de>
Date: Thu, 31 Aug 2017 18:37:21 +0200
Subject: [PATCH] remove large parts of support for reordering degrees of
freedom
---
gm/algebra.cc | 1769 -------------------------------------------------
gm/gm.h | 11 -
2 files changed, 1780 deletions(-)
diff --git a/gm/algebra.cc b/gm/algebra.cc
index da019aa68..91b300a38 100644
--- a/gm/algebra.cc
+++ b/gm/algebra.cc
@@ -113,23 +113,10 @@ USING_UGDIM_NAMESPACE
/** \brief For GetDomainPart indicating an element is meant rather than an element side */
#define NOSIDE -1
-/** \brief For LexOrderVectorsInGrid */
-#define MATTABLESIZE 32
-
/** \brief Constants for the direction of domain halfening */
#define BV_VERTICAL 0
#define BV_HORIZONTAL 1
-/** \brief For ordering of matrices */
-#define MATHPOS +1
-#define MATHNEG -1
-
-/** \brief Temp vector flag for LineOrderVectors */
-static INT ce_VCSTRONG;
-#define VCSTRONG_LEN 1
-#define VCSTRONG(p) CW_READ(p,ce_VCSTRONG)
-#define SETVCSTRONG(p,n) CW_WRITE(p,ce_VCSTRONG,n)
-
/****************************************************************************/
/* */
/* data structures used in this source file (exported data structures are */
@@ -157,8 +144,6 @@ static INT theAlgDepVarID; /* env type for ALG_DEP vars */
static INT theFindCutDirID; /* env type for FindCut dir */
static INT theFindCutVarID; /* env type for FIND_CUT vars */
-static FindCutProcPtr FindCutSet; /* pointer to find cut procedure */
-
#ifdef __TWODIM__
static MULTIGRID *GBNV_mg; /* multigrid */
static INT GBNV_MarkKey; /* key for Mark/Release */
@@ -280,8 +265,6 @@ const BV_DESC_FORMAT NS_DIM_PREFIX three_level_bvdf =
/****************************************************************************/
/* for LexOrderVectorsInGrid */
-static const INT *Order,*Sign;
-static INT SkipV,SignRad;
static DOUBLE InvMeshSize;
REP_ERR_FILE
@@ -858,44 +841,6 @@ INT NS_DIM_PREFIX CreateBlockvector_l0 (GRID *theGrid, BLOCKVECTOR **BVHandle, B
return (GM_OK);
}
-static INT CutBlockvector_l0 (GRID *theGrid, BLOCKVECTOR *theBV, INT makeVC)
-{
- if (theBV==NULL) return (GM_ERROR);
-
- switch ((theBV==GFIRSTBV(theGrid)) | ((theBV==GLASTBV(theGrid))<<1))
- {
- case 0 :
- BVSUCC(BVPRED(theBV)) = BVSUCC(theBV);
- BVPRED(BVSUCC(theBV)) = BVPRED(theBV);
-
- if (makeVC==YES)
- {
- SUCCVC(BVLASTVECTOR(BVPRED(theBV))) = BVFIRSTVECTOR(BVSUCC(theBV));
- PREDVC(BVFIRSTVECTOR(BVSUCC(theBV))) = BVLASTVECTOR(BVPRED(theBV));
- }
- break;
- case 1 :
- BVPRED(BVSUCC(theBV)) = NULL;
- GFIRSTBV(theGrid) = BVSUCC(theBV);
-
- if (makeVC==YES)
- PREDVC(BVFIRSTVECTOR(GFIRSTBV(theGrid))) = NULL;
- break;
- case 2 :
- BVSUCC(BVPRED(theBV)) = NULL;
- GLASTBV(theGrid) = BVPRED(theBV);
-
- if (makeVC==YES)
- SUCCVC(BVLASTVECTOR(BVPRED(theBV))) = NULL;
- break;
- case 3 :
- GFIRSTBV(theGrid) = GLASTBV(theGrid) = NULL;
- break;
- }
-
- return (GM_OK);
-}
-
/****************************************************************************/
/** \brief Return pointer to a new connection structure
*
@@ -4371,285 +4316,6 @@ INT NS_DIM_PREFIX MaxNextVectorClass (GRID *theGrid, ELEMENT *theElement)
return (m);
}
-/****************************************************************************/
-/* */
-/* Function: LexCompare */
-/* */
-/****************************************************************************/
-
-static INT LexCompare (VECTOR **pvec1, VECTOR **pvec2)
-{
- DOUBLE_VECTOR pv1,pv2;
- DOUBLE diff[DIM];
-
- PRINTDEBUG(gm,1,("%d: LexCompare %4d %4d",
- me,VINDEX(*pvec1),VINDEX(*pvec2)));
-
- if (SkipV)
- {
- PRINTDEBUG(gm,1,(" vecskip %4d %4d\n",
- VECSKIP(*pvec1),VECSKIP(*pvec2)));
- if (VECSKIP(*pvec1) && !VECSKIP(*pvec2))
- {
- if (SkipV==GM_PUT_AT_BEGIN)
- return (-1);
- else
- return ( 1);
- }
-
- if (VECSKIP(*pvec2) && !VECSKIP(*pvec1))
- {
- if (SkipV==GM_PUT_AT_BEGIN)
- return ( 1);
- else
- return (-1);
- }
- }
- VectorPosition(*pvec1,pv1);
- VectorPosition(*pvec2,pv2);
-
- V_DIM_SUBTRACT(pv2,pv1,diff);
- V_DIM_SCALE(InvMeshSize,diff);
-
- #ifdef __THREEDIM__
- PRINTDEBUG(gm,1,(" diff %f %f %f\n",diff[0],diff[1],diff[2]));
- #endif
-
- if (fabs(diff[Order[DIM-1]])<=ORDERRES)
- {
- #ifdef __THREEDIM__
- if (fabs(diff[Order[DIM-2]])<=ORDERRES)
- {
- if (diff[Order[DIM-3]]>0.0) return (-Sign[DIM-3]);
- else return ( Sign[DIM-3]);
- }
- else
- #endif
- if (diff[Order[DIM-2]]>0.0) return (-Sign[DIM-2]);
- else return ( Sign[DIM-2]);
- }
- else
- {
- if (diff[Order[DIM-1]]>0.0) return (-Sign[DIM-1]);
- else return ( Sign[DIM-1]);
- }
-}
-
-static INT PolarLexCompare (VECTOR **pvec1, VECTOR **pvec2)
-{
- DOUBLE_VECTOR pv1,pv2;
- DOUBLE norm1,norm2,s1,s2,c1,c2;
- DOUBLE_VECTOR polar;
-
- PRINTDEBUG(gm,1,("%d: LexCompare %4d %4d",
- me,VINDEX(*pvec1),VINDEX(*pvec2)));
-
- if (SkipV)
- {
- PRINTDEBUG(gm,1,(" vecskip %4d %4d\n",
- VECSKIP(*pvec1),VECSKIP(*pvec2)));
- if (VECSKIP(*pvec1) && !VECSKIP(*pvec2))
- {
- if (SkipV==GM_PUT_AT_BEGIN)
- return (-1);
- else
- return ( 1);
- }
-
- if (VECSKIP(*pvec2) && !VECSKIP(*pvec1))
- {
- if (SkipV==GM_PUT_AT_BEGIN)
- return ( 1);
- else
- return (-1);
- }
- }
- VectorPosition(*pvec1,pv1);
- VectorPosition(*pvec2,pv2);
-
- /* first check if one pos is (0,0) */
- V_DIM_EUKLIDNORM(pv1,norm1);
- if (fabs(norm1)<=SMALL_C)
- return (-SignRad);
- V_DIM_EUKLIDNORM(pv2,norm2);
- if (fabs(norm2)<=SMALL_C)
- return ( SignRad);
-
- /* determine radius difference polar[0] and polar angle difference polar[1] */
-
- /* radial difference between pv1 and pv2 */
- polar[0] = norm1-norm2;
-
- /* sin and cos of polar angle between pv1 and unit x-vector */
- V_DIM_SCALE(1./norm1,pv1);
- V_DIM_SCALE(1./norm2,pv2);
- V_DIM_SCALAR_PRODUCT(pv1,unit_vec[1],s1);
- V_DIM_SCALAR_PRODUCT(pv2,unit_vec[1],s2);
- if (s1*s2>=0)
- {
- V_DIM_SCALAR_PRODUCT(pv1,unit_vec[0],c1);
- V_DIM_SCALAR_PRODUCT(pv2,unit_vec[0],c2);
- polar[1] = (c1-c2);
- if (s1==0)
- polar[1] *= SIGNUM(s2);
- else
- polar[1] *= SIGNUM(s1);
- }
- else
- polar[1] = s1-s2;
-
- if (fabs(polar[Order[DIM-1]])<=ORDERRES)
- {
- if (polar[Order[DIM-2]]>0.0) return (-Sign[DIM-2]);
- else return ( Sign[DIM-2]);
- }
- else
- {
- if (polar[Order[DIM-1]]>0.0) return (-Sign[DIM-1]);
- else return ( Sign[DIM-1]);
- }
-}
-
-static int MatrixCompare (MATRIX **MatHandle1, MATRIX **MatHandle2)
-{
- INT IND1,IND2;
-
- IND1 = VINDEX(MDEST(*MatHandle1));
- IND2 = VINDEX(MDEST(*MatHandle2));
-
- if (IND1>IND2)
- return ( 1);
- else
- return (-1);
-}
-
-
-/****************************************************************************/
-/** \brief Order vectors lexicographically
-
- * @param theGrid - grid level
- * @param order - hierarchie of ordering directions (x,y[,z])
- * @param sign - signs for the directions
- * @param SpecSkipVecs - if true: GM_PUT_AT_BEGIN or GM_PUT_AT_END of vectors with skip true
- * @param AlsoOrderMatrices - if true order matrices in the same sense
-
- This function orders the vectors of one level lexicographically.
- It has the complexity of qsort which is n*log(n).
-
- * @return <ul>
- * <li> 0 if ok </li>
- * <li> 1 if error occured. </li>
- </ul>
- */
-/****************************************************************************/
-
-INT NS_DIM_PREFIX LexOrderVectorsInGrid (GRID *theGrid, INT mode, const INT *order, const INT *sign, INT which, INT SpecSkipVecs, INT AlsoOrderMatrices)
-{
- MULTIGRID *theMG;
- VECTOR **table,*theVec;
- MATRIX *theMat,**MatTable;
- INT i,entries,nm;
- HEAP *theHeap;
- INT takeSkip, takeNonSkip;
- INT MarkKey,nn;
-
- theMG = MYMG(theGrid);
- entries = NVEC(theGrid);
-
- /* calculate the diameter of the bounding rectangle of the domain */
- nn=NN(GRID_ON_LEVEL(theMG,0));
-#ifdef ModelP
- nn=UG_GlobalSumINT(nn);
-#endif
- // The following method wants the domain radius, which has been removed.
- // Dune has been setting this radius to 1.0 for years now, so I don't think it matters.
- DOUBLE BVPD_RADIUS = 1.0;
- InvMeshSize = POW2(GLEVEL(theGrid)) * pow(nn,1.0/DIM) / BVPD_RADIUS;
- assert(InvMeshSize>0.0);
-
- /* allocate memory for the node list */
- if (which==0) return (99);
- takeSkip = which & GM_TAKE_SKIP;
- takeNonSkip = which & GM_TAKE_NONSKIP;
- entries = 0;
- for (theVec=FIRSTVECTOR(theGrid); theVec!=NULL; theVec=SUCCVC(theVec))
- if ((takeSkip && VECSKIP(theVec)) || (takeNonSkip && !VECSKIP(theVec)))
- entries++;
- if (entries < 2) return(0);
- theHeap = MGHEAP(theMG);
- MarkTmpMem(theHeap,&MarkKey);
- if ((table = (VECTOR **)GetTmpMem(theHeap,entries*sizeof(VECTOR *),MarkKey))==NULL)
- {
- ReleaseTmpMem(theHeap,MarkKey);
- PrintErrorMessage('E',"LexOrderVectorsInGrid",
- "could not allocate memory from the MGHeap");
- return (2);
- }
- /* fill array of pointers to nodes */
- entries = 0;
- for (theVec=FIRSTVECTOR(theGrid); theVec!=NULL; theVec=SUCCVC(theVec))
- if ((takeSkip && VECSKIP(theVec))||(takeNonSkip && !VECSKIP(theVec))) {
- VINDEX(theVec) = entries;
- table[entries++] = theVec;
- }
- /* sort array of pointers */
- Order = order;
- Sign = sign;
- SkipV = SpecSkipVecs;
-
- if (mode==OV_POLAR)
- {
- if (Order[0]==0)
- SignRad = Sign[0];
- else
- SignRad = Sign[1];
- qsort(table,entries,sizeof(*table),
- (int (*)(const void *, const void *))PolarLexCompare);
- }
- else
- {
- qsort(table,entries,sizeof(*table),
- (int (*)(const void *, const void *))LexCompare);
- }
-
- for (i=0; i<entries; i++)
- GRID_UNLINK_VECTOR(theGrid,table[i]);
-
- for (i=0; i<entries; i++) {
- VINDEX(table[i]) = i;
- GRID_LINK_VECTOR(theGrid,table[i],PRIO(table[i]));
- }
- if (!AlsoOrderMatrices) {
- ReleaseTmpMem(theHeap,MarkKey);
- return (0);
- }
- MatTable = (MATRIX **) table;
- /* now we also order the matrices of each vector the same way
- (just using the new INDICES) */
- /* but leaving the diag at its place of course */
- for (theVec=FIRSTVECTOR(theGrid); theVec!=NULL; theVec=SUCCVC(theVec))
- {
- /* fill array for qsort */
- for (nm=0, theMat=VSTART(theVec); theMat!=NULL; theMat=MNEXT(theMat))
- MatTable[nm++] = theMat;
- if (nm < 2) continue;
- qsort(MatTable+1,nm-1,sizeof(*MatTable),
- (int (*)(const void *, const void *))MatrixCompare);
-
- /* establish pointer connections */
- MNEXT(MatTable[--nm]) = NULL;
- while (nm>0)
- {
- MNEXT(MatTable[nm-1]) = MatTable[nm];
- --nm;
- }
- VSTART(theVec) = MatTable[0];
- }
- ReleaseTmpMem(theHeap,MarkKey);
-
- return (0);
-}
-
/****************************************************************************/
/** \brief Create a new find cut procedure in environment
@@ -4745,1441 +4411,6 @@ static VECTOR *FeedbackVertexVectors (GRID *theGrid, VECTOR *LastVector, INT *nb
}
-
-#ifdef __TWODIM__
-
-static DOUBLE_VECTOR mypos;
-static INT sense;
-
-static INT SensCompare (MATRIX **MatHandle1, MATRIX **MatHandle2)
-{
- DOUBLE dx1,dy1,dx2,dy2;
- DOUBLE_VECTOR pos1,pos2;
-
- VectorPosition(MDEST(*MatHandle1),pos1);
- VectorPosition(MDEST(*MatHandle2),pos2);
- dx1 = pos1[0] - mypos[0];
- dy1 = pos1[1] - mypos[1];
- dx2 = pos2[0] - mypos[0];
- dy2 = pos2[1] - mypos[1];
-
- if (dy1>=0)
- {
- if (dy2<0)
- return (-sense);
-
- if ((dy1==0) && (dy2==0))
- {
- if (dx1>dx2)
- return (-sense);
- else
- return (sense);
- }
-
- /* both neighbours in upper half plane */
- if ((dy1*dx2-dx1*dy2)<0)
- return (-sense);
- else
- return (sense);
- }
- else
- {
- if (dy2>=0)
- return (sense);
-
- /* both neighbours in lower half plane */
- if ((dy1*dx2-dx1*dy2)<0)
- return (-sense);
- else
- return (sense);
- }
-}
-
-/****************************************************************************/
-/** \brief Reorder the matrix list of a vector in a circular order (2D only)
-
- * @param vec - order matrix list of this vector
- * @param sense - MATHPOS or MATHNEG
-
- This function reorders the matrix list of a vector in a circular order (2D only).
-
- * @return <ul>
- * <li> 0 if ok </li>
- * <li> 1 if error occured. </li>
- </ul>
- */
-/****************************************************************************/
-static INT OrderMatrices (VECTOR *vec, INT Sense)
-{
- VECTOR *nbv;
- MATRIX *mat,*MatList[MATTABLESIZE];
- INT i,flag,condition1,condition2,start,nm;
-
- /* fill matrices without diag into list */
- for (nm=0, mat=MNEXT(VSTART(vec)); mat!=NULL; mat=MNEXT(mat))
- {
- if (nm>=MATTABLESIZE)
- return (1);
-
- MatList[nm++] = mat;
- }
- VectorPosition(vec,mypos);
-
- sense = Sense;
- qsort(MatList,nm,sizeof(MATRIX*),(int (*)(const void *, const void *))SensCompare);
-
- /* find 'begin */
- flag = false;
- for (start=0; start<nm; start++)
- {
- nbv = MDEST(MatList[start]);
- condition1 = ((VCLASS(nbv)<ACTIVE_CLASS) || VCUSED(nbv));
- condition2 = (condition1 || (OBJT(MYVERTEX((NODE*)VOBJECT(nbv)))==BVOBJ));
- if (flag && !condition1)
- break;
- else if (condition2)
- flag = true;
- }
-
- /* establish pointer connections */
- for (i=start; i<nm+start-1; i++)
- MNEXT(MatList[i%nm]) = MatList[(i+1)%nm];
- MNEXT(MatList[i%nm]) = NULL;
- MNEXT(VSTART(vec)) = MatList[start%nm];
-
- return (0);
-}
-
-#endif
-
-/****************************************************************************/
-/** \brief Reorder double linked vector list by a shell algorithm
-
- * @param theGrid - pointer to grid
- * @param seed - start at this vector
-
- This function reorders double linked vector list by a shell algorithm.
-
- * @return <ul>
- * <li> 0 if ok </li>
- * <li> 1 if error occured. </li>
- </ul>
- */
-/****************************************************************************/
-
-INT NS_DIM_PREFIX ShellOrderVectors (GRID *theGrid, VECTOR *seed)
-{
- FIFO myfifo;
- void *buffer;
- VECTOR **vlist;
- VECTOR *theV;
- MATRIX *theM;
- HEAP *theHeap=MGHEAP(MYMG(theGrid));
- INT i,n;
- INT MarkKey;
-
- /* count vectors */
- n = 0;
- for (theV=FIRSTVECTOR(theGrid); theV!=NULL; theV=SUCCVC(theV))
- n++;
-
- if (n == 0)
- return(0);
-
- /* find new ordering */
- MarkTmpMem(theHeap,&MarkKey);
- buffer=(void *)GetTmpMem(theHeap,sizeof(VECTOR*)*n,MarkKey);
- vlist = (VECTOR**)GetTmpMem(theHeap,sizeof(VECTOR*)*n,MarkKey);
- fifo_init(&myfifo,buffer,sizeof(VECTOR*)*n);
- for (theV=FIRSTVECTOR(theGrid); theV!=NULL; theV=SUCCVC(theV))
- SETVCUSED(theV,0);
- fifo_in(&myfifo,(void *)seed);
- SETVCUSED(seed,1);
- i = 0;
- while(!fifo_empty(&myfifo))
- {
- theV = (VECTOR *)fifo_out(&myfifo);
- vlist[i++] = theV;
-
- #ifdef __TWODIM__
- OrderMatrices(theV,MATHPOS);
- #endif
-
- for (theM=MNEXT(VSTART(theV)); theM!=NULL; theM=MNEXT(theM))
- if (!CEXTRA(MMYCON(theM)))
- if (!VCUSED(MDEST(theM)))
- {
- fifo_in(&myfifo,(void *)MDEST(theM));
- SETVCUSED(MDEST(theM),1);
- }
- }
- assert(i==n);
-
- /* reorder vector list in grid */
- for (i=0; i<n; i++) GRID_UNLINK_VECTOR(theGrid,vlist[i]);
- for (i=0; i<n; i++) GRID_LINK_VECTOR(theGrid,vlist[i],PRIO(vlist[i]));
-
- ReleaseTmpMem(theHeap,MarkKey);
-
- return (0);
-}
-
-/****************************************************************************/
-/*
- OrderVectorAlgebraic - Reorder double linked node list by a streamline ordering
-
- SYNOPSIS:
- static INT OrderVectorAlgebraic (GRID *theGrid, INT mode, INT putSkipFirst, INT skipPat);
-
- PARAMETERS:
- * @param theGrid - pointer to grid
- * @param mode
- * @param putSkipFirst - if true put vectors with a skip pattern larger than SkipPat to begin of list
- . skipPat - s.a.
-
- DESCRIPTION:
- This function reorders double linked node list by a streamline ordering.
-
- * @return <ul>
- INT
- * <li> 0 if ok
- * <li> 1 if error occured.
- */
-/****************************************************************************/
-
-static INT OrderVectorAlgebraic (GRID *theGrid, INT mode, INT putSkipFirst, INT skipPat)
-{
- BLOCKVECTOR *theFirstBV, *theLastBV, *theCutBV, *theBV, *takeOut, **FBVList, **LBVList, **CBVList;
- VECTOR FIRST_handle,LAST_handle;
- VECTOR *FIRST_next_out,*FIRST_last_in;
- VECTOR *LAST_next_out,*LAST_last_in;
- VECTOR *CUT_begin;
- VECTOR *theVector,*theNbVector,*succVector,*predVector, *theFirstVector;
- MATRIX *theMatrix;
- DOUBLE a;
- INT i,k,cycle,nCutTot;
- INT nFIRST,nCUT,nLAST;
- INT up, down;
- HEAP *theHeap;
-
- #ifdef ModelP
- ASSERT(false); /* see TODO below */
- #endif
-
- /********************************************************************/
- /* init */
- /********************************************************************/
-
- /* cancel all BLOCKVECTORS */
- FreeAllBV(theGrid);
- SETUSED(theGrid,0);
-
- /* init USED, N_INFLOW and N_OUTFLOW */
- for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
- {
- /* reset used flag (indicating the vector hasn't yet its new place in the list) */
- SETVCUSED(theVector,0);
-
- /* count upward and downward matrices */
- up = down = 0;
- for (theMatrix=MNEXT(VSTART(theVector)); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
- {
- if (MUP(theMatrix)) up++;
- if (MDOWN(theMatrix)) down++;
- }
- SETVUP(theVector,up);
- SETVDOWN(theVector,down);
- }
-
- /* in the sequel we use (and therefore destroy) the PREDVC-list for book keeping */
-
- cycle = nCutTot = 0;
-
- /* init pointers and set the first FIRST and LAST set */
- FIRST_next_out = FIRST_last_in = &FIRST_handle;
- PREDVC(FIRST_last_in) = NULL; nFIRST = 0;
- LAST_next_out = LAST_last_in = &LAST_handle;
- PREDVC(LAST_last_in) = NULL; nLAST = 0;
- nCUT = 0;
- theFirstBV = theLastBV = theCutBV = NULL;
- for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
- {
- if (((putSkipFirst) && (VECSKIP(theVector) & skipPat) == skipPat) ||
- (VUP(theVector)==0))
- {
- /* append to FIRST list */
- nFIRST++;
- PREDVC(FIRST_last_in) = theVector;
- FIRST_last_in = theVector;
- SETVCUSED(theVector,1);
- VINDEX(theVector) = 0;
- }
- else if (VDOWN(theVector)==0)
- {
- /* append to last list */
- nLAST++;
- PREDVC(LAST_last_in) = theVector;
- LAST_last_in = theVector;
- SETVCUSED(theVector,1);
- VINDEX(theVector) = 1;
- }
- }
- PREDVC(FIRST_last_in) = PREDVC(LAST_last_in) = NULL;
-
- /* store in BLOCKVECTOR */
- if (nFIRST>0)
- {
- if (CreateBlockvector_l0(theGrid,&theBV,NULL,1)!=GM_OK) RETURN (1);
- theFirstBV = theBV;
- BVNUMBER(theFirstBV) = 3*cycle+0;
- BVFIRSTVECTOR(theFirstBV) = BVPRED(&FIRST_handle);
- }
-
- if (nLAST>0)
- {
- if (CreateBlockvector_l0(theGrid,&theBV,theFirstBV,1)!=GM_OK) RETURN (1);
- theLastBV = theBV;
- BVNUMBER(theLastBV) = 3*cycle+1;
- BVFIRSTVECTOR(theLastBV) = BVPRED(&LAST_handle);
- }
-
- do
- {
- cycle++;
-
- /****************************************************************************/
- /* find next FIRST-set in vectors not used */
- /****************************************************************************/
-
- nFIRST = nLAST = 0;
-
- theFirstVector = FIRST_last_in;
- for (FIRST_next_out=PREDVC(FIRST_next_out); FIRST_next_out!=NULL; FIRST_next_out=PREDVC(FIRST_next_out))
- for (theMatrix=MNEXT(VSTART(FIRST_next_out)); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
- {
- theNbVector = MDEST(theMatrix);
- if (MDOWN(theMatrix) && !VCUSED(theNbVector))
- {
- k = VUP(theNbVector);
- assert(k>0); /* if 0 is supposed to be VCUSED already */
- SETVUP(theNbVector,--k);
- if (k==0)
- {
- /* this vector has only matrices going down */
-
- nFIRST++;
- /* append to FIRST list */
- PREDVC(FIRST_last_in) = theNbVector;
- FIRST_last_in = theNbVector;
- PREDVC(FIRST_last_in) = NULL;
- SETVCUSED(theNbVector,1);
- VINDEX(theNbVector) = 3*cycle;
- }
- }
- if ((nCUT>0) && !VCUSED(theNbVector))
- /* we also have to push LAST vectors */
- if (MUP(theMatrix))
- {
- k = VDOWN(theNbVector);
- /*assert(k>0);*/ /* if 0 is supposed to be VCUSED already */
- if (k<=0)
- RETURN (1);
- SETVDOWN(theNbVector,--k);
- if (k==0)
- {
- /* this vector has only matrices going up */
-
- nLAST++;
- /* append to last list */
- PREDVC(LAST_last_in) = theNbVector;
- LAST_last_in = theNbVector;
- PREDVC(LAST_last_in) = NULL;
- SETVCUSED(theNbVector,1);
- VINDEX(theNbVector) = 3*cycle+1;
- }
- }
- }
- FIRST_next_out = FIRST_last_in;
-
- /* create first BLOCKVECTOR */
- if (nFIRST>0)
- {
- if (CreateBlockvector_l0(theGrid,&theBV,theFirstBV,1)!=GM_OK) RETURN (1);
- theFirstBV = theBV;
- BVNUMBER(theFirstBV) = 3*cycle+0;
- BVFIRSTVECTOR(theFirstBV) = PREDVC(theFirstVector);
- }
-
- /********************************************************************/
- /* find next LAST-set in vectors not used */
- /********************************************************************/
-
- for (LAST_next_out=PREDVC(LAST_next_out); LAST_next_out!=NULL; LAST_next_out=PREDVC(LAST_next_out))
- for (theMatrix=MNEXT(VSTART(LAST_next_out)); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
- {
- theNbVector = MDEST(theMatrix);
- if (!VCUSED(theNbVector) && MUP(theMatrix))
- {
- k = VDOWN(theNbVector);
- assert(k>0); /* if 0 is supposed to be VCUSED already */
- SETVDOWN(theNbVector,--k);
- if (k==0)
- {
- /* this vector has only matrices going up */
-
- nLAST++;
- /* append to last list */
- PREDVC(LAST_last_in) = theNbVector;
- LAST_last_in = theNbVector;
- PREDVC(LAST_last_in) = NULL;
- SETVCUSED(theNbVector,1);
- VINDEX(theNbVector) = 3*cycle+1;
- }
- }
- }
- LAST_next_out = LAST_last_in;
-
- /* create Last BLOCKVECTOR */
- if (nLAST>0)
- {
- if (CreateBlockvector_l0(theGrid,&theBV,theLastBV,0)!=GM_OK) RETURN (1);
- theLastBV = theBV;
- BVNUMBER(theLastBV) = 3*cycle+1;
- BVFIRSTVECTOR(theLastBV) = LAST_last_in;
- }
-
- /****************************************************************************/
- /* get CUT (or Feedback Vertex)-set and do what needs to be done */
- /****************************************************************************/
-
- CUT_begin = FIRST_last_in;
- FIRST_last_in = (*FindCutSet)(theGrid,FIRST_last_in,&nCUT);
- if (FIRST_last_in==NULL) nCUT = 0;
- else PREDVC(FIRST_last_in) = NULL;
-
- /* create cut BLOCKVECTOR */
- if (nCUT>0)
- {
- if (CreateBlockvector_l0(theGrid,&theBV,theFirstBV,1)!=GM_OK) RETURN (1);
- theCutBV = theBV;
- theFirstBV = theBV;
- BVNUMBER(theCutBV) = 3*(cycle+1)+2;
- BVFIRSTVECTOR(theCutBV) = PREDVC(CUT_begin);
- }
-
- /* set index */
- for (theVector=PREDVC(CUT_begin); theVector!=NULL; theVector=PREDVC(theVector))
- VINDEX(theVector) = 3*(cycle+1)+2;
-
- nCutTot += nCUT;
-
- } while (nCUT>0);
-
- UserWriteF("# %d cycles: %d cutted from %d\n",(int)cycle,(int)nCutTot,(int)NVEC(theGrid));
- a = POW((DOUBLE)NVEC(theGrid),(DOUBLE)(DIM-1)/(DOUBLE)DIM); a = (DOUBLE)nCutTot/a;
- UserWriteF("# corr. to %6.2f hyp. planes\n",(float)a);
-
- /** \todo use dlmgr macros for lists after streamwise ordering */
-
- /* insert FIRST list one-by-one to LASTVECTOR list of grid */
- LASTVECTOR(theGrid) = NULL;
- for (theVector=PREDVC(&FIRST_handle); theVector!=NULL; theVector=predVector)
- {
- predVector = PREDVC(theVector);
- PREDVC(theVector) = LASTVECTOR(theGrid);
- LASTVECTOR(theGrid) = theVector;
- }
-
- /* insert LAST list as a whole to LASTVECTOR list of grid */
- PREDVC(LAST_last_in) = LASTVECTOR(theGrid);
- LASTVECTOR(theGrid) = PREDVC(&LAST_handle);
-
- /* construct SUCCVC list */
- succVector = NULL;
- for (theVector=LASTVECTOR(theGrid); theVector!=NULL; theVector=PREDVC(theVector))
- {
- SUCCVC(theVector) = succVector;
- succVector = theVector;
- }
- SFIRSTVECTOR(theGrid) = succVector;
- PREDVC(succVector) = NULL;
-
- /* set pointers in BLOCKVECTORs */
- BVLASTVECTOR(GLASTBV(theGrid)) = LASTVECTOR(theGrid);
- for (theBV=GLASTBV(theGrid); theBV!=NULL; theBV=BVPRED(theBV))
- {
- if (BVSUCC(theBV)!=NULL && BVLASTVECTOR(theBV)==NULL)
- BVLASTVECTOR(theBV)=PREDVC(BVFIRSTVECTOR(BVSUCC(theBV)));
- if (BVFIRSTVECTOR(theBV)==NULL)
- {
- assert(0);
- BVFIRSTVECTOR(theBV)=BVENDVECTOR(theBV);
- }
- }
-
- if (mode==GM_FFLLCC)
- {
- for (theBV=GLASTBV(theGrid); BVPRED(theBV)!=NULL; theBV=BVPRED(theBV))
- {
- takeOut = BVPRED(theBV);
-
- if (BVNUMBER(takeOut)%3==2)
- {
- if (CutBlockvector_l0(theGrid,takeOut,YES))
- RETURN (1);
-
- if (InsertBlockvector_l0(theGrid,takeOut,NULL,0,YES))
- RETURN (1);
-
- theBV = BVSUCC(theBV);
- }
- }
- }
- else if (mode==GM_FFLCLC)
- {
- INT MarkKey;
-
- theHeap = MYMG(theGrid)->theHeap;
- MarkTmpMem(theHeap,&MarkKey);
- if ((FBVList=(BLOCKVECTOR **)GetTmpMem(theHeap,3*(cycle+1)*sizeof(BLOCKVECTOR*),MarkKey))==NULL) RETURN (1);
- for (i=0; i<3*cycle+3; i++) FBVList[i]=NULL;
- LBVList = FBVList + cycle + 1;
- CBVList = LBVList + cycle + 1;
-
- for (theBV=GFIRSTBV(theGrid); theBV!=NULL; theBV=BVSUCC(theBV))
- if (BVNUMBER(theBV)%3==0)
- FBVList[BVNUMBER(theBV)/3] = theBV;
- else if (BVNUMBER(theBV)%3==1)
- LBVList[BVNUMBER(theBV)/3] = theBV;
- else
- CBVList[BVNUMBER(theBV)/3] = theBV;
-
- for (i=cycle; i>0; i--)
- {
- if (CBVList[i]==NULL) continue;
- if (CutBlockvector_l0(theGrid,CBVList[i],YES)) RETURN (1);
- k=i-1; while (LBVList[k]==NULL && k>0) k--;
- if (InsertBlockvector_l0(theGrid,CBVList[i],LBVList[k],0,YES)) RETURN (1);
- }
-
- ReleaseTmpMem(theHeap,MarkKey);
- }
- else if (mode==GM_CCFFLL)
- {
- for (theBV=GFIRSTBV(theGrid); BVSUCC(theBV)!=NULL; theBV=BVSUCC(theBV))
- {
- takeOut = BVSUCC(theBV);
-
- if (BVNUMBER(takeOut)%3==2)
- {
- if (CutBlockvector_l0(theGrid,takeOut,YES))
- RETURN (1);
-
- if (InsertBlockvector_l0(theGrid,takeOut,NULL,1,YES))
- RETURN (1);
-
- theBV = BVPRED(theBV);
- }
- }
- }
-
-
- /* set VCCUT */
- for (theBV=GFIRSTBV(theGrid); theBV!=NULL; theBV=BVSUCC(theBV))
- if (BV_GEN(theBV)==BV_GEN_C)
- for (theVector=BVFIRSTVECTOR(theBV); theVector!=BVENDVECTOR(theBV); theVector=SUCCVC(theVector))
- {
- SETVCCUT(theVector,1);
- }
- else
- for (theVector=BVFIRSTVECTOR(theBV); theVector!=BVENDVECTOR(theBV); theVector=SUCCVC(theVector))
- {
- SETVCCUT(theVector,0);
- }
-
- /* check # members of succ list */
- i=0;
- for (theVector=FIRSTVECTOR(theGrid); theVector!= NULL; theVector=SUCCVC(theVector)) i++;
- if (NVEC(theGrid) != i)
- {
- UserWrite("vectorstructure corrupted\n");
- RETURN (1);
- }
-
- /* check # members of pred list */
- i=0;
- for (theVector=LASTVECTOR(theGrid); theVector!= NULL; theVector=PREDVC(theVector)) i++;
- if (NVEC(theGrid) != i)
- {
- UserWrite("vectorstructure corrupted\n");
- RETURN (1);
- }
-
- /* set index w.r.t. new order, beginning with 1 */
- i = 1;
- for (theVector=FIRSTVECTOR(theGrid); theVector!= NULL; theVector=SUCCVC(theVector))
- VINDEX(theVector) = i++;
-
-
- return (0);
-}
-
-/****************************************************************************/
-/** \brief Driver for general vector ordering
-
- * @param theMG - multigrid to order
- * @param levels - GM_ALL_LEVELS or GM_CURRENT_LEVEL
- * @param mode - GM_FCFCLL or GM_FFCCLL (see orderv command)
- * @param PutSkipFirst - if true put vectors with a skip pattern larger than SkipPat to begin of list
- * @param SkipPat - s.a.
- * @param dependency - name of user defined dependency item
- * @param dep_options - options for user dependency function
-
- This function orders VECTORs in a multigrid according to the dependency
- function provided.
-
- * @return <ul>
- * <li> GM_OK if ok
- * <li> GM_ERROR if error occured.
- </ul>
- */
-/****************************************************************************/
-
-INT NS_DIM_PREFIX OrderVectors (MULTIGRID *theMG, INT levels, INT mode, INT PutSkipFirst, INT SkipPat, const char *dependency, const char *dep_options, const char *findcut)
-{
- INT i, currlevel, baselevel;
- ALG_DEP *theAlgDep;
- FIND_CUT *theFindCut;
- GRID *theGrid;
- DependencyProcPtr DependencyProc;
-
- /* check mode */
- if ((mode!=GM_FCFCLL)&&
- (mode!=GM_FFLLCC)&&
- (mode!=GM_FFLCLC)&&
- (mode!=GM_CCFFLL)) RETURN(GM_ERROR);
-
- /* current level */
- currlevel = theMG->currentLevel;
-
- /* get algebraic dependency */
- theAlgDep=NULL;
- if (dependency!=NULL)
- {
- theAlgDep = (ALG_DEP *) SearchEnv(dependency,"/Alg Dep",theAlgDepVarID,theAlgDepDirID);
- if (theAlgDep==NULL)
- {
- UserWrite("algebraic dependency not found\n");
- RETURN(GM_ERROR);
- }
- DependencyProc = theAlgDep->DependencyProc;
- if (DependencyProc==NULL)
- {
- UserWrite("don't be stupid: implement a dependency!\n");
- RETURN(GM_ERROR);
- }
- }
-
- /* get find cut dependency */
- if (findcut==NULL)
- {
- FindCutSet = FeedbackVertexVectors;
- UserWrite("default cut set proc:\n leaving order of cyclic dependencies unchanged\n");
- }
- else
- {
- theFindCut = (FIND_CUT *) SearchEnv(findcut,"/FindCut",theFindCutVarID,theFindCutDirID);
- if (theFindCut==NULL)
- {
- UserWrite("find cut proc not found\n");
- RETURN(GM_ERROR);
- }
- FindCutSet = theFindCut->FindCutProc;
- if (FindCutSet==NULL)
- {
- UserWrite("don't be stupid: implement a find cut proc!\n");
- RETURN(GM_ERROR);
- }
- }
-
- /* go */
- if (levels==GM_ALL_LEVELS)
- baselevel = 0;
- else
- baselevel = currlevel;
- if (theAlgDep!=NULL)
- for (i=baselevel; i<=currlevel; i++)
- {
- theGrid = GRID_ON_LEVEL(theMG,i);
- if ((*DependencyProc)(theGrid,dep_options)) RETURN(GM_ERROR);
- }
- for (i=baselevel; i<=currlevel; i++)
- {
- theGrid = GRID_ON_LEVEL(theMG,i);
- if (OrderVectorAlgebraic(theGrid,mode,PutSkipFirst,SkipPat)) RETURN(GM_ERROR);
- }
-
- return (GM_OK);
-}
-
-static VECTOR *CutStrongLine (VECTOR *theVector, VECTOR *AppendVec)
-{
- MATRIX *mat;
-
- do
- {
- /* append to list */
- PREDVC(AppendVec) = theVector;
- AppendVec = theVector;
- PREDVC(AppendVec) = NULL;
- SETVCUSED(theVector,1);
-
- /* look for strong nb */
- for (mat=MNEXT(VSTART(theVector)); mat!=NULL; mat=MNEXT(mat))
- if (MSTRONG(mat))
- {
- theVector = MDEST(mat);
-
- if (!VCUSED(theVector))
- break;
- }
- }
- while (mat!=NULL);
-
- return (AppendVec);
-}
-
-static VECTOR *FindOptimalStrong (FIFO *fifo)
-{
- VECTOR *first,*vec,*nbvec,*minvec;
- MATRIX *mat;
- INT nStrong,minStrong,nUp,minUp;
-
- /* drop leading vectors used */
- while ((vec=(VECTOR*)fifo_out(fifo))!=NULL)
- if (!VCUSED(vec))
- break;
-
- if (vec==NULL)
- return (NULL);
-
- first = vec;
-
- /* determine minimal number of !used strong nbs */
- minStrong = MAX_I;
- do
- if (!VCUSED(vec))
- {
- nStrong = 0;
- for (mat=MNEXT(VSTART(vec)); mat!=NULL; mat=MNEXT(mat))
- if (MSTRONG(mat))
- {
- nbvec = MDEST(mat);
- if (!VCUSED(nbvec))
- nStrong++;
- }
- fifo_in(fifo,vec);
-
- minStrong = MIN(minStrong,nStrong);
- }
- while ((vec=(VECTOR*)fifo_out(fifo))!=first);
-
- /* determine vec with nStrong == minStrong and minimal number of up nbs */
- /* NB: fifo contains no VCUSED(vec) anymore */
- minUp = MAX_I;
- do
- {
- nUp = nStrong = 0;
- for (mat=MNEXT(VSTART(vec)); mat!=NULL; mat=MNEXT(mat))
- if (MSTRONG(mat))
- {
- nbvec = MDEST(mat);
- if (!VCUSED(nbvec))
- nStrong++;
- }
- else if (MUP(mat))
- {
- nbvec = MDEST(mat);
- if (!VCUSED(nbvec))
- nUp++;
- }
-
- fifo_in(fifo,vec);
-
- if (nStrong==minStrong)
- if (nUp<minUp)
- {
- minUp = nUp;
- minvec = vec;
- }
- }
- while ((vec=(VECTOR*)fifo_out(fifo))!=first);
-
- /* remove minvec from fifo */
- if (vec!=minvec)
- {
- do
- {
- if (vec!=minvec)
- fifo_in(fifo,vec);
- }
- while ((vec=(VECTOR*)fifo_out(fifo))!=first);
-
- /* push first again */
- fifo_in(fifo,vec);
- }
-
- return (minvec);
-}
-
-/****************************************************************************/
-/** \brief Reorder double linked vector list by a streamline ordering
-
- * @param theGrid - pointer to grid
-
- This function reorders double linked vector list by a streamline ordering.
-
- CAUTION:
- The VCUSED flag has to be initialized by the AlgDep procedure. If VCUSED==1
- then the vector is pushed to the FIRST fifo at initialization.
-
- The FindCutProc can only use the SUCCVC-list (the other one is destroyed).
-
- * @return <ul>
- * <li> 0 if ok
- * <li> 1 if error occured.
- </ul>
- */
-/****************************************************************************/
-
-static INT LineOrderVectorsAlgebraic (GRID *theGrid, INT verboselevel)
-{
- BLOCKVECTOR *theFirstBV, *theLastBV, *theCutBV, *theBV;
- FIFO FFifo,LFifo;
- VECTOR FIRST_handle,LAST_handle;
- VECTOR *FIRST_next_out,*FIRST_last_in;
- VECTOR *FIRST_line_start,*FIRST_line_end;
- VECTOR *LAST_line_end;
- VECTOR *LAST_next_out,*LAST_last_in;
- VECTOR *theVector,*theNbVector,*succVector,*predVector;
- MATRIX *theMatrix;
- DOUBLE a;
- INT i,k,init,cycle,line,nCutTot,fifosize;
- INT nCUT;
- INT up, down, strong;
- INT StrongInflow,nInflow,StrongOutflow,nOutflow,pushInflow,pushOutflow;
- INT bvn;
- char gen_label[3];
- INT MarkKey;
-
- #ifdef ModelP
- ASSERT(false); /* see TODO below */
- #endif
-
- gen_label[GM_GEN_FIRST] = 'F';
- gen_label[GM_GEN_LAST] = 'L';
- gen_label[GM_GEN_CUT] = 'C';
-
- /********************************************************************/
- /* init */
- /********************************************************************/
-
- /* dispose all BLOCKVECTORS */
- FreeAllBV(theGrid);
-
- /* init fifos */
- MarkTmpMem(MGHEAP(MYMG(theGrid)),&MarkKey);
- fifosize = 30*(INT)floor(sqrt(NVEC(theGrid)))*sizeof(VECTOR*);
- fifo_init(&FFifo,GetTmpMem(MGHEAP(MYMG(theGrid)),fifosize,MarkKey),fifosize);
- fifo_init(&LFifo,GetTmpMem(MGHEAP(MYMG(theGrid)),fifosize,MarkKey),fifosize);
-
- /* init USED, N_INFLOW and N_OUTFLOW */
- StrongInflow = StrongOutflow = 0;
- nInflow = nOutflow = 0;
- for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
- {
- SETVCSTRONG(theVector,0);
-
- /* reorder matrix lists in math pos sense */
- #ifdef __TWODIM__
- OrderMatrices(theVector,MATHNEG);
- #endif
-
- /* count upward, downward and strong matrices */
- up = down = strong = 0;
- for (theMatrix=MNEXT(VSTART(theVector)); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
- {
- if (MUP(theMatrix)) up++;
- if (MDOWN(theMatrix)) down++;
- if (MSTRONG(theMatrix)) strong++;
- }
- SETVUP(theVector,up);
- SETVDOWN(theVector,down);
- if (strong)
- SETVCSTRONG(theVector,1);
-
- if (VCUSED(theVector))
- {
- nInflow++;
- if (strong)
- StrongInflow++;
- }
- if (VCFLAG(theVector))
- {
- nOutflow++;
- if (strong)
- StrongOutflow++;
- }
- }
- pushInflow = (nInflow ==StrongInflow);
- pushOutflow = (nOutflow==StrongOutflow);
-
- /* in the sequel we use (and therefore destroy) the PREDVC-list for book keeping */
-
- cycle = nCutTot = 0;
- init = 1;
-
- /* init pointers and set the first FIRST and LAST set */
- FIRST_next_out = FIRST_last_in = &FIRST_handle;
- PREDVC(FIRST_last_in) = NULL;
- LAST_next_out = LAST_last_in = &LAST_handle;
- PREDVC(LAST_last_in) = NULL;
- nCUT = 0;
- theFirstBV = theLastBV = theCutBV = NULL;
- for (theVector=FIRSTVECTOR(theGrid); theVector!=NULL; theVector=SUCCVC(theVector))
- {
- if (VCUSED(theVector) && pushInflow)
- {
- if (fifo_in(&FFifo,theVector)!=0)
- {
- PrintErrorMessage('E',"LineOrderVectorsAlgebraic","fifo full");
- return (__LINE__);
- }
- SETVCFLAG(theVector,0);
- SETVCUSED(theVector,0);
- continue;
- }
- if (VCFLAG(theVector) && pushOutflow)
- {
- if (fifo_in(&LFifo,theVector)!=0)
- {
- PrintErrorMessage('E',"LineOrderVectorsAlgebraic","fifo full");
- return (__LINE__);
- }
- SETVCFLAG(theVector,0);
- SETVCUSED(theVector,0);
- continue;
- }
- SETVCFLAG(theVector,VCSTRONG(theVector));
- SETVCUSED(theVector,0);
- if (VUP(theVector)==0)
- {
- /* append to FIRST list */
- PREDVC(FIRST_last_in) = theVector;
- FIRST_last_in = theVector;
- SETVCUSED(theVector,1);
- }
- else if (VDOWN(theVector)==0)
- {
- /* append to LAST list */
- PREDVC(LAST_last_in) = theVector;
- LAST_last_in = theVector;
- SETVCUSED(theVector,1);
- }
- }
- PREDVC(FIRST_last_in) = PREDVC(LAST_last_in) = NULL;
-
- do
- {
- cycle++;
-
- /****************************************************************************/
- /* find next FIRST-set in vectors not used */
- /****************************************************************************/
-
- line = 0;
- do
- {
- FIRST_line_start = PREDVC(FIRST_next_out);
- FIRST_line_end = FIRST_last_in;
- for (FIRST_next_out=PREDVC(FIRST_next_out); FIRST_next_out!=NULL; FIRST_next_out=PREDVC(FIRST_next_out))
- {
- for (theMatrix=MNEXT(VSTART(FIRST_next_out)); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
- {
- theNbVector = MDEST(theMatrix);
- if (VCUSED(theNbVector)) continue;
-
- if (VCFLAG(theNbVector))
- {
- if (fifo_in(&FFifo,theNbVector)!=0)
- {
- PrintErrorMessage('E',"LineOrderVectorsAlgebraic","fifo full");
- return (__LINE__);
- }
- SETVCFLAG(theNbVector,0);
- }
- else if (MDOWN(theMatrix))
- {
- k = VUP(theNbVector);
- assert(k>0); /* if 0 is supposed to be VCUSED already */
- SETVUP(theNbVector,--k);
- if (k==0)
- {
- /* this vector has only matrices going down */
-
- /* append to FIRST list */
- PREDVC(FIRST_last_in) = theNbVector;
- FIRST_last_in = theNbVector;
- PREDVC(FIRST_last_in) = NULL;
- SETVCUSED(theNbVector,1);
- }
- }
- else if (MUP(theMatrix))
- {
- k = VDOWN(theNbVector);
- /*assert(k>0);*/ /* if 0 is supposed to be VCUSED already */
- if (k<=0)
- {
- PrintErrorMessage('E',"LineOrderVectorsAlgebraic","DOWN counter < 0");
- return (__LINE__);
- }
- SETVDOWN(theNbVector,--k);
- if (k==0)
- {
- /* this vector has only matrices going up */
-
- /* append to last list */
- PREDVC(LAST_last_in) = theNbVector;
- LAST_last_in = theNbVector;
- PREDVC(LAST_last_in) = NULL;
- SETVCUSED(theNbVector,1);
- }
- }
- }
- if (FIRST_next_out==FIRST_line_end)
- {
- /* create block line */
- if (CreateBlockvector_l0(theGrid,&theBV,theFirstBV,1)!=GM_OK)
- {
- PrintErrorMessage('E',"LineOrderVectorsAlgebraic","CreateBlockvector_l0 failed");
- return (__LINE__);
- }
- theFirstBV = theBV;
- if (nCUT>0)
- {
- SET_ORD_GCL(BVNUMBER(theFirstBV),GM_GEN_CUT,cycle,line);
- nCUT = 0;
- }
- else if (init && (line==0))
- {
- SET_ORD_GCL(BVNUMBER(theFirstBV),GM_GEN_FIRST,0,line);
- }
- else
- {
- SET_ORD_GCL(BVNUMBER(theFirstBV),GM_GEN_FIRST,cycle,line);
- }
- BVFIRSTVECTOR(theFirstBV) = FIRST_line_start;
-
- line++;
- FIRST_line_start = PREDVC(FIRST_next_out);
- FIRST_line_end = FIRST_last_in;
- }
- }
- FIRST_next_out = FIRST_last_in;
-
- /* check if there is a !used vec in the strong fifo */
- if ((theVector=FindOptimalStrong(&FFifo))!=NULL)
- FIRST_last_in = CutStrongLine(theVector,FIRST_last_in);
- }
- while (theVector!=NULL);
-
- /********************************************************************/
- /* find next LAST-set in vectors not used */
- /********************************************************************/
-
- line = 0;
- do
- {
- LAST_line_end = LAST_last_in;
- for (LAST_next_out=PREDVC(LAST_next_out); LAST_next_out!=NULL; LAST_next_out=PREDVC(LAST_next_out))
- {
- for (theMatrix=MNEXT(VSTART(LAST_next_out)); theMatrix!=NULL; theMatrix=MNEXT(theMatrix))
- {
- theNbVector = MDEST(theMatrix);
- if (VCUSED(theNbVector)) continue;
-
- if (VCFLAG(theNbVector))
- {
- /* this is a strongly coupled neighbour */
- if (fifo_in(&LFifo,theNbVector)!=0)
- return (__LINE__);
- SETVCFLAG(theNbVector,0);
- }
- else if (MUP(theMatrix))
- {
- k = VDOWN(theNbVector);
- assert(k>0); /* if 0 is supposed to be VCUSED already */
- SETVDOWN(theNbVector,--k);
- if (k==0)
- {
- /* this vector has only matrices going up */
-
- /* append to LAST list */
- PREDVC(LAST_last_in) = theNbVector;
- LAST_last_in = theNbVector;
- PREDVC(LAST_last_in) = NULL;
- SETVCUSED(theNbVector,1);
- }
- }
- else if (MDOWN(theMatrix))
- {
- k = VUP(theNbVector);
- assert(k>0); /* if 0 is supposed to be VCUSED already */
- SETVUP(theNbVector,--k);
- if (k==0)
- {
- /* this vector has only matrices going down */
-
- /* append to FIRST list */
- PREDVC(FIRST_last_in) = theNbVector;
- FIRST_last_in = theNbVector;
- PREDVC(FIRST_last_in) = NULL;
- SETVCUSED(theNbVector,1);
- }
- }
- }
- if (LAST_next_out==LAST_line_end)
- {
- if (theLastBV==NULL)
- {
- if (CreateBlockvector_l0(theGrid,&theBV,theFirstBV,1)!=GM_OK) return (__LINE__);
- }
- else
- {
- if (CreateBlockvector_l0(theGrid,&theBV,theLastBV,0)!=GM_OK) return (__LINE__);
- }
- theLastBV = theBV;
- if (init && (line==0))
- {
- SET_ORD_GCL(BVNUMBER(theLastBV),GM_GEN_LAST,0,line);
- }
- else
- {
- SET_ORD_GCL(BVNUMBER(theLastBV),GM_GEN_LAST,cycle,line);
- }
- BVFIRSTVECTOR(theLastBV) = LAST_line_end; /* end because order is reverted below */
-
- line++;
- LAST_line_end = LAST_last_in;
- }
- }
- LAST_next_out = LAST_last_in;
-
- /* check if there is a !used vec in the strong fifo */
- if ((theVector=FindOptimalStrong(&LFifo))!=NULL)
- LAST_last_in = CutStrongLine(theVector,LAST_last_in);
- }
- while (theVector!=NULL);
-
- init = 0;
-
- /****************************************************************************/
- /* get CUT (or Feedback Vertex)-set and do what needs to be done */
- /****************************************************************************/
-
- FIRST_last_in = (*FindCutSet)(theGrid,FIRST_last_in,&nCUT);
- if (FIRST_last_in==NULL) nCUT = 0;
- else PREDVC(FIRST_last_in) = NULL;
-
- nCutTot += nCUT;
-
- } while (nCUT>0);
-
- /* TODO: use dlmgr macros for lists after streamwise ordering */
-
- /* insert FIRST list one-by-one to LASTVECTOR list of grid */
- LASTVECTOR(theGrid) = NULL;
- for (theVector=PREDVC(&FIRST_handle); theVector!=NULL; theVector=predVector)
- {
- predVector = PREDVC(theVector);
- PREDVC(theVector) = LASTVECTOR(theGrid);
- LASTVECTOR(theGrid) = theVector;
- }
-
- /* insert LAST list as a whole to LASTVECTOR list of grid */
- PREDVC(LAST_last_in) = LASTVECTOR(theGrid);
- LASTVECTOR(theGrid) = PREDVC(&LAST_handle);
-
- /* construct SUCCVC list */
- succVector = NULL;
- for (theVector=LASTVECTOR(theGrid); theVector!=NULL; theVector=PREDVC(theVector))
- {
- SUCCVC(theVector) = succVector;
- succVector = theVector;
- }
- SFIRSTVECTOR(theGrid) = succVector;
- PREDVC(succVector) = NULL;
-
- /* set pointers in BLOCKVECTORs */
- BVENDVECTOR(GLASTBV(theGrid)) = NULL;
- for (theBV=GLASTBV(theGrid); theBV!=NULL; theBV=BVPRED(theBV))
- {
- if (BVSUCC(theBV)!=NULL && BVENDVECTOR(theBV)==NULL)
- BVENDVECTOR(theBV)=BVFIRSTVECTOR(BVSUCC(theBV));
- if (BVFIRSTVECTOR(theBV)==NULL)
- BVFIRSTVECTOR(theBV)=BVENDVECTOR(theBV);
- }
-
- /* set VCCUT */
- for (theBV=GFIRSTBV(theGrid); theBV!=NULL; theBV=BVSUCC(theBV))
- {
- i = 0;
- if (BV_GEN(theBV)==BV_GEN_C)
- for (theVector=BVFIRSTVECTOR(theBV); theVector!=BVENDVECTOR(theBV); theVector=SUCCVC(theVector))
- {
- SETVCCUT(theVector,1);
- i++;
- }
- else
- for (theVector=BVFIRSTVECTOR(theBV); theVector!=BVENDVECTOR(theBV); theVector=SUCCVC(theVector))
- {
- SETVCCUT(theVector,0);
- i++;
- }
- bvn = BVNUMBER(theBV);
- if (verboselevel>1)
- UserWriteF("# %d members in %c%d,%d\n",i,gen_label[ORD_GEN(bvn)],ORD_CYC(bvn),ORD_LIN(bvn));
- }
-
- if (verboselevel>0)
- {
- if (verboselevel>1)
- UserWrite("#\n# summary:\n");
- UserWriteF("# %d cycles: %d cut from %d\n",(int)cycle,(int)nCutTot,(int)NVEC(theGrid));
- a = POW((DOUBLE)NVEC(theGrid),(DOUBLE)(DIM-1)/(DOUBLE)DIM); a = (DOUBLE)nCutTot/a;
- UserWriteF("# corr. to %6.2f hyp. planes\n",(float)a);
- }
-
- /* check # members of succ list */
- i=0;
- for (theVector=FIRSTVECTOR(theGrid); theVector!= NULL; theVector=SUCCVC(theVector)) i++;
- if (NVEC(theGrid) != i)
- {
- UserWrite("vectorstructure corrupted\n");
- return (__LINE__);
- }
-
- /* check # members of pred list */
- i=0;
- for (theVector=LASTVECTOR(theGrid); theVector!= NULL; theVector=PREDVC(theVector)) i++;
- if (NVEC(theGrid) != i)
- {
- UserWrite("vectorstructure corrupted\n");
- return (__LINE__);
- }
-
- /* set index w.r.t. new order, beginning with 1 */
- i = 1;
- for (theVector=FIRSTVECTOR(theGrid); theVector!= NULL; theVector=SUCCVC(theVector))
- VINDEX(theVector) = i++;
-
- ReleaseTmpMem(MGHEAP(MYMG(theGrid)),MarkKey);
-
- return (0);
-}
-
-/****************************************************************************/
-/** \brief Driver for general vector ordering
-
- * @param theMG - multigrid to order
- * @param levels - GM_ALL_LEVELS or GM_CURRENT_LEVEL
- * @param mode - GM_FCFCLL or GM_FFCCLL (see orderv command)
- * @param dependency - name of user defined dependency item
- * @param dep_options - options for user dependency function
-
- This function orders VECTORs in a multigrid according to the dependency
- function provided.
-
- * @return <ul>
- * <li> GM_OK if ok
- * <li> GM_ERROR if error occured.
- </ul>
- */
-/****************************************************************************/
-
-INT NS_DIM_PREFIX LineOrderVectors (MULTIGRID *theMG, INT levels, const char *dependency, const char *dep_options, const char *findcut, INT verboselevel)
-{
- INT i, currlevel, baselevel;
- ALG_DEP *theAlgDep;
- FIND_CUT *theFindCut;
- GRID *theGrid;
- DependencyProcPtr DependencyProc;
- INT err;
-
- /* current level */
- currlevel = theMG->currentLevel;
-
- /* get algebraic dependency */
- theAlgDep = (ALG_DEP *) SearchEnv(dependency,"/Alg Dep",theAlgDepVarID,theAlgDepDirID);
- if (theAlgDep==NULL)
- {
- UserWrite("algebraic dependency not found\n");
- return(GM_ERROR);
- }
- DependencyProc = theAlgDep->DependencyProc;
- if (DependencyProc==NULL)
- {
- UserWrite("don't be stupid: implement a dependency!\n");
- return(GM_ERROR);
- }
-
- /* get find cut dependency */
- if (findcut==NULL)
- {
- FindCutSet = FeedbackVertexVectors;
- UserWrite("default cut set proc:\n leaving order of cyclic dependencies unchanged\n");
- }
- else
- {
- theFindCut = (FIND_CUT *) SearchEnv(findcut,"/FindCut",theFindCutVarID,theFindCutDirID);
- if (theFindCut==NULL)
- {
- UserWrite("find cut proc not found\n");
- return(GM_ERROR);
- }
- FindCutSet = theFindCut->FindCutProc;
- if (FindCutSet==NULL)
- {
- UserWrite("don't be stupid: implement a find cut proc!\n");
- return(GM_ERROR);
- }
- }
- if (AllocateControlEntry(VECTOR_CW,VCSTRONG_LEN,&ce_VCSTRONG) != GM_OK)
- return (GM_ERROR);
-
- /* go */
- if (levels==GM_ALL_LEVELS)
- baselevel = 0;
- else
- baselevel = currlevel;
- for (i=baselevel; i<=currlevel; i++)
- {
- theGrid = GRID_ON_LEVEL(theMG,i);
- if ((*DependencyProc)(theGrid,dep_options))
- {
- PrintErrorMessage('E',"LineOrderVectors","DependencyProc failed");
- return (GM_ERROR);
- }
- if ((err=LineOrderVectorsAlgebraic(theGrid,verboselevel))!=0)
- {
- PrintErrorMessage('E',"LineOrderVectors","LineOrderVectorsAlgebraic failed");
- return (GM_ERROR);
- }
- }
-
- FreeControlEntry(ce_VCSTRONG);
-
- return (GM_OK);
-}
-
-/****************************************************************************/
-/** \brief Prepare for lineorderv
-
- * @param theGrid - pointer to grid
-
- This function resets flags in vectors for the line ordering algorithm.
-
- * @return <ul>
- * <li> 0 if ok
- * <li> 1 if error occured.
- </ul>
- */
-/****************************************************************************/
-
-INT NS_DIM_PREFIX PrepareForLineorderVectors (GRID *theGrid)
-{
- VECTOR *vec;
-
- for (vec=FIRSTVECTOR(theGrid); vec!=NULL; vec=SUCCVC(vec))
- {
- SETVCUSED(vec,false);
- SETVCFLAG(vec,false);
- }
-
- return (0);
-}
-
-/****************************************************************************/
-/** \brief Flag begin or end vectors for lineorderv
-
- * @param elem - set flag in marked vectors of the element
- * @param dt - data types of vectors needed
- * @param ot - object types of vectors needed
- * @param mark - set GM_LOV_BEGIN or GM_LOV_END or no flag
-
- This function sets flags in vectors for the line ordering algorithm.
-
- * @return <ul>
- * <li> 0 if ok
- * <li> 1 if error occured.
- </ul>
- */
-/****************************************************************************/
-
-INT NS_DIM_PREFIX MarkBeginEndForLineorderVectors (ELEMENT *elem, INT dt, INT ot, const INT *mark)
-{
- VECTOR *vList[MAX_ELEM_VECTORS];
- INT i,cnt;
-
- if (GetVectorsOfDataTypesInObjects(elem,dt,ot,&cnt,vList)!=GM_OK)
- REP_ERR_RETURN (1);
-
- for (i=0; i<cnt; i++)
- switch (mark[i])
- {
- case GM_LOV_BEGIN : SETVCUSED(vList[i],true);
- case GM_LOV_END : SETVCFLAG(vList[i],true);
- }
-
- return (0);
-}
-
-/****************************************************************************/
-/** \brief Revert order in vector list
-
- * @param theGrid - pointer to grid
-
- This function revertes the order in the vector list.
-
- * @return <ul>
- * <li> 0 if ok
- * <li> 1 if error occured.
- </ul>
- */
-/****************************************************************************/
-
-INT NS_DIM_PREFIX RevertVecOrder (GRID *theGrid)
-{
- VECTOR *vec,*tmp;
- BLOCKVECTOR *bv;
-
- for (vec=FIRSTVECTOR(theGrid); vec!=NULL; vec=PREDVC(vec))
- {
- SWAP(PREDVC(vec),SUCCVC(vec),tmp);
- }
- SWAP(SFIRSTVECTOR(theGrid),LASTVECTOR(theGrid),tmp);
-
- /* also change the blockvectors */
- for (bv=GFIRSTBV(theGrid); bv!=NULL; bv=BVSUCC(bv))
- {
- tmp = BVFIRSTVECTOR(bv);
- BVFIRSTVECTOR(bv) = (BVENDVECTOR(bv)==NULL) ? FIRSTVECTOR(theGrid) : SUCCVC(BVENDVECTOR(bv));
- BVENDVECTOR(bv) = SUCCVC(tmp);
- }
-
- return (0);
-}
-
/****************************************************************************/
/** \brief Dependency function for lexicographic ordering
diff --git a/gm/gm.h b/gm/gm.h
index 492e59e62..52c691f8c 100644
--- a/gm/gm.h
+++ b/gm/gm.h
@@ -348,10 +348,6 @@ enum NodeType {CORNER_NODE,
#define TOPNODE(p) ((p)->iv.topnode)
*/
-/** \brief Modes for LexOrderVectorsInGrid */
-enum {OV_CARTES,
- OV_POLAR};
-
/****************************************************************************/
/* */
/* format definition data structure */
@@ -3549,13 +3545,6 @@ NS_PREFIX BLOCK_DESC *GetMGUDBlockDescriptor (NS_PREFIX BLOCK_ID id);
/* ordering of degrees of freedom */
ALG_DEP *CreateAlgebraicDependency (const char *name, DependencyProcPtr DependencyProc);
FIND_CUT *CreateFindCutProc (const char *name, FindCutProcPtr FindCutProc);
-INT LexOrderVectorsInGrid (GRID *theGrid, INT mode, const INT *order, const INT *sign, INT which, INT SpecSkipVecs, INT AlsoOrderMatrices);
-INT OrderVectors (MULTIGRID *theMG, INT levels, INT mode, INT PutSkipFirst, INT SkipPat, const char *dependency, const char *dep_options, const char *findcut);
-INT ShellOrderVectors (GRID *theGrid, VECTOR *seed);
-INT PrepareForLineorderVectors (GRID *theGrid);
-INT MarkBeginEndForLineorderVectors (ELEMENT *elem, INT dt, INT ot, const INT *mark);
-INT LineOrderVectors (MULTIGRID *theMG, INT levels, const char *dependency, const char *dep_options, const char *findcut, INT verboselevel);
-INT RevertVecOrder (GRID *theGrid);
/* functions for evaluation-fct management */
INT InitEvalProc (void);
--
GitLab