Annotation of OpenXM_contrib/PHC/Ada/Continuation/vlprs_tables.adb, Revision 1.1
1.1 ! maekawa 1: package body vLpRs_Tables is
! 2:
! 3: -- I. The v-table :
! 4:
! 5: procedure v_pipe ( v : in out Vector; p : in Vector;
! 6: vr : in double_float ) is
! 7:
! 8: tmp0,tmp1 : double_float;
! 9:
! 10: begin
! 11: tmp0 := v(0);
! 12: v(0) := vr;
! 13: for i in 1..v'last loop
! 14: tmp1 := v(i);
! 15: v(i) := tmp0 - p(i-1)*v(i-1);
! 16: tmp0 := tmp1;
! 17: end loop;
! 18: end v_pipe;
! 19:
! 20: procedure v_pipe ( v,p : in Vector; vr : in double_float;
! 21: vrp : in out Vector ) is
! 22: begin
! 23: vrp(0) := vr;
! 24: for i in 1..v'last loop
! 25: vrp(i) := v(i-1) - p(i-1)*vrp(i-1);
! 26: end loop;
! 27: end v_pipe;
! 28:
! 29: -- II. The L-table :
! 30:
! 31: procedure L_pipe ( l : in out Vector; p : in Vector;
! 32: lr : in double_float ) is
! 33:
! 34: tmp0,tmp1 : double_float;
! 35:
! 36: begin
! 37: tmp0 := l(0);
! 38: l(0) := lr;
! 39: for i in 1..l'last loop
! 40: tmp1 := l(i);
! 41: l(i) := tmp0 - p(i-1)*l(i-1);
! 42: tmp0 := tmp1;
! 43: end loop;
! 44: end L_pipe;
! 45:
! 46: procedure L_pipe ( l,p : in Vector; lr : in double_float;
! 47: lrp : in out Vector ) is
! 48: begin
! 49: lrp(0) := lr;
! 50: for i in 1..l'last loop
! 51: lrp(i) := l(i-1) - p(i-1)*lrp(i-1);
! 52: end loop;
! 53: end L_pipe;
! 54:
! 55: -- The combined full version for v- and L-table :
! 56:
! 57: procedure vL_full ( s,l,v : in Vector; srp,dsp,p,lrp,vrp : out Vector;
! 58: rt1,rt2 : in out Matrix ) is
! 59:
! 60: srm1,tmpsrp,tmpdsp : Vector(1..s'last-1);
! 61: tmpp : Vector(0..s'last-1);
! 62: prev_lrp,new_lrp,prev_vrp,new_vrp : Vector(s'range);
! 63: k_last : integer;
! 64:
! 65: begin
! 66: srm1(srm1'first) := s(0);
! 67: for i in srm1'first+1..srm1'last loop
! 68: srm1(i) := s(0)*srm1(i-1);
! 69: end loop;
! 70: s_pipe(srm1,s(1),tmpsrp,tmpdsp); srm1 := tmpsrp;
! 71: for j in rt1'range(2) loop
! 72: rt1(1,j) := tmpdsp(j);
! 73: end loop;
! 74: tmpp(0) := 1.0;
! 75: prev_lrp(0) := l(0); new_lrp(0) := l(1);
! 76: prev_vrp(0) := v(0); new_vrp(0) := v(1);
! 77: new_vrp(1) := prev_vrp(0) - new_vrp(0);
! 78: new_lrp(1) := prev_lrp(0) - new_lrp(0);
! 79: prev_lrp(0..1) := new_lrp(0..1);
! 80: prev_vrp(0..1) := new_vrp(0..1);
! 81: for i in 2..s'last loop
! 82: s_pipe(srm1,s(i),tmpsrp,tmpdsp); srm1 := tmpsrp;
! 83: for j in rt2'range(2) loop
! 84: rt2(1,j) := tmpdsp(j);
! 85: end loop;
! 86: if i < s'last
! 87: then k_last := i; -- compute one additional row
! 88: else k_last := i-1;
! 89: end if;
! 90: for k in 1..k_last loop
! 91: tmpp(k) := rt1(k,k)/rt2(k,k);
! 92: for j in k+1..rt2'last(2) loop
! 93: rt2(k+1,j) := rt1(k,j) - tmpp(k)*rt2(k,j);
! 94: end loop;
! 95: end loop;
! 96: rt1 := rt2;
! 97: new_vrp(0) := v(i);
! 98: new_lrp(0) := l(i);
! 99: for k in 1..i loop
! 100: new_vrp(k) := prev_vrp(k-1) - tmpp(k-1)*new_vrp(k-1);
! 101: new_lrp(k) := prev_lrp(k-1) - tmpp(k-1)*new_lrp(k-1);
! 102: end loop;
! 103: prev_vrp(0..i) := new_vrp(0..i);
! 104: prev_lrp(0..i) := new_lrp(0..i);
! 105: end loop;
! 106: srp := tmpsrp;
! 107: dsp := tmpdsp;
! 108: p := tmpp;
! 109: lrp := new_lrp;
! 110: vrp := new_vrp;
! 111: end vL_full;
! 112:
! 113: -- III. The p-table :
! 114:
! 115: procedure p_full ( s : in Vector; srp,dsp,p : out Vector;
! 116: rt1,rt2 : in out Matrix ) is
! 117: begin
! 118: RR_full(s,srp,dsp,p,rt1,rt2);
! 119: end p_full;
! 120:
! 121: procedure p_pipe ( rt1,rt2 : in Matrix; p : out Vector ) is
! 122: begin
! 123: p(0) := 1.0;
! 124: for i in p'first+1..p'last loop
! 125: p(i) := rt1(i,i)/rt2(i,i);
! 126: end loop;
! 127: end p_pipe;
! 128:
! 129: -- IV. The R-table :
! 130:
! 131: procedure R_full ( s : in Vector; srp,dsp,p : out Vector;
! 132: rt1,rt2 : in out Matrix ) is
! 133:
! 134: srm1,tmpsrp,tmpdsp : Vector(1..s'last-1);
! 135: tmpp : Vector(0..s'last-1);
! 136:
! 137: begin
! 138: srm1(srm1'first) := s(0);
! 139: for i in srm1'first+1..srm1'last loop
! 140: srm1(i) := s(0)*srm1(i-1);
! 141: end loop;
! 142: s_pipe(srm1,s(1),tmpsrp,tmpdsp); srm1 := tmpsrp;
! 143: for j in tmpdsp'range loop
! 144: rt1(1,j) := tmpdsp(j);
! 145: end loop;
! 146: tmpp(0) := 1.0;
! 147: for i in 2..s'last loop
! 148: s_pipe(srm1,s(i),tmpsrp,tmpdsp); srm1 := tmpsrp;
! 149: for j in tmpdsp'range loop
! 150: rt2(1,j) := tmpdsp(j);
! 151: end loop;
! 152: for k in 1..i-1 loop
! 153: tmpp(k) := rt1(k,k)/rt2(k,k);
! 154: for j in k+1..rt2'last(2) loop
! 155: rt2(k+1,j) := rt1(k,j) - tmpp(k)*rt2(k,j);
! 156: end loop;
! 157: end loop;
! 158: rt1 := rt2;
! 159: end loop;
! 160: srp := tmpsrp;
! 161: dsp := tmpdsp;
! 162: p := tmpp;
! 163: end R_full;
! 164:
! 165: procedure RR_full ( s : in Vector; srp,dsp,p : out Vector;
! 166: rt1,rt2 : in out Matrix ) is
! 167:
! 168: srm1,tmpsrp,tmpdsp : Vector(1..s'last-1);
! 169: tmpp : Vector(0..s'last-1);
! 170: k_last,j_first : integer;
! 171:
! 172: begin
! 173: srm1(srm1'first) := s(0);
! 174: for i in srm1'first+1..srm1'last loop
! 175: srm1(i) := s(0)*srm1(i-1);
! 176: end loop;
! 177: s_pipe(srm1,s(1),tmpsrp,tmpdsp); srm1 := tmpsrp;
! 178: for j in tmpdsp'range loop
! 179: rt1(1,j) := tmpdsp(j);
! 180: end loop;
! 181: tmpp(0) := 1.0;
! 182: for i in 2..s'last loop
! 183: s_pipe(srm1,s(i),tmpsrp,tmpdsp); srm1 := tmpsrp;
! 184: for j in tmpdsp'range loop
! 185: rt2(1,j) := tmpdsp(j);
! 186: end loop;
! 187: if i < s'last
! 188: then k_last := i; -- compute one additional row
! 189: else k_last := i-1;
! 190: end if;
! 191: for k in 1..k_last loop
! 192: tmpp(k) := rt1(k,k)/rt2(k,k);
! 193: if k < rt2'last(2)
! 194: then j_first := k;
! 195: else j_first := k+1;
! 196: end if;
! 197: for j in j_first..rt2'last(2) loop -- start earlier with columns
! 198: rt2(k+1,j) := rt1(k,j) - tmpp(k)*rt2(k,j);
! 199: end loop;
! 200: end loop;
! 201: rt1 := rt2;
! 202: end loop;
! 203: srp := tmpsrp;
! 204: dsp := tmpdsp;
! 205: p := tmpp;
! 206: end RR_full;
! 207:
! 208: procedure R_pipe ( rt1 : in Matrix; s,p : in Vector; rt2 : in out Matrix ) is
! 209: begin
! 210: rt2(1,1) := s(1);
! 211: for j in 2..s'last loop
! 212: rt2(1,j) := s(j);
! 213: for i in 2..j loop
! 214: rt2(i,j) := rt1(i-1,j) - p(i-1)*rt2(i-1,j);
! 215: end loop;
! 216: end loop;
! 217: end R_pipe;
! 218:
! 219: procedure RR_pipe ( rt1 : in Matrix; s,p : in Vector; rt2 : in out Matrix ) is
! 220:
! 221: i_last : integer;
! 222:
! 223: begin
! 224: rt2(1,1) := s(1);
! 225: for j in 2..s'last loop
! 226: rt2(1,j) := s(j);
! 227: if j < rt2'last(1)
! 228: then i_last := j+1; -- compute one additional row
! 229: else i_last := j;
! 230: end if;
! 231: for i in 2..i_last loop
! 232: rt2(i,j) := rt1(i-1,j) - p(i-1)*rt2(i-1,j);
! 233: end loop;
! 234: end loop;
! 235: end RR_pipe;
! 236:
! 237: -- V. The s-table :
! 238:
! 239: procedure s_full ( s : in Vector; srp,dsp : out Vector ) is
! 240:
! 241: srm1,tmpsrp : Vector(1..s'last-1);
! 242:
! 243: begin
! 244: srm1(srm1'first) := s(0);
! 245: for i in srm1'first+1..srm1'last loop
! 246: srm1(i) := s(0)*srm1(i-1);
! 247: end loop;
! 248: for i in 1..s'last loop
! 249: s_pipe(srm1,s(i),tmpsrp,dsp);
! 250: srm1 := tmpsrp;
! 251: end loop;
! 252: srp := tmpsrp;
! 253: end s_full;
! 254:
! 255: procedure s_pipe ( srp : in out Vector; sr : in double_float;
! 256: dsp : out Vector ) is
! 257:
! 258: tmp : double_float := srp(1);
! 259: srpow : double_float := sr;
! 260:
! 261: begin
! 262: srp(1) := srpow;
! 263: dsp(1) := srpow - tmp;
! 264: for i in 2..srp'last loop
! 265: srpow := srpow*sr;
! 266: tmp := srp(i);
! 267: srp(i) := srpow;
! 268: dsp(i) := srpow - tmp;
! 269: end loop;
! 270: end s_pipe;
! 271:
! 272: procedure s_pipe ( sr1 : in Vector; sr : in double_float;
! 273: srp,dsp : out Vector ) is
! 274:
! 275: srpow : double_float := sr;
! 276:
! 277: begin
! 278: srp(1) := srpow;
! 279: dsp(1) := srpow - sr1(1);
! 280: for i in sr1'first+1..sr1'last loop
! 281: srpow := srpow*sr;
! 282: srp(i) := srpow;
! 283: dsp(i) := srpow - sr1(i);
! 284: end loop;
! 285: end s_pipe;
! 286:
! 287: end vLpRs_Tables;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to vlprs_tables.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