Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/standard_floating_qr_decomposition.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 2: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
! 3:
! 4: package body Standard_Floating_QR_Decomposition 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: function dmax1 ( a,b : double_float ) return double_float is
! 20:
! 21: -- DESCRIPTION : returns the maximum of a and b.
! 22:
! 23: begin
! 24: if a >= b
! 25: then return a;
! 26: else return b;
! 27: end if;
! 28: end dmax1;
! 29:
! 30: function dsign ( a,b : double_float ) return double_float is
! 31:
! 32: -- DESCRIPTION : returns the absolute value of a, set to the same sign as b.
! 33:
! 34: begin
! 35: if b >= 0.0
! 36: then return abs(a);
! 37: else return -abs(a);
! 38: end if;
! 39: end dsign;
! 40:
! 41: procedure dswap ( a : in out Standard_Floating_Matrices.Matrix;
! 42: k1,k2 : in natural ) is
! 43:
! 44: -- DESCRIPTION :
! 45: -- Swaps the columns k1 and k2 in the matrix a.
! 46:
! 47: tmp : double_float;
! 48:
! 49: begin
! 50: for i in a'range(1) loop
! 51: tmp := a(i,k1); a(i,k1) := a(i,k2); a(i,k2) := tmp;
! 52: end loop;
! 53: end dswap;
! 54:
! 55: function dnrm2 ( a : Standard_Floating_Matrices.Matrix; row,col : natural )
! 56: return double_float is
! 57:
! 58: -- DESCRIPTION :
! 59: -- Computes the 2-norm of the vector in the column col of the matrix,
! 60: -- starting at the given row.
! 61:
! 62: sum : double_float := 0.0;
! 63:
! 64: begin
! 65: for i in row..a'last(1) loop
! 66: sum := sum + a(i,col)*a(i,col);
! 67: end loop;
! 68: return SQRT(sum);
! 69: end dnrm2;
! 70:
! 71: function ddot ( a : Standard_Floating_Matrices.Matrix; row,k1,k2 : natural )
! 72: return double_float is
! 73:
! 74: -- DESCRIPTION :
! 75: -- Returns the inner product of the vectors in the columns k1 and k2,
! 76: -- starting at the given row.
! 77:
! 78: res : double_float := 0.0;
! 79:
! 80: begin
! 81: for i in row..a'last(1) loop
! 82: res := res + a(i,k1)*a(i,k2);
! 83: end loop;
! 84: return res;
! 85: end ddot;
! 86:
! 87: procedure daxpy ( a : in out Standard_Floating_Matrices.Matrix;
! 88: f : in double_float; row,k1,k2 : in integer ) is
! 89:
! 90: -- DESCRIPTION :
! 91: -- The column k2 is added with f times the column k1, starting at row.
! 92:
! 93: begin
! 94: for i in row..a'last(1) loop
! 95: a(i,k2) := a(i,k2) + f*a(i,k1);
! 96: end loop;
! 97: end daxpy;
! 98:
! 99: procedure dscal ( a : in out Standard_Floating_Matrices.Matrix;
! 100: f : in double_float; row,col : in integer ) is
! 101:
! 102: -- DESCRIPTION :
! 103: -- Multiplies the column col of the matrix with f, starting at row.
! 104:
! 105: begin
! 106: for i in row..a'last(1) loop
! 107: a(i,col) := f*a(i,col);
! 108: end loop;
! 109: end dscal;
! 110:
! 111: procedure QRD ( x : in out Standard_Floating_Matrices.Matrix;
! 112: qraux : in out Standard_Floating_Vectors.Vector;
! 113: jpvt : in out Standard_Natural_Vectors.Vector;
! 114: piv : in boolean ) is
! 115:
! 116: n : constant natural := x'length(1); -- number of rows
! 117: p : constant natural := x'length(2); -- number of columns
! 118: work : Standard_Floating_Vectors.Vector(x'range(2));
! 119: jj,jp,lp1,lup,maxj,pl,pu : integer;
! 120: maxnrm,tt,nrmxl,t : double_float;
! 121: negj,swapj : boolean;
! 122:
! 123: begin
! 124: pl := 1; pu := 0;
! 125: if piv
! 126: then for j in x'range(2) loop -- rearrange columns according to jpvt
! 127: swapj := (jpvt(j) > 0);
! 128: negj := (jpvt(j) < 0);
! 129: jpvt(j) := j;
! 130: if negj
! 131: then jpvt(j) := -j;
! 132: end if;
! 133: if (swapj and then (j /= pl))
! 134: then dswap(x,pl,j);
! 135: jpvt(j) := jpvt(pl);
! 136: jpvt(pl) := j;
! 137: pl := pl + 1;
! 138: end if;
! 139: end loop;
! 140: pu := p;
! 141: for j in 1..p loop
! 142: jj := p - j + 1;
! 143: if jpvt(jj) < 0
! 144: then jpvt(jj) := -jpvt(jj);
! 145: if jj /= pu
! 146: then dswap(x,pu,jj);
! 147: jp := jpvt(pu);
! 148: jpvt(pu) := jpvt(jj);
! 149: jpvt(jj) := jp;
! 150: end if;
! 151: pu := pu - 1;
! 152: end if;
! 153: end loop;
! 154: end if;
! 155: for j in pl..pu loop -- compute norms of the free columns
! 156: qraux(j) := dnrm2(x,1,j);
! 157: work(j) := qraux(j);
! 158: end loop;
! 159: lup := min0(n,p); -- perform the householder reduction of x
! 160: for l in 1..lup loop
! 161: if (l >= pl and l < pu)
! 162: then maxnrm := 0.0; -- locate column with largest norm and
! 163: maxj := l; -- bring it in the pivot position
! 164: for j in l..pu loop
! 165: if qraux(j) > maxnrm
! 166: then maxnrm := qraux(j);
! 167: maxj := j;
! 168: end if;
! 169: end loop;
! 170: if maxj /= l
! 171: then dswap(x,l,maxj);
! 172: qraux(maxj) := qraux(l);
! 173: work(maxj) := work(l);
! 174: jp := jpvt(maxj);
! 175: jpvt(maxj) := jpvt(l);
! 176: jpvt(l) := jp;
! 177: end if;
! 178: end if;
! 179: qraux(l) := 0.0;
! 180: if l /= n
! 181: then nrmxl := dnrm2(x,l,l); -- householder transformation for column l
! 182: if nrmxl /= 0.0
! 183: then if x(l,l) /= 0.0
! 184: then nrmxl := dsign(nrmxl,x(l,l));
! 185: end if;
! 186: dscal(x,1.0/nrmxl,l,l);
! 187: x(l,l) := 1.0 + x(l,l);
! 188: lp1 := l + 1; -- apply the transformation to the remaining
! 189: for j in lp1..p loop -- columns, updating the norms
! 190: t := -ddot(x,l,l,j)/x(l,l);
! 191: daxpy(x,t,l,l,j);
! 192: if (j >= pl) and (j <= pu) and (qraux(j) /= 0.0)
! 193: then tt := 1.0 - (abs(x(l,j))/qraux(j))**2;
! 194: tt := dmax1(tt,0.0);
! 195: t := tt;
! 196: tt := 1.0 + 0.05*tt*(qraux(j)/work(j))**2;
! 197: if tt /= 1.0
! 198: then qraux(j) := qraux(j)*SQRT(t);
! 199: else qraux(j) := dnrm2(x,l+1,j);
! 200: work(j) := qraux(j);
! 201: end if;
! 202: end if;
! 203: end loop;
! 204: qraux(l) := x(l,l); -- save the transformation
! 205: x(l,l) := -nrmxl;
! 206: end if;
! 207: end if;
! 208: end loop;
! 209: end QRD;
! 210:
! 211: procedure Permute_Columns ( x : in out Standard_Floating_Matrices.Matrix;
! 212: jpvt : in Standard_Natural_Vectors.Vector ) is
! 213:
! 214: res : Standard_Floating_Matrices.Matrix(x'range(1),x'range(2));
! 215:
! 216: begin
! 217: for k in jpvt'range loop
! 218: for i in res'range(1) loop
! 219: res(i,k) := x(i,jpvt(k));
! 220: end loop;
! 221: end loop;
! 222: x := res;
! 223: end Permute_Columns;
! 224:
! 225: procedure Permute ( x : in out Standard_Floating_Vectors.Vector;
! 226: jpvt : in Standard_Natural_Vectors.Vector ) is
! 227:
! 228: res : Standard_Floating_Vectors.Vector(x'range);
! 229:
! 230: begin
! 231: for k in jpvt'range loop
! 232: res(k) := x(jpvt(k));
! 233: end loop;
! 234: x := res;
! 235: end Permute;
! 236:
! 237: procedure Basis ( qr : in out Standard_Floating_Matrices.Matrix;
! 238: x : in Standard_Floating_Matrices.Matrix ) is
! 239:
! 240: sum : double_float;
! 241: wrk : Standard_Floating_Vectors.Vector(qr'range(2));
! 242:
! 243: begin
! 244: for j in x'range(2) loop -- compute jth column of q
! 245: for i in qr'range(1) loop
! 246: sum := x(i,j);
! 247: for k in qr'first(2)..(j-1) loop
! 248: sum := sum - qr(i,k)*qr(k,j);
! 249: end loop;
! 250: wrk(i) := sum/qr(j,j);
! 251: end loop;
! 252: for i in qr'range(1) loop
! 253: qr(i,j) := wrk(i);
! 254: end loop;
! 255: end loop;
! 256: end Basis;
! 257:
! 258: end Standard_Floating_QR_Decomposition;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to standard_floating_qr_decomposition.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