Annotation of OpenXM_contrib/PHC/Ada/Continuation/correctors.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Natural_Vectors;
! 2: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
! 3: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
! 4: with Standard_Complex_QR_Decomposition; use Standard_Complex_QR_Decomposition;
! 5: with Standard_Complex_Least_Squares; use Standard_Complex_Least_Squares;
! 6: with Process_io; use Process_io;
! 7:
! 8: package body Correctors is
! 9:
! 10: -- AUXILIARIES FOR IMPLEMENTING MULTIPLE CORRECTORS :
! 11:
! 12: procedure Equals ( s : in Solu_Info_Array; x : in Vector; i : in natural;
! 13: d : in double_float; j : in out natural ) is
! 14: eq : boolean;
! 15:
! 16: begin
! 17: while j < i loop
! 18: eq := true;
! 19: for k in x'range loop
! 20: eq := Equal(s(j).sol.v(k),x(k),d);
! 21: exit when not eq;
! 22: end loop;
! 23: exit when eq;
! 24: j := j + 1;
! 25: end loop;
! 26: end Equals;
! 27:
! 28: generic
! 29: with procedure Corrector ( s : in out Solu_Info; c : in Corr_Pars );
! 30: procedure Multiple_Silent_Corrector
! 31: ( s : in out Solu_Info_Array;
! 32: pivot : in out natural; dist_sols : in double_float;
! 33: c : in Corr_Pars; fail : out boolean );
! 34:
! 35: -- DESCRIPTION :
! 36: -- Allows to build any multiple silent corrector,
! 37: -- depending on the corrector supplied as generic parameter.
! 38:
! 39: generic
! 40: with procedure Corrector ( file : in file_type;
! 41: s : in out Solu_Info; c : in Corr_Pars );
! 42: procedure Multiple_Reporting_Corrector
! 43: ( file : in file_type; s : in out Solu_Info_Array;
! 44: pivot : in out natural; dist_sols : in double_float;
! 45: c : in Corr_Pars; fail : out boolean );
! 46:
! 47: -- DESCRIPTION :
! 48: -- Allows to build any multiple reporting corrector,
! 49: -- depending on the corrector supplied as generic parameter.
! 50:
! 51: procedure Multiple_Silent_Corrector
! 52: ( s : in out Solu_Info_Array;
! 53: pivot : in out natural; dist_sols : in double_float;
! 54: c : in Corr_Pars; fail : out boolean ) is
! 55:
! 56: i,j : natural;
! 57: ffail : boolean := false;
! 58:
! 59: begin
! 60: i := pivot;
! 61: loop
! 62: Corrector(s(i),c);
! 63: ffail := (((s(i).resa > c.epsaf) and then (s(i).cora > c.epsax))
! 64: or else ((s(i).resr > c.epsrf) and then (s(i).corr > c.epsrx)));
! 65: if not ffail
! 66: then j := 1;
! 67: Equals(s,s(i).sol.v,i,dist_sols,j);
! 68: if j /= i
! 69: then ffail := true;
! 70: end if;
! 71: end if;
! 72: exit when ffail;
! 73: i := i + 1;
! 74: if i > s'last
! 75: then i := s'first;
! 76: end if;
! 77: exit when (i = pivot);
! 78: end loop;
! 79: if ffail
! 80: then pivot := i;
! 81: end if;
! 82: fail := ffail;
! 83: end Multiple_Silent_Corrector;
! 84:
! 85: procedure Multiple_Reporting_Corrector
! 86: ( file : in file_type; s : in out Solu_Info_Array;
! 87: pivot : in out natural; dist_sols : in double_float;
! 88: c : in Corr_Pars; fail : out boolean ) is
! 89:
! 90: i,j : natural;
! 91: ffail : boolean := false;
! 92:
! 93: begin
! 94: i := pivot;
! 95: loop
! 96: Write_path(file,i);
! 97: Corrector(file,s(i),c);
! 98: sWrite(file,s(i).sol.all);
! 99: ffail := (((s(i).resa > c.epsaf) and then (s(i).cora > c.epsax))
! 100: or else ((s(i).resr > c.epsrf) and then (s(i).corr > c.epsrx)));
! 101: if not ffail
! 102: then j := 1;
! 103: Equals(s,s(i).sol.v,i,dist_sols,j);
! 104: if j /= i
! 105: then ffail := true;
! 106: end if;
! 107: end if;
! 108: exit when ffail;
! 109: i := i + 1;
! 110: if i > s'last
! 111: then i := s'first;
! 112: end if;
! 113: exit when (i = pivot);
! 114: end loop;
! 115: if ffail
! 116: then pivot := i;
! 117: end if;
! 118: fail := ffail;
! 119: end Multiple_Reporting_Corrector;
! 120:
! 121: -- TARGET PROCEDURES :
! 122:
! 123: procedure Affine_Single_Loose_Normal_Silent_Corrector
! 124: ( s : in out Solu_Info; c : in Corr_Pars ) is
! 125:
! 126: info : integer;
! 127: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 128: j : Matrix(y'range,y'range);
! 129: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 130: nit : natural := 0;
! 131: normv : double_float;
! 132:
! 133: begin
! 134: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 135: normv := Norm(s.sol.v);
! 136: if normv + 1.0 /= 1.0
! 137: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 138: end if;
! 139: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 140: j := dH(s.sol.v,s.sol.t);
! 141: lufac(j,y'last,ipvt,info);
! 142: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
! 143: Min(y);
! 144: lusolve(j,y'last,ipvt,y); s.cora := Norm(y);
! 145: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 146: y := H(s.sol.v,s.sol.t); s.resa := Norm(y);
! 147: normv := Norm(s.sol.v);
! 148: if normv + 1.0 /= 1.0
! 149: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 150: end if;
! 151: nit := nit + 1;
! 152: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 153: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 154: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 155: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 156: exception
! 157: when numeric_error => return;
! 158: end Affine_Single_Loose_Normal_Silent_Corrector;
! 159:
! 160: procedure Affine_Single_Loose_Normal_Reporting_Corrector
! 161: ( file : in file_type;
! 162: s : in out Solu_Info; c : in Corr_Pars ) is
! 163:
! 164: info : integer;
! 165: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 166: j : Matrix(y'range,y'range);
! 167: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 168: nit : natural := 0;
! 169: normv : double_float;
! 170:
! 171: begin
! 172: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 173: normv := Norm(s.sol.v);
! 174: if normv + 1.0 /= 1.0
! 175: then s.corr := s.cora / normv;
! 176: s.resr := s.resa / normv;
! 177: end if;
! 178: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 179: j := dH(s.sol.v,s.sol.t);
! 180: lufac(j,y'last,ipvt,info);
! 181: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
! 182: Min(y);
! 183: lusolve(j,y'last,ipvt,y); s.cora := Norm(y);
! 184: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 185: y := H(s.sol.v,s.sol.t); s.resa := Norm(y);
! 186: normv := Norm(s.sol.v);
! 187: if normv + 1.0 /= 1.0
! 188: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 189: end if;
! 190: cWrite(file,s.cora,s.corr,s.resa,s.resr); -- WRITE STATISTICS
! 191: nit := nit + 1;
! 192: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 193: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 194: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 195: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 196: exception
! 197: when numeric_error => return;
! 198: end Affine_Single_Loose_Normal_Reporting_Corrector;
! 199:
! 200: procedure Affine_Single_Loose_Conditioned_Silent_Corrector
! 201: ( s : in out Solu_Info; c : in Corr_Pars ) is
! 202:
! 203: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 204: n : constant natural := y'last;
! 205: j : Matrix(y'range,y'range);
! 206: qraux : Standard_Complex_Vectors.Vector(y'range)
! 207: := (y'range => Create(0.0));
! 208: dum : Standard_Complex_Vectors.Vector(y'range);
! 209: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 210: info : integer;
! 211: nit : natural := 0;
! 212: normv : double_float;
! 213:
! 214: begin
! 215: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 216: normv := Norm(s.sol.v);
! 217: if normv + 1.0 /= 1.0
! 218: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 219: end if;
! 220: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 221: j := dH(s.sol.v,s.sol.t);
! 222: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
! 223: Min(y); -- FOLLOWED BY LEAST SQUARES
! 224: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
! 225: s.cora := Norm(y);
! 226: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 227: y := H(s.sol.v,s.sol.t); s.resa := Norm(y);
! 228: normv := Norm(s.sol.v);
! 229: if normv + 1.0 /= 1.0
! 230: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 231: end if;
! 232: nit := nit + 1;
! 233: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 234: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 235: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 236: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 237: j := dH(s.sol.v,s.sol.t);
! 238: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
! 239: exception
! 240: when numeric_error => return;
! 241: end Affine_Single_Loose_Conditioned_Silent_Corrector;
! 242:
! 243: procedure Affine_Single_Loose_Conditioned_Reporting_Corrector
! 244: ( file : in file_type;
! 245: s : in out Solu_Info; c : in Corr_Pars ) is
! 246:
! 247: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 248: n : constant natural := y'last;
! 249: j : Matrix(y'range,y'range);
! 250: qraux : Standard_Complex_Vectors.Vector(y'range)
! 251: := (y'range => Create(0.0));
! 252: dum : Standard_Complex_Vectors.Vector(y'range);
! 253: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 254: info : integer;
! 255: nit : natural := 0;
! 256: normv : double_float;
! 257:
! 258: begin
! 259: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 260: normv := Norm(s.sol.v);
! 261: if normv + 1.0 /= 1.0
! 262: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 263: end if;
! 264: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 265: j := dH(s.sol.v,s.sol.t);
! 266: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
! 267: Min(y); -- FOLLOWED BY LEAST SQUARES
! 268: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
! 269: s.cora := Norm(y);
! 270: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 271: y := H(s.sol.v,s.sol.t); s.resa := Norm(y);
! 272: normv := Norm(s.sol.v);
! 273: if normv + 1.0 /= 1.0
! 274: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 275: end if;
! 276: cWrite(file,s.cora,s.corr,s.resa,s.resr); -- WRITE STATISTICS
! 277: cWrite(file,s.rcond,s.sol.m);
! 278: nit := nit + 1;
! 279: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 280: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 281: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 282: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 283: j := dH(s.sol.v,s.sol.t);
! 284: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
! 285: exception
! 286: when numeric_error => return;
! 287: end Affine_Single_Loose_Conditioned_Reporting_Corrector;
! 288:
! 289: procedure Affine_Single_Severe_Normal_Silent_Corrector
! 290: ( s : in out Solu_Info; c : in Corr_Pars ) is
! 291:
! 292: info : integer;
! 293: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 294: j : Matrix(y'range,y'range);
! 295: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 296: nit : natural := 0;
! 297: normv,ncora,nresa,ncorr,nresr : double_float;
! 298: stop : boolean;
! 299:
! 300: begin
! 301: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 302: normv := Norm(s.sol.v);
! 303: if normv + 1.0 /= 1.0
! 304: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 305: end if;
! 306: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 307: j := dH(s.sol.v,s.sol.t);
! 308: lufac(j,y'last,ipvt,info);
! 309: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
! 310: Min(y);
! 311: lusolve(j,y'last,ipvt,y); ncora := Norm(y);
! 312: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 313: y := H(s.sol.v,s.sol.t); nresa := Norm(y);
! 314: normv := Norm(s.sol.v);
! 315: if normv + 1.0 /= 1.0
! 316: then ncorr := ncora / normv; nresr := nresa / normv;
! 317: end if;
! 318: nit := nit + 1;
! 319: stop := (((ncora > s.cora) and then (nresa > s.resa))
! 320: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
! 321: -- STOP WHEN DIVERGENCE
! 322: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
! 323: exit when stop;
! 324: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 325: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 326: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 327: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 328: exception
! 329: when numeric_error => return;
! 330: end Affine_Single_Severe_Normal_Silent_Corrector;
! 331:
! 332: procedure Affine_Single_Severe_Normal_Reporting_Corrector
! 333: ( file : in file_type;
! 334: s : in out Solu_Info; c : in Corr_Pars ) is
! 335:
! 336: info : integer;
! 337: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 338: j : Matrix(y'range,y'range);
! 339: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 340: nit : natural := 0;
! 341: normv,ncora,nresa,ncorr,nresr : double_float;
! 342: stop : boolean;
! 343:
! 344: begin
! 345: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 346: normv := Norm(s.sol.v);
! 347: if normv + 1.0 /= 1.0
! 348: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 349: end if;
! 350: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 351: j := dH(s.sol.v,s.sol.t);
! 352: lufac(j,y'last,ipvt,info);
! 353: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
! 354: Min(y);
! 355: lusolve(j,y'last,ipvt,y); ncora := Norm(y);
! 356: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 357: y := H(s.sol.v,s.sol.t); nresa := Norm(y);
! 358: normv := Norm(s.sol.v);
! 359: if normv + 1.0 /= 1.0
! 360: then ncorr := ncora / normv; nresr := nresa / normv;
! 361: end if;
! 362: cWrite(file,ncora,ncorr,nresa,nresr); -- WRITE STATISTICS
! 363: nit := nit + 1;
! 364: stop := (((ncora > s.cora) and then (nresa > s.resa))
! 365: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
! 366: -- STOP WHEN DIVERGENCE
! 367: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
! 368: exit when stop;
! 369: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 370: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 371: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 372: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 373: exception
! 374: when numeric_error => return;
! 375: end Affine_Single_Severe_Normal_Reporting_Corrector;
! 376:
! 377: procedure Affine_Single_Severe_Conditioned_Silent_Corrector
! 378: ( s : in out Solu_Info; c : in Corr_Pars ) is
! 379:
! 380: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 381: n : constant natural := y'last;
! 382: j : Matrix(y'range,y'range);
! 383: qraux : Standard_Complex_Vectors.Vector(y'range)
! 384: := (y'range => Create(0.0));
! 385: dum : Standard_Complex_Vectors.Vector(y'range);
! 386: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 387: info : integer;
! 388: nit : natural := 0;
! 389: normv,ncora,nresa,ncorr,nresr : double_float;
! 390: stop : boolean;
! 391:
! 392: begin
! 393: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 394: normv := Norm(s.sol.v);
! 395: if normv + 1.0 /= 1.0
! 396: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 397: end if;
! 398: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 399: j := dH(s.sol.v,s.sol.t);
! 400: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
! 401: Min(y); -- FOLLOWED BY LEAST SQUARES
! 402: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
! 403: ncora := Norm(y);
! 404: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 405: y := H(s.sol.v,s.sol.t); nresa := Norm(y);
! 406: normv := Norm(s.sol.v);
! 407: if normv + 1.0 /= 1.0
! 408: then ncorr := ncora / normv; nresr := nresa / normv;
! 409: end if;
! 410: nit := nit + 1;
! 411: stop := (((ncora > s.cora) and then (nresa > s.resa))
! 412: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
! 413: -- STOP WHEN DIVERGENCE
! 414: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
! 415: exit when stop;
! 416: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 417: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 418: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 419: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 420: j := dH(s.sol.v,s.sol.t);
! 421: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
! 422: exception
! 423: when numeric_error => return;
! 424: end Affine_Single_Severe_Conditioned_Silent_Corrector;
! 425:
! 426: procedure Affine_Single_Severe_Conditioned_Reporting_Corrector
! 427: ( file : in file_type;
! 428: s : in out Solu_Info; c : in Corr_Pars ) is
! 429:
! 430: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 431: n : constant natural := y'last;
! 432: j : Matrix(y'range,y'range);
! 433: qraux : Standard_Complex_Vectors.Vector(y'range)
! 434: := (y'range => Create(0.0));
! 435: dum : Standard_Complex_Vectors.Vector(y'range);
! 436: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 437: info : integer;
! 438: nit : natural := 0;
! 439: normv,ncora,nresa,ncorr,nresr : double_float;
! 440: stop : boolean;
! 441:
! 442: begin
! 443: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 444: normv := Norm(s.sol.v);
! 445: if normv + 1.0 /= 1.0
! 446: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 447: end if;
! 448: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 449: j := dH(s.sol.v,s.sol.t);
! 450: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
! 451: Min(y); -- FOLLOWED BY LEAST SQUARES
! 452: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
! 453: ncora := Norm(y);
! 454: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 455: y := H(s.sol.v,s.sol.t); nresa := Norm(y);
! 456: normv := Norm(s.sol.v);
! 457: if normv + 1.0 /= 1.0
! 458: then ncorr := ncora / normv; nresr := nresa / normv;
! 459: end if;
! 460: cWrite(file,ncora,ncorr,nresa,nresr); -- WRITE STATISTICS
! 461: cWrite(file,s.rcond,s.sol.m);
! 462: nit := nit + 1;
! 463: stop := (((ncora > s.cora) and then (nresa > s.resa))
! 464: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
! 465: -- STOP WHEN DIVERGENCE
! 466: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
! 467: exit when stop;
! 468: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 469: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 470: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 471: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 472: j := dH(s.sol.v,s.sol.t);
! 473: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
! 474: exception
! 475: when numeric_error => return;
! 476: end Affine_Single_Severe_Conditioned_Reporting_Corrector;
! 477:
! 478: procedure Affine_Multiple_Loose_Normal_Silent_Corrector
! 479: ( s : in out Solu_Info_Array;
! 480: pivot : in out natural; dist_sols : in double_float;
! 481: c : in Corr_Pars; fail : out boolean ) is
! 482:
! 483: procedure Single_Corrector is
! 484: new Affine_Single_Loose_Normal_Silent_Corrector(Norm,H,dH);
! 485: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
! 486:
! 487: begin
! 488: Corrector(s,pivot,dist_sols,c,fail);
! 489: end Affine_Multiple_Loose_Normal_Silent_Corrector;
! 490:
! 491: procedure Affine_Multiple_Loose_Normal_Reporting_Corrector
! 492: ( file : in file_type; s : in out Solu_Info_Array;
! 493: pivot : in out natural; dist_sols : in double_float;
! 494: c : in Corr_Pars; fail : out boolean ) is
! 495:
! 496: procedure Single_Corrector is
! 497: new Affine_Single_Loose_Normal_Reporting_Corrector(Norm,H,dH);
! 498: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
! 499:
! 500: begin
! 501: Corrector(file,s,pivot,dist_sols,c,fail);
! 502: end Affine_Multiple_Loose_Normal_Reporting_Corrector;
! 503:
! 504: procedure Affine_Multiple_Loose_Conditioned_Silent_Corrector
! 505: ( s : in out Solu_Info_Array;
! 506: pivot : in out natural; dist_sols : in double_float;
! 507: c : in Corr_Pars; fail : out boolean ) is
! 508:
! 509: procedure Single_Corrector is
! 510: new Affine_Single_Loose_Conditioned_Silent_Corrector(Norm,H,dH);
! 511: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
! 512:
! 513: begin
! 514: Corrector(s,pivot,dist_sols,c,fail);
! 515: end Affine_Multiple_Loose_Conditioned_Silent_Corrector;
! 516:
! 517: procedure Affine_Multiple_Loose_Conditioned_Reporting_Corrector
! 518: ( file : in file_type; s : in out Solu_Info_Array;
! 519: pivot : in out natural; dist_sols : in double_float;
! 520: c : in Corr_Pars; fail : out boolean ) is
! 521:
! 522: procedure Single_Corrector is
! 523: new Affine_Single_Loose_Conditioned_Reporting_Corrector(Norm,H,dH);
! 524: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
! 525:
! 526: begin
! 527: Corrector(file,s,pivot,dist_sols,c,fail);
! 528: end Affine_Multiple_Loose_Conditioned_Reporting_Corrector;
! 529:
! 530: procedure Affine_Multiple_Severe_Normal_Silent_Corrector
! 531: ( s : in out Solu_Info_Array;
! 532: pivot : in out natural; dist_sols : in double_float;
! 533: c : in Corr_Pars; fail : out boolean ) is
! 534:
! 535: procedure Single_Corrector is
! 536: new Affine_Single_Severe_Normal_Silent_Corrector(Norm,H,dH);
! 537: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
! 538:
! 539: begin
! 540: Corrector(s,pivot,dist_sols,c,fail);
! 541: end Affine_Multiple_Severe_Normal_Silent_Corrector;
! 542:
! 543: procedure Affine_Multiple_Severe_Normal_Reporting_Corrector
! 544: ( file : in file_type; s : in out Solu_Info_Array;
! 545: pivot : in out natural; dist_sols : in double_float;
! 546: c : in Corr_Pars; fail : out boolean ) is
! 547:
! 548: procedure Single_Corrector is
! 549: new Affine_Single_Severe_Normal_Reporting_Corrector(Norm,H,dH);
! 550: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
! 551:
! 552: begin
! 553: Corrector(file,s,pivot,dist_sols,c,fail);
! 554: end Affine_Multiple_Severe_Normal_Reporting_Corrector;
! 555:
! 556: procedure Affine_Multiple_Severe_Conditioned_Silent_Corrector
! 557: ( s : in out Solu_Info_Array;
! 558: pivot : in out natural; dist_sols : in double_float;
! 559: c : in Corr_Pars; fail : out boolean ) is
! 560:
! 561: procedure Single_Corrector is
! 562: new Affine_Single_Severe_Conditioned_Silent_Corrector(Norm,H,dH);
! 563: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
! 564:
! 565: begin
! 566: Corrector(s,pivot,dist_sols,c,fail);
! 567: end Affine_Multiple_Severe_Conditioned_Silent_Corrector;
! 568:
! 569: procedure Affine_Multiple_Severe_Conditioned_Reporting_Corrector
! 570: ( file : in file_type; s : in out Solu_Info_Array;
! 571: pivot : in out natural; dist_sols : in double_float;
! 572: c : in Corr_Pars; fail : out boolean ) is
! 573:
! 574: procedure Single_Corrector is
! 575: new Affine_Single_Severe_Conditioned_Reporting_Corrector(Norm,H,dH);
! 576: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
! 577:
! 578: begin
! 579: Corrector(file,s,pivot,dist_sols,c,fail);
! 580: end Affine_Multiple_Severe_Conditioned_Reporting_Corrector;
! 581:
! 582: procedure Projective_Single_Loose_Normal_Silent_Corrector
! 583: ( s : in out Solu_Info; c : in Corr_Pars ) is
! 584:
! 585: info : integer;
! 586: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 587: j : Matrix(y'range,y'range);
! 588: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 589: nit : natural := 0;
! 590: normv : double_float;
! 591:
! 592: begin
! 593: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 594: normv := Norm(s.sol.v);
! 595: if normv + 1.0 /= 1.0
! 596: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 597: end if;
! 598: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 599: j := dH(s.sol.v,s.sol.t);
! 600: for jj in y'range loop
! 601: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
! 602: end loop;
! 603: lufac(j,y'last,ipvt,info);
! 604: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
! 605: Min(y);
! 606: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
! 607: lusolve(j,y'last,ipvt,y); s.cora := Norm(y);
! 608: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 609: y := H(s.sol.v,s.sol.t);
! 610: y(y'last) := Create(0.0); s.resa := Norm(y);
! 611: normv := Norm(s.sol.v);
! 612: if normv + 1.0 /= 1.0
! 613: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 614: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
! 615: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
! 616: end loop;
! 617: end if;
! 618: nit := nit + 1;
! 619: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 620: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 621: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 622: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 623: exception
! 624: when numeric_error => return;
! 625: end Projective_Single_Loose_Normal_Silent_Corrector;
! 626:
! 627: procedure Projective_Single_Loose_Normal_Reporting_Corrector
! 628: ( file : in file_type;
! 629: s : in out Solu_Info; c : in Corr_Pars ) is
! 630:
! 631: info : integer;
! 632: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 633: j : Matrix(y'range,y'range);
! 634: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 635: nit : natural := 0;
! 636: normv : double_float;
! 637:
! 638: begin
! 639: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 640: normv := Norm(s.sol.v);
! 641: if normv + 1.0 /= 1.0
! 642: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 643: end if;
! 644: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 645: j := dH(s.sol.v,s.sol.t);
! 646: for jj in y'range loop
! 647: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
! 648: end loop;
! 649: lufac(j,y'last,ipvt,info);
! 650: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
! 651: Min(y);
! 652: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
! 653: lusolve(j,y'last,ipvt,y); s.cora := Norm(y);
! 654: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 655: y := H(s.sol.v,s.sol.t);
! 656: y(y'last) := Create(0.0); s.resa := Norm(y);
! 657: normv := Norm(s.sol.v);
! 658: if normv + 1.0 /= 1.0
! 659: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 660: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
! 661: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
! 662: end loop;
! 663: end if;
! 664: cWrite(file,s.cora,s.corr,s.resa,s.resr); -- WRITE STATISTICS
! 665: nit := nit + 1;
! 666: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 667: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 668: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 669: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 670: exception
! 671: when numeric_error => return;
! 672: end Projective_Single_Loose_Normal_Reporting_Corrector;
! 673:
! 674: procedure Projective_Single_Loose_Conditioned_Silent_Corrector
! 675: ( s : in out Solu_Info; c : in Corr_Pars ) is
! 676:
! 677: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 678: n : constant natural := y'last;
! 679: j : Matrix(y'range,y'range);
! 680: qraux : Standard_Complex_Vectors.Vector(y'range)
! 681: := (y'range => Create(0.0));
! 682: dum : Standard_Complex_Vectors.Vector(y'range);
! 683: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 684: info : integer;
! 685: nit : natural := 0;
! 686: normv : double_float;
! 687:
! 688: begin
! 689: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 690: normv := Norm(s.sol.v);
! 691: if normv + 1.0 /= 1.0
! 692: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 693: end if;
! 694: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 695: j := dH(s.sol.v,s.sol.t);
! 696: for jj in y'range loop
! 697: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
! 698: end loop;
! 699: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
! 700: Min(y); -- FOLLOWED BY LEAST SQUARES
! 701: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
! 702: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
! 703: s.cora := Norm(y);
! 704: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 705: y := H(s.sol.v,s.sol.t);
! 706: y(y'last) := Create(0.0); s.resa := Norm(y);
! 707: normv := Norm(s.sol.v);
! 708: if normv + 1.0 /= 1.0
! 709: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 710: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
! 711: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
! 712: end loop;
! 713: end if;
! 714: nit := nit + 1;
! 715: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 716: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 717: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 718: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 719: j := dH(s.sol.v,s.sol.t);
! 720: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
! 721: exception
! 722: when numeric_error => return;
! 723: end Projective_Single_Loose_Conditioned_Silent_Corrector;
! 724:
! 725: procedure Projective_Single_Loose_Conditioned_Reporting_Corrector
! 726: ( file : in file_type;
! 727: s : in out Solu_Info; c : in Corr_Pars ) is
! 728:
! 729: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 730: n : constant natural := y'last;
! 731: j : Matrix(y'range,y'range);
! 732: qraux : Standard_Complex_Vectors.Vector(y'range)
! 733: := (y'range => Create(0.0));
! 734: dum : Standard_Complex_Vectors.Vector(y'range);
! 735: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 736: info : integer;
! 737: nit : natural := 0;
! 738: normv : double_float;
! 739:
! 740: begin
! 741: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 742: normv := Norm(s.sol.v);
! 743: if normv + 1.0 /= 1.0
! 744: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 745: end if;
! 746: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 747: j := dH(s.sol.v,s.sol.t);
! 748: for jj in y'range loop
! 749: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
! 750: end loop;
! 751: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
! 752: Min(y); -- FOLLOWED BY LEAST SQUARES
! 753: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
! 754: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
! 755: s.cora := Norm(y);
! 756: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 757: y := H(s.sol.v,s.sol.t);
! 758: y(y'last) := Create(0.0); s.resa := Norm(y);
! 759: normv := Norm(s.sol.v);
! 760: if normv + 1.0 /= 1.0
! 761: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 762: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
! 763: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
! 764: end loop;
! 765: end if;
! 766: cWrite(file,s.cora,s.corr,s.resa,s.resr); -- WRITE STATISTICS
! 767: cWrite(file,s.rcond,s.sol.m);
! 768: nit := nit + 1;
! 769: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 770: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 771: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 772: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 773: j := dH(s.sol.v,s.sol.t);
! 774: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
! 775: exception
! 776: when numeric_error => return;
! 777: end Projective_Single_Loose_Conditioned_Reporting_Corrector;
! 778:
! 779: procedure Projective_Single_Severe_Normal_Silent_Corrector
! 780: ( s : in out Solu_Info; c : in Corr_Pars ) is
! 781:
! 782: info : integer;
! 783: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 784: j : Matrix(y'range,y'range);
! 785: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 786: nit : natural := 0;
! 787: normv,ncora,nresa,ncorr,nresr : double_float;
! 788: stop : boolean;
! 789:
! 790: begin
! 791: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 792: normv := Norm(s.sol.v);
! 793: if normv + 1.0 /= 1.0
! 794: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 795: end if;
! 796: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 797: j := dH(s.sol.v,s.sol.t);
! 798: for jj in y'range loop
! 799: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
! 800: end loop;
! 801: lufac(j,y'last,ipvt,info);
! 802: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
! 803: Min(y);
! 804: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
! 805: lusolve(j,y'last,ipvt,y); ncora := Norm(y);
! 806: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 807: y := H(s.sol.v,s.sol.t);
! 808: y(y'last) := Create(0.0); nresa := Norm(y);
! 809: normv := Norm(s.sol.v);
! 810: if normv + 1.0 /= 1.0
! 811: then ncorr := ncora / normv; nresr := nresa / normv;
! 812: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
! 813: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
! 814: end loop;
! 815: end if;
! 816: nit := nit + 1;
! 817: stop := (((ncora > s.cora) and then (nresa > s.resa))
! 818: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
! 819: -- STOP WHEN DIVERGENCE
! 820: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
! 821: exit when stop;
! 822: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 823: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 824: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 825: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 826: exception
! 827: when numeric_error => return;
! 828: end Projective_Single_Severe_Normal_Silent_Corrector;
! 829:
! 830: procedure Projective_Single_Severe_Normal_Reporting_Corrector
! 831: ( file : in file_type;
! 832: s : in out Solu_Info; c : in Corr_Pars ) is
! 833:
! 834: info : integer;
! 835: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 836: j : Matrix(y'range,y'range);
! 837: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 838: nit : natural := 0;
! 839: normv,ncora,nresa,ncorr,nresr : double_float;
! 840: stop : boolean;
! 841:
! 842: begin
! 843: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 844: normv := Norm(s.sol.v);
! 845: if normv + 1.0 /= 1.0
! 846: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 847: end if;
! 848: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 849: j := dH(s.sol.v,s.sol.t);
! 850: for jj in y'range loop
! 851: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
! 852: end loop;
! 853: lufac(j,y'last,ipvt,info);
! 854: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
! 855: Min(y);
! 856: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
! 857: lusolve(j,y'last,ipvt,y); ncora := Norm(y);
! 858: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 859: y := H(s.sol.v,s.sol.t);
! 860: y(y'last) := Create(0.0); nresa := Norm(y);
! 861: normv := Norm(s.sol.v);
! 862: if normv + 1.0 /= 1.0
! 863: then ncorr := ncora / normv; nresr := nresa / normv;
! 864: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
! 865: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
! 866: end loop;
! 867: end if;
! 868: cWrite(file,ncora,ncorr,nresa,nresr); -- WRITE STATISTICS
! 869: nit := nit + 1;
! 870: stop := (((ncora > s.cora) and then (nresa > s.resa))
! 871: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
! 872: -- STOP WHEN DIVERGENCE
! 873: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
! 874: exit when stop;
! 875: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 876: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 877: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 878: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 879: exception
! 880: when numeric_error => return;
! 881: end Projective_Single_Severe_Normal_Reporting_Corrector;
! 882:
! 883: procedure Projective_Single_Severe_Conditioned_Silent_Corrector
! 884: ( s : in out Solu_Info; c : in Corr_Pars ) is
! 885:
! 886: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 887: n : constant natural := y'last;
! 888: j : Matrix(y'range,y'range);
! 889: qraux : Standard_Complex_Vectors.Vector(y'range)
! 890: := (y'range => Create(0.0));
! 891: dum : Standard_Complex_Vectors.Vector(y'range);
! 892: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 893: info : integer;
! 894: nit : natural := 0;
! 895: normv,ncora,nresa,ncorr,nresr : double_float;
! 896: stop : boolean;
! 897:
! 898: begin
! 899: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 900: normv := Norm(s.sol.v);
! 901: if normv + 1.0 /= 1.0
! 902: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 903: end if;
! 904: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 905: j := dH(s.sol.v,s.sol.t);
! 906: for jj in y'range loop
! 907: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
! 908: end loop;
! 909: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
! 910: Min(y); -- FOLLOWED BY LEAST SQUARES
! 911: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
! 912: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
! 913: ncora := Norm(y);
! 914: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 915: y := H(s.sol.v,s.sol.t);
! 916: y(y'last) := Create(0.0); nresa := Norm(y);
! 917: normv := Norm(s.sol.v);
! 918: if normv + 1.0 /= 1.0
! 919: then ncorr := ncora / normv; nresr := nresa / normv;
! 920: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
! 921: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
! 922: end loop;
! 923: end if;
! 924: nit := nit + 1;
! 925: stop := (((ncora > s.cora) and then (nresa > s.resa))
! 926: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
! 927: -- STOP WHEN DIVERGENCE
! 928: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
! 929: exit when stop;
! 930: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 931: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 932: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 933: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 934: j := dH(s.sol.v,s.sol.t);
! 935: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
! 936: exception
! 937: when numeric_error => return;
! 938: end Projective_Single_Severe_Conditioned_Silent_Corrector;
! 939:
! 940: procedure Projective_Single_Severe_Conditioned_Reporting_Corrector
! 941: ( file : in file_type;
! 942: s : in out Solu_Info; c : in Corr_Pars ) is
! 943:
! 944: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
! 945: n : constant natural := y'last;
! 946: j : Matrix(y'range,y'range);
! 947: qraux : Standard_Complex_Vectors.Vector(y'range)
! 948: := (y'range => Create(0.0));
! 949: dum : Standard_Complex_Vectors.Vector(y'range);
! 950: ipvt : Standard_Natural_Vectors.Vector(y'range);
! 951: info : integer;
! 952: nit : natural := 0;
! 953: normv,ncora,nresa,ncorr,nresr : double_float;
! 954: stop : boolean;
! 955:
! 956: begin
! 957: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
! 958: normv := Norm(s.sol.v);
! 959: if normv + 1.0 /= 1.0
! 960: then s.corr := s.cora / normv; s.resr := s.resa / normv;
! 961: end if;
! 962: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
! 963: j := dH(s.sol.v,s.sol.t);
! 964: for jj in y'range loop
! 965: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
! 966: end loop;
! 967: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
! 968: Min(y); -- FOLLOWED BY LEAST SQUARES
! 969: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
! 970: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
! 971: ncora := Norm(y);
! 972: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
! 973: y := H(s.sol.v,s.sol.t);
! 974: y(y'last) := Create(0.0); nresa := Norm(y);
! 975: normv := Norm(s.sol.v);
! 976: if normv + 1.0 /= 1.0
! 977: then ncorr := ncora / normv; nresr := nresa / normv;
! 978: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
! 979: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
! 980: end loop;
! 981: end if;
! 982: cWrite(file,ncora,ncorr,nresa,nresr); -- WRITE STATISTICS
! 983: cWrite(file,s.rcond,s.sol.m);
! 984: nit := nit + 1;
! 985: stop := (((ncora > s.cora) and then (nresa > s.resa))
! 986: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
! 987: -- STOP WHEN DIVERGENCE
! 988: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
! 989: exit when stop;
! 990: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
! 991: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
! 992: end loop; -- STOP WHEN DESIRED PRECISION REACHED
! 993: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
! 994: j := dH(s.sol.v,s.sol.t);
! 995: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
! 996: exception
! 997: when numeric_error => return;
! 998: end Projective_Single_Severe_Conditioned_Reporting_Corrector;
! 999:
! 1000: procedure Projective_Multiple_Loose_Normal_Silent_Corrector
! 1001: ( s : in out Solu_Info_Array;
! 1002: pivot : in out natural; dist_sols : in double_float;
! 1003: c : in Corr_Pars; fail : out boolean ) is
! 1004:
! 1005: procedure Single_Corrector is
! 1006: new Projective_Single_Loose_Normal_Silent_Corrector(Norm,H,dH);
! 1007: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
! 1008:
! 1009: begin
! 1010: Corrector(s,pivot,dist_sols,c,fail);
! 1011: end Projective_Multiple_Loose_Normal_Silent_Corrector;
! 1012:
! 1013: procedure Projective_Multiple_Loose_Normal_Reporting_Corrector
! 1014: ( file : in file_type; s : in out Solu_Info_Array;
! 1015: pivot : in out natural; dist_sols : in double_float;
! 1016: c : in Corr_Pars; fail : out boolean ) is
! 1017:
! 1018: procedure Single_Corrector is
! 1019: new Projective_Single_Loose_Normal_Reporting_Corrector(Norm,H,dH);
! 1020: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
! 1021:
! 1022: begin
! 1023: Corrector(file,s,pivot,dist_sols,c,fail);
! 1024: end Projective_Multiple_Loose_Normal_Reporting_Corrector;
! 1025:
! 1026: procedure Projective_Multiple_Loose_Conditioned_Silent_Corrector
! 1027: ( s : in out Solu_Info_Array;
! 1028: pivot : in out natural; dist_sols : in double_float;
! 1029: c : in Corr_Pars; fail : out boolean ) is
! 1030:
! 1031: procedure Single_Corrector is
! 1032: new Projective_Single_Loose_Conditioned_Silent_Corrector(Norm,H,dH);
! 1033: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
! 1034:
! 1035: begin
! 1036: Corrector(s,pivot,dist_sols,c,fail);
! 1037: end Projective_Multiple_Loose_Conditioned_Silent_Corrector;
! 1038:
! 1039: procedure Projective_Multiple_Loose_Conditioned_Reporting_Corrector
! 1040: ( file : in file_type; s : in out Solu_Info_Array;
! 1041: pivot : in out natural; dist_sols : in double_float;
! 1042: c : in Corr_Pars; fail : out boolean ) is
! 1043:
! 1044: procedure Single_Corrector is
! 1045: new Projective_Single_Loose_Normal_Reporting_Corrector(Norm,H,dH);
! 1046: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
! 1047:
! 1048: begin
! 1049: Corrector(file,s,pivot,dist_sols,c,fail);
! 1050: end Projective_Multiple_Loose_Conditioned_Reporting_Corrector;
! 1051:
! 1052: procedure Projective_Multiple_Severe_Normal_Silent_Corrector
! 1053: ( s : in out Solu_Info_Array;
! 1054: pivot : in out natural; dist_sols : in double_float;
! 1055: c : in Corr_Pars; fail : out boolean ) is
! 1056:
! 1057: procedure Single_Corrector is
! 1058: new Projective_Single_Severe_Normal_Silent_Corrector(Norm,H,dH);
! 1059: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
! 1060:
! 1061: begin
! 1062: Corrector(s,pivot,dist_sols,c,fail);
! 1063: end Projective_Multiple_Severe_Normal_Silent_Corrector;
! 1064:
! 1065: procedure Projective_Multiple_Severe_Normal_Reporting_Corrector
! 1066: ( file : in file_type; s : in out Solu_Info_Array;
! 1067: pivot : in out natural; dist_sols : in double_float;
! 1068: c : in Corr_Pars; fail : out boolean ) is
! 1069:
! 1070: procedure Single_Corrector is
! 1071: new Projective_Single_Loose_Normal_Reporting_Corrector(Norm,H,dH);
! 1072: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
! 1073:
! 1074: begin
! 1075: Corrector(file,s,pivot,dist_sols,c,fail);
! 1076: end Projective_Multiple_Severe_Normal_Reporting_Corrector;
! 1077:
! 1078: procedure Projective_Multiple_Severe_Conditioned_Silent_Corrector
! 1079: ( s : in out Solu_Info_Array;
! 1080: pivot : in out natural; dist_sols : in double_float;
! 1081: c : in Corr_Pars; fail : out boolean ) is
! 1082:
! 1083: procedure Single_Corrector is
! 1084: new Projective_Single_Severe_Conditioned_Silent_Corrector(Norm,H,dH);
! 1085: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
! 1086:
! 1087: begin
! 1088: Corrector(s,pivot,dist_sols,c,fail);
! 1089: end Projective_Multiple_Severe_Conditioned_Silent_Corrector;
! 1090:
! 1091: procedure Projective_Multiple_Severe_Conditioned_Reporting_Corrector
! 1092: ( file : in file_type; s : in out Solu_Info_Array;
! 1093: pivot : in out natural; dist_sols : in double_float;
! 1094: c : in Corr_Pars; fail : out boolean ) is
! 1095:
! 1096: procedure Single_Corrector is
! 1097: new Projective_Single_Severe_Conditioned_Reporting_Corrector(Norm,H,dH);
! 1098: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
! 1099:
! 1100: begin
! 1101: Corrector(file,s,pivot,dist_sols,c,fail);
! 1102: end Projective_Multiple_Severe_Conditioned_Reporting_Corrector;
! 1103:
! 1104: end Correctors;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to correctors.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Continuation