1 // -*- Mode : c++ -*-
2 //
3 // SUMMARY :
4 // USAGE :
5 // ORG :
6 // AUTHOR : Frederic Hecht
7 // E-MAIL : hecht@ann.jussieu.fr
8 //
9
10 /*
11
12 This file is part of Freefem++
13
14 Freefem++ is free software; you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation; either version 2.1 of the License, or
17 (at your option) any later version.
18
19 Freefem++ is distributed in the hope that it will be useful,
20 but WITHOUT ANY WARRANTY; without even the implied warranty of
21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 GNU Lesser General Public License for more details.
23
24 You should have received a copy of the GNU Lesser General Public License
25 along with Freefem++; if not, write to the Free Software
26 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
27 */
28 #include <cmath>
29 #include "error.hpp"
30 #include <iostream>
31 #include <fstream>
32 #include "RNM.hpp"
33 #include "rgraph.hpp"
34 #include "fem.hpp"
35 using namespace Fem2D;
36 #include "FESpacen.hpp"
37 #include "FESpace.hpp"
38
39 #define mmax(a,b)(a>b?a:b)
40 #define mmin(a,b)(a<b?a:b)
41 #define ffalse 0
42 #define ttrue 1
43 extern long verbosity;
44 /* -- translated by f2c (version of 23 May 1992 14:18:33).
45 You must link the resulting object file with the libraries:
46 -lF77 -lI77 -lm -lc (in that order)
47 */
48
49
50 #define integer long
51 #define logical long
52
53 int gibbs1_(integer* n,integer* record,integer* ptvois);
54 int gibbs2_(integer* n,integer* record,integer* criter);
55 int gibbsa_(integer* n,integer* ptvois,integer* vois,integer* r,integer* m,
56 integer* nv,integer* nx,integer* ny,integer* nn,integer* w1,integer* w2,
57 integer* pfold,integer* pfnew,integer* impre,integer* nfout);
58 int gibbsb_(integer* x,integer* y,integer* n,integer* ptvois,
59 integer* vois,integer* nx,integer* ny,integer* nv,integer* nn,integer* m,
60 integer* wh,integer* wl,integer* r, integer* impre, integer* nfout);
61 int gibbsc_(integer* nz,integer* nv,integer* niveau,integer* n,integer* );
62 int gibbsd_(integer* racine,integer* n,integer* ptvois,integer*
63 vois,integer* nv,integer* r,integer* niveau);
64 int gibbst_(integer* n,integer* p,integer* nv,integer* nn,integer* ptvois,integer* vois,
65 integer* m,integer* r,integer* new_,integer* option,
66 integer* pfnew,integer* impre,integer* nfout);
67
gibbs1_(integer * n,integer * record,integer * ptvois)68 /* Subroutine */ int gibbs1_(integer* n,integer* record,integer* ptvois)
69 {
70 /* System generated locals */
71 integer i__1;
72
73 /* Local variables */
74 static integer crit, i, j, l, r, rec;
75
76 /* -----------------------------------------------------------------------
77 */
78 /* routine appele par gibbs0 */
79 /* -----------------------------------------------------------------------
80 */
81 /* but: trie record (ensemble de n sommet) telle que l'ordre des somm
82 */
83 /* soit croissant (ordre du sommet i est ptvois(i+1)-ptvois(i)) */
84 /* -----------------------------------------------------------------------
85 */
86
87 /* Parameter adjustments */
88 --ptvois;
89 --record;
90
91 /* Function Body */
92 if (*n <= 1) {
93 return 0;
94 }
95 l = *n / 2 + 1;
96 r = *n;
97 L2:
98 if (l <= 1) {
99 goto L20;
100 }
101 --l;
102 rec = record[l];
103 crit = ptvois[record[l] + 1] - ptvois[record[l]];
104 goto L3;
105 L20:
106 rec = record[r];
107 crit = ptvois[record[r] + 1] - ptvois[record[r]];
108 record[r] = record[1];
109 --r;
110 if (r == 1) {
111 goto L999;
112 }
113 L3:
114 j = l;
115 L4:
116 i = j;
117 j <<= 1;
118 if ((i__1 = j - r) < 0) {
119 goto L5;
120 } else if (i__1 == 0) {
121 goto L6;
122 } else {
123 goto L8;
124 }
125 L5:
126 if (ptvois[record[j] + 1] - ptvois[record[j]] < ptvois[record[j + 1] + 1]
127 - ptvois[record[j + 1]]) {
128 ++j;
129 }
130 L6:
131 if (crit >= ptvois[record[j] + 1] - ptvois[record[j]]) {
132 goto L8;
133 }
134 record[i] = record[j];
135 goto L4;
136 L8:
137 record[i] = rec;
138 goto L2;
139 L999:
140 record[1] = rec;
141 return 0;
142 } /* gibbs1_ */
143
gibbs2_(integer * n,integer * record,integer * criter)144 /* Subroutine */ int gibbs2_(integer* n,integer* record,integer* criter)
145 {
146 static integer crit, i, j, l, r, rec;
147
148
149 /* trie record selon les valeurs de criter(record(.)) croissantes */
150
151
152 /* Parameter adjustments */
153 --criter;
154 --record;
155
156 /* Function Body */
157 if (*n <= 1) {
158 return 0;
159 }
160 l = *n / 2 + 1;
161 r = *n;
162 L2:
163 if (l <= 1) {
164 goto L20;
165 }
166 --l;
167 rec = record[l];
168 crit = criter[rec];
169 goto L3;
170 L20:
171 rec = record[r];
172 crit = criter[rec];
173 record[r] = record[1];
174 --r;
175 if (r == 1) {
176 goto L999;
177 }
178 L3:
179 j = l;
180 L4:
181 i = j;
182 j <<= 1;
183 if (j - r < 0) {
184 goto L5;
185 } else if (j == r) {
186 goto L6;
187 } else {
188 goto L8;
189 }
190 L5:
191 if (criter[record[j]] < criter[record[j + 1]]) {
192 ++j;
193 }
194 L6:
195 if (crit >= criter[record[j]]) {
196 goto L8;
197 }
198 record[i] = record[j];
199 goto L4;
200 L8:
201 record[i] = rec;
202 goto L2;
203 L999:
204 record[1] = rec;
205 return 0;
206 } /* gibbs2_ */
207
208 /* Subroutine */
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)209 int gibbsa_(integer* n,integer* ptvois,integer* vois,integer* r,integer* m,
210 integer* nv,integer* nx,integer* ny,integer* nn,integer* w1,integer* w2,
211 integer* pfold,integer* pfnew,integer* impre,integer* nfout)
212 {
213 /* System generated locals */
214 integer i__1, i__2, i__3, i__4;
215
216 /* Builtin functions */
217
218 /* Local variables */
219 static integer nbcc, degi, bold, bnew, i, j, k, p, degre, x, y, p1, p2;
220 /* extern Subroutine int gibbs1_();*/
221 static integer pf;
222 /* extern Subroutine int gibbsb_(), gibbsd_(), gibbst_();*/
223 static integer nbpass, niveau, pf1, option, old, new_, opt, new1;
224 long pfn1,pfo1,m1e9;
225 m1e9 = 1000000000;
226
227 /* -----------------------------------------------------------------------
228 */
229 /* but: calculer une renumerotation des sommets d'un graphe defini par:
230 */
231 /* par la methode de gibbs */
232 /* -----------------------------------------------------------------------
233 */
234 /* entree */
235 /* -------- */
236 /* n = nb de sommet du graphe */
237 /* les voisins d'un sommet i ont pour numero : */
238 /* ( vois(j) , j=ptvois(i),ptvois(i+1)-1 ) */
239
240 /* impre parametre d'impression */
241 /* nfout numero du fichier pour impression */
242
243 /* sortie */
244 /* ------ */
245 /* r(1:n) tableau donnant la nouvelle numerotation: */
246 /* r(i) = nouveau numero du sommet i */
247 /* pfolf = ancien profile */
248 /* pfnew = nouveau profile */
249
250 /* tableau de travail : */
251 /* -------------------- */
252 /* m(n) */
253 /* nv(0:n+n) */
254 /* nx(n) */
255 /* ny(n) */
256 /* nn(0:n) */
257 /* w1(n) */
258 /* w2(n) */
259
260 /* -----------------------------------------------------------------------
261 */
262 /* programmeur f. hecht le 3/02/1987 */
263 /* -----------------------------------------------------------------------
264 */
265
266 /* tri des voisins d'un sommet du graphe par degre croissant */
267 /* --------------------------------------------------------------- */
268 /* Parameter adjustments */
269 --w2;
270 --w1;
271 --ny;
272 --nx;
273 --m;
274 --r;
275 --vois;
276 --ptvois;
277
278 /* Function Body */
279 p2 = ptvois[1] - 1;
280 i__1 = *n;
281 for (i = 1; i <= i__1; ++i) {
282 p1 = p2 + 1;
283 p2 = ptvois[i + 1] - 1;
284 i__2 = p2 - p1 + 1;
285 gibbs1_(&i__2, &vois[p1], &ptvois[1]);
286 /* if(impre.le.-9) then */
287 /* write (nfout,*) 'les voisin de ',i,'sont: ', (vois(j),j=p1,p
288 2) */
289 /* endif */
290 /* L10: */
291 }
292 i__1 = *n;
293 for (i = 1; i <= i__1; ++i) {
294 r[i] = 0;
295 /* L20: */
296 }
297 /* boucle sur les composante connexe du graphe */
298 new_ = 0;
299 nbcc = 0;
300 L30:
301 if (new_ < *n) {
302 ++nbcc;
303 /* recherche d'une racine y (un sommet non numerote) de degree m
304 ini */
305 y = 0;
306 degre = *n + 1;
307 i__1 = *n;
308 for (i = 1; i <= i__1; ++i) {
309 if (r[i] <= 0) {
310 degi = ptvois[i + 1] - ptvois[i];
311 if (degi < degre) {
312 degre = degi;
313 y = i;
314 }
315 }
316 /* L40: */
317 }
318 if (y == 0) {
319 return -3;/* s_stop("fatal erreur gibbs 2 : pb racine", 33L); */
320 }
321 gibbsd_(&y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
322 nbpass = 0;
323 L50:
324 ++nbpass;
325 x = y;
326 p = niveau;
327 k = 0;
328 i__1 = nv[p + 1];
329 for (i = nv[p] + 1; i <= i__1; ++i) {
330 ++k;
331 m[k] = nv[i];
332 /* L60: */
333 }
334 gibbs1_(&k, &m[1], &ptvois[1]);
335 i__1 = k;
336 for (i = 1; i <= i__1; ++i) {
337 y = m[i];
338 gibbsd_(&y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
339 if (niveau > p) {
340 goto L50;
341 }
342 /* L70: */
343 }
344 y = m[1];
345 /* if(impre.lt.0) then */
346 /* write(nfout,*) */
347 /* + ' nb de pass pour trouver le pseudo diametre',nbpass */
348 /* + ,' x=',x,',y=',y,' de la composante connexe ',nbcc */
349
350 /* write (nfout,*) ('-',i=1,78) */
351 /* endif */
352 /* optimisation de la descendance de la numerotation */
353 /* ------------------------------------------------- */
354 gibbsb_(&x, &y, n, &ptvois[1], &vois[1], &nx[1], &ny[1], nv, nn, &m[1]
355 , &w1[1], &w2[1], &r[1], impre, nfout);
356
357 /* renumerotation de cuthill mac kee avec la meilleur des 4 option
358 s */
359 /* --------------------------------------------------------------
360 --- */
361 pf = 1073741824;
362 option = -2;
363 new1 = new_;
364 for (opt = -2; opt <= 2; ++opt) {
365 new_ = new1;
366 if (opt != 0) {
367 gibbst_(n, &p, nv, nn, &ptvois[1], &vois[1], &m[1], &r[1], &
368 new_, &opt, &pf1, impre, nfout);
369 if (pf1 < pf) {
370 pf = pf1;
371 option = opt;
372 }
373 }
374 /* L80: */
375 }
376 /* if(impre.ne.0) write (nfout,*) ' on a choisi l''option ',
377 */
378 /* + option,', new =',new */
379 new_ = new1;
380 gibbst_(n, &p, nv, nn, &ptvois[1], &vois[1], &m[1], &r[1], &new_, &
381 option, &pf1, impre, nfout);
382 goto L30;
383 }
384 /* if(impre.ne.0) write(nfout,*) */
385 /* + ' nb de composante connexe du graphe =',nbcc */
386 /* calcul du profile */
387 *pfold = 0;
388 *pfnew = 0;
389 pfo1=0;
390 pfn1=0;
391 bnew = 0;
392 bold = 0;
393 i__1 = *n;
394 for (i = 1; i <= i__1; ++i) {
395 old = i;
396 new_ = r[i];
397 i__2 = ptvois[i + 1] - 1;
398 for (j = ptvois[i]; j <= i__2; ++j) {
399 /* Computing MIN */
400 i__3 = old, i__4 = vois[j];
401 old = mmin(i__3,i__4);
402 /* Computing MIN */
403 i__3 = new_, i__4 = r[vois[j]];
404 new_ = mmin(i__3,i__4);
405 /* L100: */
406 }
407 *pfold = *pfold + i - old + 1;
408 /* Computing MAX */
409 i__2 = bold, i__3 = i - old + 1;
410 bold = mmax(i__2,i__3);
411 *pfnew = *pfnew + r[i] - new_ + 1;
412 /* Computing MAX */
413 i__2 = bnew, i__3 = r[i] - new_ + 1;
414 bnew = mmax(i__2,i__3);
415 if(*pfold>m1e9) { pfo1+=*pfold/m1e9; *pfold = *pfold%m1e9;}
416 if(*pfnew>m1e9) { pfn1+=*pfnew/m1e9; *pfnew = *pfnew%m1e9;}
417
418 /* L110: */
419 }
420 if(pfo1 || pfn1)
421 { // change unit of pf
422 if(pfo1==pfn1)
423 {
424 *pfnew= pfn1*10 + (*pfnew >= *pfold) ;
425 *pfold= pfo1*10 + (*pfnew <= *pfold) ;
426 }
427 else
428 {
429 *pfnew= pfn1 ;
430 *pfold= pfo1;
431 }
432
433 }
434 /* if(impre.ne.0) then */
435 /* write(nfout,*)'profile old = ',pfold,', profile new = ',pfnew
436 */
437 /* write(nfout,*)'1/2 bande old = ',bold ,', 1/2 band new = ',bnew
438 */
439 /* endif */
440 return 0;
441 } /* gibbsa_ */
442
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 *,integer *)443 /* Subroutine */ int gibbsb_(integer* x,integer* y,integer* n,integer* ptvois,
444 integer* vois,integer* nx,integer* ny,integer* nv,integer* nn,integer* m,
445 integer* wh,integer* wl,integer* r, integer* , integer* )
446 {
447 /* System generated locals */
448 integer i__1, i__2;
449
450 /* Local variables */
451 static int flag_;
452 static integer i, j, k, p, s, h0, i1, l0, i2;
453 /* extern Subroutine int gibbs1_(); */
454 static integer lg;
455 /* extern Subroutine int gibbsd_(), gibbsc_();*/
456 static integer niveau, mxcanx, mxcany, nbc;
457
458 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
459 */
460 /* ......................................................................
461 */
462 /* attention on met la descente optimiser dans r <0 ou nulle */
463 /* .......................................................................
464 */
465 /* Parameter adjustments */
466 --r;
467 --m;
468 --ny;
469 --nx;
470 --vois;
471 --ptvois;
472
473 /* Function Body */
474 gibbsd_(y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
475 gibbsc_(&ny[1], nv, &niveau, n, &mxcany);
476 gibbsd_(x, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
477 p = niveau;
478 gibbsc_(&nx[1], nv, &niveau, n, &mxcanx);
479 flag_ = ffalse;
480 i__1 = *n;
481 for (i = 1; i <= i__1; ++i) {
482 if (nx[i] + ny[i] == p) {
483 r[i] = -nx[i];
484 } else if (nx[i] >= 0) {
485 flag_ = ttrue;
486 r[i] = -1073741824;
487 } else {
488 if (r[i] <= 0) {
489 r[i] = -1073741822;
490 }
491 }
492 /* L20: */
493 }
494 if (flag_) {
495 /* calcul des composantes connexe du graphe sans les sommets de
496 nn */
497 /* ------------------------------------------------------------
498 --- */
499 j = *n;
500 k = 0;
501 nbc = 0;
502 nv[nbc] = j;
503 L30:
504 ++k;
505 if (k <= *n) {
506 if (r[k] == -1073741824) {
507 /* recherche de la fermeture transitive partant de k
508 */
509 ++nbc;
510 i = -1;
511 s = k;
512 L40:
513 ++i;
514 wl[i] = ptvois[s];
515 wh[i] = ptvois[s + 1];
516 ++j;
517 nv[j] = s;
518 r[s] = -1073741823;
519 L50:
520 if (i >= 0) {
521 if (wl[i] < wh[i]) {
522 s = vois[wl[i]];
523 ++wl[i];
524 if (r[s] == -1073741824) {
525 goto L40;
526 }
527 goto L50;
528 }
529 --i;
530 goto L50;
531 }
532 nv[nbc] = j;
533 m[nbc] = nbc;
534 }
535 goto L30;
536 }
537 /* if(impre.lt.0) write(nfout,*) */
538 /* + ' nb de composante connexe du graphe reduit =',nbc */
539
540 /* --------------- fin de construction des composantes connexes------
541 --- */
542 /* nv(0)=n */
543 /* if(impre.le.-10) write(nfout,5555)'nv(0:n+n) = ',(nv(i),i=0,
544 n+n) */
545 gibbs1_(&nbc, &m[1], nv);
546 /* if(impre.le.-10)write(nfout,5555)'trie m =',(m(i),i=1,nbc)
547 */
548 i__1 = p;
549 for (i = 0; i <= i__1; ++i) {
550 nn[i] = 0;
551 /* L60: */
552 }
553 i__1 = *n;
554 for (i = 1; i <= i__1; ++i) {
555 j = -r[i];
556 if (j >= 0 && j <= p) {
557 ++nn[j];
558 }
559 /* L70: */
560 }
561
562 /* boucle sur les composante connexes par ordre croissantes */
563 /* -------------------------------------------------------- */
564 for (k = nbc; k >= 1; --k) {
565 i = m[k];
566 i1 = nv[i - 1] + 1;
567 i2 = nv[i];
568 lg = i2 - i1 + 1;
569 /* if(impre.le.-7) */
570 /* + write(nfout,*) k,' composante ',i,',lg=',lg,',i1,i2
571 =',i1,i2 */
572 /* if(impre.le.-8) */
573 /* + write (nfout,5555)' ',(nv(i),i=i1,i2) */
574 h0 = 0;
575 l0 = 0;
576 i__1 = p;
577 for (j = 0; j <= i__1; ++j) {
578 wh[j] = nn[j];
579 wl[j] = nn[j];
580 /* L90: */
581 }
582 i__1 = i2;
583 for (i = i1; i <= i__1; ++i) {
584 s = nv[i];
585 ++wh[nx[s]];
586 ++wl[p - ny[s]];
587 /* L100: */
588 }
589 i__1 = p;
590 for (j = 0; j <= i__1; ++j) {
591 if (wh[j] != nn[j]) {
592 /* Computing MAX */
593 i__2 = wh[j];
594 h0 = mmax(i__2,h0);
595 }
596 if (wl[j] != nn[j]) {
597 /* Computing MAX */
598 i__2 = wl[j];
599 l0 = mmax(i__2,l0);
600 }
601 /* L110: */
602 }
603 if (h0 < l0 || (h0 == l0 && mxcanx <= mxcany)) {
604 /* if(impre.le.-2) write(nfout,*) */
605 /* + ' h0 = ',h0,',l0 = ',l0,' ------- XXXX
606 --------' */
607 i__1 = i2;
608 for (i = i1; i <= i__1; ++i) {
609 s = nv[i];
610 r[s] = -nx[s];
611 ++nn[-r[s]];
612 /* L120: */
613 }
614 } else {
615 /* if (impre.le.-2) write(nfout,*) */
616 /* + ' h0 = ',h0,',l0 = ',l0,' ------- YYYY
617 --------' */
618 i__1 = i2;
619 for (i = i1; i <= i__1; ++i) {
620 s = nv[i];
621 r[s] = -p + ny[s];
622 ++nn[-r[s]];
623 /* L130: */
624 }
625 }
626 /* L140: */
627 }
628 }
629 /* on met les nouveaux niveaux de la descendance optimiser dans nn */
630 /* -----------------------------------------------------------------
631 */
632 i__1 = *n;
633 for (i = 1; i <= i__1; ++i) {
634 if (r[i] > 0) {
635 nn[i] = -1;
636 } else if (r[i] == -1073741822) {
637 nn[i] = -2;
638 } else {
639 nn[i] = -r[i];
640 }
641 /* L150: */
642 }
643 /* if(impre.le.-10)write (nfout,5555)' nn(i)=',(nn(i),i=1,n) */
644 /* 5555 format(' -------- ',a,/,5(15x,10(i5)/)) */
645 return 0;} /* gibbsb_ */
646
gibbsc_(integer * nz,integer * nv,integer * niveau,integer * n,integer * mxz)647 /* Subroutine */ int gibbsc_(integer* nz,integer* nv,integer* niveau,integer* n,integer* mxz)
648 {
649 /* System generated locals */
650 integer i__1, i__2, i__3;
651
652 /* Local variables */
653 static integer i, j;
654
655 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
656 */
657 /* Parameter adjustments */
658 --nz;
659
660 /* Function Body */
661 i__1 = *n;
662 for (i = 1; i <= i__1; ++i) {
663 nz[i] = -1;
664 /* L10: */
665 }
666 *mxz = 0;
667 i__1 = *niveau;
668 for (i = 0; i <= i__1; ++i) {
669 /* Computing MAX */
670 i__2 = *mxz, i__3 = nv[i + 1] - nv[i];
671 *mxz = mmax(i__2,i__3);
672 i__2 = nv[i + 1];
673 for (j = nv[i] + 1; j <= i__2; ++j) {
674 if(nv[j] > *n)
675 printf(" bug in gibbsc_ ##### %ld %ld %ld %ld \n",j,nv[j],i,*niveau );
676 if(nv[j] <= *n)
677 nz[nv[j]] = i;
678 /* L20: */
679 }
680 }
681 return 0;} /* gibbsc_ */
682
gibbsd_(integer * racine,integer * n,integer * ptvois,integer * vois,integer * nv,integer * r,integer * niveau)683 /* Subroutine */ int gibbsd_(integer* racine,integer* n,integer* ptvois,integer* vois,integer* nv,integer* r,integer* niveau)
684 {
685 /* System generated locals */
686 integer i__1, i__2;
687
688 /* Local variables */
689 static integer i, k, s, sv, stk, stk1, stk2, nvni=-1;
690
691 /* -----------------------------------------------------------------------
692 */
693 /* but construire la structure des descendant de racine du graphe */
694 /* -----------------------------------------------------------------------
695 */
696 /* sortie : */
697 /* -------- */
698 /* nv est la structure des niveaux */
699 /* les sommets du niveau (i =0,niveau_ sont defini par : */
700 /* (nv(j),j=nv(i),nv(i+1)-1) */
701
702 /* le tableau r(i) n'est modifier que sur les sommets */
703 /* de la composante connexe du graphe contenant la racine */
704
705 /* -----------------------------------------------------------------------
706 */
707
708 /* on demark tout les sommets non remuneroter */
709 /* -------------------------------------------------- */
710 /* Parameter adjustments */
711 --r;
712 --vois;
713 --ptvois;
714
715 /* Function Body */
716 i__1 = *n;
717 for (i = 1; i <= i__1; ++i) {
718 if (r[i] < 0) {
719 r[i] = 0;
720 }
721 /* L10: */
722 }
723
724 /* initialisation */
725
726 stk = *n ;// correct bug FH june 2011 ....
727 nv[0] = stk;
728 stk2 = stk;
729 *niveau = 0;
730 ++stk;
731 nv[stk] = *racine;
732 r[*racine] = -1;
733 L20:
734 if (stk2 < stk) {
735 ++(*niveau);
736 stk1 = stk2 + 1;
737 nvni=nv[*niveau];/* save value */
738 nv[*niveau] = stk;
739 stk2 = stk;
740 /* print *,' ------- niveau =',niveau,' stk=',stk1,stk2 */
741 i__1 = stk2;
742 for (k = stk1; k <= i__1; ++k) {
743 s = nv[k];
744 /* print *,'----------------- s=',s */
745 i__2 = ptvois[s + 1] - 1;
746 for (i = ptvois[s]; i <= i__2; ++i) {
747 /* pour tout les sommets (sv) voisin */
748 /* d'un sommet (s) du niveau precedent */
749 sv = vois[i];
750 /* print *,' voisin =',sv */
751 /* si le sommet n'est pas marque on le marque et
752 on l'ajout */
753 if (r[sv] == 0) {
754 ++stk;
755 nv[stk] = sv;
756 r[sv] = -1;
757 }
758 /* L30: */
759 }
760 /* L40: */
761 }
762 goto L20;
763 }
764 // if(nvni>0) nv[*niveau]=nvni;
765 --(*niveau);
766 /* call pnv(' gibbsd ',n,nv,niveau) */
767 return 0;} /* gibbsd_ */
768
769
gibbst_(integer * n,integer * p,integer * nv,integer * nn,integer * ptvois,integer * vois,integer * m,integer * r,integer * new_,integer * option,integer * pfnew,integer * impre,integer *)770 /* Subroutine */ int gibbst_(integer* n,integer* p,integer* nv,integer* nn,integer* ptvois,integer* vois,integer* m,integer* r,integer* new_,integer* option,
771 integer* pfnew,integer* impre,integer* )
772 {
773 /* System generated locals */
774 integer i__1, i__2, i__3, i__4, i__5;
775
776 /* Local variables */
777 static integer nbsc, bnew, knew, step, plus, i, j, k, s, debut, i1, i2;
778 /* extern Subroutine int gibbs2_();*/
779 static integer fin;
780
781
782 /* construction de la stucture de niveau dans nv a partir de nn */
783 /* ------------------------------------------------------------ */
784 /* Parameter adjustments */
785 --r;
786 --m;
787 --vois;
788 --ptvois;
789
790 /* Function Body */
791 nv[0] = *n;
792 i__1 = *p + 1;
793 for (i = 1; i <= i__1; ++i) {
794 nv[i] = 0;
795 /* L150: */
796 }
797 i__1 = *n;
798 for (i = 1; i <= i__1; ++i) {
799 if (nn[i] >= 0) {
800 ++nv[nn[i] + 1];
801 }
802 /* L160: */
803 }
804 i__1 = *p;
805 for (i = 0; i <= i__1; ++i) {
806 nv[i + 1] += nv[i];
807 /* L170: */
808 }
809 i__1 = *n;
810 for (i = 1; i <= i__1; ++i) {
811 if (nn[i] >= 0) {
812 j = nn[i];
813 ++nv[j];
814 nv[nv[j]] = i;
815 }
816 /* L180: */
817 }
818 for (i = *p; i >= 0; --i) {
819 nv[i + 1] = nv[i];
820 /* L190: */
821 }
822 nv[0] = *n;
823 nbsc = nv[*p + 1] - nv[0];
824 /* --- fin de la construction ------------------------------------ */
825 if (*option == -2) {
826 i__1 = *impre - 1;
827 }
828 i__1 = *n;
829 for (i = 1; i <= i__1; ++i) {
830 m[i] = *n * 3 + ptvois[i + 1] - ptvois[i];
831 /* L10: */
832 }
833 if ((((int)*option) == 1)||(((int)*option) == -1)) {
834 debut = 0;
835 fin = *p;
836 step = 1;
837 } else {
838 debut = *p;
839 fin = 0;
840 step = -1;
841 }
842 i__1 = fin;
843 i__2 = step;
844 for (i = debut; i__2 < 0 ? i >= i__1 : i <= i__1; i += i__2) {
845 i1 = nv[i] + 1;
846 i2 = nv[i + 1];
847 i__3 = i2 - i1 + 1;
848 gibbs2_(&i__3, &nv[i1], &m[1]);
849 i__3 = i2;
850 for (j = i1; j <= i__3; ++j) {
851 s = nv[j];
852 i__4 = ptvois[s + 1] - 1;
853 for (k = ptvois[s]; k <= i__4; ++k) {
854 /* Computing MIN */
855 i__5 = m[vois[k]];
856 m[vois[k]] = mmin(i__5,j);
857 /* L20: */
858 }
859 /* L30: */
860 }
861 /* L40: */
862 }
863 if (*option > 0) {
864 knew = *new_;
865 plus = 1;
866 } else {
867 knew = *new_ + nbsc + 1;
868 plus = -1;
869 }
870 *new_ += nbsc;
871 /* if(option.gt.0) then */
872 /* do 60 k = debut , fin , step */
873 /* do 60 j = nv(k+1),nv(k)+1,-1 */
874 /* knew = knew + plus */
875 /* r(nv(j)) = knew */
876 /* 60 continue */
877 /* else */
878 i__2 = fin;
879 i__1 = step;
880 for (k = debut; i__1 < 0 ? k >= i__2 : k <= i__2; k += i__1) {
881 i__3 = nv[k + 1];
882 for (j = nv[k] + 1; j <= i__3; ++j) {
883 knew += plus;
884 r[nv[j]] = knew;
885 /* L70: */
886 }
887 }
888 /* endif */
889 *pfnew = 0;
890 bnew = 0;
891 i__3 = *n;
892 for (i = 1; i <= i__3; ++i) {
893 k = r[i];
894 if (k > 0) {
895 i__1 = ptvois[i + 1] - 1;
896 for (j = ptvois[i]; j <= i__1; ++j) {
897 if (r[vois[j]] > 0) {
898 /* Computing MIN */
899 i__2 = k, i__4 = r[vois[j]];
900 k = mmin(i__2,i__4);
901 }
902 /* L100: */
903 }
904 *pfnew = *pfnew + r[i] - k + 1;
905 /* Computing MAX */
906 i__1 = bnew, i__2 = r[i] - k + 1;
907 bnew = mmax(i__1,i__2);
908 }
909 /* L110: */
910 }
911 /* if(impre.lt.0.or.impre.gt.2) then */
912 /* write(nfout,*) ' option =',option,', profile =',pfnew */
913 /* + ,', 1/2 bande =',bnew,', new=',new,', nbss composante=',nbsc
914 */
915 /* endif */
916 return 0;} /* gibbst_ */
917
gibbsv(integer * ptvoi,integer * vois,integer * lvois,integer * w,integer * v)918 /* function */ int Mesh::gibbsv (integer* ptvoi,
919 integer* vois,integer* lvois,integer* w,integer* v)
920 {
921 /* System generated locals */
922 integer i__2;
923
924 /* Local variables */
925 integer i, j, k, T, ss, iii, ptv, ptv1;
926 integer nbss = nv , nbt = nt;
927 /*--- Prepare les donees pour gibbsa en construisant ptvoi, vois, lvois -
928 ------------*/
929 /* in */
930 /* --- nbnt =3 pour des triangles 2D, */
931 /* nbt = nb de triangle */
932 /* nbss = nb de sommets */
933 /* nsea = numeros de 3 sommets de chaque triangle (me) */
934 /* out */
935 /* --- ptvoi, vois, lvois, err */
936 /* tableaux de travail w, v */
937 /*-----------------------------------------------------------------------
938 ----------*/
939 /* Parameter adjustments */
940 --v;
941 --w;
942 --vois;
943 --ptvoi;
944
945 /* Function Body */
946 for (i = 1; i <= nbss; ++i) {
947 w[i] = -1;
948 ptvoi[i] = 0; }
949 ptvoi[nbss + 1] = 0;
950 for (i = 0; i < nbt; ++i) {
951 for (j = 0; j < 3; ++j) {
952 ss = number(triangles[i][j])+1;
953 ++ptvoi[ss + 1];
954 w[ss] = 0;}
955 }
956
957 for (i = 1; i <= nbss; ++i)
958 ptvoi[i + 1] += ptvoi[i];
959
960 for (i = 0; i < nbt; ++i) {
961 for (j = 0; j < 3; ++j) {
962 ss = number(triangles[i][j])+1;
963 ++ptvoi[ss];
964 v[ptvoi[ss]] = i;
965 }
966 }
967 ptv1 = 0;
968 iii = 1;
969 for (i = 1; i <= nbss; ++i) {
970 ptv = ptv1 + 1;
971 ptv1 = ptvoi[i];
972 ptvoi[i] = iii;
973 i__2 = ptv1;
974 for (j = ptv; j <= i__2; ++j) {
975 T = v[j];
976 for (k = 0; k < 3; ++k) {
977 ss = number(triangles[T][k])+1; /* nsea[k + T * nsea_dim1]; */
978 if (w[ss] != i) {
979 w[ss] = i;
980 if (iii > *lvois) return 2 ;
981 /* print*,'pas assez de place memoire' */
982
983 vois[iii] = ss;
984 ++iii;}
985 }
986 }
987 }
988 ptvoi[nbss + 1] = iii;
989 *lvois = iii - 1;
990 return 0; /* OK */
991 return 0;} /* gibbsv_ */
992
renum()993 int Mesh::renum()
994 /* --------
995 renumber vertices by gibbs method; updates triangle and edge array
996 in: mesh
997 out: mesh
998 auxiliary arrays: ptvois,vois,r,m,nv,nx,ny,nn,w1,w2,f
999 all of size nv+1 except vois (10(nv+1)) and nv (2(nv+1))
1000 err = -1 : memory alloc pb; err = -3: fatal erreur gibbs 2 : pb racine
1001 */
1002 {
1003 long pfold=0, pfnew=0;
1004 long* ptvois=NULL;
1005 long* vois=NULL;
1006 long* nn =NULL;
1007 long* r =NULL;
1008 long* m =NULL;
1009 long* nnv =NULL;
1010 long* nx =NULL;
1011 long* ny =NULL;
1012 long* w1 =NULL;
1013 long* w2=NULL;
1014 long nbvoisin = 10*nv;
1015 long printint=0, iodev=6;
1016 int err=0;
1017 ptvois = new long[nv+1]; //(long*)calloc((long)(nv + 1) , sizeof(long));
1018 nn = new long[3*nt]; //(long*)calloc(3 * nt ,sizeof(long));
1019 vois = new long[nbvoisin+10]; //(long*)calloc((long)(nbvoisin + 10) , sizeof(long));
1020 r = new long[nv+1]; //(long*)calloc((long)(nv + 1) , sizeof(long));
1021 if((!ptvois)||(!nn)||(!vois)||(!r)) return -1;
1022 err = gibbsv(ptvois,vois,&nbvoisin,r,nn) ;
1023 delete [] nn; // free(nn);
1024 if(err==0)
1025 {
1026 m = new long[nv+1];
1027 nn = new long[nv+1];
1028 nnv = new long[(nv+1)<<1];
1029 nx = new long[nv+1];
1030 ny = new long[nv+1];
1031 w1 = new long[nv+1];
1032 w2 = new long[nv+1];
1033 long lnv = nv;
1034 err = gibbsa_ (&lnv, ptvois, vois, r, m, nnv, nx, ny, nn, w1, w2, &pfold, &pfnew,
1035 &printint, &iodev);
1036 delete [] m;
1037 delete [] nnv;
1038 delete [] nn;
1039 delete [] nx;
1040 delete [] ny;
1041 delete [] w1;
1042 delete [] w2;
1043 }
1044
1045 delete [] vois;
1046 delete [] ptvois;
1047 if(verbosity>1)
1048 cout << " -- Mesh: Gibbs: old skyline = " << pfold << " new skyline = " << pfnew << endl;
1049 if (err == 0 && (pfnew <= pfold))
1050 {
1051 int i,j;
1052 for ( i=0;i<nv;i++)
1053 r[i] -= 1;
1054 Vertex * f= new Vertex[nv];
1055 for (i = 0; i < nv; ++i)
1056 f[i] = vertices[i];
1057 for (i = 0; i < nv; ++i)
1058 vertices[r[i]] = f[i];
1059 delete [] f;
1060
1061 for (j = 0; j < nt; ++j) // updates triangle array
1062 triangles[j].Renum(vertices,r);
1063
1064 for (j = 0; j < neb; ++j) // updates edge array
1065 bedges[j].Renum(vertices,r);
1066
1067 if (BoundaryAdjacencesHead )
1068 {
1069 for (int i=0;i<nv;i++)
1070 BoundaryAdjacencesHead[i]=-1;
1071 int j2=0;
1072 for (int j=0;j<neb;j++)
1073 for (int k=0;k<2;k++,j2++)
1074 {
1075 int v = number(bedges[j][k]);
1076 ffassert(v >=0 && v < nv);
1077 BoundaryAdjacencesLink[j2]=BoundaryAdjacencesHead[v];
1078 BoundaryAdjacencesHead[v]=j2;
1079 }
1080 }
1081 for (int it=0;it<nt;it++)
1082 for (int j=0;j<3;j++)
1083 TriangleConteningVertex[(*this)(it,j)]=it;
1084
1085 if (quadtree) {
1086 delete quadtree;quadtree=0;
1087 MakeQuadTree(); }
1088
1089
1090 }
1091 delete [] r;
1092 return err;
1093 }
renum()1094 int FESpace::renum()
1095 /* --------
1096 renumber vertices by gibbs method; updates triangle and edge array
1097 in: mesh
1098 out: mesh
1099 auxiliary arrays: ptvois,vois,r,m,nv,nx,ny,nn,w1,w2,f
1100 all of size nv+1 except vois (10(nv+1)) and nv (2(nv+1))
1101 err = -1 : memory alloc pb; err = -3: fatal erreur gibbs 2 : pb racine
1102 */
1103 {
1104 if (cdef==0) return -2;
1105 if (cdef->NodesOfElement==0) return -2;
1106 int nv = NbOfNodes;
1107 long pfold =0, pfnew=0;
1108 long* ptvois=NULL;
1109 long* vois=NULL;
1110 long* nn =NULL;
1111 long* r =NULL;
1112 long* m =NULL;
1113 long* nnv =NULL;
1114 long* nx =NULL;
1115 long* ny =NULL;
1116 long* w1 =NULL;
1117 long* w2=NULL;
1118 long nbvperelem = 20; // pour le P1
1119 // cout << "gibbs: nbvperelem =" << nbvperelem << endl;
1120 long nbvoisin = (nbvperelem)*nv;
1121 long printint=0, iodev=6;
1122 int err=0;
1123 int nnx= SizeToStoreAllNodeofElement();
1124 ptvois = new long[nv+1];
1125 nn = new long[nnx];
1126 vois = new long[nbvoisin+100];
1127 r = new long[nv+1];
1128 if((!ptvois)||(!nn)||(!vois)||(!r)) return -1;
1129 err = gibbsv(ptvois,vois,&nbvoisin,r,nn) ;
1130 delete [] nn;
1131 if(nbvoisin <= nv) // no renumbering if no enought neighbour FH. nov 2017
1132 // thank to beniamin.bogosel@cmap.polytechnique.fr
1133 {
1134 err=1;
1135 if(verbosity>1) cout << " No renumbering ( nbvois " << nbvoisin << " <= " <<nv <<" nb nodes )" << endl;
1136 }
1137 else if(err==0 )
1138 {
1139 m = new long[nv+1];
1140 nn = new long[nv+1];
1141 nnv = new long[(nv+1)<<1];
1142 nx = new long[nv+1];
1143 ny = new long[nv+1];
1144 w1 = new long[nv+1];
1145 w2 = new long[nv+1];
1146 long lnv = nv;
1147 err = gibbsa_ (&lnv, ptvois, vois, r, m, nnv, nx, ny, nn, w1, w2, &pfold, &pfnew,
1148 &printint, &iodev);
1149 delete [] m;
1150 delete [] nnv;
1151 delete [] nn;
1152 delete [] nx;
1153 delete [] ny;
1154 delete [] w1;
1155 delete [] w2;
1156 }
1157 else
1158 cerr << " Not enought memory (bug)" << nbvoisin << " " << nbvoisin/NbOfElements << endl;
1159 delete [] vois;
1160 delete [] ptvois;
1161 if (err==0 && verbosity>1)
1162 cout << " FESpace:Gibbs: old skyline = " << pfold << " new skyline = " << pfnew << endl;
1163 if (err == 0 && (pfnew <= pfold))
1164 {
1165 int i;
1166 // cout << *this << endl;
1167 for ( i=0;i<nv;i++)
1168 r[i] -= 1;
1169 cdef->renum(r,nnx);
1170 }
1171
1172 delete [] r;
1173 return err;
1174 }
gibbsv(integer * ptvoi,integer * vois,integer * lvois,integer * w,integer * v)1175 /* function */ int FESpace::gibbsv (integer* ptvoi,
1176 integer* vois,integer* lvois,integer* w,integer* v)
1177 {
1178 /* System generated locals */
1179 integer i__2;
1180
1181 /* Local variables */
1182
1183 integer i, j, k, T, ss, iii, ptv, ptv1;
1184 integer nbss = NbOfNodes , nbt = NbOfElements;
1185 integer jj,kk;
1186 /*--- Prepare les donees pour gibbsa en construisant ptvoi, vois, lvois -
1187 ------------*/
1188 /* in */
1189 /* --- nbnt =3 pour des triangles 2D, */
1190 /* nbt = nb de triangle */
1191 /* nbss = nb de sommets */
1192 /* nsea = numeros de 3 sommets de chaque triangle (me) */
1193 /* out */
1194 /* --- ptvoi, vois, lvois, err */
1195 /* tableaux de travail w, v */
1196 /*-----------------------------------------------------------------------
1197 ----------*/
1198 /* Parameter adjustments */
1199 --v;
1200 --w;
1201 --vois;
1202 --ptvoi;
1203
1204 /* Function Body */
1205 for (i = 1; i <= nbss; ++i) {
1206 w[i] = -1;
1207 ptvoi[i] = 0; }
1208 ptvoi[nbss + 1] = 0;
1209 for (i = 0; i < nbt; ++i) {
1210 for (j = 0; j < NbOfNodesInElement(i) ; ++j) {
1211 ss = (*this)(i,j) +1;
1212 ++ptvoi[ss + 1];
1213 w[ss] = 0;}
1214 }
1215
1216 for (i = 1; i <= nbss; ++i)
1217 ptvoi[i + 1] += ptvoi[i];
1218
1219 for (i = 0; i < nbt; ++i) {
1220 for (j = 0,jj= NbOfNodesInElement(i); j <jj ; ++j) {
1221 ss = (*this)(i,j) +1;//number(triangles[i][j])+1;
1222 if ( ! ( ss>0 && ss <= nbss))
1223 {
1224 cout << "bug " << ss << " " << i << " " << j << endl;
1225 exit(1);
1226 }
1227 ++ptvoi[ss];
1228 v[ptvoi[ss]] = i;
1229 }
1230 }
1231 ptv1 = 0;
1232 iii = 1;
1233 for (i = 1; i <= nbss; ++i) {
1234 ptv = ptv1 + 1;
1235 ptv1 = ptvoi[i];
1236 ptvoi[i] = iii;
1237 i__2 = ptv1;
1238 for (j = ptv; j <= i__2; ++j) {
1239 T = v[j];
1240 for (k = 0,kk=NbOfNodesInElement(T); k < kk; ++k) {
1241 ss = (*this)(T,k) +1;//number(triangles[T][k])+1; /* nsea[k + T * nsea_dim1]; */
1242 if (w[ss] != i) {
1243 w[ss] = i;
1244 if (iii > *lvois) return 2 ;
1245 /* print*,'pas assez de place memoire' */
1246
1247 vois[iii] = ss;
1248 ++iii;}
1249 }
1250 }
1251 }
1252 ptvoi[nbss + 1] = iii;
1253 *lvois = iii - 1;
1254 return 0; /* OK */
1255 return 0;} /* gibbsv_ */
1256 /* message d'erreur: *err = 2; print*,'pas assez de place memoire' */
1257
1258 /* message d'erreur: *err = 2; print*,'pas assez de place memoire' */
1259