source: git/src/network.c @ 815b25e

stereo-2025
Last change on this file since 815b25e was 66be513, checked in by Olly Betts <olly@…>, 8 months ago

Fix to find more delta-star opportunities

This allows more complicated networks to be reduced further so
fewer and/or smaller matrices need solving, which is faster.

  • Property mode set to 100644
File size: 21.0 KB
Line 
1/* network.c
2 * Survex network reduction - find patterns and apply network reductions
3 * Copyright (C) 1991-2002,2005 Olly Betts
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
18 */
19
20#if 0
21#define DEBUG_INVALID 1
22#define VALIDATE 1
23#define DUMP_NETWORK 1
24#endif
25
26#include <config.h>
27
28#include "validate.h"
29#include "debug.h"
30#include "cavern.h"
31#include "message.h"
32#include "netbits.h"
33#include "network.h"
34#include "out.h"
35
36/* type field isn't vital - join3 is unused except for deltastar, so
37 * we can set its value to indicate which type this is:
38 * join3 == NULL for noose, join3 == join1 for ||, otherwise D* */
39#ifdef EXPLICIT_STACKRED_TYPE
40#define SET_NOOSE(SR) (SR)->type = 1
41#define IS_NOOSE(SR) ((SR)->type == 1)
42#define SET_PARALLEL(SR) (SR)->type = 0
43#define IS_PARALLEL(SR) ((SR)->type == 0)
44#define SET_DELTASTAR(SR) (SR)->type = 2
45#define IS_DELTASTAR(SR) ((SR)->type == 2)
46#else
47#define IS_NOOSE(SR) ((SR)->join3 == NULL)
48#define SET_NOOSE(SR) (SR)->join3 = NULL
49#define IS_PARALLEL(SR) ((SR)->join3 == (SR)->join1)
50#define SET_PARALLEL(SR) (SR)->join3 = (SR)->join1
51#define IS_DELTASTAR(SR) (!IS_NOOSE(SR) && !IS_PARALLEL(SR))
52#define SET_DELTASTAR(SR) NOP
53#endif
54
55typedef struct StackRed {
56   struct Link *join1, *join2, *join3;
57#ifdef EXPLICIT_STACKRED_TYPE
58   int type; /* 1 => noose, 0 => parallel legs, 2 => delta-star */
59#endif
60   struct StackRed *next;
61} stackRed;
62
63static stackRed *ptrRed; /* Ptr to TRaverse linked list for C*-*< , -*=*- */
64
65/* can be altered by -z<letters> on command line */
66unsigned long optimize = BITA('l') | BITA('p') | BITA('d');
67/* Lollipops, Parallel legs, Iterate mx, Delta* */
68
69extern void
70remove_subnets(void)
71{
72   node *stn, *stn2, *stn3, *stn4;
73   int dirn, dirn2, dirn3, dirn4;
74   stackRed *trav;
75   linkfor *newleg, *newleg2;
76   bool fMore = true;
77
78   ptrRed = NULL;
79
80   out_current_action(msg(/*Simplifying network*/129));
81
82   while (fMore) {
83      fMore = false;
84      if (optimize & BITA('l')) {
85#if PRINT_NETBITS
86         printf("replacing lollipops\n");
87#endif
88         /*        _
89          *       ( )
90          *        * stn
91          *        |
92          *        * stn2
93          *       / \
94          * stn4 *   * stn3  -->  stn4 *---* stn3
95          *      :   :                 :   :
96          */
97         /* NB can have non-fixed 0 nodes */
98         FOR_EACH_STN(stn, stnlist) {
99            if (three_node(stn)) {
100               dirn = -1;
101               if (stn->leg[1]->l.to == stn) dirn++;
102               if (stn->leg[0]->l.to == stn) dirn += 2;
103               if (dirn < 0) continue;
104
105               stn2 = stn->leg[dirn]->l.to;
106               if (fixed(stn2)) continue;
107
108               SVX_ASSERT(three_node(stn2));
109
110               dirn2 = reverse_leg_dirn(stn->leg[dirn]);
111               dirn2 = (dirn2 + 1) % 3;
112               stn3 = stn2->leg[dirn2]->l.to;
113               if (stn2 == stn3) continue; /* dumb-bell - leave alone */
114
115               dirn3 = reverse_leg_dirn(stn2->leg[dirn2]);
116
117               trav = osnew(stackRed);
118               newleg2 = (linkfor*)osnew(linkrev);
119
120               newleg = copy_link(stn3->leg[dirn3]);
121
122               dirn2 = (dirn2 + 1) % 3;
123               stn4 = stn2->leg[dirn2]->l.to;
124               dirn4 = reverse_leg_dirn(stn2->leg[dirn2]);
125#if 0
126               printf("Noose found with stn...stn4 = \n");
127               print_prefix(stn->name); putnl();
128               print_prefix(stn2->name); putnl();
129               print_prefix(stn3->name); putnl();
130               print_prefix(stn4->name); putnl();
131#endif
132
133               addto_link(newleg, stn2->leg[dirn2]);
134
135               /* remove stn and stn2 */
136               remove_stn_from_list(&stnlist, stn);
137               remove_stn_from_list(&stnlist, stn2);
138
139               /* stack noose and replace with a leg between stn3 and stn4 */
140               trav->join1 = stn3->leg[dirn3];
141               newleg->l.to = stn4;
142               newleg->l.reverse = dirn4 | FLAG_DATAHERE | FLAG_REPLACEMENTLEG;
143
144               trav->join2 = stn4->leg[dirn4];
145               newleg2->l.to = stn3;
146               newleg2->l.reverse = dirn3 | FLAG_REPLACEMENTLEG;
147
148               stn3->leg[dirn3] = newleg;
149               stn4->leg[dirn4] = newleg2;
150
151               trav->next = ptrRed;
152               SET_NOOSE(trav);
153#if PRINT_NETBITS
154               printf("remove noose\n");
155#endif
156               ptrRed = trav;
157               fMore = true;
158            }
159         }
160      }
161
162      if (optimize & BITA('p')) {
163#if PRINT_NETBITS
164         printf("replacing parallel legs\n");
165#endif
166         FOR_EACH_STN(stn, stnlist) {
167            /*
168             *  :            :
169             *  * stn3       * stn3
170             *  |            |
171             *  * stn        |
172             * ( )      -->  |
173             *  * stn2       |
174             *  |            |
175             *  * stn4       * stn4
176             *  :            :
177             */
178            if (three_node(stn)) {
179               stn2 = stn->leg[0]->l.to;
180               if (stn2 == stn->leg[1]->l.to) {
181                  dirn = 2;
182               } else if (stn2 == stn->leg[2]->l.to) {
183                  dirn = 1;
184               } else {
185                  if (stn->leg[1]->l.to != stn->leg[2]->l.to) continue;
186                  stn2 = stn->leg[1]->l.to;
187                  dirn = 0;
188               }
189
190               /* stn == stn2 => noose */
191               if (stn == stn2 || fixed(stn2)) continue;
192
193               SVX_ASSERT(three_node(stn2));
194
195               stn3 = stn->leg[dirn]->l.to;
196               /* 3 parallel legs (=> nothing else) so leave */
197               if (stn3 == stn2) continue;
198
199               dirn3 = reverse_leg_dirn(stn->leg[dirn]);
200               dirn2 = (0 + 1 + 2 - reverse_leg_dirn(stn->leg[(dirn + 1) % 3])
201                        - reverse_leg_dirn(stn->leg[(dirn + 2) % 3]));
202
203               stn4 = stn2->leg[dirn2]->l.to;
204               dirn4 = reverse_leg_dirn(stn2->leg[dirn2]);
205
206               trav = osnew(stackRed);
207
208               newleg = copy_link(stn->leg[(dirn + 1) % 3]);
209               /* use newleg2 for scratch */
210               newleg2 = copy_link(stn->leg[(dirn + 2) % 3]);
211                 {
212#ifdef NO_COVARIANCES
213                    vars sum;
214                    var prod;
215                    delta temp, temp2;
216                    addss(&sum, &newleg->v, &newleg2->v);
217                    SVX_ASSERT2(!fZeros(&sum), "loop of zero variance found");
218                    mulss(&prod, &newleg->v, &newleg2->v);
219                    mulsd(&temp, &newleg2->v, &newleg->d);
220                    mulsd(&temp2, &newleg->v, &newleg2->d);
221                    adddd(&temp, &temp, &temp2);
222                    divds(&newleg->d, &temp, &sum);
223                    sdivvs(&newleg->v, &prod, &sum);
224#else
225                    svar inv1, inv2, sum;
226                    delta temp, temp2;
227                    /* if leg one is an equate, we can just ignore leg two
228                     * whatever it is */
229                    if (invert_svar(&inv1, &newleg->v)) {
230                       if (invert_svar(&inv2, &newleg2->v)) {
231                          addss(&sum, &inv1, &inv2);
232                          if (!invert_svar(&newleg->v, &sum)) {
233                             BUG("matrix singular in parallel legs replacement");
234                          }
235
236                          mulsd(&temp, &inv1, &newleg->d);
237                          mulsd(&temp2, &inv2, &newleg2->d);
238                          adddd(&temp, &temp, &temp2);
239                          mulsd(&newleg->d, &newleg->v, &temp);
240                       } else {
241                          /* leg two is an equate, so just ignore leg 1 */
242                          linkfor *tmpleg;
243                          tmpleg = newleg;
244                          newleg = newleg2;
245                          newleg2 = tmpleg;
246                       }
247                    }
248#endif
249                 }
250               osfree(newleg2);
251               newleg2 = (linkfor*)osnew(linkrev);
252
253               addto_link(newleg, stn2->leg[dirn2]);
254               addto_link(newleg, stn3->leg[dirn3]);
255
256#if 0
257               printf("Parallel found with stn...stn4 = \n");
258               (dump_node)(stn); (dump_node)(stn2); (dump_node)(stn3); (dump_node)(stn4);
259               printf("dirns = %d %d %d %d\n", dirn, dirn2, dirn3, dirn4);
260#endif
261               SVX_ASSERT2(stn3->leg[dirn3]->l.to == stn, "stn3 end of || doesn't recip");
262               SVX_ASSERT2(stn4->leg[dirn4]->l.to == stn2, "stn4 end of || doesn't recip");
263               SVX_ASSERT2(stn->leg[(dirn+1)%3]->l.to == stn2 && stn->leg[(dirn + 2) % 3]->l.to == stn2, "|| legs aren't");
264
265               /* remove stn and stn2 (already discarded triple parallel) */
266               /* so stn!=stn4 <=> stn2!=stn3 */
267               remove_stn_from_list(&stnlist, stn);
268               remove_stn_from_list(&stnlist, stn2);
269
270               /* stack parallel and replace with a leg between stn3 and stn4 */
271               trav->join1 = stn3->leg[dirn3];
272               newleg->l.to = stn4;
273               newleg->l.reverse = dirn4 | FLAG_DATAHERE | FLAG_REPLACEMENTLEG;
274
275               trav->join2 = stn4->leg[dirn4];
276               newleg2->l.to = stn3;
277               newleg2->l.reverse = dirn3 | FLAG_REPLACEMENTLEG;
278
279               stn3->leg[dirn3] = newleg;
280               stn4->leg[dirn4] = newleg2;
281
282               trav->next = ptrRed;
283               SET_PARALLEL(trav);
284#if PRINT_NETBITS
285               printf("remove parallel\n");
286#endif
287               ptrRed = trav;
288               fMore = true;
289            }
290         }
291      }
292
293      if (optimize & BITA('d')) {
294         node *stn5, *stn6;
295         int dirn5, dirn6;
296         linkfor *legAB, *legBC, *legCA;
297#if PRINT_NETBITS
298         printf("replacing deltas with stars\n");
299#endif
300         FOR_EACH_STN(stn, stnlist) {
301            /*    printf("*");*/
302            /*
303             *          :                     :
304             *          * stn5                * stn5
305             *          |                     |
306             *          * stn2                |
307             *         / \        -->         O stnZ
308             *        |   |                  / \
309             *    stn *---* stn3            /   \
310             *       /     \               /     \
311             * stn4 *       * stn6   stn4 *       * stn6
312             *      :       :             :       :
313             */
314            if (three_node(stn)) {
315               for (int dirn12 = 0; dirn12 <= 2; dirn12++) {
316                  stn2 = stn->leg[dirn12]->l.to;
317                  if (stn2 == stn || fixed(stn2)) continue;
318                  SVX_ASSERT(three_node(stn2));
319                  int dirn13 = (dirn12 + 1) % 3;
320                  stn3 = stn->leg[dirn13]->l.to;
321                  if (stn3 == stn || stn3 == stn2 || fixed(stn3)) continue;
322                  SVX_ASSERT(three_node(stn3));
323                  int dirn23 = reverse_leg_dirn(stn->leg[dirn12]);
324                  dirn23 = (dirn23 + 1) % 3;
325                  if (stn2->leg[dirn23]->l.to != stn3) {
326                      dirn23 = (dirn23 + 1) % 3;
327                      if (stn2->leg[dirn23]->l.to != stn3) {
328                          continue;
329                      }
330                  }
331                  legAB = copy_link(stn->leg[dirn12]);
332                  legBC = copy_link(stn2->leg[dirn23]);
333                  legCA = copy_link(stn3->leg[reverse_leg_dirn(stn->leg[dirn13])]);
334                  dirn = (0 + 1 + 2) - dirn12 - dirn13;
335                  dirn2 = (0 + 1 + 2) - dirn23 - reverse_leg_dirn(stn->leg[dirn12]);
336                  dirn3 = (0 + 1 + 2) - reverse_leg_dirn(stn->leg[dirn13]) - reverse_leg_dirn(stn2->leg[dirn23]);
337                  stn4 = stn->leg[dirn]->l.to;
338                  stn5 = stn2->leg[dirn2]->l.to;
339                  stn6 = stn3->leg[dirn3]->l.to;
340                  if (stn4 == stn2 || stn4 == stn3 || stn5 == stn3) continue;
341                  dirn4 = reverse_leg_dirn(stn->leg[dirn]);
342                  dirn5 = reverse_leg_dirn(stn2->leg[dirn2]);
343                  dirn6 = reverse_leg_dirn(stn3->leg[dirn3]);
344#if 0
345                  printf("delta-star, stn ... stn6 are:\n");
346                  (dump_node)(stn);
347                  (dump_node)(stn2);
348                  (dump_node)(stn3);
349                  (dump_node)(stn4);
350                  (dump_node)(stn5);
351                  (dump_node)(stn6);
352#endif
353                  SVX_ASSERT(stn4->leg[dirn4]->l.to == stn);
354                  SVX_ASSERT(stn5->leg[dirn5]->l.to == stn2);
355                  SVX_ASSERT(stn6->leg[dirn6]->l.to == stn3);
356
357                  trav = osnew(stackRed);
358                  {
359                    linkfor *legAZ, *legBZ, *legCZ;
360                    node *stnZ;
361                    prefix *nameZ;
362                    svar invAB, invBC, invCA, tmp, sum, inv;
363                    var vtmp;
364                    svar sumAZBZ, sumBZCZ, sumCZAZ;
365                    delta temp, temp2;
366
367                    /* FIXME: ought to handle cases when some legs are
368                     * equates, but handle as a special case maybe? */
369                    if (!invert_svar(&invAB, &legAB->v)) break;
370                    if (!invert_svar(&invBC, &legBC->v)) break;
371                    if (!invert_svar(&invCA, &legCA->v)) break;
372
373                    addss(&sum, &legBC->v, &legCA->v);
374                    addss(&tmp, &sum, &legAB->v);
375                    if (!invert_svar(&inv, &tmp)) {
376                       /* impossible - loop of zero variance */
377                       BUG("loop of zero variance found");
378                    }
379
380                    legAZ = osnew(linkfor);
381                    legBZ = osnew(linkfor);
382                    legCZ = osnew(linkfor);
383
384                    /* AZBZ */
385                    /* done above: addvv(&sum, &legBC->v, &legCA->v); */
386                    mulss(&vtmp, &sum, &inv);
387                    smulvs(&sumAZBZ, &vtmp, &legAB->v);
388
389                    adddd(&temp, &legBC->d, &legCA->d);
390                    divds(&temp2, &temp, &sum);
391                    mulsd(&temp, &invAB, &legAB->d);
392                    subdd(&temp, &temp2, &temp);
393                    mulsd(&legBZ->d, &sumAZBZ, &temp);
394
395                    /* leg vectors after transform are determined up to
396                     * a constant addition, so arbitrarily fix AZ = 0 */
397                    legAZ->d[2] = legAZ->d[1] = legAZ->d[0] = 0;
398
399                    /* BZCZ */
400                    addss(&sum, &legCA->v, &legAB->v);
401                    mulss(&vtmp, &sum, &inv);
402                    smulvs(&sumBZCZ, &vtmp, &legBC->v);
403
404                    /* CZAZ */
405                    addss(&sum, &legAB->v, &legBC->v);
406                    mulss(&vtmp, &sum, &inv);
407                    smulvs(&sumCZAZ, &vtmp, &legCA->v);
408
409                    adddd(&temp, &legAB->d, &legBC->d);
410                    divds(&temp2, &temp, &sum);
411                    mulsd(&temp, &invCA, &legCA->d);
412                    /* NB: swapped arguments to negate answer for legCZ->d */
413                    subdd(&temp, &temp, &temp2);
414                    mulsd(&legCZ->d, &sumCZAZ, &temp);
415
416                    osfree(legAB);
417                    osfree(legBC);
418                    osfree(legCA);
419
420                    /* Now add two, subtract third, and scale by 0.5 */
421                    addss(&sum, &sumAZBZ, &sumCZAZ);
422                    subss(&sum, &sum, &sumBZCZ);
423                    mulsc(&legAZ->v, &sum, 0.5);
424
425                    addss(&sum, &sumBZCZ, &sumAZBZ);
426                    subss(&sum, &sum, &sumCZAZ);
427                    mulsc(&legBZ->v, &sum, 0.5);
428
429                    addss(&sum, &sumCZAZ, &sumBZCZ);
430                    subss(&sum, &sum, &sumAZBZ);
431                    mulsc(&legCZ->v, &sum, 0.5);
432
433                    nameZ = osnew(prefix);
434                    nameZ->pos = osnew(pos);
435                    nameZ->ident.p = NULL;
436                    nameZ->shape = 3;
437                    stnZ = osnew(node);
438                    stnZ->name = nameZ;
439                    nameZ->stn = stnZ;
440                    nameZ->up = NULL;
441                    nameZ->min_export = nameZ->max_export = 0;
442                    unfix(stnZ);
443                    add_stn_to_list(&stnlist, stnZ);
444                    legAZ->l.to = stnZ;
445                    legAZ->l.reverse = 0 | FLAG_DATAHERE | FLAG_REPLACEMENTLEG;
446                    legBZ->l.to = stnZ;
447                    legBZ->l.reverse = 1 | FLAG_DATAHERE | FLAG_REPLACEMENTLEG;
448                    legCZ->l.to = stnZ;
449                    legCZ->l.reverse = 2 | FLAG_DATAHERE | FLAG_REPLACEMENTLEG;
450                    stnZ->leg[0] = (linkfor*)osnew(linkrev);
451                    stnZ->leg[1] = (linkfor*)osnew(linkrev);
452                    stnZ->leg[2] = (linkfor*)osnew(linkrev);
453                    stnZ->leg[0]->l.to = stn4;
454                    stnZ->leg[0]->l.reverse = dirn4;
455                    stnZ->leg[1]->l.to = stn5;
456                    stnZ->leg[1]->l.reverse = dirn5;
457                    stnZ->leg[2]->l.to = stn6;
458                    stnZ->leg[2]->l.reverse = dirn6;
459                    addto_link(legAZ, stn4->leg[dirn4]);
460                    addto_link(legBZ, stn5->leg[dirn5]);
461                    addto_link(legCZ, stn6->leg[dirn6]);
462                    /* stack stuff */
463                    trav->join1 = stn4->leg[dirn4];
464                    trav->join2 = stn5->leg[dirn5];
465                    trav->join3 = stn6->leg[dirn6];
466                    trav->next = ptrRed;
467                    SET_DELTASTAR(trav);
468#if PRINT_NETBITS
469                    printf("remove delta*\n");
470#endif
471                    ptrRed = trav;
472                    fMore = true;
473
474                    remove_stn_from_list(&stnlist, stn);
475                    remove_stn_from_list(&stnlist, stn2);
476                    remove_stn_from_list(&stnlist, stn3);
477                    stn4->leg[dirn4] = legAZ;
478                    stn5->leg[dirn5] = legBZ;
479                    stn6->leg[dirn6] = legCZ;
480                  }
481                  break;
482               }
483            }
484         }
485      }
486
487   }
488}
489
490extern void
491replace_subnets(void)
492{
493   stackRed *ptrOld;
494   node *stn2, *stn3, *stn4;
495   int dirn2, dirn3, dirn4;
496
497   /* help to catch bad accesses */
498   stn2 = stn3 = stn4 = NULL;
499   dirn2 = dirn3 = dirn4 = 0;
500
501   out_current_action(msg(/*Calculating network*/130));
502
503   while (ptrRed != NULL) {
504      /*  printf("replace_subnets() type %d\n", ptrRed->type);*/
505
506#if PRINT_NETBITS
507      printf("replace_subnets\n");
508      if (IS_NOOSE(ptrRed)) printf("isnoose\n");
509      if (IS_PARALLEL(ptrRed)) printf("isparallel\n");
510      if (IS_DELTASTAR(ptrRed)) printf("isdelta*\n");
511#endif
512
513      if (!IS_DELTASTAR(ptrRed)) {
514         linkfor *leg;
515         leg = ptrRed->join1; leg = reverse_leg(leg);
516         stn3 = leg->l.to; dirn3 = reverse_leg_dirn(leg);
517         leg = ptrRed->join2; leg = reverse_leg(leg);
518         stn4 = leg->l.to; dirn4 = reverse_leg_dirn(leg);
519
520         if (!fixed(stn3) || !fixed(stn4)) {
521             SVX_ASSERT(!fixed(stn3) && !fixed(stn4));
522             goto skip;
523         }
524         SVX_ASSERT(data_here(stn3->leg[dirn3]));
525      }
526
527      if (IS_NOOSE(ptrRed)) {
528         /* noose (hanging-loop) */
529         node *stn;
530         delta e;
531         linkfor *leg;
532         int zero;
533
534         leg = stn3->leg[dirn3];
535         stn2 = ptrRed->join1->l.to;
536         dirn2 = reverse_leg_dirn(ptrRed->join1);
537
538         zero = fZeros(&leg->v);
539         if (!zero) {
540            delta tmp;
541            subdd(&e, &POSD(stn4), &POSD(stn3));
542            subdd(&tmp, &e, &leg->d);
543            divds(&e, &tmp, &leg->v);
544         }
545         if (data_here(ptrRed->join1)) {
546            adddd(&POSD(stn2), &POSD(stn3), &ptrRed->join1->d);
547            if (!zero) {
548               delta tmp;
549               mulsd(&tmp, &ptrRed->join1->v, &e);
550               adddd(&POSD(stn2), &POSD(stn2), &tmp);
551            }
552         } else {
553            subdd(&POSD(stn2), &POSD(stn3), &stn2->leg[dirn2]->d);
554            if (!zero) {
555               delta tmp;
556               mulsd(&tmp, &stn2->leg[dirn2]->v, &e);
557               adddd(&POSD(stn2), &POSD(stn2), &tmp);
558            }
559         }
560         dirn2 = (dirn2 + 2) % 3; /* point back at stn again */
561         stn = stn2->leg[dirn2]->l.to;
562#if 0
563         printf("Replacing noose with stn...stn4 = \n");
564         print_prefix(stn->name); putnl();
565         print_prefix(stn2->name); putnl();
566         print_prefix(stn3->name); putnl();
567         print_prefix(stn4->name); putnl();
568#endif
569         if (data_here(stn2->leg[dirn2]))
570            adddd(&POSD(stn), &POSD(stn2), &stn2->leg[dirn2]->d);
571         else
572            subdd(&POSD(stn), &POSD(stn2), &reverse_leg(stn2->leg[dirn2])->d);
573
574         /* the "rope" of the noose is a new articulation */
575         stn2->leg[dirn2]->l.reverse |= FLAG_ARTICULATION;
576         reverse_leg(stn2->leg[dirn2])->l.reverse |= FLAG_ARTICULATION;
577
578         add_stn_to_list(&fixedlist, stn);
579         add_stn_to_list(&fixedlist, stn2);
580
581         osfree(stn3->leg[dirn3]);
582         stn3->leg[dirn3] = ptrRed->join1;
583         osfree(stn4->leg[dirn4]);
584         stn4->leg[dirn4] = ptrRed->join2;
585      } else if (IS_PARALLEL(ptrRed)) {
586         /* parallel legs */
587         node *stn;
588         delta e, e2;
589         linkfor *leg;
590         int dirn;
591
592         stn = ptrRed->join1->l.to;
593         stn2 = ptrRed->join2->l.to;
594
595         dirn = reverse_leg_dirn(ptrRed->join1);
596         dirn2 = reverse_leg_dirn(ptrRed->join2);
597
598         leg = stn3->leg[dirn3];
599
600         if (leg->l.reverse & FLAG_ARTICULATION) {
601            ptrRed->join1->l.reverse |= FLAG_ARTICULATION;
602            stn->leg[dirn]->l.reverse |= FLAG_ARTICULATION;
603            ptrRed->join2->l.reverse |= FLAG_ARTICULATION;
604            stn2->leg[dirn2]->l.reverse |= FLAG_ARTICULATION;
605         }
606
607         if (fZeros(&leg->v))
608            e[0] = e[1] = e[2] = 0.0;
609         else {
610            delta tmp;
611            subdd(&e, &POSD(stn4), &POSD(stn3));
612            subdd(&tmp, &e, &leg->d);
613            divds(&e, &tmp, &leg->v);
614         }
615
616         if (data_here(ptrRed->join1)) {
617            leg = ptrRed->join1;
618            adddd(&POSD(stn), &POSD(stn3), &leg->d);
619         } else {
620            leg = stn->leg[dirn];
621            subdd(&POSD(stn), &POSD(stn3), &leg->d);
622         }
623         mulsd(&e2, &leg->v, &e);
624         adddd(&POSD(stn), &POSD(stn), &e2);
625
626         if (data_here(ptrRed->join2)) {
627            leg = ptrRed->join2;
628            adddd(&POSD(stn2), &POSD(stn4), &leg->d);
629         } else {
630            leg = stn2->leg[dirn2];
631            subdd(&POSD(stn2), &POSD(stn4), &leg->d);
632         }
633         mulsd(&e2, &leg->v, &e);
634         subdd(&POSD(stn2), &POSD(stn2), &e2);
635#if 0
636         printf("Replacing parallel with stn...stn4 = \n");
637         print_prefix(stn->name); putnl();
638         print_prefix(stn2->name); putnl();
639         print_prefix(stn3->name); putnl();
640         print_prefix(stn4->name); putnl();
641#endif
642
643         add_stn_to_list(&fixedlist, stn);
644         add_stn_to_list(&fixedlist, stn2);
645
646         osfree(stn3->leg[dirn3]);
647         stn3->leg[dirn3] = ptrRed->join1;
648         osfree(stn4->leg[dirn4]);
649         stn4->leg[dirn4] = ptrRed->join2;
650      } else if (IS_DELTASTAR(ptrRed)) {
651         node *stnZ;
652         node *stn[3];
653         int dirn[3];
654         linkfor *legs[3];
655         int i;
656         linkfor *leg;
657
658         legs[0] = ptrRed->join1;
659         legs[1] = ptrRed->join2;
660         legs[2] = ptrRed->join3;
661
662         /* work out ends as we don't bother stacking them */
663         leg = reverse_leg(legs[0]);
664         stn[0] = leg->l.to;
665         dirn[0] = reverse_leg_dirn(leg);
666         stnZ = stn[0]->leg[dirn[0]]->l.to;
667         SVX_ASSERT(fixed(stnZ));
668         if (!fixed(stnZ)) {
669            SVX_ASSERT(!fixed(stn[0]));
670            goto skip;
671         }
672         stn[1] = stnZ->leg[1]->l.to;
673         dirn[1] = reverse_leg_dirn(stnZ->leg[1]);
674         stn[2] = stnZ->leg[2]->l.to;
675         dirn[2] = reverse_leg_dirn(stnZ->leg[2]);
676         /*print_prefix(stnZ->name);printf(" %p\n",(void*)stnZ);*/
677
678         for (i = 0; i < 3; i++) {
679            SVX_ASSERT2(fixed(stn[i]), "stn not fixed for D*");
680
681            leg = stn[i]->leg[dirn[i]];
682
683            SVX_ASSERT2(data_here(leg), "data not on leg for D*");
684            SVX_ASSERT2(leg->l.to == stnZ, "bad sub-network for D*");
685
686            stn2 = legs[i]->l.to;
687
688            if (data_here(legs[i])) {
689               adddd(&POSD(stn2), &POSD(stn[i]), &legs[i]->d);
690            } else {
691               subdd(&POSD(stn2), &POSD(stn[i]), &reverse_leg(legs[i])->d);
692            }
693
694            if (!fZeros(&leg->v)) {
695               delta e, tmp;
696               subdd(&e, &POSD(stnZ), &POSD(stn[i]));
697               subdd(&e, &e, &leg->d);
698               divds(&tmp, &e, &leg->v);
699               if (data_here(legs[i])) {
700                  mulsd(&e, &legs[i]->v, &tmp);
701               } else {
702                  mulsd(&e, &reverse_leg(legs[i])->v, &tmp);
703               }
704               adddd(&POSD(stn2), &POSD(stn2), &e);
705            }
706            add_stn_to_list(&fixedlist, stn2);
707            osfree(leg);
708            stn[i]->leg[dirn[i]] = legs[i];
709            /* transfer the articulation status of the radial legs */
710            if (stnZ->leg[i]->l.reverse & FLAG_ARTICULATION) {
711               legs[i]->l.reverse |= FLAG_ARTICULATION;
712               reverse_leg(legs[i])->l.reverse |= FLAG_ARTICULATION;
713            }
714            osfree(stnZ->leg[i]);
715            stnZ->leg[i] = NULL;
716         }
717/*printf("---%f %f %f\n",POS(stnZ, 0), POS(stnZ, 1), POS(stnZ, 2));*/
718         remove_stn_from_list(&fixedlist, stnZ);
719         osfree(stnZ->name);
720         osfree(stnZ);
721      } else {
722         BUG("ptrRed has unknown type");
723      }
724
725skip:
726      ptrOld = ptrRed;
727      ptrRed = ptrRed->next;
728      osfree(ptrOld);
729   }
730}
Note: See TracBrowser for help on using the repository browser.