![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to standard_complex_least_squares.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 standard_complex_least_squares.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices |
1.1 ! maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 2: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 3:
! 4: package body Standard_Complex_Least_Squares is
! 5:
! 6: -- AUXILIARIES :
! 7:
! 8: function min0 ( a,b : integer ) return integer is
! 9:
! 10: -- DESCRIPTION : returns the minimum of a and b.
! 11:
! 12: begin
! 13: if a <= b
! 14: then return a;
! 15: else return b;
! 16: end if;
! 17: end min0;
! 18:
! 19: procedure zcopy ( n,start : in natural;
! 20: x : in Standard_Complex_Vectors.Vector;
! 21: y : out Standard_Complex_Vectors.Vector ) is
! 22:
! 23: -- DESCRIPTION :
! 24: -- Copies the n entries of x to y, beginning at start.
! 25:
! 26: begin
! 27: for i in start..start+n-1 loop
! 28: y(i) := x(i);
! 29: end loop;
! 30: end zcopy;
! 31:
! 32: function zdotc ( row : natural;
! 33: x : Standard_Complex_Matrices.Matrix;
! 34: y : Standard_Complex_Vectors.Vector )
! 35: return Complex_Number is
! 36:
! 37: -- DESCRIPTION :
! 38: -- Dot product of two vectors : x(row..x'last(1),row)*y(row..y'last).
! 39:
! 40: res : Complex_Number := Create(0.0);
! 41:
! 42: begin
! 43: for i in row..y'last loop
! 44: res := res + Conjugate(x(i,row))*y(i);
! 45: end loop;
! 46: return res;
! 47: end zdotc;
! 48:
! 49: procedure zaxpy ( n,row,col : in natural; f : in Complex_Number;
! 50: x : in Standard_Complex_Matrices.Matrix;
! 51: y : in out Standard_Complex_Vectors.Vector ) is
! 52:
! 53: -- DESCRIPTION :
! 54: -- y(i) := y(i) + f*x(i,col) for i in row..row+n-1.
! 55:
! 56: begin
! 57: for i in row..row+n-1 loop
! 58: y(i) := y(i) + f*x(i,col);
! 59: end loop;
! 60: end zaxpy;
! 61:
! 62: procedure QRLS ( x : in out Standard_Complex_Matrices.Matrix;
! 63: ldx,n,k : in natural;
! 64: qraux,y : in Standard_Complex_Vectors.Vector;
! 65: qy,qty,b,rsd,xb : out Standard_Complex_Vectors.Vector;
! 66: job : in natural; info : out natural ) is
! 67:
! 68: cb,cqy,cqty,cr,cxb : boolean;
! 69: jj,ju,kp1 : integer;
! 70: t,temp : Complex_Number;
! 71:
! 72: begin
! 73: info := 0; -- set info flag
! 74: cqy := (job/10000 /= 0); -- determine what to compute
! 75: cqty := (job mod 10000 /= 0);
! 76: cb := ((job mod 1000)/100 /= 0);
! 77: cr := ((job mod 100)/10 /= 0);
! 78: cxb := ((job mod 10) /= 0);
! 79: ju := min0(k,n-1);
! 80: if ju = 0 -- special action when n = 1
! 81: then if cqy then qy(1) := y(1); end if;
! 82: if cqty then qty(1) := y(1); end if;
! 83: if cxb then xb(1) := y(1); end if;
! 84: if cb
! 85: then if AbsVal(x(1,1)) = 0.0
! 86: then info := 1;
! 87: else b(1) := y(1)/x(1,1);
! 88: end if;
! 89: end if;
! 90: if cr then rsd(1) := Create(0.0); end if;
! 91: return;
! 92: end if;
! 93: if cqy -- set up to compute qy or qty
! 94: then zcopy(n,y'first,y,qy);
! 95: end if;
! 96: if cqty
! 97: then zcopy(n,y'first,y,qty);
! 98: end if;
! 99: if cqy -- compute qy
! 100: then for j in 1..ju loop
! 101: jj := ju - j + 1;
! 102: if AbsVal(qraux(jj)) /= 0.0
! 103: then temp := x(jj,jj);
! 104: x(jj,jj) := qraux(jj);
! 105: t := -zdotc(jj,x,qy)/x(jj,jj);
! 106: zaxpy(n-jj+1,jj,jj,t,x,qy);
! 107: x(jj,jj) := temp;
! 108: end if;
! 109: end loop;
! 110: end if;
! 111: if cqty -- compute trans(q)*y
! 112: then for j in 1..ju loop
! 113: if AbsVal(qraux(j)) /= 0.0
! 114: then temp := x(j,j);
! 115: x(j,j) := qraux(j);
! 116: t := -zdotc(j,x,qty)/x(j,j);
! 117: zaxpy(n-j+1,j,j,t,x,qty);
! 118: x(j,j) := temp;
! 119: end if;
! 120: end loop;
! 121: end if;
! 122: if cb -- set up to compute b,rsd, or xb
! 123: then zcopy(k,qty'first,qty,b);
! 124: end if;
! 125: kp1 := k + 1;
! 126: if cxb then zcopy(k,qty'first,qty,xb); end if;
! 127: if (cr and (k < n))
! 128: then zcopy(n-k,kp1,qty,rsd);
! 129: end if;
! 130: if (cxb and (kp1 <= n))
! 131: then for i in kp1..n loop
! 132: xb(i) := Create(0.0);
! 133: end loop;
! 134: end if;
! 135: if cr
! 136: then for i in 1..k loop
! 137: rsd(i) := Create(0.0);
! 138: end loop;
! 139: end if;
! 140: if cb -- compute b
! 141: then for j in 1..k loop
! 142: jj := k - j + 1;
! 143: if AbsVal(x(jj,jj)) = 0.0
! 144: then info := jj; exit;
! 145: end if;
! 146: b(jj) := b(jj)/x(jj,jj);
! 147: if jj /= 1
! 148: then t := -b(jj);
! 149: zaxpy(jj-1,1,jj,t,x,b);
! 150: end if;
! 151: end loop;
! 152: end if;
! 153: if cr or cxb -- compute rsd or xb as requested
! 154: then for j in 1..ju loop
! 155: jj := ju - j + 1;
! 156: if AbsVal(qraux(jj)) /= 0.0
! 157: then temp := x(jj,jj);
! 158: x(jj,jj) := qraux(jj);
! 159: if cr
! 160: then t := -zdotc(jj,x,rsd)/x(jj,jj);
! 161: zaxpy(n-jj+1,jj,jj,t,x,rsd);
! 162: end if;
! 163: if cxb
! 164: then t := -zdotc(jj,x,xb)/x(jj,jj);
! 165: zaxpy(n-jj+1,jj,jj,t,x,xb);
! 166: end if;
! 167: x(jj,jj) := temp;
! 168: end if;
! 169: end loop;
! 170: end if;
! 171: end QRLS;
! 172:
! 173: end Standard_Complex_Least_Squares;