![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to multprec_complex_linear_solvers.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices|
Return to multprec_complex_linear_solvers.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices |
1.1 ! maekawa 1: with Multprec_Complex_Numbers; use Multprec_Complex_Numbers;
! 2:
! 3: package body Multprec_Complex_Linear_Solvers is
! 4:
! 5: -- AUXLILIARIES :
! 6:
! 7: function dconjg ( x : Complex_Number ) return Complex_Number is
! 8:
! 9: -- DESCRIPTION :
! 10: -- Returns the complex conjugate of x.
! 11:
! 12: res : Complex_Number;
! 13: re : Floating_Number := REAL_PART(x);
! 14: im : Floating_Number := IMAG_PART(x);
! 15:
! 16: begin
! 17: Min(im);
! 18: res := Create(re,im);
! 19: Clear(re); Clear(im);
! 20: return res;
! 21: end dconjg;
! 22:
! 23: function csign ( x,y : Complex_Number ) return Complex_Number is
! 24:
! 25: -- DESCRIPTION :
! 26: -- Returns |x|*y/|y|.
! 27:
! 28: res : Complex_Number;
! 29: fltacc : Floating_Number;
! 30: cmpacc : Complex_Number;
! 31:
! 32: begin
! 33: fltacc := AbsVal(x);
! 34: res := Create(fltacc);
! 35: Clear(fltacc);
! 36: Mul(res,y);
! 37: fltacc := AbsVal(y);
! 38: cmpacc := Create(fltacc);
! 39: Div(res,cmpacc);
! 40: Clear(fltacc);
! 41: Clear(cmpacc);
! 42: return res;
! 43: end csign;
! 44:
! 45: function Inverse_Abs_Sum ( z : Vector ) return Floating_Number is
! 46:
! 47: -- DESCRIPTION :
! 48: -- Returns the reciprocal of the sum of the absolute values in z.
! 49:
! 50: res,sum,acc : Floating_Number;
! 51:
! 52: begin
! 53: sum := Create(0);
! 54: for i in z'range loop
! 55: acc := AbsVal(z(i));
! 56: Add(sum,acc);
! 57: Clear(acc);
! 58: end loop;
! 59: res := Create(1);
! 60: Div(res,sum);
! 61: Clear(sum);
! 62: return res;
! 63: end Inverse_Abs_Sum;
! 64:
! 65: -- TARGET ROUTINES :
! 66:
! 67: procedure Scale ( a : in out Matrix; b : in out Vector ) is
! 68:
! 69: fac : Complex_Number;
! 70:
! 71: function Maximum ( a : in Matrix; i : in integer ) return Complex_Number is
! 72:
! 73: res : integer := a'first(2);
! 74: max : Floating_Number := AbsVal(a(i,res));
! 75: tmp : Floating_Number;
! 76:
! 77: begin
! 78: for j in a'first(2)+1..a'last(2) loop
! 79: tmp := AbsVal(a(i,j));
! 80: if tmp > max
! 81: then max := tmp; res := j;
! 82: end if;
! 83: end loop;
! 84: return a(i,res);
! 85: end Maximum;
! 86:
! 87: procedure Divide ( a : in out Matrix; b : in out Vector;
! 88: i : in integer; fac : in Complex_Number ) is
! 89: begin
! 90: for j in a'range(2) loop
! 91: a(i,j) := a(i,j)/fac;
! 92: end loop;
! 93: b(i) := b(i)/fac;
! 94: end Divide;
! 95:
! 96: begin
! 97: for i in a'range(1) loop
! 98: fac := Maximum(a,i);
! 99: Divide(a,b,i,fac);
! 100: end loop;
! 101: end Scale;
! 102:
! 103: -- TARGET ROUTINES :
! 104:
! 105: procedure lufac ( a : in out Matrix; n : in integer;
! 106: ipvt : out Standard_Natural_Vectors.Vector;
! 107: info : out natural ) is
! 108:
! 109: kp1,l,nm1 : integer;
! 110: smax,fltacc : Floating_Number;
! 111: temp,cmpacc : Complex_Number;
! 112:
! 113: begin
! 114: info := 0;
! 115: nm1 := n - 1;
! 116: if nm1 >= 1
! 117: then for k in 1..nm1 loop
! 118: kp1 := k + 1;
! 119: l := k;
! 120: smax := AbsVal(a(k,k)); -- find the pivot index l
! 121: for i in kp1..n loop
! 122: fltacc := AbsVal(a(i,k));
! 123: if fltacc > smax
! 124: then l := i;
! 125: Copy(fltacc,smax);
! 126: end if;
! 127: Clear(fltacc);
! 128: end loop;
! 129: ipvt(k) := l;
! 130: if Equal(smax,0.0) -- this column is already triangularized
! 131: then info := k;
! 132: else if l /= k -- interchange if necessary
! 133: then Copy(a(l,k),temp);
! 134: Copy(a(k,k),a(l,k));
! 135: Copy(temp,a(k,k)); Clear(temp);
! 136: end if;
! 137: cmpacc := Create(1); -- compute multipliers
! 138: Div(cmpacc,a(k,k));
! 139: Min(cmpacc);
! 140: for i in kp1..n loop
! 141: Mul(a(i,k),cmpacc);
! 142: end loop;
! 143: Clear(cmpacc);
! 144: for j in kp1..n loop -- row elimination
! 145: Copy(a(l,j),temp);
! 146: if l /= k
! 147: then Copy(a(k,j),a(l,j));
! 148: Copy(temp,a(k,j));
! 149: end if;
! 150: for i in kp1..n loop
! 151: cmpacc := temp*a(i,k);
! 152: Add(a(i,j),cmpacc);
! 153: Clear(cmpacc);
! 154: end loop;
! 155: end loop;
! 156: Clear(temp);
! 157: end if;
! 158: Clear(smax);
! 159: end loop;
! 160: end if;
! 161: ipvt(n) := n;
! 162: fltacc := AbsVal(a(n,n));
! 163: if Equal(fltacc,0.0)
! 164: then info := n;
! 165: end if;
! 166: Clear(fltacc);
! 167: end lufac;
! 168:
! 169: procedure lufco ( a : in out Matrix; n : in integer;
! 170: ipvt : out Standard_Natural_Vectors.Vector;
! 171: rcond : out Floating_Number ) is
! 172:
! 173: -- NOTE :
! 174: -- rcond = 1/(norm(a)*(estimate of norm(inverse(a))))
! 175: -- estimate = norm(z)/norm(y) where a*z = y and ctrans(a)*y = e.
! 176: -- ctrans(a) is the conjugate transpose of a.
! 177: -- The components of e are chosen to cause maximum local
! 178: -- growth in teh elements of w where ctrans(u)*w = e.
! 179: -- The vectors are frequently rescaled to avoid overflow.
! 180:
! 181: z : Multprec_Complex_Vectors.Vector(1..n);
! 182: info,kb,kp1,l : integer;
! 183: s,sm,sum,anorm,ynorm,fltacc1,fltacc2 : Floating_Number;
! 184: ek,t,wk,wkm,cmpacc1,cmpacc2 : Complex_Number;
! 185: ipvtt : Standard_Natural_Vectors.Vector(1..n);
! 186:
! 187: begin
! 188: anorm := Create(0); -- compute 1-norm of a
! 189: for j in 1..n loop
! 190: sum := Create(0);
! 191: for i in 1..n loop
! 192: fltacc1 := AbsVal(a(i,j));
! 193: Add(sum,fltacc1);
! 194: Clear(fltacc1);
! 195: end loop;
! 196: if sum > anorm
! 197: then Copy(sum,anorm);
! 198: end if;
! 199: Clear(sum);
! 200: end loop;
! 201: lufac(a,n,ipvtt,info); -- factor
! 202: ipvt := ipvtt;
! 203: ek := Create(1); -- solve ctrans(u)*w = e
! 204: for j in 1..n loop
! 205: z(j) := Create(0);
! 206: end loop;
! 207: for k in 1..n loop
! 208: fltacc1 := AbsVal(z(k));
! 209: if not Equal(fltacc1,0.0)
! 210: then cmpacc1 := -z(k);
! 211: Clear(ek);
! 212: ek := csign(ek,cmpacc1);
! 213: Clear(cmpacc1);
! 214: end if;
! 215: Clear(fltacc1);
! 216: cmpacc1 := ek - z(k);
! 217: fltacc1 := AbsVal(cmpacc1);
! 218: Clear(cmpacc1);
! 219: fltacc2 := AbsVal(a(k,k));
! 220: if fltacc1 > fltacc2
! 221: then s := fltacc2/fltacc1;
! 222: cmpacc1 := Create(s);
! 223: Mul(z,cmpacc1);
! 224: Mul(ek,cmpacc1);
! 225: Clear(cmpacc1);
! 226: Clear(s);
! 227: end if;
! 228: Clear(fltacc1); Clear(fltacc2);
! 229: Clear(wk); wk := ek - z(k);
! 230: Clear(wkm); wkm := ek + z(k);
! 231: Min(wkm);
! 232: s := AbsVal(wk);
! 233: sm := AbsVal(wkm);
! 234: fltacc1 := AbsVal(a(k,k));
! 235: if Equal(fltacc1,0.0)
! 236: then Clear(wk); wk := Create(1);
! 237: Clear(wkm); wkm := Create(1);
! 238: else cmpacc1 := dconjg(a(k,k));
! 239: Div(wk,cmpacc1);
! 240: Div(wkm,cmpacc1);
! 241: Clear(cmpacc1);
! 242: end if;
! 243: Clear(fltacc1);
! 244: kp1 := k + 1;
! 245: if kp1 <= n
! 246: then for j in kp1..n loop
! 247: cmpacc2 := dconjg(a(k,j));
! 248: cmpacc1 := wkm*cmpacc2;
! 249: Add(cmpacc1,z(j));
! 250: fltacc1 := AbsVal(cmpacc1);
! 251: Add(sm,fltacc1);
! 252: Clear(fltacc1);
! 253: Clear(cmpacc1);
! 254: cmpacc1 := wk*cmpacc2;
! 255: Add(z(j),cmpacc1);
! 256: Clear(cmpacc1);
! 257: Clear(cmpacc2);
! 258: fltacc1 := AbsVal(z(j));
! 259: Add(s,fltacc1);
! 260: Clear(fltacc1);
! 261: end loop;
! 262: if s < sm
! 263: then t := wkm - wk;
! 264: Copy(wkm,wk);
! 265: for j in kp1..n loop
! 266: cmpacc1 := t*dconjg(a(k,j));
! 267: Add(z(j),cmpacc1);
! 268: Clear(cmpacc1);
! 269: end loop;
! 270: Clear(t);
! 271: end if;
! 272: end if;
! 273: Copy(wk,z(k));
! 274: end loop;
! 275: s := Inverse_Abs_Sum(z);
! 276: cmpacc1 := Create(s);
! 277: Mul(z,cmpacc1);
! 278: Clear(cmpacc1);
! 279: Clear(s);
! 280: for k in 1..n loop -- solve ctrans(l)*y = w
! 281: kb := n+1-k;
! 282: if kb < n
! 283: then t := Create(0);
! 284: for i in (kb+1)..n loop
! 285: cmpacc1 := dconjg(a(i,kb))*z(i);
! 286: Add(t,cmpacc1);
! 287: Clear(cmpacc1);
! 288: end loop;
! 289: Add(z(kb),t);
! 290: end if;
! 291: fltacc1 := AbsVal(z(kb));
! 292: if fltacc1 > 1.0
! 293: then s := Create(1);
! 294: Div(s,fltacc1);
! 295: cmpacc1 := Create(s);
! 296: Mul(z,cmpacc1);
! 297: Clear(cmpacc1);
! 298: Clear(s);
! 299: end if;
! 300: Clear(fltacc1);
! 301: l := ipvtt(kb);
! 302: if l /= kb
! 303: then Copy(z(l),t);
! 304: Copy(z(kb),z(l));
! 305: Copy(t,z(kb));
! 306: Clear(t);
! 307: end if;
! 308: end loop;
! 309: s := Inverse_Abs_Sum(z);
! 310: cmpacc1 := Create(s);
! 311: Clear(s);
! 312: Mul(z,cmpacc1);
! 313: Clear(cmpacc1);
! 314: ynorm := Create(1);
! 315: for k in 1..n loop -- solve l*v = y
! 316: l := ipvtt(k);
! 317: if l /= k
! 318: then Copy(z(l),t);
! 319: Copy(z(k),z(l));
! 320: Copy(t,z(k));
! 321: else Copy(z(l),t);
! 322: end if;
! 323: if k < n
! 324: then for i in (k+1)..n loop
! 325: cmpacc1 := t*a(i,k);
! 326: Add(z(i),cmpacc1);
! 327: Clear(cmpacc1);
! 328: end loop;
! 329: end if;
! 330: Clear(t);
! 331: fltacc1 := AbsVal(z(k));
! 332: if fltacc1 > 1.0
! 333: then s := Create(1);
! 334: Div(s,fltacc1);
! 335: cmpacc1 := Create(s);
! 336: Mul(z,cmpacc1);
! 337: Clear(cmpacc1);
! 338: Mul(ynorm,s);
! 339: Clear(s);
! 340: end if;
! 341: Clear(fltacc1);
! 342: end loop;
! 343: s := Inverse_Abs_Sum(z);
! 344: cmpacc1 := Create(s);
! 345: Mul(z,cmpacc1);
! 346: Clear(cmpacc1);
! 347: Mul(ynorm,s);
! 348: Clear(s);
! 349: for k in 1..n loop -- solve u*z = v
! 350: kb := n+1-k;
! 351: fltacc1 := AbsVal(z(kb));
! 352: fltacc2 := AbsVal(a(kb,kb));
! 353: if fltacc1 > fltacc2
! 354: then s := fltacc2/fltacc1;
! 355: cmpacc1 := Create(s);
! 356: Mul(z,cmpacc1);
! 357: Mul(ynorm,s);
! 358: Clear(s);
! 359: end if;
! 360: if Equal(fltacc2,0.0)
! 361: then Clear(z(kb));
! 362: z(kb) := Create(1);
! 363: else Div(z(kb),a(kb,kb));
! 364: end if;
! 365: Clear(fltacc1);
! 366: Clear(fltacc2);
! 367: t := -z(kb);
! 368: for i in 1..(kb-1) loop
! 369: cmpacc1 := t*a(i,kb);
! 370: Add(z(i),cmpacc1);
! 371: Clear(cmpacc1);
! 372: end loop;
! 373: Clear(t);
! 374: end loop;
! 375: s := Inverse_Abs_Sum(z); -- make znorm = 1.0
! 376: cmpacc1 := Create(s);
! 377: Mul(z,cmpacc1);
! 378: Clear(cmpacc1);
! 379: Mul(ynorm,s);
! 380: if Equal(anorm,0.0)
! 381: then rcond := Create(0);
! 382: else rcond := ynorm/anorm;
! 383: end if;
! 384: Clear(ynorm); Clear(s); Clear(sm);
! 385: Clear(ek); Clear(t); Clear(wk); Clear(wkm);
! 386: Clear(z);
! 387: end lufco;
! 388:
! 389: procedure lusolve ( a : in Matrix; n : in integer;
! 390: ipvt : in Standard_Natural_Vectors.Vector;
! 391: b : in out Vector ) is
! 392:
! 393: l,nm1,kb : integer;
! 394: temp,acc : Complex_Number;
! 395:
! 396: begin
! 397: nm1 := n-1;
! 398: if nm1 >= 1 -- solve l*y = b
! 399: then for k in 1..nm1 loop
! 400: l := ipvt(k);
! 401: Copy(b(l),temp);
! 402: if l /= k
! 403: then Copy(b(k),b(l));
! 404: Copy(temp,b(k));
! 405: end if;
! 406: for i in (k+1)..n loop
! 407: acc := temp*a(i,k);
! 408: Add(b(i),acc);
! 409: Clear(acc);
! 410: end loop;
! 411: Clear(temp);
! 412: end loop;
! 413: end if;
! 414: for k in 1..n loop -- solve u*x = y
! 415: kb := n+1-k;
! 416: Div(b(kb),a(kb,kb));
! 417: temp := -b(kb);
! 418: for j in 1..(kb-1) loop
! 419: acc := temp*a(j,kb);
! 420: Add(b(j),acc);
! 421: Clear(acc);
! 422: end loop;
! 423: Clear(temp);
! 424: end loop;
! 425: end lusolve;
! 426:
! 427: procedure Triangulate ( a : in out Matrix; n,m : in integer ) is
! 428:
! 429: max,cbs : Floating_Number;
! 430: temp : Complex_Number;
! 431: pivot,k,kcolumn : integer;
! 432:
! 433: begin
! 434: k := 1;
! 435: kcolumn := 1;
! 436: while (k <= n) and (kcolumn <= m) loop
! 437: max := Create(0); -- find pivot
! 438: pivot := 0;
! 439: for l in k..n loop
! 440: cbs := AbsVal(a(l,kcolumn));
! 441: if cbs > max
! 442: then max := cbs;
! 443: pivot := l;
! 444: end if;
! 445: end loop;
! 446: if pivot = 0
! 447: then kcolumn := kcolumn + 1;
! 448: else if pivot /= k -- interchange if necessary
! 449: then for i in 1..m loop
! 450: temp := a(pivot,i);
! 451: a(pivot,i) := a(k,i);
! 452: a(k,i) := temp;
! 453: end loop;
! 454: end if;
! 455: for j in (kcolumn+1)..m loop -- triangulate a
! 456: a(k,j) := a(k,j) / a(k,kcolumn);
! 457: end loop;
! 458: a(k,kcolumn) := Create(1);
! 459: for i in (k+1)..n loop
! 460: for j in (kcolumn+1)..m loop
! 461: a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
! 462: end loop;
! 463: a(i,kcolumn) := Create(0);
! 464: end loop;
! 465: k := k + 1;
! 466: kcolumn := kcolumn + 1;
! 467: end if;
! 468: end loop;
! 469: end Triangulate;
! 470:
! 471: procedure Diagonalize ( a : in out Matrix; n,m : in integer ) is
! 472:
! 473: max : Floating_Number;
! 474: temp : Complex_Number;
! 475: pivot,k,kcolumn : integer;
! 476:
! 477: begin
! 478: k := 1;
! 479: kcolumn := 1;
! 480: while (k <= n) and (kcolumn <= m) loop
! 481: max := Create(0); -- find pivot
! 482: for l in k..n loop
! 483: if AbsVal(a(l,kcolumn)) > max
! 484: then max := AbsVal(a(l,kcolumn));
! 485: pivot := l;
! 486: end if;
! 487: end loop;
! 488: if Equal(max,0.0)
! 489: then kcolumn := kcolumn + 1;
! 490: else if pivot /= k -- interchange if necessary
! 491: then for i in 1..m loop
! 492: temp := a(pivot,i);
! 493: a(pivot,i) := a(k,i);
! 494: a(k,i) := temp;
! 495: end loop;
! 496: end if;
! 497: for j in (kcolumn+1)..m loop -- diagonalize a
! 498: a(k,j) := a(k,j) / a(k,kcolumn);
! 499: end loop;
! 500: a(k,kcolumn) := Create(1);
! 501: for i in 1..(k-1) loop
! 502: for j in (kcolumn+1)..m loop
! 503: a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
! 504: end loop;
! 505: end loop;
! 506: for i in (k+1)..n loop
! 507: for j in (kcolumn+1)..m loop
! 508: a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
! 509: end loop;
! 510: end loop;
! 511: for j in 1..(k-1) loop
! 512: a(j,kcolumn) := Create(0);
! 513: end loop;
! 514: for j in (k+1)..n loop
! 515: a(j,kcolumn) := Create(0);
! 516: end loop;
! 517: k := k + 1;
! 518: kcolumn := kcolumn + 1;
! 519: end if;
! 520: end loop;
! 521: end Diagonalize;
! 522:
! 523: end Multprec_Complex_Linear_Solvers;