← Back to branch summary

~ubuntu-branches/debian/lenny/fpc/lenny

« back to all changes in this revision

Viewing changes to fpcsrc/packages/extra/numlib/mdt.pas

  • Committer: Bazaar Package Importer
  • Author(s): Mazen Neifer, Torsten Werner, Mazen Neifer
  • Date: 2008-05-17 17:12:11 UTC
  • mfrom: (3.1.9 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080517171211-9qi33xhd9evfa0kg
Tags: 2.2.0-dfsg1-9
http://bugs.debian.org/477743
[ Torsten Werner ]
* Add Mazen Neifer to Uploaders field.

[ Mazen Neifer ]
* Moved FPC sources into a version dependent directory from /usr/share/fpcsrc
  to /usr/share/fpcsrc/${FPCVERSION}. This allow installing more than on FPC
  release.
* Fixed far call issue in compiler preventing building huge binearies.
  (closes: #477743)
* Updated building dependencies, recomennded and suggested packages.
* Moved fppkg to fp-utils as it is just a helper tool and is not required by
  compiler.

Show diffs side-by-side

added added

removed removed

Lines of Context:
fpcsrc/packages/extra/numlib/mdt.pas
 
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, k, 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, ii,
 
318
   imin1, jmin1, indexpivot, iplus1, indi, indj, indk, indp       : ArbInt;
 
319
   ra, h, absh, maxim, pivot, ct, norma, sumrowi, normt, normr, s : 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, m, lbj, lbi, ubi, ls,
 
619
             ii, jj, ll, s, js, jl, ubj       : ArbInt;
 
620
  ra, 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
  randseed := 12345
 
950
End.