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