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
29 #include <math.h>
30 #include <stdio.h>
31 #include <stdlib.h>
32 #include <iostream>
33 #include <fstream>
34 using namespace std;
35 using namespace bamg;
36 #include "Mesh2.h"
37
38 #define mmax(a,b)(a>b?a:b)
39 #define mmin(a,b)(a<b?a:b)
40 #define ffalse 0
41 #define ttrue 1
42
43 /* -- translated by f2c (version of 23 May 1992 14:18:33).
44 You must link the resulting object file with the libraries:
45 -lF77 -lI77 -lm -lc (in that order)
46 */
47
48
49 #define integer long
50 #define logical long
51
52 int gibbs1_(integer* n,integer* record,integer* ptvois);
53 int gibbs2_(integer* n,integer* record,integer* criter);
54 int gibbsa_(integer* n,integer* ptvois,integer* vois,integer* r,integer* m,
55 integer* nv,integer* nx,integer* ny,integer* nn,integer* w1,integer* w2,
56 integer* pfold,integer* pfnew,integer* impre,integer* nfout);
57 int gibbsb_(integer* x,integer* y,integer* n,integer* ptvois,
58 integer* vois,integer* nx,integer* ny,integer* nv,integer* nn,integer* m,
59 integer* wh,integer* wl,integer* r, integer* impre, integer* nfout);
60 int gibbsc_(integer* nz,integer* nv,integer* niveau,integer* n,integer* );
61 int gibbsd_(integer* racine,integer* n,integer* ptvois,integer*
62 vois,integer* nv,integer* r,integer* niveau);
63 int gibbst_(integer* n,integer* p,integer* nv,integer* nn,integer* ptvois,integer* vois,
64 integer* m,integer* r,integer* new_,integer* option,
65 integer* pfnew,integer* impre,integer* nfout);
66
gibbs1_(integer * n,integer * record,integer * ptvois)67 /* Subroutine */ int gibbs1_(integer* n,integer* record,integer* ptvois)
68 {
69 /* System generated locals */
70 integer i__1;
71
72 /* Local variables */
73 static integer crit, i, j, l, r, rec;
74
75 /* -----------------------------------------------------------------------
76 */
77 /* routine appele par gibbs0 */
78 /* -----------------------------------------------------------------------
79 */
80 /* but: trie record (ensemble de n sommet) telle que l'ordre des somm
81 */
82 /* soit croissant (ordre du sommet i est ptvois(i+1)-ptvois(i)) */
83 /* -----------------------------------------------------------------------
84 */
85
86 /* Parameter adjustments */
87 --ptvois;
88 --record;
89
90 /* Function Body */
91 if (*n <= 1) {
92 return 0;
93 }
94 l = *n / 2 + 1;
95 r = *n;
96 L2:
97 if (l <= 1) {
98 goto L20;
99 }
100 --l;
101 rec = record[l];
102 crit = ptvois[record[l] + 1] - ptvois[record[l]];
103 goto L3;
104 L20:
105 rec = record[r];
106 crit = ptvois[record[r] + 1] - ptvois[record[r]];
107 record[r] = record[1];
108 --r;
109 if (r == 1) {
110 goto L999;
111 }
112 L3:
113 j = l;
114 L4:
115 i = j;
116 j <<= 1;
117 if ((i__1 = j - r) < 0) {
118 goto L5;
119 } else if (i__1 == 0) {
120 goto L6;
121 } else {
122 goto L8;
123 }
124 L5:
125 if (ptvois[record[j] + 1] - ptvois[record[j]] < ptvois[record[j + 1] + 1]
126 - ptvois[record[j + 1]]) {
127 ++j;
128 }
129 L6:
130 if (crit >= ptvois[record[j] + 1] - ptvois[record[j]]) {
131 goto L8;
132 }
133 record[i] = record[j];
134 goto L4;
135 L8:
136 record[i] = rec;
137 goto L2;
138 L999:
139 record[1] = rec;
140 return 0;
141 } /* gibbs1_ */
142
gibbs2_(integer * n,integer * record,integer * criter)143 /* Subroutine */ int gibbs2_(integer* n,integer* record,integer* criter)
144 {
145 static integer crit, i, j, l, r, rec;
146
147
148 /* trie record selon les valeurs de criter(record(.)) croissantes */
149
150
151 /* Parameter adjustments */
152 --criter;
153 --record;
154
155 /* Function Body */
156 if (*n <= 1) {
157 return 0;
158 }
159 l = *n / 2 + 1;
160 r = *n;
161 L2:
162 if (l <= 1) {
163 goto L20;
164 }
165 --l;
166 rec = record[l];
167 crit = criter[rec];
168 goto L3;
169 L20:
170 rec = record[r];
171 crit = criter[rec];
172 record[r] = record[1];
173 --r;
174 if (r == 1) {
175 goto L999;
176 }
177 L3:
178 j = l;
179 L4:
180 i = j;
181 j <<= 1;
182 if (j - r < 0) {
183 goto L5;
184 } else if (j == r) {
185 goto L6;
186 } else {
187 goto L8;
188 }
189 L5:
190 if (criter[record[j]] < criter[record[j + 1]]) {
191 ++j;
192 }
193 L6:
194 if (crit >= criter[record[j]]) {
195 goto L8;
196 }
197 record[i] = record[j];
198 goto L4;
199 L8:
200 record[i] = rec;
201 goto L2;
202 L999:
203 record[1] = rec;
204 return 0;
205 } /* gibbs2_ */
206
207 /* 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)208 int gibbsa_(integer* n,integer* ptvois,integer* vois,integer* r,integer* m,
209 integer* nv,integer* nx,integer* ny,integer* nn,integer* w1,integer* w2,
210 integer* pfold,integer* pfnew,integer* impre,integer* nfout)
211 {
212 /* System generated locals */
213 integer i__1, i__2, i__3, i__4;
214
215 /* Builtin functions */
216 /* Subroutine */ int s_stop();
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
225 /* -----------------------------------------------------------------------
226 */
227 /* but: calculer une renumerotation des sommets d'un graphe defini par:
228 */
229 /* par la methode de gibbs */
230 /* -----------------------------------------------------------------------
231 */
232 /* entree */
233 /* -------- */
234 /* n = nb de sommet du graphe */
235 /* les voisins d'un sommet i ont pour numero : */
236 /* ( vois(j) , j=ptvois(i),ptvois(i+1)-1 ) */
237
238 /* impre parametre d'impression */
239 /* nfout numero du fichier pour impression */
240
241 /* sortie */
242 /* ------ */
243 /* r(1:n) tableau donnant la nouvelle numerotation: */
244 /* r(i) = nouveau numero du sommet i */
245 /* pfolf = ancien profile */
246 /* pfnew = nouveau profile */
247
248 /* tableau de travail : */
249 /* -------------------- */
250 /* m(n) */
251 /* nv(0:n+n) */
252 /* nx(n) */
253 /* ny(n) */
254 /* nn(0:n) */
255 /* w1(n) */
256 /* w2(n) */
257
258 /* -----------------------------------------------------------------------
259 */
260 /* programmeur f. hecht le 3/02/1987 */
261 /* -----------------------------------------------------------------------
262 */
263
264 /* tri des voisins d'un sommet du graphe par degre croissant */
265 /* --------------------------------------------------------------- */
266 /* Parameter adjustments */
267 --w2;
268 --w1;
269 --ny;
270 --nx;
271 --m;
272 --r;
273 --vois;
274 --ptvois;
275
276 /* Function Body */
277 p2 = ptvois[1] - 1;
278 i__1 = *n;
279 for (i = 1; i <= i__1; ++i) {
280 p1 = p2 + 1;
281 p2 = ptvois[i + 1] - 1;
282 i__2 = p2 - p1 + 1;
283 gibbs1_(&i__2, &vois[p1], &ptvois[1]);
284 /* if(impre.le.-9) then */
285 /* write (nfout,*) 'les voisin de ',i,'sont: ', (vois(j),j=p1,p
286 2) */
287 /* endif */
288 /* L10: */
289 }
290 i__1 = *n;
291 for (i = 1; i <= i__1; ++i) {
292 r[i] = 0;
293 /* L20: */
294 }
295 /* boucle sur les composante connexe du graphe */
296 new_ = 0;
297 nbcc = 0;
298 L30:
299 if (new_ < *n) {
300 ++nbcc;
301 /* recherche d'une racine y (un sommet non numerote) de degree m
302 ini */
303 y = 0;
304 degre = *n + 1;
305 i__1 = *n;
306 for (i = 1; i <= i__1; ++i) {
307 if (r[i] <= 0) {
308 degi = ptvois[i + 1] - ptvois[i];
309 if (degi < degre) {
310 degre = degi;
311 y = i;
312 }
313 }
314 /* L40: */
315 }
316 if (y == 0) {
317 return -3;/* s_stop("fatal erreur gibbs 2 : pb racine", 33L); */
318 }
319 gibbsd_(&y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
320 nbpass = 0;
321 L50:
322 ++nbpass;
323 x = y;
324 p = niveau;
325 k = 0;
326 i__1 = nv[p + 1];
327 for (i = nv[p] + 1; i <= i__1; ++i) {
328 ++k;
329 m[k] = nv[i];
330 /* L60: */
331 }
332 gibbs1_(&k, &m[1], &ptvois[1]);
333 i__1 = k;
334 for (i = 1; i <= i__1; ++i) {
335 y = m[i];
336 gibbsd_(&y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
337 if (niveau > p) {
338 goto L50;
339 }
340 /* L70: */
341 }
342 y = m[1];
343 /* if(impre.lt.0) then */
344 /* write(nfout,*) */
345 /* + ' nb de pass pour trouver le pseudo diametre',nbpass */
346 /* + ,' x=',x,',y=',y,' de la composante connexe ',nbcc */
347
348 /* write (nfout,*) ('-',i=1,78) */
349 /* endif */
350 /* optimisation de la descendance de la numerotation */
351 /* ------------------------------------------------- */
352 gibbsb_(&x, &y, n, &ptvois[1], &vois[1], &nx[1], &ny[1], nv, nn, &m[1]
353 , &w1[1], &w2[1], &r[1], impre, nfout);
354
355 /* renumerotation de cuthill mac kee avec la meilleur des 4 option
356 s */
357 /* --------------------------------------------------------------
358 --- */
359 pf = 1073741824;
360 option = -2;
361 new1 = new_;
362 for (opt = -2; opt <= 2; ++opt) {
363 new_ = new1;
364 if (opt != 0) {
365 gibbst_(n, &p, nv, nn, &ptvois[1], &vois[1], &m[1], &r[1], &
366 new_, &opt, &pf1, impre, nfout);
367 if (pf1 < pf) {
368 pf = pf1;
369 option = opt;
370 }
371 }
372 /* L80: */
373 }
374 /* if(impre.ne.0) write (nfout,*) ' on a choisi l''option ',
375 */
376 /* + option,', new =',new */
377 new_ = new1;
378 gibbst_(n, &p, nv, nn, &ptvois[1], &vois[1], &m[1], &r[1], &new_, &
379 option, &pf1, impre, nfout);
380 goto L30;
381 }
382 /* if(impre.ne.0) write(nfout,*) */
383 /* + ' nb de composante connexe du graphe =',nbcc */
384 /* calcul du profile */
385 *pfold = 0;
386 *pfnew = 0;
387 bnew = 0;
388 bold = 0;
389 i__1 = *n;
390 for (i = 1; i <= i__1; ++i) {
391 old = i;
392 new_ = r[i];
393 i__2 = ptvois[i + 1] - 1;
394 for (j = ptvois[i]; j <= i__2; ++j) {
395 /* Computing MIN */
396 i__3 = old, i__4 = vois[j];
397 old = mmin(i__3,i__4);
398 /* Computing MIN */
399 i__3 = new_, i__4 = r[vois[j]];
400 new_ = mmin(i__3,i__4);
401 /* L100: */
402 }
403 *pfold = *pfold + i - old + 1;
404 /* Computing MAX */
405 i__2 = bold, i__3 = i - old + 1;
406 bold = mmax(i__2,i__3);
407 *pfnew = *pfnew + r[i] - new_ + 1;
408 /* Computing MAX */
409 i__2 = bnew, i__3 = r[i] - new_ + 1;
410 bnew = mmax(i__2,i__3);
411 /* L110: */
412 }
413 /* if(impre.ne.0) then */
414 /* write(nfout,*)'profile old = ',pfold,', profile new = ',pfnew
415 */
416 /* write(nfout,*)'1/2 bande old = ',bold ,', 1/2 band new = ',bnew
417 */
418 /* endif */
419 return 0;
420 } /* gibbsa_ */
421
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)422 /* Subroutine */ int gibbsb_(integer* x,integer* y,integer* n,integer* ptvois,
423 integer* vois,integer* nx,integer* ny,integer* nv,integer* nn,integer* m,
424 integer* wh,integer* wl,integer* r, integer* impre, integer* nfout)
425 {
426 /* System generated locals */
427 integer i__1, i__2;
428
429 /* Local variables */
430 static flag_;
431 static integer i, j, k, p, s, h0, i1, l0, i2;
432 /* extern Subroutine int gibbs1_(); */
433 static integer lg;
434 /* extern Subroutine int gibbsd_(), gibbsc_();*/
435 static integer niveau, mxcanx, mxcany, nbc;
436
437 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
438 */
439 /* ......................................................................
440 */
441 /* attention on met la descente optimiser dans r <0 ou nulle */
442 /* .......................................................................
443 */
444 /* Parameter adjustments */
445 --r;
446 --m;
447 --ny;
448 --nx;
449 --vois;
450 --ptvois;
451
452 /* Function Body */
453 gibbsd_(y, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
454 gibbsc_(&ny[1], nv, &niveau, n, &mxcany);
455 gibbsd_(x, n, &ptvois[1], &vois[1], nv, &r[1], &niveau);
456 p = niveau;
457 gibbsc_(&nx[1], nv, &niveau, n, &mxcanx);
458 flag_ = ffalse;
459 i__1 = *n;
460 for (i = 1; i <= i__1; ++i) {
461 if (nx[i] + ny[i] == p) {
462 r[i] = -nx[i];
463 } else if (nx[i] >= 0) {
464 flag_ = ttrue;
465 r[i] = -1073741824;
466 } else {
467 if (r[i] <= 0) {
468 r[i] = -1073741822;
469 }
470 }
471 /* L20: */
472 }
473 if (flag_) {
474 /* calcul des composantes connexe du graphe sans les sommets de
475 nn */
476 /* ------------------------------------------------------------
477 --- */
478 j = *n;
479 k = 0;
480 nbc = 0;
481 nv[nbc] = j;
482 L30:
483 ++k;
484 if (k <= *n) {
485 if (r[k] == -1073741824) {
486 /* recherche de la fermeture transitive partant de k
487 */
488 ++nbc;
489 i = -1;
490 s = k;
491 L40:
492 ++i;
493 wl[i] = ptvois[s];
494 wh[i] = ptvois[s + 1];
495 ++j;
496 nv[j] = s;
497 r[s] = -1073741823;
498 L50:
499 if (i >= 0) {
500 if (wl[i] < wh[i]) {
501 s = vois[wl[i]];
502 ++wl[i];
503 if (r[s] == -1073741824) {
504 goto L40;
505 }
506 goto L50;
507 }
508 --i;
509 goto L50;
510 }
511 nv[nbc] = j;
512 m[nbc] = nbc;
513 }
514 goto L30;
515 }
516 /* if(impre.lt.0) write(nfout,*) */
517 /* + ' nb de composante connexe du graphe reduit =',nbc */
518
519 /* --------------- fin de construction des composantes connexes------
520 --- */
521 /* nv(0)=n */
522 /* if(impre.le.-10) write(nfout,5555)'nv(0:n+n) = ',(nv(i),i=0,
523 n+n) */
524 gibbs1_(&nbc, &m[1], nv);
525 /* if(impre.le.-10)write(nfout,5555)'trie m =',(m(i),i=1,nbc)
526 */
527 i__1 = p;
528 for (i = 0; i <= i__1; ++i) {
529 nn[i] = 0;
530 /* L60: */
531 }
532 i__1 = *n;
533 for (i = 1; i <= i__1; ++i) {
534 j = -r[i];
535 if (j >= 0 && j <= p) {
536 ++nn[j];
537 }
538 /* L70: */
539 }
540
541 /* boucle sur les composante connexes par ordre croissantes */
542 /* -------------------------------------------------------- */
543 for (k = nbc; k >= 1; --k) {
544 i = m[k];
545 i1 = nv[i - 1] + 1;
546 i2 = nv[i];
547 lg = i2 - i1 + 1;
548 /* if(impre.le.-7) */
549 /* + write(nfout,*) k,' composante ',i,',lg=',lg,',i1,i2
550 =',i1,i2 */
551 /* if(impre.le.-8) */
552 /* + write (nfout,5555)' ',(nv(i),i=i1,i2) */
553 h0 = 0;
554 l0 = 0;
555 i__1 = p;
556 for (j = 0; j <= i__1; ++j) {
557 wh[j] = nn[j];
558 wl[j] = nn[j];
559 /* L90: */
560 }
561 i__1 = i2;
562 for (i = i1; i <= i__1; ++i) {
563 s = nv[i];
564 ++wh[nx[s]];
565 ++wl[p - ny[s]];
566 /* L100: */
567 }
568 i__1 = p;
569 for (j = 0; j <= i__1; ++j) {
570 if (wh[j] != nn[j]) {
571 /* Computing MAX */
572 i__2 = wh[j];
573 h0 = mmax(i__2,h0);
574 }
575 if (wl[j] != nn[j]) {
576 /* Computing MAX */
577 i__2 = wl[j];
578 l0 = mmax(i__2,l0);
579 }
580 /* L110: */
581 }
582 if (h0 < l0 || h0 == l0 && mxcanx <= mxcany) {
583 /* if(impre.le.-2) write(nfout,*) */
584 /* + ' h0 = ',h0,',l0 = ',l0,' ------- XXXX
585 --------' */
586 i__1 = i2;
587 for (i = i1; i <= i__1; ++i) {
588 s = nv[i];
589 r[s] = -nx[s];
590 ++nn[-r[s]];
591 /* L120: */
592 }
593 } else {
594 /* if (impre.le.-2) write(nfout,*) */
595 /* + ' h0 = ',h0,',l0 = ',l0,' ------- YYYY
596 --------' */
597 i__1 = i2;
598 for (i = i1; i <= i__1; ++i) {
599 s = nv[i];
600 r[s] = -p + ny[s];
601 ++nn[-r[s]];
602 /* L130: */
603 }
604 }
605 /* L140: */
606 }
607 }
608 /* on met les nouveaux niveaux de la descendance optimiser dans nn */
609 /* -----------------------------------------------------------------
610 */
611 i__1 = *n;
612 for (i = 1; i <= i__1; ++i) {
613 if (r[i] > 0) {
614 nn[i] = -1;
615 } else if (r[i] == -1073741822) {
616 nn[i] = -2;
617 } else {
618 nn[i] = -r[i];
619 }
620 /* L150: */
621 }
622 /* if(impre.le.-10)write (nfout,5555)' nn(i)=',(nn(i),i=1,n) */
623 /* 5555 format(' -------- ',a,/,5(15x,10(i5)/)) */
624 return 0;} /* gibbsb_ */
625
gibbsc_(integer * nz,integer * nv,integer * niveau,integer * n,integer * mxz)626 /* Subroutine */ int gibbsc_(integer* nz,integer* nv,integer* niveau,integer* n,integer* mxz)
627 {
628 /* System generated locals */
629 integer i__1, i__2, i__3;
630
631 /* Local variables */
632 static integer i, j;
633
634 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
635 */
636 /* Parameter adjustments */
637 --nz;
638
639 /* Function Body */
640 i__1 = *n;
641 for (i = 1; i <= i__1; ++i) {
642 nz[i] = -1;
643 /* L10: */
644 }
645 *mxz = 0;
646 i__1 = *niveau;
647 for (i = 0; i <= i__1; ++i) {
648 /* Computing MAX */
649 i__2 = *mxz, i__3 = nv[i + 1] - nv[i];
650 *mxz = mmax(i__2,i__3);
651 i__2 = nv[i + 1];
652 for (j = nv[i] + 1; j <= i__2; ++j) {
653 nz[nv[j]] = i;
654 /* L20: */
655 }
656 }
657 return 0;} /* gibbsc_ */
658
gibbsd_(integer * racine,integer * n,integer * ptvois,integer * vois,integer * nv,integer * r,integer * niveau)659 /* Subroutine */ int gibbsd_(integer* racine,integer* n,integer* ptvois,integer* vois,integer* nv,integer* r,integer* niveau)
660 {
661 /* System generated locals */
662 integer i__1, i__2;
663
664 /* Local variables */
665 static integer i, k, s, sv, stk, stk1, stk2;
666
667 /* -----------------------------------------------------------------------
668 */
669 /* but construire la structure des descendant de racine du graphe */
670 /* -----------------------------------------------------------------------
671 */
672 /* sortie : */
673 /* -------- */
674 /* nv est la structure des niveaux */
675 /* les sommets du niveau (i =0,niveau_ sont defini par : */
676 /* (nv(j),j=nv(i),nv(i+1)-1) */
677
678 /* le tableau r(i) n'est modifier que sur les sommets */
679 /* de la composante connexe du graphe contenant la racine */
680
681 /* -----------------------------------------------------------------------
682 */
683
684 /* on demark tout les sommets non remuneroter */
685 /* -------------------------------------------------- */
686 /* Parameter adjustments */
687 --r;
688 --vois;
689 --ptvois;
690
691 /* Function Body */
692 i__1 = *n;
693 for (i = 1; i <= i__1; ++i) {
694 if (r[i] < 0) {
695 r[i] = 0;
696 }
697 /* L10: */
698 }
699
700 /* initialisation */
701
702 stk = *n - 1;
703 nv[0] = stk;
704 stk2 = stk;
705 *niveau = 0;
706 ++stk;
707 nv[stk] = *racine;
708 r[*racine] = -1;
709 L20:
710 if (stk2 < stk) {
711 ++(*niveau);
712 stk1 = stk2 + 1;
713 nv[*niveau] = stk;
714 stk2 = stk;
715 /* print *,' ------- niveau =',niveau,' stk=',stk1,stk2 */
716 i__1 = stk2;
717 for (k = stk1; k <= i__1; ++k) {
718 s = nv[k];
719 /* print *,'----------------- s=',s */
720 i__2 = ptvois[s + 1] - 1;
721 for (i = ptvois[s]; i <= i__2; ++i) {
722 /* pour tout les sommets (sv) voisin */
723 /* d'un sommet (s) du niveau precedent */
724 sv = vois[i];
725 /* print *,' voisin =',sv */
726 /* si le sommet n'est pas marque on le marque et
727 on l'ajout */
728 if (r[sv] == 0) {
729 ++stk;
730 nv[stk] = sv;
731 r[sv] = -1;
732 }
733 /* L30: */
734 }
735 /* L40: */
736 }
737 goto L20;
738 }
739 --(*niveau);
740 /* call pnv(' gibbsd ',n,nv,niveau) */
741 return 0;} /* gibbsd_ */
742
743
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 * nfout)744 /* Subroutine */ int gibbst_(integer* n,integer* p,integer* nv,integer* nn,integer* ptvois,integer* vois,integer* m,integer* r,integer* new_,integer* option,
745 integer* pfnew,integer* impre,integer* nfout)
746 {
747 /* System generated locals */
748 integer i__1, i__2, i__3, i__4, i__5;
749
750 /* Local variables */
751 static integer nbsc, bnew, knew, step, plus, i, j, k, s, debut, i1, i2;
752 /* extern Subroutine int gibbs2_();*/
753 static integer fin;
754
755
756 /* construction de la stucture de niveau dans nv a partir de nn */
757 /* ------------------------------------------------------------ */
758 /* Parameter adjustments */
759 --r;
760 --m;
761 --vois;
762 --ptvois;
763
764 /* Function Body */
765 nv[0] = *n;
766 i__1 = *p + 1;
767 for (i = 1; i <= i__1; ++i) {
768 nv[i] = 0;
769 /* L150: */
770 }
771 i__1 = *n;
772 for (i = 1; i <= i__1; ++i) {
773 if (nn[i] >= 0) {
774 ++nv[nn[i] + 1];
775 }
776 /* L160: */
777 }
778 i__1 = *p;
779 for (i = 0; i <= i__1; ++i) {
780 nv[i + 1] += nv[i];
781 /* L170: */
782 }
783 i__1 = *n;
784 for (i = 1; i <= i__1; ++i) {
785 if (nn[i] >= 0) {
786 j = nn[i];
787 ++nv[j];
788 nv[nv[j]] = i;
789 }
790 /* L180: */
791 }
792 for (i = *p; i >= 0; --i) {
793 nv[i + 1] = nv[i];
794 /* L190: */
795 }
796 nv[0] = *n;
797 nbsc = nv[*p + 1] - nv[0];
798 /* --- fin de la construction ------------------------------------ */
799 if (*option == -2) {
800 i__1 = *impre - 1;
801 }
802 i__1 = *n;
803 for (i = 1; i <= i__1; ++i) {
804 m[i] = *n * 3 + ptvois[i + 1] - ptvois[i];
805 /* L10: */
806 }
807 if ((((int)*option) == 1)||(((int)*option) == -1)) {
808 debut = 0;
809 fin = *p;
810 step = 1;
811 } else {
812 debut = *p;
813 fin = 0;
814 step = -1;
815 }
816 i__1 = fin;
817 i__2 = step;
818 for (i = debut; i__2 < 0 ? i >= i__1 : i <= i__1; i += i__2) {
819 i1 = nv[i] + 1;
820 i2 = nv[i + 1];
821 i__3 = i2 - i1 + 1;
822 gibbs2_(&i__3, &nv[i1], &m[1]);
823 i__3 = i2;
824 for (j = i1; j <= i__3; ++j) {
825 s = nv[j];
826 i__4 = ptvois[s + 1] - 1;
827 for (k = ptvois[s]; k <= i__4; ++k) {
828 /* Computing MIN */
829 i__5 = m[vois[k]];
830 m[vois[k]] = mmin(i__5,j);
831 /* L20: */
832 }
833 /* L30: */
834 }
835 /* L40: */
836 }
837 if (*option > 0) {
838 knew = *new_;
839 plus = 1;
840 } else {
841 knew = *new_ + nbsc + 1;
842 plus = -1;
843 }
844 *new_ += nbsc;
845 /* if(option.gt.0) then */
846 /* do 60 k = debut , fin , step */
847 /* do 60 j = nv(k+1),nv(k)+1,-1 */
848 /* knew = knew + plus */
849 /* r(nv(j)) = knew */
850 /* 60 continue */
851 /* else */
852 i__2 = fin;
853 i__1 = step;
854 for (k = debut; i__1 < 0 ? k >= i__2 : k <= i__2; k += i__1) {
855 i__3 = nv[k + 1];
856 for (j = nv[k] + 1; j <= i__3; ++j) {
857 knew += plus;
858 r[nv[j]] = knew;
859 /* L70: */
860 }
861 }
862 /* endif */
863 *pfnew = 0;
864 bnew = 0;
865 i__3 = *n;
866 for (i = 1; i <= i__3; ++i) {
867 k = r[i];
868 if (k > 0) {
869 i__1 = ptvois[i + 1] - 1;
870 for (j = ptvois[i]; j <= i__1; ++j) {
871 if (r[vois[j]] > 0) {
872 /* Computing MIN */
873 i__2 = k, i__4 = r[vois[j]];
874 k = mmin(i__2,i__4);
875 }
876 /* L100: */
877 }
878 *pfnew = *pfnew + r[i] - k + 1;
879 /* Computing MAX */
880 i__1 = bnew, i__2 = r[i] - k + 1;
881 bnew = mmax(i__1,i__2);
882 }
883 /* L110: */
884 }
885 /* if(impre.lt.0.or.impre.gt.2) then */
886 /* write(nfout,*) ' option =',option,', profile =',pfnew */
887 /* + ,', 1/2 bande =',bnew,', new=',new,', nbss composante=',nbsc
888 */
889 /* endif */
890 return 0;} /* gibbst_ */
891
892 /* function */
gibbsv(integer * ptvoi,integer * vois,integer * lvois,integer * w,integer * v)893 int Triangles::gibbsv (integer* ptvoi,
894 integer* vois,integer* lvois,integer* w,integer* v)
895 {
896 /* System generated locals */
897 integer i__2;
898
899 /* Local variables */
900 integer i, j, k, T, ss, iii, ptv, ptv1;
901 integer nbss = nbv;
902 /*--- Prepare les donees pour gibbsa en construisant ptvoi, vois, lvois -
903 ------------*/
904 /* in */
905 /* --- nbnt =3 pour des triangles 2D, */
906 /* nbt = nb de triangle */
907 /* nbss = nb de sommets */
908 /* nsea = numeros de 3 sommets de chaque triangle (me) */
909 /* out */
910 /* --- ptvoi, vois, lvois, err */
911 /* tableaux de travail w, v */
912 /*-----------------------------------------------------------------------
913 ----------*/
914 /* Parameter adjustments */
915 --v;
916 --w;
917 --vois;
918 --ptvoi;
919 long nt = nbt-NbOutT;
920 /* Function Body */
921 for (i = 1; i <= nbss; ++i) {
922 w[i] = -1;
923 ptvoi[i] = 0; }
924 ptvoi[nbss + 1] = 0;
925 for (i = 0; i < nt; ++i)
926 {
927 assert(triangles[i].link);
928 for (j = 0; j < 3; ++j)
929 {
930 ss = Number(triangles[i][j])+1;
931 ++ptvoi[ss + 1];
932 w[ss] = 0;
933 }
934 }
935
936 for (i = 1; i <= nbss; ++i)
937 ptvoi[i + 1] += ptvoi[i];
938
939 for (i = 0; i < nt; ++i)
940 if (triangles[i].link)
941 for (j = 0; j < 3; ++j)
942 {
943 ss = Number(triangles[i][j])+1;
944 ++ptvoi[ss];
945 v[ptvoi[ss]] = i;
946 }
947
948 ptv1 = 0;
949 iii = 1;
950 for (i = 1; i <= nbss; ++i) {
951 ptv = ptv1 + 1;
952 ptv1 = ptvoi[i];
953 ptvoi[i] = iii;
954 i__2 = ptv1;
955 for (j = ptv; j <= i__2; ++j) {
956 T = v[j];
957 for (k = 0; k < 3; ++k) {
958 ss = Number(triangles[T][k])+1; /* nsea[k + T * nsea_dim1]; */
959 if (w[ss] != i) {
960 w[ss] = i;
961 if (iii > *lvois) return 2 ;
962 /* print*,'pas assez de place memoire' */
963
964 vois[iii] = ss;
965 ++iii;}
966 }
967 }
968 }
969 ptvoi[nbss + 1] = iii;
970 *lvois = iii - 1;
971 return 0; /* OK */
972 return 0;} /* gibbsv_ */
973
gibbs()974 int Triangles::gibbs()
975 /* --------
976 renumber vertices by gibbs method; updates triangle and edge array
977 in: mesh
978 out: mesh
979 auxiliary arrays: ptvois,vois,r,m,nv,nx,ny,nn,w1,w2,f
980 all of size nv+1 except vois (10(nv+1)) and nv (2(nv+1))
981 err = -1 : memory alloc pb; err = -3: fatal erreur gibbs 2 : pb racine
982 */
983 {
984 long nv = nbv;
985 long nt = nbt-NbOutT;
986 long i, j, pfold, pfnew;
987 long* ptvois=NULL;
988 long* vois=NULL;
989 long* nn =NULL;
990 long* r =NULL;
991 long* m =NULL;
992 long* nnv =NULL;
993 long* nx =NULL;
994 long* ny =NULL;
995 long* w1 =NULL;
996 long* w2=NULL;
997 long nbvoisin = 10*nv;
998 long printint=0, iodev=6;
999 int err=0;
1000 ptvois = new long[nv+1]; //(long*)calloc((long)(nv + 1) , sizeof(long));
1001 nn = new long[3*nt]; //(long*)calloc(3 * nt ,sizeof(long));
1002 vois = new long[nbvoisin+10]; //(long*)calloc((long)(nbvoisin + 10) , sizeof(long));
1003 r = new long[nv+1]; //(long*)calloc((long)(nv + 1) , sizeof(long));
1004 if((!ptvois)||(!nn)||(!vois)||(!r)) return -1;
1005 err = gibbsv(ptvois,vois,&nbvoisin,r,nn) ;
1006 delete [] nn; // free(nn);
1007 if(err==0)
1008 {
1009 m = new long[nv+1];
1010 nn = new long[nv+1];
1011 nnv = new long[(nv+1)<<1];
1012 nx = new long[nv+1];
1013 ny = new long[nv+1];
1014 w1 = new long[nv+1];
1015 w2 = new long[nv+1];
1016 long lnv = nv;
1017 err = gibbsa_ (&lnv, ptvois, vois, r, m, nnv, nx, ny, nn, w1, w2, &pfold, &pfnew,
1018 &printint, &iodev);
1019 delete [] m;
1020 delete [] nnv;
1021 delete [] nn;
1022 delete [] nx;
1023 delete [] ny;
1024 delete [] w1;
1025 delete [] w2;
1026 }
1027
1028 delete [] vois;
1029 delete [] ptvois;
1030 /*
1031 if (err == 0 && (pfnew <= pfold))
1032 {
1033 A<bVertex> f(nv);
1034 for (i = 0; i < nv; ++i)
1035 { f[i].x = v[i].x;
1036 f[i].y = v[i].y;
1037 f[i].where = v[i].where;
1038 }
1039 for (i = 0; i < nv; ++i)
1040 { v[r[i] - 1].x = f[i].x;
1041 v[r[i] - 1].y = f[i].y;
1042 v[r[i] - 1].where = f[i].where;
1043 }
1044
1045 for (j = 0; j < nt; ++j) // updates triangle array
1046 for (i = 0; i < 3; i++)
1047 t[j].v[i] = &v[r[no(t[j].v[i])] - 1];
1048
1049 for (j = 0; j < ne; ++j) // updates edge array
1050 {
1051 e[j].in = &v[r[no(e[j].in)] - 1];
1052 e[j].out = &v[r[no(e[j].out)] - 1];
1053 }
1054 f.destroy();
1055 if (!NumThinGrid)
1056 { NumThinGrid= new int [nv];
1057 for (i=0;i<nv;i++) NumThinGrid[i]=i;// Same numbering
1058 }
1059 for (i=0;i<nv;i++) NumThinGrid[i]=r[NumThinGrid[i]]-1;
1060
1061 }
1062 */
1063 delete [] r;
1064 return err;
1065 }
1066
1067 /* message d'erreur: *err = 2; print*,'pas assez de place memoire' */
1068