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 Unit was originally undocumented, but is probably an variant of DET. 11 Det accepts vectors as arguments, while MDT calculates determinants for 12 matrices. 13 14 Contrary to the other undocumented units, this unit is exported in the 15 DLL. 16 17 See the file COPYING.FPC, included in this distribution, 18 for details about the copyright. 19 20 This program is distributed in the hope that it will be useful, 21 but WITHOUT ANY WARRANTY; without even the implied warranty of 22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 23 24 **********************************************************************} 25 Unit mdt; 26 27 interface 28 {$I DIRECT.INC} 29 30 uses typ, dsl, omv; 31 32 Procedure mdtgen(n, rwidth: ArbInt; Var alu: ArbFloat; Var p: ArbInt; 33 Var ca:ArbFloat; Var term: ArbInt); 34 35 Procedure mdtgtr(n: ArbInt; Var l, d, u, l1, d1, u1, u2: ArbFloat; 36 Var p: boolean; Var ca: ArbFloat; Var term: ArbInt); 37 38 Procedure mdtgsy(n, rwidth: ArbInt; Var a: ArbFloat; Var pp:ArbInt; 39 Var qq:boolean; Var ca:ArbFloat; Var term:ArbInt); 40 41 Procedure mdtgpd(n, rwidth: ArbInt; Var al, ca: ArbFloat; Var term: ArbInt); 42 43 Procedure mdtgba(n, lb, rb, rwa: ArbInt; Var a: ArbFloat; rwl: ArbInt; 44 Var l:ArbFloat; Var p: ArbInt; Var ca: ArbFloat; Var term:ArbInt); 45 46 Procedure mdtgpb(n, lb, rwidth: ArbInt; Var al, ca: ArbFloat; 47 Var term: ArbInt); 48 49 Procedure mdtdtr(n: ArbInt; Var l, d, u, l1, d1, u1: ArbFloat; 50 Var term:ArbInt); 51 52 implementation 53 54 Procedure mdtgen(n, rwidth: ArbInt; Var alu: ArbFloat; Var p: ArbInt; 55 Var ca:ArbFloat; Var term: ArbInt); 56 57 Var 58 indi, indk, nsr, ind, i, j, k, indexpivot : ArbInt; 59 normr, sumrowi, pivot, l, normt, maxim, h, s : ArbFloat; 60 palu, sumrow, t : ^arfloat1; 61 pp : ^arint1; 62 singular : boolean; 63 Begin 64 If (n<1) Or (rwidth<1) Then 65 Begin 66 term := 3; 67 exit 68 End; {wrong input} 69 palu := @alu; 70 pp := @p; 71 nsr := n*sizeof(ArbFloat); 72 getmem(sumrow, nsr); 73 getmem(t, nsr); 74 normr := 0; 75 singular := false ; 76 For i:=1 To n Do 77 Begin 78 ind := (i-1)*rwidth; 79 pp^[i] := i; 80 sumrowi := 0; 81 For j:=1 To n Do 82 sumrowi := sumrowi+abs(palu^[ind+j]); 83 sumrow^[i] := sumrowi; 84 h := 2*random-1; 85 t^[i] := sumrowi*h; 86 h := abs(h); 87 If normr<h Then normr := h; 88 If sumrowi=0 Then 89 singular := true 90 End; {i} 91 For k:=1 To n Do 92 Begin 93 maxim := 0; 94 indexpivot := k; 95 For i:=k To n Do 96 Begin 97 ind := (i-1)*rwidth; 98 sumrowi := sumrow^[i]; 99 If sumrowi <> 0 Then 100 Begin 101 h := abs(palu^[ind+k])/sumrowi; 102 If maxim<h Then 103 Begin 104 maxim := h; 105 indexpivot := i 106 End {maxim<h} 107 End {sumrow <> 0} 108 End; {i} 109 If maxim=0 Then 110 singular := true 111 Else 112 Begin 113 If indexpivot <> k Then 114 Begin 115 ind := (indexpivot-1)*rwidth; 116 indk := (k-1)*rwidth; 117 For j:=1 To n Do 118 Begin 119 h := palu^[ind+j]; 120 palu^[ind+j] := palu^[indk+j]; 121 palu^[indk+j] := h 122 End; {j} 123 h := t^[indexpivot]; 124 t^[indexpivot] := t^[k]; 125 t^[k] := h; 126 pp^[k] := indexpivot; 127 sumrow^[indexpivot] := sumrow^[k] 128 End; {indexpivot <> k} 129 pivot := palu^[(k-1)*rwidth+k]; 130 For i:=k+1 To n Do 131 Begin 132 ind := (i-1)*rwidth; 133 l := palu^[ind+k]/pivot; 134 palu^[ind+k] := l; 135 If l <> 0 Then 136 Begin 137 For j:=k+1 To n Do 138 palu^[ind+j] := palu^[ind+j]-l*palu^[(k-1)*rwidth+j]; 139 If Not singular Then 140 t^[i] := t^[i]-l*t^[k] 141 End {l <> 0} 142 End {i} 143 End {maxim <> 0} 144 End; {k} 145 If Not singular Then 146 Begin 147 normt := 0; 148 For i:=n Downto 1 Do 149 Begin 150 indi := (i-1)*rwidth; 151 s := 0; 152 For j:=i+1 To n Do 153 s := s+t^[j]*palu^[indi+j]; 154 t^[i] := (t^[i]-s)/palu^[indi+i]; 155 h := abs(t^[i]); 156 If normt<h Then 157 normt := h 158 End; {i} 159 ca := normt/normr 160 End; {not singular} 161 If singular Then 162 Begin 163 term := 4; 164 ca := giant 165 End 166 Else 167 term := 1; 168 freemem(sumrow, nsr); 169 freemem(t, nsr) 170 End; {mdtgen} 171 172 Procedure mdtgtr(n: ArbInt; Var l, d, u, l1, d1, u1, u2: ArbFloat; 173 Var p: boolean; Var ca: ArbFloat; Var term: ArbInt); 174 175 Var 176 i, j, nmin1, sr : ArbInt; 177 normr, normt, sumrowi, h, lj, di, ui, ll : ArbFloat; 178 sing : boolean; 179 pd, pu, pd1, pu1, pu2, t, sumrow : ^arfloat1; 180 pl, pl1 : ^arfloat2; 181 pp : ^arbool1; 182 Begin 183 If n<1 Then 184 Begin 185 term := 3; 186 exit 187 End; {wrong input} 188 pl := @l; 189 pd := @d; 190 pu := @u; 191 pl1 := @l1; 192 pd1 := @d1; 193 pu1 := @u1; 194 pu2 := @u2; 195 pp := @p; 196 sr := sizeof(ArbFloat); 197 move(pl^, pl1^, (n-1)*sr); 198 move(pd^, pd1^, n*sr); 199 move(pu^, pu1^, (n-1)*sr); 200 getmem(t, n*sr); 201 getmem(sumrow, n*sr); 202 normr := 0; 203 sing := false; 204 nmin1 := n-1; 205 For i:=1 To n Do 206 Begin 207 pp^[i] := false; 208 If i=1 Then 209 sumrowi := abs(pd^[1])+abs(pu^[1]) 210 Else 211 If i=n Then 212 sumrowi := abs(pl^[n])+abs(pd^[n]) 213 Else 214 sumrowi := abs(pl^[i])+abs(pd^[i])+abs(pu^[i]); 215 sumrow^[i] := sumrowi; 216 h := 2*random-1; 217 t^[i] := sumrowi*h; 218 h := abs(h); 219 If normr<h Then 220 normr := h; 221 If sumrowi=0 Then 222 sing := true 223 End; {i} 224 j := 1; 225 while (j <> n) Do 226 Begin 227 i := j; 228 j := j+1; 229 lj := pl1^[j]; 230 If lj <> 0 Then 231 Begin 232 di := pd1^[i]; 233 If di=0 Then 234 pp^[i] := true 235 Else 236 pp^[i] := abs(di/sumrow^[i])<abs(lj/sumrow^[j]); 237 If pp^[i] Then 238 Begin 239 ui := pu1^[i]; 240 pd1^[i] := lj; 241 pu1^[i] := pd1^[j]; 242 pl1^[j] := di/lj; 243 ll := pl1^[j]; 244 pd1^[j] := ui-ll*pd1^[j]; 245 If i<nmin1 Then 246 Begin 247 pu2^[i] := pu1^[j]; 248 pu1^[j] := -ll*pu2^[i] 249 End; {i<nmin1} 250 sumrow^[j] := sumrow^[i]; 251 If (Not sing) Then 252 Begin 253 h := t^[i]; 254 t^[i] := t^[j]; 255 t^[j] := h-ll*t^[i] 256 End {not sing} 257 End {pp^[i]} 258 Else 259 Begin 260 pl1^[j] := lj/di; 261 ll := pl1^[j]; 262 pd1^[j] := pd1^[j]-ll*pu1^[i]; 263 If i<nmin1 Then 264 pu2^[i] := 0; 265 If (Not sing) Then 266 t^[j] := t^[j]-ll*t^[i] 267 End {not pp^[i]} 268 End {lj<>0} 269 Else 270 Begin 271 If i<nmin1 Then 272 pu2^[i] := 0; 273 If pd1^[i]=0 Then 274 sing := true 275 End {lj=0} 276 End; {j} 277 If pd1^[n]=0 Then 278 sing := true; 279 If (Not sing) Then 280 Begin 281 normt := 0; 282 t^[n] := t^[n]/pd1^[n]; 283 h := abs(t^[n]); 284 If normt<h Then 285 normt := h; 286 If n > 1 Then 287 Begin 288 t^[nmin1] := (t^[nmin1]-pu1^[nmin1]*t^[n])/pd1^[nmin1]; 289 h := abs(t^[nmin1]) 290 End; {n > 1} 291 If normt<h Then 292 normt := h; 293 For i:=n-2 Downto 1 Do 294 Begin 295 t^[i] := (t^[i]-pu1^[i]*t^[i+1]-pu2^[i]*t^[i+2])/pd1^[i]; 296 h := abs(t^[i]); 297 If normt<h Then 298 normt := h 299 End; {i} 300 ca := normt/normr 301 End; {not sing} 302 If (sing) Then 303 Begin 304 term := 4; 305 ca := giant 306 End {sing} 307 Else 308 term := 1; 309 freemem(t, n*sr); 310 freemem(sumrow, n*sr) 311 End; {mdtgtr} 312 313 Procedure mdtgsy(n, rwidth: ArbInt; Var a: ArbFloat; Var pp:ArbInt; 314 Var qq:boolean; Var ca:ArbFloat; Var term:ArbInt); 315 316 Var 317 i, j, kmin1, k, kplus1, kmin2, imin2, nsr, nsi, nsb, 318 imin1, jmin1, indexpivot, iplus1, indi, indj, indk, indp : ArbInt; 319 h, absh, maxim, pivot, ct, norma, sumrowi, normt, normr : ArbFloat; 320 alt, l, d, t, u, v, l1, d1, u1, t1 : ^arfloat1; 321 p : ^arint1; 322 q : ^arbool1; 323 Begin 324 If (n<1) Or (rwidth<1) Then 325 Begin 326 term := 3; 327 exit 328 End; {if} 329 alt := @a; 330 p := @pp; 331 q := @qq; 332 nsr := n*sizeof(ArbFloat); 333 nsi := n*sizeof(ArbInt); 334 nsb := n*sizeof(boolean); 335 getmem(l, nsr); 336 getmem(d, nsr); 337 getmem(t, nsr); 338 getmem(u, nsr); 339 getmem(v, nsr); 340 getmem(l1, nsr); 341 getmem(d1, nsr); 342 getmem(u1, nsr); 343 getmem(t1, nsr); 344 norma := 0; 345 For i:=1 To n Do 346 Begin 347 indi := (i-1)*rwidth; 348 p^[i] := i; 349 sumrowi := 0; 350 For j:=1 To i Do 351 sumrowi := sumrowi+abs(alt^[indi+j]); 352 For j:=i+1 To n Do 353 sumrowi := sumrowi+abs(alt^[(j-1)*rwidth+i]); 354 If norma<sumrowi Then 355 norma := sumrowi 356 End; {i} 357 kmin1 := -1; 358 k := 0; 359 kplus1 := 1; 360 while k<n Do 361 Begin 362 kmin2 := kmin1; 363 kmin1 := k; 364 k := kplus1; 365 kplus1 := kplus1+1; 366 indk := kmin1*rwidth; 367 If k>3 Then 368 Begin 369 t^[2] := alt^[rwidth+2]*alt^[indk+1]+alt^[2*rwidth+2]*alt^[indk+2]; 370 For i:=3 To kmin2 Do 371 Begin 372 indi := (i-1)*rwidth; 373 t^[i] := alt^[indi+i-1]*alt^[indk+i-2]+alt^[indi+i] 374 *alt^[indk+i-1]+alt^[indi+rwidth+i]*alt^[indk+i] 375 End; {i} 376 t^[kmin1] := alt^[indk-rwidth+kmin2]*alt^[indk+k-3] 377 +alt^[indk-rwidth+kmin1]*alt^[indk+kmin2] 378 +alt^[indk+kmin1]; 379 h := alt^[indk+k]; 380 For j:=2 To kmin1 Do 381 h := h-t^[j]*alt^[indk+j-1]; 382 t^[k] := h; 383 alt^[indk+k] := h-alt^[indk+kmin1]*alt^[indk+kmin2] 384 End {k>3} 385 Else 386 If k=3 Then 387 Begin 388 t^[2] := alt^[rwidth+2]*alt^[2*rwidth+1]+alt^[2*rwidth+2]; 389 h := alt^[2*rwidth+3]-t^[2]*alt^[2*rwidth+1]; 390 t^[3] := h; 391 alt^[2*rwidth+3] := h-alt^[2*rwidth+2]*alt^[2*rwidth+1] 392 End {k=3} 393 Else 394 If k=2 Then 395 t^[2] := alt^[rwidth+2]; 396 maxim := 0; 397 For i:=kplus1 To n Do 398 Begin 399 indi := (i-1)*rwidth; 400 h := alt^[indi+k]; 401 For j:=2 To k Do 402 h := h-t^[j]*alt^[indi+j-1]; 403 absh := abs(h); 404 If maxim<absh Then 405 Begin 406 maxim := absh; 407 indexpivot := i 408 End; {if} 409 alt^[indi+k] := h 410 End; {i} 411 If maxim <> 0 Then 412 Begin 413 If indexpivot>kplus1 Then 414 Begin 415 indp := (indexpivot-1)*rwidth; 416 indk := k*rwidth; 417 p^[kplus1] := indexpivot; 418 For j:=1 To k Do 419 Begin 420 h := alt^[indk+j]; 421 alt^[indk+j] := alt^[indp+j]; 422 alt^[indp+j] := h 423 End; {j} 424 For j:=indexpivot Downto kplus1 Do 425 Begin 426 indj := (j-1)*rwidth; 427 h := alt^[indj+kplus1]; 428 alt^[indj+kplus1] := alt^[indp+j]; 429 alt^[indp+j] := h 430 End; {j} 431 For i:=indexpivot To n Do 432 Begin 433 indi := (i-1)*rwidth; 434 h := alt^[indi+kplus1]; 435 alt^[indi+kplus1] := alt^[indi+indexpivot]; 436 alt^[indi+indexpivot] := h 437 End {i} 438 End; {if} 439 pivot := alt^[k*rwidth+k]; 440 For i:=k+2 To n Do 441 alt^[(i-1)*rwidth+k] := alt^[(i-1)*rwidth+k]/pivot 442 End {maxim <> 0} 443 End; {k} 444 d^[1] := alt^[1]; 445 i := 1; 446 while i<n Do 447 Begin 448 imin1 := i; 449 i := i+1; 450 u^[imin1] := alt^[(i-1)*rwidth+imin1]; 451 l^[imin1] := u^[imin1]; 452 d^[i] := alt^[(i-1)*rwidth+i] 453 End; {i} 454 mdtgtr(n, l^[1], d^[1], u^[1], l1^[1], d1^[1], u1^[1], v^[1], 455 q^[1], ct, term); 456 alt^[1] := d1^[1]; 457 alt^[rwidth+1] := l1^[1]; 458 alt^[rwidth+2] := d1^[2]; 459 alt^[2] := u1^[1]; 460 imin1 := 1; 461 i := 2; 462 while i<n Do 463 Begin 464 imin2 := imin1; 465 imin1 := i; 466 i := i+1; 467 indi := imin1*rwidth; 468 alt^[indi+imin1] := l1^[imin1]; 469 alt^[indi+i] := d1^[i]; 470 alt^[(imin1-1)*rwidth+i] := u1^[imin1]; 471 alt^[(imin2-1)*rwidth+i] := v^[imin2] 472 End; {i} 473 If term=1 Then 474 Begin 475 normr := 0; 476 For i:=1 To n Do 477 Begin 478 t^[i] := 2*random-1; 479 h := t^[i]; 480 h := abs(h); 481 If normr<h Then 482 normr := h 483 End; {i} 484 i := 0; 485 while i<n Do 486 Begin 487 imin1 := i; 488 i := i+1; 489 j := 1; 490 h := t^[i]; 491 while j<imin1 Do 492 Begin 493 jmin1 := j; 494 j := j+1; 495 h := h-alt^[(i-1)*rwidth+jmin1]*t^[j] 496 End; {j} 497 t^[i] := h 498 End; {i} 499 dslgtr(n, l1^[1], d1^[1], u1^[1], v^[1], q^[1], t^[1], t1^[1], term); 500 i := n+1; 501 imin1 := n; 502 normt := 0; 503 while i>2 Do 504 Begin 505 iplus1 := i; 506 i := imin1; 507 imin1 := imin1-1; 508 h := t1^[i]; 509 For j:=iplus1 To n Do 510 h := h-alt^[(j-1)*rwidth+imin1]*t1^[j]; 511 t1^[i] := h; 512 h := abs(h); 513 If normt<h Then 514 normt := h 515 End; {i} 516 ca := norma*normt/normr 517 End {term=1} 518 Else ca := giant; 519 freemem(l, nsr); 520 freemem(d, nsr); 521 freemem(t, nsr); 522 freemem(u, nsr); 523 freemem(v, nsr); 524 freemem(l1, nsr); 525 freemem(d1, nsr); 526 freemem(u1, nsr); 527 freemem(t1, nsr) 528 End; {mdtgsy} 529 530 Procedure mdtgpd(n, rwidth: ArbInt; Var al, ca: ArbFloat; Var term: ArbInt); 531 532 Var 533 posdef : boolean; 534 i, j, k, kmin1, indk, indi : ArbInt; 535 h, lkk, normr, normt, sumrowi, norma : ArbFloat; 536 pal, t : ^arfloat1; 537 Begin 538 If (n<1) Or (rwidth<1) Then 539 Begin 540 term := 3; 541 exit 542 End; {wrong input} 543 getmem(t, sizeof(ArbFloat)*n); 544 pal := @al; 545 normr := 0; 546 posdef := true; 547 norma := 0; 548 For i:=1 To n Do 549 Begin 550 sumrowi := 0; 551 For j:=1 To i Do 552 sumrowi := sumrowi+abs(pal^[(i-1)*rwidth+j]); 553 For j:=i+1 To n Do 554 sumrowi := sumrowi+abs(pal^[(j-1)*rwidth+i]); 555 If norma<sumrowi Then 556 norma := sumrowi; 557 t^[i] := 2*random-1; 558 h := t^[i]; 559 h := abs(h); 560 If normr<h Then 561 normr := h 562 End; {i} 563 k := 0; 564 while (k<n) and posdef Do 565 Begin 566 kmin1 := k; 567 k := k+1; 568 indk := (k-1)*rwidth; 569 lkk := pal^[indk+k]; 570 For j:=1 To kmin1 Do 571 lkk := lkk-sqr(pal^[indk+j]); 572 If lkk <= 0 Then 573 posdef := false 574 Else 575 Begin 576 pal^[indk+k] := sqrt(lkk); 577 lkk := pal^[indk+k]; 578 For i:=k+1 To n Do 579 Begin 580 indi := (i-1)*rwidth; 581 h := pal^[indi+k]; 582 For j:=1 To kmin1 Do 583 h := h-pal^[indk+j]*pal^[indi+j]; 584 pal^[indi+k] := h/lkk 585 End; {i} 586 h := t^[k]; 587 For j:=1 To kmin1 Do 588 h := h-pal^[indk+j]*t^[j]; 589 t^[k] := h/lkk 590 End {posdef} 591 End; {k} 592 If posdef Then 593 Begin 594 normt := 0; 595 For i:=n Downto 1 Do 596 Begin 597 h := t^[i]; 598 For j:=i+1 To n Do 599 h := h-pal^[(j-1)*rwidth+i]*t^[j]; 600 t^[i] := h/pal^[(i-1)*rwidth+i]; 601 h := abs(t^[i]); 602 If normt<h Then 603 normt := h 604 End; {i} 605 ca := norma*normt/normr 606 End; {posdef} 607 If posdef Then 608 term := 1 609 Else 610 term := 2; 611 freemem(t, sizeof(ArbFloat)*n); 612 End; {mdtgpd} 613 614 Procedure mdtgba(n, lb, rb, rwa: ArbInt; Var a: ArbFloat; rwl: ArbInt; 615 Var l:ArbFloat; Var p: ArbInt; Var ca: ArbFloat; Var term:ArbInt); 616 617 Var 618 sr, i, j, k, ipivot, lbj, lbi, ubi, ls, 619 ii, jj, ll, jl, ubj : ArbInt; 620 normr, sumrowi, pivot, normt, maxim, h : ArbFloat; 621 pl, au, sumrow, t, row : ^arfloat1; 622 pp : ^arint1; 623 624 Begin 625 If (n<1) Or (lb<0) Or (rb<0) Or (lb>n-1) Or (rb>n-1) Or (rwl<0) Or (rwa<1) Then 626 Begin 627 term := 3; 628 exit 629 End; {term=3} 630 sr := sizeof(ArbFloat); 631 au := @a; 632 pl := @l; 633 pp := @p; 634 ll := lb+rb+1; 635 ls := ll*sr; 636 getmem(sumrow, n*sr); 637 getmem(t, n*sr); 638 getmem(row, ls); 639 lbi := n-rb+1; 640 lbj := 0; 641 jj := 1; 642 For i:=lb Downto 1 Do 643 Begin 644 move(au^[i+jj], au^[jj], (ll-i)*sr); 645 fillchar(au^[jj+ll-i], i*sr, 0); 646 jj := jj+rwa 647 End; {i} 648 jj := (n-rb)*rwa+ll; 649 For i:=1 To rb Do 650 Begin 651 fillchar(au^[jj], i*sr, 0); 652 jj := jj+rwa-1 653 End; {i} 654 normr := 0; 655 term := 1; 656 ii := 1; 657 For i:=1 To n Do 658 Begin 659 pp^[i] := i; 660 sumrowi := omvn1v(au^[ii], ll); 661 ii := ii+rwa; 662 sumrow^[i] := sumrowi; 663 h := 2*random-1; 664 t^[i] := sumrowi*h; 665 h := abs(h); 666 If normr<h Then 667 normr := h; 668 If sumrowi=0 Then 669 term := 4 670 End; {i} 671 ubi := lb; 672 jj := 1; 673 For k:=1 To n Do 674 Begin 675 maxim := 0; 676 ipivot := k; 677 ii := jj; 678 If ubi<n Then 679 ubi := ubi+1; 680 For i:=k To ubi Do 681 Begin 682 sumrowi := sumrow^[i]; 683 If sumrowi <> 0 Then 684 Begin 685 h := abs(au^[ii])/sumrowi; 686 ii := ii+rwa; 687 If maxim<h Then 688 Begin 689 maxim := h; 690 ipivot := i 691 End {maxim<h} 692 End {sumrowi <> 0} 693 End; {i} 694 If maxim=0 Then 695 Begin 696 lbj := 1; 697 ubj := ubi-k; 698 For j:=lbj To ubj Do 699 pl^[(k-1)*rwl+j] := 0; 700 For i:=k+1 To ubi Do 701 Begin 702 ii := (i-1)*rwa; 703 For j:=2 To ll Do 704 au^[ii+j-1] := au^[ii+j]; 705 au^[ii+ll] := 0 706 End; {i} 707 term := 4 708 End {maxim=0} 709 Else 710 Begin 711 If ipivot <> k Then 712 Begin 713 ii := (ipivot-1)*rwa+1; 714 move(au^[ii], row^, ls); 715 move(au^[jj], au^[ii], ls); 716 move(row^, au^[jj], ls); 717 h := t^[ipivot]; 718 t^[ipivot] := t^[k]; 719 t^[k] := h; 720 pp^[k] := ipivot; 721 sumrow^[ipivot] := sumrow^[k] 722 End; {ipivot <> k} 723 pivot := au^[jj]; 724 jl := 0; 725 ii := jj; 726 For i:=k+1 To ubi Do 727 Begin 728 jl := jl+1; 729 ii := ii+rwa; 730 h := au^[ii]/pivot; 731 pl^[(k-1)*rwl+jl] := h; 732 For j:=0 To ll-2 Do 733 au^[ii+j] := au^[ii+j+1]-h*au^[jj+j+1]; 734 au^[ii+ll-1] := 0; 735 If term=1 Then 736 t^[i] := t^[i]-h*t^[k] 737 End {i} 738 End; {maxim <> 0} 739 jj := jj+rwa 740 End; {k} 741 If term=1 Then 742 Begin 743 normt := 0; 744 ubj := -lb-1; 745 jj := n*rwa+1; 746 For i:=n Downto 1 Do 747 Begin 748 jj := jj-rwa; 749 If ubj<rb Then 750 ubj := ubj+1; 751 h := t^[i]; 752 For j:=1 To ubj+lb Do 753 h := h-au^[jj+j]*t^[i+j]; 754 t^[i] := h/au^[jj]; 755 h := abs(t^[i]); 756 If normt<h Then 757 normt := h 758 End; {i} 759 ca := normt/normr 760 End {term=1} 761 Else 762 ca := giant; 763 freemem(sumrow, n*sr); 764 freemem(t, n*sr); 765 freemem(row, ls) 766 End; {mdtgba} 767 768 Procedure mdtgpb(n, lb, rwidth: ArbInt; Var al, ca: ArbFloat; 769 Var term: ArbInt); 770 771 Var 772 posdef : boolean; 773 i, j, k, r, p, q, ll, llmin1, jmin1, indi : ArbInt; 774 h, normr, normt, sumrowi, alim, norma : ArbFloat; 775 pal, t : ^arfloat1; 776 777 Procedure decomp(i, r: ArbInt); 778 779 Var 780 k, ii, ir : ArbInt; 781 Begin 782 ii := (i-1)*rwidth; 783 ir := (r-1)*rwidth; 784 h := pal^[ii+j]; 785 q := ll-j+p; 786 For k:=p To jmin1 Do 787 Begin 788 h := h-pal^[ii+k]*pal^[ir+q]; 789 q := q+1 790 End; {k} 791 If j<ll Then 792 pal^[ii+j] := h/pal^[ir+ll] 793 End; {decomp} 794 795 Procedure lmin1t(i: ArbInt); 796 797 Var 798 k:ArbInt; 799 Begin 800 h := t^[i]; 801 q := i; 802 For k:=llmin1 Downto p Do 803 Begin 804 q := q-1; 805 h := h-pal^[indi+k]*t^[q] 806 End; {k} 807 t^[i] := h/alim 808 End; {lmin1t} 809 810 Begin 811 If (lb<0) Or (lb>n-1) Or (n<1) Or (rwidth<1) Then 812 Begin 813 term := 3; 814 exit 815 End; {wrong input} 816 pal := @al; 817 getmem(t, n*sizeof(ArbFloat)); 818 ll := lb+1; 819 normr := 0; 820 p := ll+1; 821 norma := 0; 822 For i:=1 To n Do 823 Begin 824 If p>1 Then 825 p := p-1; 826 indi := (i-1)*rwidth+p; 827 sumrowi := omvn1v(pal^[indi], ll-p+1); 828 r := i; 829 j := ll; 830 while (r<n) and (j>1) Do 831 Begin 832 r := r+1; 833 j := j-1; 834 sumrowi := sumrowi+abs(pal^[(r-1)*rwidth+j]) 835 End; {r,j} 836 If norma<sumrowi Then 837 norma := sumrowi; 838 h := 2*random-1; 839 t^[i] := h; 840 h := abs(h); 841 If normr<h Then 842 normr := h 843 End; {i} 844 llmin1 := ll-1; 845 p := ll+1; 846 i := 0; 847 posdef := true ; 848 while (i<n) and posdef Do 849 Begin 850 i := i+1; 851 indi := (i-1)*rwidth; 852 If p>1 Then 853 p := p-1; 854 r := i-ll+p; 855 j := p-1; 856 while j<llmin1 Do 857 Begin 858 jmin1 := j; 859 j := j+1; 860 decomp(i, r); 861 r := r+1 862 End; {j} 863 jmin1 := llmin1; 864 j := ll; 865 decomp(i, i); 866 If h <= 0 Then 867 posdef := false 868 Else 869 Begin 870 alim := sqrt(h); 871 pal^[indi+ll] := alim; 872 lmin1t(i) 873 End 874 End; {i} 875 If posdef Then 876 Begin 877 normt := 0; 878 p := ll+1; 879 For i:=n Downto 1 Do 880 Begin 881 If p>1 Then 882 p := p-1; 883 q := i; 884 h := t^[i]; 885 For k:=llmin1 Downto p Do 886 Begin 887 q := q+1; 888 h := h-pal^[(q-1)*rwidth+k]*t^[q] 889 End; {k} 890 t^[i] := h/pal^[(i-1)*rwidth+ll]; 891 h := abs(t^[i]); 892 If normt<h Then 893 normt := h 894 End; {i} 895 ca := norma*normt/normr 896 End; {posdef} 897 If posdef Then 898 term := 1 899 Else 900 term := 2; 901 freemem(t, n*sizeof(ArbFloat)); 902 End; {mdtgpb} 903 904 Procedure mdtdtr(n: ArbInt; Var l, d, u, l1, d1, u1: ArbFloat; 905 Var term:ArbInt); 906 907 Var 908 i, j, s : ArbInt; 909 lj, di : ArbFloat; 910 pd, pu, pd1, pu1 : ^arfloat1; 911 pl, pl1 : ^arfloat2; 912 913 Begin 914 If n<1 Then 915 Begin 916 term := 3; 917 exit 918 End; {wrong input} 919 pl := @l; 920 pd := @d; 921 pu := @u; 922 pl1 := @l1; 923 pd1 := @d1; 924 pu1 := @u1; 925 s := sizeof(ArbFloat); 926 move(pl^, pl1^, (n-1)*s); 927 move(pd^, pd1^, n*s); 928 move(pu^, pu1^, (n-1)*s); 929 j := 1; 930 di := pd1^[j]; 931 If di=0 Then 932 term := 2 933 Else 934 term := 1; 935 while (term=1) and (j <> n) Do 936 Begin 937 i := j; 938 j := j+1; 939 lj := pl1^[j]/di; 940 pl1^[j] := lj; 941 di := pd1^[j]-lj*pu1^[i]; 942 pd1^[j] := di; 943 If di=0 Then 944 term := 2 945 End {j} 946 End; {mdtdtr} 947 948 Begin 949 {$ifdef fixate_random} 950 randseed := 12345 951 {$endif} 952 End. 953