femgibbs.cpp

FreeFEM is an implementation of the GFEM language dedicated to the finite element method. It provid

// Emacs will be in -*- Mode: c++ -*-
//
// ********** DO NOT REMOVE THIS BANNER **********
//
// SUMMARY: Language for a Finite Element Method
//
// AUTHORS:  C. Prud'homme
// ORG    :  
// E-MAIL : [e-mail]
//
// ORIG-DATE:     June-94
// LAST-MOD: 12-Aug-02 at 08:42:57 by Christophe Prud'homme
//
// DESCRIPTION:  
/*
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.
  
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

*/

// DESCRIP-END.
//

//
// system includes
//
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

//
// Freefem includes
//
#include <femMesh.hpp>
#include <femGraphic.hpp>

#define mmax(a,b)(a>b?a:b)
#define mmin(a,b)(a<b?a:b)
#define ffalse 0
#define ttrue 1

namespace fem
{
/*  -- translated by f2c (version of 23 May 1992  14:18:33).
   You must link the resulting object file with the libraries:
   -lF77 -lI77 -lm -lc   (in that order)
 */

/* #include "f2c.h" */


/* Subroutine */
int 
femMesh::gibbs1_ (integer * n, integer * record, integer * ptvois)
{
  /* System generated locals */
  integer         i__1;

  /* Local variables */
  static integer  crit, i, j, l, r, rec;

/* -----------------------------------------------------------------------
 */
/*     routine appele par gibbs0 */
/* -----------------------------------------------------------------------
 */
/*     but: trie record (ensemble de n sommet) telle que l'ordre des somm 
 */
/*     soit croissant (ordre du sommet i est ptvois(i+1)-ptvois(i)) */
/* -----------------------------------------------------------------------
 */

  /* Parameter adjustments */
  --ptvois;
  --record;

  /* Function Body */
  if (*n <= 1)
     {
       return 0;
     }
  l = *n / 2 + 1;
  r = *n;
L2:
  if (l <= 1)
     {
       goto L20;
     }
  --l;
  rec = record[l];
  crit = ptvois[record[l] + 1] - ptvois[record[l]];
  goto L3;
L20:
  rec = record[r];
  crit = ptvois[record[r] + 1] - ptvois[record[r]];
  record[r] = record[1];
  --r;
  if (r == 1)
     {
       goto L999;
     }
L3:
  j = l;
L4:
  i = j;
  j <<= 1;
  if ((i__1 = j - r) < 0)
     {
       goto L5;
     }
  else if (i__1 == 0)
     {
       goto L6;
     }
  else
     {
       goto L8;
     }
L5:
  if (ptvois[record[j] + 1] - ptvois[record[j]] < ptvois[record[j + 1] + 1]
      - ptvois[record[j + 1]])
     {
       ++j;
     }
L6:
  if (crit >= ptvois[record[j] + 1] - ptvois[record[j]])
     {
       goto L8;
     }
  record[i] = record[j];
  goto L4;
L8:
  record[i] = rec;
  goto L2;
L999:
  record[1] = rec;
  return 0;
}				/* gibbs1_ */

/* Subroutine */
int 
femMesh::gibbs2_ (integer * n, integer * record, integer * criter)
{
  static integer  crit, i, j, l, r, rec;


/*     trie record selon les valeurs de criter(record(.)) croissantes */


  /* Parameter adjustments */
  --criter;
  --record;

  /* Function Body */
  if (*n <= 1)
     {
       return 0;
     }
  l = *n / 2 + 1;
  r = *n;
L2:
  if (l <= 1)
     {
       goto L20;
     }
  --l;
  rec = record[l];
  crit = criter[rec];
  goto L3;
L20:
  rec = record[r];
  crit = criter[rec];
  record[r] = record[1];
  --r;
  if (r == 1)
     {
       goto L999;
     }
L3:
  j = l;
L4:
  i = j;
  j <<= 1;
  if (j - r < 0)
     {
       goto L5;
     }
  else if (j == r)
     {
       goto L6;
     }
  else
     {
       goto L8;
     }
L5:
  if (criter[record[j]] < criter[record[j + 1]])
     {
       ++j;
     }
L6:
  if (crit >= criter[record[j]])
     {
       goto L8;
     }
  record[i] = record[j];
  goto L4;
L8:
  record[i] = rec;
  goto L2;
L999:
  record[1] = rec;
  return 0;
}				/* gibbs2_ */

/* Subroutine */
int 
femMesh::gibbsa_ (integer * n, integer * ptvois, integer * vois, integer * r, integer * m,
                        integer * nv, integer * nx, integer * ny, integer * nn, integer * w1, integer * w2,
                        integer * pfold, integer * pfnew, integer * impre, integer * nfout)
{
  /* System generated locals */
  integer         i__1, i__2, i__3, i__4;

  /* Builtin functions */
  /* Subroutine */ int s_stop ();

  /* Local variables */
  static integer  nbcc, degi, bold, bnew, i, j, k, p, degre, x, y, p1,
                  p2;

/*    extern  Subroutine  int gibbs1_(); */
  static integer  pf;

/*    extern  Subroutine int gibbsb_(), gibbsd_(), gibbst_(); */
  static integer  nbpass, niveau, pf1, option, old, new_, opt, new1;

/* -----------------------------------------------------------------------
 */
/*  but: calculer une renumerotation des sommets d'un graphe defini par: 
 */
/*     par la methode de gibbs */
/* -----------------------------------------------------------------------
 */
/*  entree */
/* -------- */
/*     n = nb de sommet du graphe */
/*      les voisins d'un sommet i ont pour numero : */
/*     ( vois(j) , j=ptvois(i),ptvois(i+1)-1 ) */

/*     impre   parametre d'impression */
/*     nfout   numero du fichier pour impression */

/*  sortie */
/*  ------ */
/*     r(1:n) tableau donnant la nouvelle numerotation: */
/*       r(i) = nouveau numero du sommet i */
/*     pfolf = ancien  profile */
/*     pfnew = nouveau profile */

/*  tableau de travail : */
/*  -------------------- */
/*     m(n) */
/*     nv(0:n+n) */
/*     nx(n) */
/*     ny(n) */
/*     nn(0:n) */
/*     w1(n) */
/*     w2(n) */

/* -----------------------------------------------------------------------
 */
/*     programmeur f. hecht  le 3/02/1987 */
/* -----------------------------------------------------------------------
 */

/*     tri des voisins d'un sommet du graphe par degre croissant */
/* --------------------------------------------------------------- */
  /* Parameter adjustments */
  --w2;
  --w1;
  --ny;
  --nx;
  --m;
  --r;
  --vois;
  --ptvois;

  /* Function Body */
  p2 = ptvois[1] - 1;
  i__1 = *n;
  for (i = 1; i <= i__1; ++i)
     {
       p1 = p2 + 1;
       p2 = ptvois[i + 1] - 1;
       i__2 = p2 - p1 + 1;
       gibbs1_ (&i__2, &vois[p1], &ptvois[1]);
/*       if(impre.le.-9) then */
/*        write (nfout,*) 'les voisin de ',i,'sont: ', (vois(j),j=p1,p
   2) */
/*       endif */
/* L10: */
     }
  i__1 = *n;
  for (i = 1; i <= i__1; ++i)
     {
       r[i] = 0;
/* L20: */
     }
/*     boucle sur les composante connexe du graphe */
  new_ = 0;
  nbcc = 0;
L30:
  if (new_ < *n)
     {
       ++nbcc;
/*       recherche d'une racine y (un sommet non numerote) de degree m
   ini */
       y = 0;
       degre = *n + 1;
       i__1 = *n;
       for (i = 1; i <= i__1; ++i)
	  {
	    if (r[i] <= 0)
	       {
		 degi = ptvois[i + 1] - ptvois[i];
		 if (degi < degre)
		    {
		      degre = degi;
		      y = i;
		    }
	       }
/* L40: */
	  }
       if (y == 0)
	  {
	    return -3;		/*  s_stop("fatal erreur  gibbs 2 : pb racine", 33L); */
	  }
       gibbsd_ (&y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
       nbpass = 0;
     L50:
       ++nbpass;
       x = y;
       p = niveau;
       k = 0;
       i__1 = nv[p + 1];
       for (i = nv[p] + 1; i <= i__1; ++i)
	  {
	    ++k;
	    m[k] = nv[i];
/* L60: */
	  }
       gibbs1_ (&k, &m[1], &ptvois[1]);
       i__1 = k;
       for (i = 1; i <= i__1; ++i)
	  {
	    y = m[i];
	    gibbsd_ (&y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
	    if (niveau > p)
	       {
		 goto L50;
	       }
/* L70: */
	  }
       y = m[1];
/*        if(impre.lt.0) then */
/*          write(nfout,*) */
/*     +  '    nb de pass pour trouver le pseudo diametre',nbpass */
/*     +         ,' x=',x,',y=',y,' de la composante connexe ',nbcc */

/*          write (nfout,*) ('-',i=1,78) */
/*        endif */
/*       optimisation de la descendance de la numerotation */
/*       ------------------------------------------------- */
       gibbsb_ (&x, &y, n, &ptvois[1], &vois[1], &nx[1], &ny[1], nv, nn, &m[1]
		,&w1[1], &w2[1], &r[1], impre, nfout);

/*     renumerotation de cuthill mac kee avec la meilleur des 4 option
   s */
/*     --------------------------------------------------------------
   --- */
       pf = 1073741824;
       option = -2;
       new1 = new_;
       for (opt = -2; opt <= 2; ++opt)
	  {
	    new_ = new1;
	    if (opt != 0)
	       {
		 gibbst_ (n, &p, nv, nn, &ptvois[1], &vois[1], &m[1], &r[1], &
			  new_, &opt, &pf1, impre, nfout);
		 if (pf1 < pf)
		    {
		      pf = pf1;
		      option = opt;
		    }
	       }
/* L80: */
	  }
/*        if(impre.ne.0) write (nfout,*) '    on a choisi l''option ',
 */
/*     +             option,', new =',new */
       new_ = new1;
       gibbst_ (n, &p, nv, nn, &ptvois[1], &vois[1], &m[1], &r[1], &new_, &
		option, &pf1, impre, nfout);
       goto L30;
     }
/*      if(impre.ne.0) write(nfout,*) */
/*     +       '   nb de composante connexe du graphe =',nbcc */
/*     calcul du profile */
  *pfold = 0;
  *pfnew = 0;
  bnew = 0;
  bold = 0;
  i__1 = *n;
  for (i = 1; i <= i__1; ++i)
     {
       old = i;
       new_ = r[i];
       i__2 = ptvois[i + 1] - 1;
       for (j = ptvois[i]; j <= i__2; ++j)
	  {
/* Computing MIN */
	    i__3 = old, i__4 = vois[j];
	    old = mmin (i__3, i__4);
/* Computing MIN */
	    i__3 = new_, i__4 = r[vois[j]];
	    new_ = mmin (i__3, i__4);
/* L100: */
	  }
       *pfold = *pfold + i - old + 1;
/* Computing MAX */
       i__2 = bold, i__3 = i - old + 1;
       bold = mmax (i__2, i__3);
       *pfnew = *pfnew + r[i] - new_ + 1;
/* Computing MAX */
       i__2 = bnew, i__3 = r[i] - new_ + 1;
       bnew = mmax (i__2, i__3);
/* L110: */
     }
/*      if(impre.ne.0) then */
/*        write(nfout,*)'profile  old  = ',pfold,', profile  new = ',pfnew
 */
/*        write(nfout,*)'1/2 bande old = ',bold ,', 1/2 band new = ',bnew 
 */
/*      endif */
  return 0;
}				/* gibbsa_ */

/* Subroutine */
int 
femMesh::gibbsb_ (integer * x, integer * y, integer * n, integer * ptvois,
                        integer * vois, integer * nx, integer * ny, integer * nv, integer * nn, integer * m,
                        integer * wh, integer * wl, integer * r, integer * impre, integer * nfout)
{
  /* System generated locals */
  integer         i__1, i__2;

  /* Local variables */
  static int flag_;
  static integer  i, j, k, p, s, h0, i1, l0, i2;

/*    extern  Subroutine  int gibbs1_(); */
  static integer  lg;

/*    extern  Subroutine  int gibbsd_(), gibbsc_(); */
  static integer  niveau, mxcanx, mxcany, nbc;

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 */
/* ...................................................................... 
 */
/*     attention on met la descente optimiser dans r <0 ou nulle */
/* .......................................................................
 */
  /* Parameter adjustments */
  --r;
  --m;
  --ny;
  --nx;
  --vois;
  --ptvois;

  /* Function Body */
  gibbsd_ (y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
  gibbsc_ (&ny[1], nv, &niveau, n, &mxcany);
  gibbsd_ (x, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
  p = niveau;
  gibbsc_ (&nx[1], nv, &niveau, n, &mxcanx);
  flag_ = ffalse;
  i__1 = *n;
  for (i = 1; i <= i__1; ++i)
     {
       if (nx[i] + ny[i] == p)
	  {
	    r[i] = -nx[i];
	  }
       else if (nx[i] >= 0)
	  {
	    flag_ = ttrue;
	    r[i] = -1073741824;
	  }
       else
	  {
	    if (r[i] <= 0)
	       {
		 r[i] = -1073741822;
	       }
	  }
/* L20: */
     }
  if (flag_)
     {
/*       calcul des composantes connexe du graphe sans les sommets de 
   nn */
/*       ------------------------------------------------------------
   --- */
       j = *n;
       k = 0;
       nbc = 0;
       nv[nbc] = j;
     L30:
       ++k;
       if (k <= *n)
	  {
	    if (r[k] == -1073741824)
	       {
/*           recherche de la fermeture transitive partant de k
 */
		 ++nbc;
		 i = -1;
		 s = k;
	       L40:
		 ++i;
		 wl[i] = ptvois[s];
		 wh[i] = ptvois[s + 1];
		 ++j;
		 nv[j] = s;
		 r[s] = -1073741823;
	       L50:
		 if (i >= 0)
		    {
		      if (wl[i] < wh[i])
			 {
			   s = vois[wl[i]];
			   ++wl[i];
			   if (r[s] == -1073741824)
			      {
				goto L40;
			      }
			   goto L50;
			 }
		      --i;
		      goto L50;
		    }
		 nv[nbc] = j;
		 m[nbc] = nbc;
	       }
	    goto L30;
	  }
/*        if(impre.lt.0) write(nfout,*) */
/*     +         ' nb de composante connexe du graphe reduit =',nbc */

/* --------------- fin de construction des composantes connexes------
   --- */
/*        nv(0)=n */
/*        if(impre.le.-10) write(nfout,5555)'nv(0:n+n) = ',(nv(i),i=0,
   n+n) */
       gibbs1_ (&nbc, &m[1], nv);
/*        if(impre.le.-10)write(nfout,5555)'trie m =',(m(i),i=1,nbc) 
 */
       i__1 = p;
       for (i = 0; i <= i__1; ++i)
	  {
	    nn[i] = 0;
/* L60: */
	  }
       i__1 = *n;
       for (i = 1; i <= i__1; ++i)
	  {
	    j = -r[i];
	    if (j >= 0 && j <= p)
	       {
		 ++nn[j];
	       }
/* L70: */
	  }