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