1 { 2 This file is part of the Numlib package. 3 Copyright (c) 1986-2000 by 4 Kees van Ginneken, Wil Kortsmit and Loek van Reij of the 5 Computational centre of the Eindhoven University of Technology 6 7 FPC port Code by Marco van de Voort (marco@freepascal.org) 8 documentation by Michael van Canneyt (Michael@freepascal.org) 9 10 Determinants for different kinds of matrices (different with respect 11 to symmetry) 12 13 See the file COPYING.FPC, included in this distribution, 14 for details about the copyright. 15 16 This program is distributed in the hope that it will be useful, 17 but WITHOUT ANY WARRANTY; without even the implied warranty of 18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 19 20 **********************************************************************} 21 unit det; 22 23 interface 24 {$I DIRECT.INC} 25 26 uses typ; 27 28 {Generic determinant} 29 procedure detgen(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt); 30 31 {determinant symmetrical matrix} 32 procedure detgsy(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt); 33 34 {determinant of a symmetrical positive definitive matrix} 35 procedure detgpd(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt); 36 37 {determinant of a generic bandmatrix} 38 procedure detgba(n, l, r: ArbInt; var a, f: ArbFloat; var k, term:ArbInt); 39 40 {determinant of a symmetrical positive definitive bandmatrix} 41 procedure detgpb(n, l: ArbInt; var a, f: ArbFloat; var k, term:ArbInt); 42 43 {determinant of a tridiagonal matrix} 44 procedure detgtr(n: ArbInt; var l, d, u, f: ArbFloat; var k, term:ArbInt); 45 46 { moved to the TYP unit because of a bug in FPC 1.0.x FK 47 var og : ArbFloat absolute ogx; 48 bg : ArbFloat absolute bgx; 49 MaxExp : ArbInt absolute maxexpx; 50 } 51 52 implementation 53 54 uses mdt; 55 56 procedure detgen(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt); 57 58 var 59 kk, ind, ind1, ns, i : ArbInt; 60 u, ca : ArbFloat; 61 pa, acopy : ^arfloat1; 62 p : ^arint1; 63 begin 64 if (n<1) or (rwidth<1) then 65 begin 66 term:=3; exit 67 end; {wrong input} 68 pa:=@a; 69 ns:=n*sizeof(ArbFloat); 70 getmem(p, n*sizeof(ArbInt)); 71 getmem(acopy, n*ns); 72 ind:=1; ind1:=1; 73 for i:=1 to n do 74 begin 75 move(pa^[ind1], acopy^[ind], ns); 76 ind1:=ind1+rwidth; ind:=ind+n 77 end; {i} 78 mdtgen(n, n, acopy^[1], p^[1], ca, term); 79 if term=1 then 80 begin 81 f:=1; k:=0; kk:=1; ind:=1; 82 while (kk<=n) and (f<>0) do 83 begin 84 u:=acopy^[ind]; 85 while (u<>0) and (abs(u)<og) do 86 begin 87 u:=u/og; k:=k-maxexp 88 end; {underflow control} 89 while abs(u)>bg do 90 begin 91 u:=u/bg; k:=k+maxexp 92 end; {overflow control} 93 f:=f*u; 94 if p^[kk]<>kk then f:=-f; 95 while (f<>0) and (abs(f)<og) do 96 begin 97 f:=f/og; k:=k-maxexp 98 end; {underflow control} 99 while abs(f)>bg do 100 begin 101 f:=f/bg; k:=k+maxexp 102 end; {overflow control} 103 kk:=kk+1; ind:=ind+n+1 104 end; {kk} 105 end {term=1} 106 else {term=4} 107 begin 108 f:=0; k:=0; term:=1 109 end; 110 freemem(p, n*sizeof(ArbInt)); 111 freemem(acopy, n*ns) 112 end; {detgen} 113 114 procedure detgsy(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt); 115 116 var i, kk, ind, ind1, s : ArbInt; 117 u, ca : ArbFloat; 118 pa, acopy : ^arfloat1; 119 p : ^arint1; 120 q : ^arbool1; 121 begin 122 if (n<1) or (rwidth<1) then 123 begin 124 term:=3; exit 125 end; {wrong input} 126 pa:=@a; 127 getmem(p, n*sizeof(ArbInt)); 128 getmem(q, n*sizeof(boolean)); 129 s:=sizeof(ArbFloat); 130 getmem(acopy, n*n*s); 131 ind:=1; ind1:=1; 132 for i:=1 to n do 133 begin 134 move(pa^[ind1], acopy^[ind], i*s); 135 ind1:=ind1+rwidth; ind:=ind+n 136 end; {i} 137 mdtgsy(n, n, acopy^[1], p^[1], q^[1], ca, term); 138 if term=1 then 139 begin 140 f:=1; k:=0; kk:=1; ind:=1; 141 while (kk<=n) and (f<>0) do 142 begin 143 u:=acopy^[ind]; 144 while (u<>0) and (abs(u)<og) do 145 begin 146 u:=u/og; k:=k-maxexp 147 end; {underflow control} 148 while abs(u)>bg do 149 begin 150 u:=u/bg; k:=k+maxexp 151 end; {overflow control} 152 f:=f*u; 153 if q^[kk] then f:=-f; 154 while (f<>0) and (abs(f)<og) do 155 begin 156 f:=f/og; k:=k-maxexp 157 end; {underflow control} 158 while abs(f)>bg do 159 begin 160 f:=f/bg; k:=k+maxexp 161 end; {overflow control} 162 kk:=kk+1; ind:=ind+n+1 163 end; {kk} 164 end {term=1} 165 else {term=4} 166 begin 167 term:=1; f:=0; k:=0 168 end; 169 freemem(p, n*sizeof(ArbInt)); 170 freemem(q, n*sizeof(boolean)); 171 freemem(acopy, n*n*s) 172 end; {detgsy} 173 174 procedure detgpd(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt); 175 176 var 177 i, kk, ind, ind1, s : ArbInt; 178 u, ca : ArbFloat; 179 pa, acopy : ^arfloat1; 180 begin 181 if (n<1) or (rwidth<1) then 182 begin 183 term:=3; exit 184 end; {wrong input} 185 pa:=@a; 186 s:=sizeof(ArbFloat); 187 getmem(acopy, n*n*s); 188 ind:=1; ind1:=1; 189 for i:=1 to n do 190 begin 191 move(pa^[ind1], acopy^[ind], i*s); 192 ind1:=ind1+rwidth; ind:=ind+n 193 end; {i} 194 mdtgpd(n, n, acopy^[1], ca, term); 195 if term=1 then 196 begin 197 f:=1; k:=0; kk:=1; ind:=1; 198 while kk<=n do 199 begin 200 u:=sqr(acopy^[ind]); 201 while u < og do 202 begin 203 u:=u/og; k:=k-maxexp 204 end; {underflow control} 205 while u > bg do 206 begin 207 u:=u/bg; k:=k+maxexp 208 end; {overflow control} 209 f:=f*u; 210 while f < og do 211 begin 212 f:=f/og; k:=k-maxexp 213 end; {underflow control} 214 while f > bg do 215 begin 216 f:=f/bg; k:=k+maxexp 217 end; {overflow control} 218 kk:=kk+1; ind:=ind+n+1 219 end; {kk} 220 end; {term=1} 221 freemem(acopy, n*n*s) 222 end; {detgpd} 223 224 procedure detgba(n, l, r: ArbInt; var a, f: ArbFloat; 225 var k, term:ArbInt); 226 var 227 rwidth, s, ns, kk, ii, i, jj, ll : ArbInt; 228 u, ca : ArbFloat; 229 pa, l1, acopy : ^arfloat1; 230 p : ^arint1; 231 begin 232 if (n<1) or (l<0) or (r<0) or (l>n-1) or (r>n-1) then 233 begin 234 term:=3; exit 235 end; {wrong input} 236 pa:=@a; 237 s:=sizeof(ArbFloat); ns:=n*s; 238 ll:=l+r+1; 239 getmem(acopy, ll*ns); 240 getmem(l1, l*ns); 241 getmem(p, n*sizeof(ArbInt)); 242 jj:=1; ii:=1; 243 for i:=1 to n do 244 begin 245 if i <= l+1 then 246 begin 247 if i <= (n-r) then rwidth:=r+i else rwidth:=n 248 end else 249 if i <= (n-r) then rwidth:=ll else rwidth:=n-i+l+1; 250 if i > l then kk:=ii else kk:=ii+l-i+1; 251 move(pa^[jj], acopy^[kk], rwidth*s); 252 jj:=jj+rwidth; ii:=ii+ll; 253 end; {i} 254 mdtgba(n, l, r, ll, acopy^[1], l, l1^[1], p^[1], ca, term); 255 if term=1 then 256 begin 257 f:=1; k:=0; kk:=1; ii:=1; 258 while (kk<=n) and (f<>0) do 259 begin 260 u:=acopy^[ii]; 261 while (u<>0) and (abs(u)<og) do 262 begin 263 u:=u/og; k:=k-maxexp 264 end; {underflow control} 265 while abs(u)>bg do 266 begin 267 u:=u/bg; k:=k+maxexp 268 end; {overflow control} 269 f:=f*u; 270 if p^[kk]<>kk then f:=-f; 271 while (f<>0) and (abs(f)<og) do 272 begin 273 f:=f/og; k:=k-maxexp 274 end; {underflow control} 275 while abs(f)>bg do 276 begin 277 f:=f/bg; k:=k+maxexp 278 end; {overflow control} 279 kk:=kk+1; ii:=ii+ll 280 end; {kk} 281 end {term=1} 282 else {term=4} 283 begin 284 term:=1; f:=0; k:=0 285 end; 286 freemem(acopy, ll*ns); 287 freemem(l1, l*ns); 288 freemem(p, n*sizeof(ArbInt)) 289 end; {detgba} 290 291 procedure detgpb(n, l: ArbInt; var a, f: ArbFloat; var k, term: ArbInt); 292 293 var 294 rwidth, kk, ii, ns, ll, jj, i, s : ArbInt; 295 u, ca : ArbFloat; 296 pa, acopy : ^arfloat1; 297 begin 298 if (n<1) or (l<0) or (l>n-1) then 299 begin 300 term:=3; exit 301 end; {wrong input} 302 pa:=@a; 303 ll:=l+1; 304 s:=sizeof(ArbFloat); ns:=s*n; 305 getmem(acopy, ll*ns); 306 jj:=1; ii:=1; 307 for i:=1 to n do 308 begin 309 if i>l then rwidth:=ll else rwidth:=i; 310 move(pa^[jj], acopy^[ii+ll-rwidth], rwidth*s); 311 jj:=jj+rwidth; ii:=ii+ll 312 end; {i} 313 mdtgpb(n, l, ll, acopy^[1], ca, term); 314 if term=1 then 315 begin 316 f:=1; k:=0; kk:=1; ii:=ll; 317 while (kk<=n) do 318 begin 319 u:=sqr(acopy^[ii]); 320 while u < og do 321 begin 322 u:=u/og; k:=k-maxexp 323 end; {underflow control} 324 while u > bg do 325 begin 326 u:=u/bg; k:=k+maxexp 327 end; {overflow control} 328 f:=f*u; 329 while f < og do 330 begin 331 f:=f/og; k:=k-maxexp 332 end; {underflow control} 333 while f > bg do 334 begin 335 f:=f/bg; k:=k+maxexp 336 end; {overflow control} 337 kk:=kk+1; ii:=ii+ll 338 end; {kk} 339 end; {term=1} 340 freemem(acopy, ll*ns); 341 end; {detgpb} 342 343 procedure detgtr(n: ArbInt; var l, d, u, f: ArbFloat; var k, term:ArbInt); 344 345 var 346 ns, kk : ArbInt; 347 uu, ca : ArbFloat; 348 pl, pd, pu, l1, d1, u1, u2 : ^arfloat1; 349 p : ^arbool1; 350 begin 351 if n<1 then 352 begin 353 term:=3; exit 354 end; {wrong input} 355 pl:=@l; pd:=@d; pu:=@u; 356 ns:=n*sizeof(ArbFloat); 357 getmem(l1, ns); 358 getmem(d1, ns); 359 getmem(u1, ns); 360 getmem(u2, ns); 361 getmem(p, n*sizeof(boolean)); 362 mdtgtr(n, pl^[1], pd^[1], pu^[1], l1^[1], d1^[1], u1^[1], u2^[1], 363 p^[1], ca, term); 364 if term=1 then 365 begin 366 f:=1; k:=0; kk:=1; 367 while (kk<=n) and (f<>0) do 368 begin 369 if p^[kk] then f:=-f; 370 uu:=d1^[kk]; 371 while (uu<>0) and (abs(uu)<og) do 372 begin 373 uu:=uu/og; k:=k-maxexp 374 end; {underflow control} 375 while abs(uu)>bg do 376 begin 377 uu:=uu/bg; k:=k+maxexp 378 end; {overflow control} 379 f:=f*uu; 380 while (f<>0) and (abs(f)<og) do 381 begin 382 f:=f/og; k:=k-maxexp 383 end; {underflow control} 384 while abs(f)>bg do 385 begin 386 f:=f/bg; k:=k+maxexp 387 end; {overflow control} 388 kk:=kk+1 389 end; {kk} 390 end {term=1} 391 else {term=4} 392 begin 393 term:=1; f:=0; k:=0 394 end; 395 freemem(l1, ns); 396 freemem(d1, ns); 397 freemem(u1, ns); 398 freemem(u2, ns); 399 freemem(p, n*sizeof(boolean)); 400 end; {detgtr} 401 402 end. 403