[BACK]Return to standard_root_refiners.adb CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / PHC / Ada / Continuation

Annotation of OpenXM_contrib/PHC/Ada/Continuation/standard_root_refiners.adb, Revision 1.1

1.1     ! maekawa     1: with integer_io;                         use integer_io;
        !             2: with Standard_Floating_Numbers_io;       use Standard_Floating_Numbers_io;
        !             3: with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
        !             4: with Standard_Complex_Numbers_io;        use Standard_Complex_Numbers_io;
        !             5: with Standard_Complex_Vectors;           use Standard_Complex_Vectors;
        !             6: with Standard_Natural_Vectors;
        !             7: with Standard_Complex_Matrices;          use Standard_Complex_Matrices;
        !             8: with Standard_Complex_Linear_Solvers;    use Standard_Complex_Linear_Solvers;
        !             9: with Standard_Complex_Norms_Equals;      use Standard_Complex_Norms_Equals;
        !            10: with Standard_Complex_Solutions_io;      use Standard_Complex_Solutions_io;
        !            11:
        !            12: package body Standard_Root_Refiners is
        !            13:
        !            14: -- AUXILIARIES :
        !            15:
        !            16:   function Is_Real ( sol : Solution; tol : double_float ) return boolean is
        !            17:   begin
        !            18:     for i in sol.v'range loop
        !            19:       if ABS(IMAG_PART(sol.v(i))) > tol
        !            20:        then return false;
        !            21:       end if;
        !            22:     end loop;
        !            23:     return true;
        !            24:   end Is_Real;
        !            25:
        !            26:   function Equal ( s1,s2 : Solution; tol : double_float ) return boolean is
        !            27:   begin
        !            28:     for i in s1.v'range loop
        !            29:       if not Equal(s1.v(i),s2.v(i),tol)
        !            30:        then return false;
        !            31:       end if;
        !            32:     end loop;
        !            33:     return true;
        !            34:   end Equal;
        !            35:
        !            36:   function Complex_Conjugate ( s : Solution ) return Solution is
        !            37:
        !            38:     res : Solution(s.n) := s;
        !            39:
        !            40:   begin
        !            41:     for i in res.v'range loop
        !            42:       res.v(i) := Create(REAL_PART(s.v(i)),-IMAG_PART(s.v(i)));
        !            43:     end loop;
        !            44:     return res;
        !            45:   end Complex_Conjugate;
        !            46:
        !            47:   function Is_Clustered ( sol : Solution; nb : natural; sols : Solution_List;
        !            48:                           tol : double_float ) return natural is
        !            49:
        !            50:     tmp : Solution_List := sols;
        !            51:     cnt : natural := 0;
        !            52:
        !            53:   begin
        !            54:     while not Is_Null(tmp) loop
        !            55:       cnt := cnt + 1;
        !            56:       if cnt /= nb
        !            57:        then if Equal(sol,Head_Of(tmp).all,tol)
        !            58:             then return cnt;
        !            59:             end if;
        !            60:       end if;
        !            61:       tmp := Tail_Of(tmp);
        !            62:     end loop;
        !            63:     return nb;
        !            64:   end Is_Clustered;
        !            65:
        !            66:   function Is_Clustered ( sol : Solution; nb : natural; sols : Solution_Array;
        !            67:                           tol : double_float ) return natural is
        !            68:   begin
        !            69:     for i in sols'range loop
        !            70:       if i /= nb
        !            71:        then if Equal(sol,sols(i).all,tol)
        !            72:              then return i;
        !            73:             end if;
        !            74:       end if;
        !            75:     end loop;
        !            76:     return nb;
        !            77:   end Is_Clustered;
        !            78:
        !            79:   function Multiplicity ( sol : Solution; sols : Solution_List;
        !            80:                           tol : double_float ) return natural is
        !            81:
        !            82:     tmp : Solution_List := sols;
        !            83:     cnt : natural := 0;
        !            84:
        !            85:   begin
        !            86:     while not Is_Null(tmp) loop
        !            87:       if Equal(sol,Head_Of(tmp).all,tol)
        !            88:        then cnt := cnt + 1;
        !            89:       end if;
        !            90:       tmp := Tail_Of(tmp);
        !            91:     end loop;
        !            92:     return cnt;
        !            93:   end Multiplicity;
        !            94:
        !            95:   function Multiplicity ( sol : Solution; sols : Solution_Array;
        !            96:                           tol : double_float ) return natural is
        !            97:     cnt : natural := 0;
        !            98:
        !            99:   begin
        !           100:     for i in sols'range loop
        !           101:       if Equal(sol,sols(i).all,tol)
        !           102:        then cnt := cnt + 1;
        !           103:       end if;
        !           104:     end loop;
        !           105:     return cnt;
        !           106:   end Multiplicity;
        !           107:
        !           108:   procedure Write_Bar ( file : in file_type ) is
        !           109:   begin
        !           110:     for i in 1..75 loop
        !           111:       put(file,'=');
        !           112:     end loop;
        !           113:     new_line(file);
        !           114:   end Write_Bar;
        !           115:
        !           116:   procedure Write_Info ( file : in file_type; zero : in Solution;
        !           117:                          i,n,numb : in natural; fail : in boolean ) is
        !           118:
        !           119:   -- DESCRIPTION :
        !           120:   --   The information concerning the zero is written
        !           121:
        !           122:   begin
        !           123:    -- Write_Bar(file);
        !           124:     put(file,"solution "); put(file,i,1); put(file," :        ");
        !           125:     put(file," start residual : "); put(file,zero.res,2,3,3); new_line(file);
        !           126:     put(file,zero);
        !           127:   end Write_Info;
        !           128:
        !           129:   procedure Root_Accounting
        !           130:                ( file : in file_type; ls : in out Link_to_Solution;
        !           131:                  nb : in natural; sa : in out Solution_Array;
        !           132:                  fail : in boolean; tolsing,tolclus : in double_float;
        !           133:                  nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus : in out natural ) is
        !           134:
        !           135:   -- DESCRIPTION :
        !           136:   --   This procedure does root accounting of the solution sol, w.r.t. the
        !           137:   --   solution list sols.  Information will be provided concerning the type
        !           138:   --   of solution.
        !           139:
        !           140:   begin
        !           141:     if fail
        !           142:      then put_line(file," no solution ==");
        !           143:           nbfail := nbfail + 1;
        !           144:           ls.m := 0;
        !           145:      else if Is_Real(ls.all,10.0**(-13))
        !           146:           then put(file," real ");
        !           147:                nbreal := nbreal + 1;
        !           148:           else put(file," complex ");
        !           149:                nbcomp := nbcomp + 1;
        !           150:           end if;
        !           151:           if sa(nb).rco < tolsing
        !           152:           then declare
        !           153:                   m : natural := Multiplicity(ls.all,sa,tolclus);
        !           154:                 begin
        !           155:                   if m = 1
        !           156:                    then m := 0;
        !           157:                  end if;
        !           158:                  ls.m := m;
        !           159:                end;
        !           160:                put_line(file,"singular ==");
        !           161:                nbsing := nbsing + 1;
        !           162:            else declare
        !           163:                   nb2 : natural := Is_Clustered(ls.all,nb,sa,tolclus);
        !           164:                 begin
        !           165:                   if nb2 = nb
        !           166:                   then put_line(file,"regular ==");
        !           167:                        nbreg := nbreg + 1;
        !           168:                    else put(file,"clustered : ");
        !           169:                        put(file,nb2,1);
        !           170:                        put_line(file," ==");
        !           171:                        nbclus := nbclus + 1;
        !           172:                   end if;
        !           173:                  ls.m := 1;
        !           174:                end;
        !           175:           end if;
        !           176:     end if;
        !           177:   end Root_Accounting;
        !           178:
        !           179:   procedure Root_Accounting
        !           180:                 ( ls : in out Link_to_Solution; nb : in natural;
        !           181:                   sa : in out Solution_Array; fail : in boolean;
        !           182:                   tolsing,tolclus : in double_float ) is
        !           183:
        !           184:   -- DESCRIPTION :
        !           185:   --   This procedure does root accounting of the solution sol, w.r.t. the
        !           186:   --   solution list sols.  Information will be provided concerning the type
        !           187:   --   of solution.
        !           188:
        !           189:   begin
        !           190:     if fail
        !           191:      then ls.m := 0;
        !           192:      elsif sa(nb).rco < tolsing
        !           193:          then declare
        !           194:                 m : natural := Multiplicity(ls.all,sa,tolclus);
        !           195:               begin
        !           196:                 if m = 1
        !           197:                  then ls.m := 0;
        !           198:                  else ls.m := m;
        !           199:                 end if;
        !           200:               end;
        !           201:          else ls.m := 1;
        !           202:     end if;
        !           203:   end Root_Accounting;
        !           204:
        !           205:   procedure Write_Global_Info
        !           206:              ( file : in file_type;
        !           207:                tot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus : in natural ) is
        !           208:
        !           209:   begin
        !           210:     Write_Bar(file);
        !           211:     put(file,"A list of "); put(file,tot,1);
        !           212:     put_line(file," solutions has been refined :");
        !           213:     put(file,"Number of regular solutions   : "); put(file,nbreg,1);
        !           214:     put_line(file,".");
        !           215:     put(file,"Number of singular solutions  : "); put(file,nbsing,1);
        !           216:     put_line(file,".");
        !           217:     put(file,"Number of real solutions      : "); put(file,nbreal,1);
        !           218:     put_line(file,".");
        !           219:     put(file,"Number of complex solutions   : "); put(file,nbcomp,1);
        !           220:     put_line(file,".");
        !           221:     put(file,"Number of clustered solutions : "); put(file,nbclus,1);
        !           222:     put_line(file,".");
        !           223:     put(file,"Number of failures            : "); put(file,nbfail,1);
        !           224:     put_line(file,".");
        !           225:     Write_Bar(file);
        !           226:   end Write_Global_Info;
        !           227:
        !           228: -- TARGET ROUTINES :
        !           229:
        !           230:   procedure Silent_Newton
        !           231:                ( p_eval : in Eval_Poly_Sys; j_eval : in Eval_Jaco_Mat;
        !           232:                  zero : in out Solution; epsxa,epsfa : in double_float;
        !           233:                  numit : in out natural; max : in natural;
        !           234:                  fail : out boolean ) is
        !           235:
        !           236:     n : natural := zero.n;
        !           237:     jacobian : matrix(1..n,1..n);
        !           238:     ipvt : Standard_Natural_Vectors.Vector(1..n);
        !           239:     y,deltax : Vector(1..n);
        !           240:
        !           241:   begin
        !           242:     y := eval(p_eval,zero.v);               -- y = f(zero)
        !           243:     for i in 1..max loop
        !           244:       jacobian := eval(j_eval,zero.v);      -- solve jacobian*deltax = -f(zero)
        !           245:       lufco(jacobian,n,ipvt,zero.rco);
        !           246:       if 1.0 + zero.rco = 1.0
        !           247:        then fail := (Sum_Norm(y) > epsfa);
        !           248:             return;                         -- singular Jacobian matrix
        !           249:       end if;
        !           250:       deltax := -y;
        !           251:       lusolve(jacobian,n,ipvt,deltax);
        !           252:       Add(zero.v,deltax);                   -- make the updates
        !           253:       y := eval(p_eval,zero.v);
        !           254:       zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
        !           255:       numit := numit + 1;
        !           256:       if ( zero.err < epsxa ) or else ( zero.res < epsfa )
        !           257:                                             -- stopping criteria
        !           258:        then fail := false; exit;
        !           259:        elsif numit >= max
        !           260:            then fail := true; exit;
        !           261:       end if;
        !           262:     end loop;
        !           263:     jacobian := eval(j_eval,zero.v);        -- compute condition number
        !           264:     lufco(jacobian,n,ipvt,zero.rco);
        !           265:   exception
        !           266:     when numeric_error | constraint_error => fail := true; return;
        !           267:   end Silent_Newton;
        !           268:
        !           269:   procedure Silent_Newton
        !           270:                ( p_eval : in Standard_Complex_Poly_SysFun.Evaluator;
        !           271:                  j_eval : in Standard_Complex_Jaco_Matrices.Evaluator;
        !           272:                  zero : in out Solution; epsxa,epsfa : in double_float;
        !           273:                  numit : in out natural; max : in natural;
        !           274:                  fail : out boolean ) is
        !           275:
        !           276:     n : natural := zero.n;
        !           277:     jacobian : matrix(1..n,1..n);
        !           278:     ipvt : Standard_Natural_Vectors.Vector(1..n);
        !           279:     y,deltax : Vector(1..n);
        !           280:
        !           281:   begin
        !           282:     y := p_eval(zero.v);                    -- y = f(zero)
        !           283:     for i in 1..max loop
        !           284:       jacobian := j_eval(zero.v);           -- solve jacobian*deltax = -f(zero)
        !           285:       lufco(jacobian,n,ipvt,zero.rco);
        !           286:       if 1.0 + zero.rco = 1.0
        !           287:        then fail := (Sum_Norm(y) > epsfa);
        !           288:             return;                         -- singular Jacobian matrix
        !           289:       end if;
        !           290:       deltax := -y;
        !           291:       lusolve(jacobian,n,ipvt,deltax);
        !           292:       Add(zero.v,deltax);                   -- make the updates
        !           293:       y := p_eval(zero.v);
        !           294:       zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
        !           295:       numit := numit + 1;
        !           296:       if ( zero.err < epsxa ) or else ( zero.res < epsfa )
        !           297:                                             -- stopping criteria
        !           298:        then fail := false; exit;
        !           299:        elsif numit >= max
        !           300:            then fail := true; exit;
        !           301:       end if;
        !           302:     end loop;
        !           303:     jacobian := j_eval(zero.v);            -- compute condition number
        !           304:     lufco(jacobian,n,ipvt,zero.rco);
        !           305:   exception
        !           306:     when numeric_error | constraint_error => fail := true; return;
        !           307:   end Silent_Newton;
        !           308:
        !           309:   procedure Reporting_Newton
        !           310:                ( file : in file_type;
        !           311:                  p_eval : in Eval_Poly_Sys; j_eval : in Eval_Jaco_Mat;
        !           312:                  zero : in out Solution; epsxa,epsfa : in double_float;
        !           313:                  numit : in out natural; max : in natural;
        !           314:                  fail : out boolean ) is
        !           315:
        !           316:     n : natural := zero.n;
        !           317:     jacobian : matrix(1..n,1..n);
        !           318:     ipvt : Standard_Natural_Vectors.Vector(1..n);
        !           319:     y,deltax : Vector(1..n);
        !           320:
        !           321:   begin
        !           322:     y := Eval(p_eval,zero.v);              -- y = f(zero)
        !           323:     for i in 1..max loop
        !           324:       jacobian := eval(j_eval,zero.v);     -- solve jacobian*deltax = -f(zero)
        !           325:       lufco(jacobian,n,ipvt,zero.rco);
        !           326:       if 1.0 + zero.rco = 1.0              -- singular jacobian matrix
        !           327:        then fail := (Sum_Norm(y) > epsfa); -- accuracy not reached yet
        !           328:             return;
        !           329:       end if;
        !           330:       deltax := -y;
        !           331:       lusolve(jacobian,n,ipvt,deltax);
        !           332:       Add(zero.v,deltax);                  -- make the updates
        !           333:       y := eval(p_eval,zero.v);
        !           334:       zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
        !           335:       numit := numit + 1;
        !           336:       put(file,"Step "); put(file,numit,4); new_line(file);      -- output
        !           337:       put(file," |errxa| : "); put(file,zero.err); new_line(file);
        !           338:       put(file," |errfa| : "); put(file,zero.res); new_line(file);
        !           339:       if ( zero.err < epsxa ) or else ( zero.res < epsfa )
        !           340:                                                   -- stopping criteria
        !           341:        then fail := false; exit;
        !           342:        elsif numit >= max
        !           343:            then fail := true; exit;
        !           344:       end if;
        !           345:     end loop;
        !           346:     jacobian := eval(j_eval,zero.v);              -- compute condition number
        !           347:     lufco(jacobian,n,ipvt,zero.rco);
        !           348:   exception
        !           349:     when numeric_error | constraint_error => fail := true; return;
        !           350:   end Reporting_Newton;
        !           351:
        !           352:   procedure Reporting_Newton
        !           353:                ( file : in file_type;
        !           354:                  p_eval : in Standard_Complex_Poly_SysFun.Evaluator;
        !           355:                  j_eval : in Standard_Complex_Jaco_Matrices.Evaluator;
        !           356:                  zero : in out Solution; epsxa,epsfa : in double_float;
        !           357:                  numit : in out natural; max : in natural;
        !           358:                  fail : out boolean ) is
        !           359:
        !           360:     n : natural := zero.n;
        !           361:     jacobian : matrix(1..n,1..n);
        !           362:     ipvt : Standard_Natural_Vectors.Vector(1..n);
        !           363:     y,deltax : Vector(1..n);
        !           364:
        !           365:   begin
        !           366:     y := p_eval(zero.v);                   -- y = f(zero)
        !           367:     for i in 1..max loop
        !           368:       jacobian := j_eval(zero.v);          -- solve jacobian*deltax = -f(zero)
        !           369:       lufco(jacobian,n,ipvt,zero.rco);
        !           370:       if 1.0 + zero.rco = 1.0              -- singular jacobian matrix
        !           371:        then fail := (Sum_Norm(y) > epsfa); -- accuracy not reached yet
        !           372:             return;
        !           373:       end if;
        !           374:       deltax := -y;
        !           375:       lusolve(jacobian,n,ipvt,deltax);
        !           376:       Add(zero.v,deltax);                  -- make the updates
        !           377:       y := p_eval(zero.v);
        !           378:       zero.err := Sum_Norm(deltax); zero.res := Sum_Norm(y);
        !           379:       numit := numit + 1;
        !           380:       put(file,"Step "); put(file,numit,4); new_line(file);      -- output
        !           381:       put(file," |errxa| : "); put(file,zero.err); new_line(file);
        !           382:       put(file," |errfa| : "); put(file,zero.res); new_line(file);
        !           383:       if ( zero.err < epsxa ) or else ( zero.res < epsfa )
        !           384:                                                   -- stopping criteria
        !           385:        then fail := false; exit;
        !           386:        elsif numit >= max
        !           387:            then fail := true; exit;
        !           388:       end if;
        !           389:     end loop;
        !           390:     jacobian := j_eval(zero.v);                   -- compute condition number
        !           391:     lufco(jacobian,n,ipvt,zero.rco);
        !           392:   exception
        !           393:     when numeric_error | constraint_error => fail := true; return;
        !           394:   end Reporting_Newton;
        !           395:
        !           396:   procedure Silent_Root_Refiner
        !           397:                ( p : in Poly_Sys; sols : in out Solution_List;
        !           398:                  epsxa,epsfa,tolsing : in double_float;
        !           399:                  numit : in out natural; max : in natural ) is
        !           400:
        !           401:     n : natural := p'length;
        !           402:     p_eval : Eval_Poly_Sys(1..n) := Create(p);
        !           403:     jac : Jaco_Mat(1..n,1..n) := Create(p);
        !           404:     jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
        !           405:     numb : natural;
        !           406:     fail : boolean;
        !           407:     sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
        !           408:
        !           409:   begin
        !           410:     for i in sa'range loop
        !           411:       numb := 0;
        !           412:       sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
        !           413:       if sa(i).res < 1.0
        !           414:        then Silent_Newton(p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           415:        else fail := true;
        !           416:       end if;
        !           417:       Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
        !           418:       numit := numit + numb;
        !           419:     end loop;
        !           420:     clear(p_eval); clear(jac); clear(jac_eval);
        !           421:     Deep_Clear(sols); sols := Create(sa); Clear(sa);
        !           422:   end Silent_Root_Refiner;
        !           423:
        !           424:   procedure Silent_Root_Refiner
        !           425:                ( p : in Standard_Complex_Poly_SysFun.Evaluator;
        !           426:                  j : in Standard_Complex_Jaco_Matrices.Evaluator;
        !           427:                  sols : in out Solution_List;
        !           428:                  epsxa,epsfa,tolsing : in double_float;
        !           429:                  numit : in out natural; max : in natural ) is
        !           430:
        !           431:     n : constant natural := Head_Of(sols).n;
        !           432:     numb : natural;
        !           433:     fail : boolean;
        !           434:     sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
        !           435:
        !           436:   begin
        !           437:     for i in sa'range loop
        !           438:       numb := 0;
        !           439:       sa(i).res := Sum_Norm(p(sa(i).v));
        !           440:       if sa(i).res < 1.0
        !           441:        then Silent_Newton(p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           442:        else fail := true;
        !           443:       end if;
        !           444:       Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
        !           445:       numit := numit + numb;
        !           446:     end loop;
        !           447:     Deep_Clear(sols); sols := Create(sa); Clear(sa);
        !           448:   end Silent_Root_Refiner;
        !           449:
        !           450:   procedure Silent_Root_Refiner
        !           451:                ( p : in Poly_Sys; sols,refsols : in out Solution_List;
        !           452:                  epsxa,epsfa,tolsing : in double_float;
        !           453:                  numit : in out natural; max : in natural ) is
        !           454:
        !           455:     n : natural := p'length;
        !           456:     p_eval : Eval_Poly_Sys(1..n) := Create(p);
        !           457:     jac : Jaco_Mat(1..n,1..n) := Create(p);
        !           458:     jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
        !           459:     numb : natural;
        !           460:     fail : boolean;
        !           461:     sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
        !           462:     refsols_last : Solution_List;
        !           463:
        !           464:   begin
        !           465:     for i in sa'range loop
        !           466:       numb := 0;
        !           467:       sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
        !           468:       if sa(i).res < 1.0
        !           469:        then Silent_Newton(p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           470:             if not fail
        !           471:              then Append(refsols,refsols_last,sa(i).all);
        !           472:             end if;
        !           473:        else fail := true;
        !           474:       end if;
        !           475:       Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
        !           476:       numit := numit + numb;
        !           477:     end loop;
        !           478:     clear(p_eval); clear(jac); clear(jac_eval);
        !           479:     Deep_Clear(sols); sols := Create(sa); Clear(sa);
        !           480:   end Silent_Root_Refiner;
        !           481:
        !           482:   procedure Silent_Root_Refiner
        !           483:                ( p : in Standard_Complex_Poly_SysFun.Evaluator;
        !           484:                  j : in Standard_Complex_Jaco_Matrices.Evaluator;
        !           485:                  sols,refsols : in out Solution_List;
        !           486:                  epsxa,epsfa,tolsing : in double_float;
        !           487:                  numit : in out natural; max : in natural ) is
        !           488:
        !           489:     n : constant natural := Head_Of(sols).n;
        !           490:     numb : natural;
        !           491:     fail : boolean;
        !           492:     sa : Solution_Array(1..Length_Of(sols)) := Create(sols);
        !           493:     refsols_last : Solution_List;
        !           494:
        !           495:   begin
        !           496:     for i in sa'range loop
        !           497:       numb := 0;
        !           498:       sa(i).res := Sum_Norm(p(sa(i).v));
        !           499:       if sa(i).res < 1.0
        !           500:        then Silent_Newton(p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           501:             if not fail
        !           502:              then Append(refsols,refsols_last,sa(i).all);
        !           503:             end if;
        !           504:        else fail := true;
        !           505:       end if;
        !           506:       Root_Accounting(sa(i),i,sa(sa'first..i),fail,tolsing,epsxa);
        !           507:       numit := numit + numb;
        !           508:     end loop;
        !           509:     Deep_Clear(sols); sols := Create(sa); Clear(sa);
        !           510:   end Silent_Root_Refiner;
        !           511:
        !           512:   procedure Reporting_Root_Refiner
        !           513:                ( file : in file_type;
        !           514:                  p : in Poly_Sys; sols : in out Solution_List;
        !           515:                  epsxa,epsfa,tolsing : in double_float;
        !           516:                  numit : in out natural; max : in natural;
        !           517:                  wout : in boolean ) is
        !           518:
        !           519:     n : natural := p'length;
        !           520:     p_eval : Eval_Poly_Sys(1..n) := Create(p);
        !           521:     jac : Jaco_Mat(1..n,1..n) := Create(p);
        !           522:     jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
        !           523:     numb : natural;
        !           524:     nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
        !           525:     nbtot : constant natural := Length_Of(sols);
        !           526:     fail : boolean;
        !           527:     sa : Solution_Array(1..nbtot) := Create(sols);
        !           528:
        !           529:   begin
        !           530:     new_line(file);
        !           531:     put_line(file,"THE SOLUTIONS :"); new_line(file);
        !           532:     put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
        !           533:     Write_Bar(file);
        !           534:     for i in sa'range loop
        !           535:       numb := 0;
        !           536:       sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
        !           537:       if sa(i).res < 1.0
        !           538:        then
        !           539:          if wout
        !           540:           then Reporting_Newton
        !           541:                  (file,p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           542:           else Silent_Newton
        !           543:                  (p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           544:          end if;
        !           545:        else
        !           546:          fail := true;
        !           547:       end if;
        !           548:       Write_Info(file,sa(i).all,i,n,numb,fail);
        !           549:       Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
        !           550:                       nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
        !           551:       numit := numit + numb;
        !           552:     end loop;
        !           553:     Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
        !           554:     clear(p_eval); clear(jac); clear(jac_eval);
        !           555:     Deep_Clear(sols); sols := Create(sa); Clear(sa);
        !           556:   end Reporting_Root_Refiner;
        !           557:
        !           558:   procedure Reporting_Root_Refiner
        !           559:                ( file : in file_type;
        !           560:                  p : in Standard_Complex_Poly_SysFun.Evaluator;
        !           561:                  j : in Standard_Complex_Jaco_Matrices.Evaluator;
        !           562:                  sols : in out Solution_List;
        !           563:                  epsxa,epsfa,tolsing : in double_float;
        !           564:                  numit : in out natural; max : in natural;
        !           565:                  wout : in boolean ) is
        !           566:
        !           567:     n : constant natural := Head_Of(sols).n;
        !           568:     numb : natural;
        !           569:     nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
        !           570:     nbtot : constant natural := Length_Of(sols);
        !           571:     fail : boolean;
        !           572:     sa : Solution_Array(1..nbtot) := Create(sols);
        !           573:
        !           574:   begin
        !           575:     new_line(file);
        !           576:     put_line(file,"THE SOLUTIONS :"); new_line(file);
        !           577:     put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
        !           578:     Write_Bar(file);
        !           579:     for i in sa'range loop
        !           580:       numb := 0;
        !           581:       sa(i).res := Sum_Norm(p(sa(i).v));
        !           582:       if sa(i).res < 1.0
        !           583:        then
        !           584:          if wout
        !           585:           then Reporting_Newton(file,p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           586:           else Silent_Newton(p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           587:          end if;
        !           588:        else
        !           589:          fail := true;
        !           590:       end if;
        !           591:       Write_Info(file,sa(i).all,i,n,numb,fail);
        !           592:       Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
        !           593:                       nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
        !           594:       numit := numit + numb;
        !           595:     end loop;
        !           596:     Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
        !           597:     Deep_Clear(sols); sols := Create(sa); Clear(sa);
        !           598:   end Reporting_Root_Refiner;
        !           599:
        !           600:   procedure Reporting_Root_Refiner
        !           601:                ( file : in file_type;
        !           602:                  p : in Poly_Sys; sols,refsols : in out Solution_List;
        !           603:                  epsxa,epsfa,tolsing : in double_float;
        !           604:                  numit : in out natural; max : in natural;
        !           605:                  wout : in boolean ) is
        !           606:
        !           607:     n : natural := p'length;
        !           608:     p_eval : Eval_Poly_Sys(1..n) := Create(p);
        !           609:     jac : Jaco_Mat(1..n,1..n) := Create(p);
        !           610:     jac_eval : Eval_Jaco_Mat(1..n,1..n) := Create(jac);
        !           611:     numb : natural;
        !           612:     nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
        !           613:     nbtot : constant natural := Length_Of(sols);
        !           614:     fail : boolean;
        !           615:     sa : Solution_Array(1..nbtot) := Create(sols);
        !           616:     refsols_last : Solution_List;
        !           617:
        !           618:   begin
        !           619:     new_line(file);
        !           620:     put_line(file,"THE SOLUTIONS :"); new_line(file);
        !           621:     put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
        !           622:     Write_Bar(file);
        !           623:     for i in sa'range loop
        !           624:       numb := 0;
        !           625:       sa(i).res := Sum_Norm(Eval(p_eval,sa(i).v));
        !           626:       if sa(i).res < 1.0
        !           627:        then
        !           628:          if wout
        !           629:           then Reporting_Newton
        !           630:                  (file,p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           631:           else Silent_Newton
        !           632:                  (p_eval,jac_eval,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           633:          end if;
        !           634:          if not fail
        !           635:           then Append(refsols,refsols_last,sa(i).all);
        !           636:          end if;
        !           637:        else
        !           638:          fail := true;
        !           639:       end if;
        !           640:       Write_Info(file,sa(i).all,i,n,numb,fail);
        !           641:       Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
        !           642:                       nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
        !           643:       numit := numit + numb;
        !           644:     end loop;
        !           645:     Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
        !           646:     clear(p_eval); clear(jac); clear(jac_eval);
        !           647:     Deep_Clear(sols); sols := Create(sa); Clear(sa);
        !           648:   end Reporting_Root_Refiner;
        !           649:
        !           650:   procedure Reporting_Root_Refiner
        !           651:                ( file : in file_type;
        !           652:                  p : in Standard_Complex_Poly_SysFun.Evaluator;
        !           653:                  j : in Standard_Complex_Jaco_Matrices.Evaluator;
        !           654:                  sols,refsols : in out Solution_List;
        !           655:                  epsxa,epsfa,tolsing : in double_float;
        !           656:                  numit : in out natural; max : in natural;
        !           657:                  wout : in boolean ) is
        !           658:
        !           659:     n : constant natural := Head_Of(sols).n;
        !           660:     numb : natural;
        !           661:     nbfail,nbreg,nbsing,nbclus,nbreal,nbcomp : natural := 0;
        !           662:     nbtot : constant natural := Length_Of(sols);
        !           663:     fail : boolean;
        !           664:     sa : Solution_Array(1..nbtot) := Create(sols);
        !           665:     refsols_last : Solution_List;
        !           666:
        !           667:   begin
        !           668:     new_line(file);
        !           669:     put_line(file,"THE SOLUTIONS :"); new_line(file);
        !           670:     put(file,nbtot,1); put(file," "); put(file,n,1); new_line(file);
        !           671:     Write_Bar(file);
        !           672:     for i in sa'range loop
        !           673:       numb := 0;
        !           674:       sa(i).res := Sum_Norm(p(sa(i).v));
        !           675:       if sa(i).res < 1.0
        !           676:        then
        !           677:          if wout
        !           678:           then Reporting_Newton(file,p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           679:           else Silent_Newton(p,j,sa(i).all,epsxa,epsfa,numb,max,fail);
        !           680:          end if;
        !           681:          if not fail
        !           682:           then Append(refsols,refsols_last,sa(i).all);
        !           683:          end if;
        !           684:        else
        !           685:          fail := true;
        !           686:       end if;
        !           687:       Write_Info(file,sa(i).all,i,n,numb,fail);
        !           688:       Root_Accounting(file,sa(i),i,sa(sa'first..i),fail,tolsing,epsxa,
        !           689:                       nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
        !           690:       numit := numit + numb;
        !           691:     end loop;
        !           692:     Write_Global_Info(file,nbtot,nbfail,nbreal,nbcomp,nbreg,nbsing,nbclus);
        !           693:     Deep_Clear(sols); sols := Create(sa); Clear(sa);
        !           694:   end Reporting_Root_Refiner;
        !           695:
        !           696: end Standard_Root_Refiners;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>