← Back to branch summary

~ubuntu-branches/ubuntu/saucy/lazarus/saucy

« back to all changes in this revision

Viewing changes to components/tachart/numlib_fix/sle.pas

  • Committer: Package Import Robot
  • Author(s): Paul Gevers, Abou Al Montacir, Bart Martens, Paul Gevers
  • Date: 2013-06-08 14:12:17 UTC
  • mfrom: (1.1.9)
  • Revision ID: package-import@ubuntu.com-20130608141217-7k0cy9id8ifcnutc
Tags: 1.0.8+dfsg-1
http://bugs.debian.org/708695
[ Abou Al Montacir ]
* New upstream major release and multiple maintenace release offering many
  fixes and new features marking a new milestone for the Lazarus development
  and its stability level.
  - The detailed list of changes can be found here:
    http://wiki.lazarus.freepascal.org/Lazarus_1.0_release_notes
    http://wiki.lazarus.freepascal.org/Lazarus_1.0_fixes_branch
* LCL changes:
  - LCL is now a normal package.
      + Platform independent parts of the LCL are now in the package LCLBase
      + LCL is automatically recompiled when switching the target platform,
        unless pre-compiled binaries for this target are already installed.
      + No impact on existing projects.
      + Linker options needed by LCL are no more added to projects that do
        not use the LCL package.
  - Minor changes in LCL basic classes behaviour
      + TCustomForm.Create raises an exception if a form resource is not
        found.
      + TNotebook and TPage: a new implementation of these classes was added.
      + TDBNavigator: It is now possible to have focusable buttons by setting
        Options = [navFocusableButtons] and TabStop = True, useful for
        accessibility and for devices with neither mouse nor touch screen.
      + Names of TControlBorderSpacing.GetSideSpace and GetSpace were swapped
        and are now consistent. GetSideSpace = Around + GetSpace.
      + TForm.WindowState=wsFullscreen was added
      + TCanvas.TextFitInfo was added to calculate how many characters will
        fit into a specified Width. Useful for word-wrapping calculations.
      + TControl.GetColorResolvingParent and
        TControl.GetRGBColorResolvingParent were added, simplifying the work
        to obtain the final color of the control while resolving clDefault
        and the ParentColor.
      + LCLIntf.GetTextExtentExPoint now has a good default implementation
        which works in any platform not providing a specific implementation.
        However, Widgetset specific implementation is better, when available.
      + TTabControl was reorganized. Now it has the correct class hierarchy
        and inherits from TCustomTabControl as it should.
  - New unit in the LCL:
      + lazdialogs.pas: adds non-native versions of various native dialogs,
        for example TLazOpenDialog, TLazSaveDialog, TLazSelectDirectoryDialog.
        It is used by widgetsets which either do not have a native dialog, or
        do not wish to use it because it is limited. These dialogs can also be
        used by user applications directly.
      + lazdeviceapis.pas: offers an interface to more hardware devices such
        as the accelerometer, GPS, etc. See LazDeviceAPIs
      + lazcanvas.pas: provides a TFPImageCanvas descendent implementing
        drawing in a LCL-compatible way, but 100% in Pascal.
      + lazregions.pas. LazRegions is a wholly Pascal implementation of
        regions for canvas clipping, event clipping, finding in which control
        of a region tree one an event should reach, for drawing polygons, etc.
      + customdrawncontrols.pas, customdrawndrawers.pas,
        customdrawn_common.pas, customdrawn_android.pas and
        customdrawn_winxp.pas: are the Lazarus Custom Drawn Controls -controls
        which imitate the standard LCL ones, but with the difference that they
        are non-native and support skinning.
  - New APIs added to the LCL to improve support of accessibility software
    such as screen readers.
* IDE changes:
  - Many improvments.
  - The detailed list of changes can be found here:
    http://wiki.lazarus.freepascal.org/New_IDE_features_since#v1.0_.282012-08-29.29
    http://wiki.lazarus.freepascal.org/Lazarus_1.0_release_notes#IDE_Changes
* Debugger / Editor changes:
  - Added pascal sources and breakpoints to the disassembler
  - Added threads dialog.
* Components changes:
  - TAChart: many fixes and new features
  - CodeTool: support Delphi style generics and new syntax extensions.
  - AggPas: removed to honor free licencing. (Closes: Bug#708695)
[Bart Martens]
* New debian/watch file fixing issues with upstream RC release.
[Abou Al Montacir]
* Avoid changing files in .pc hidden directory, these are used by quilt for
  internal purpose and could lead to surprises during build.
[Paul Gevers]
* Updated get-orig-source target and it compinion script orig-tar.sh so that they
  repack the source file, allowing bug 708695 to be fixed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
components/tachart/numlib_fix/sle.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
    !! modifies randseed, might not exactly work as TP version!!!
 
11
 
 
12
    Solve set of linear equations of the type Ax=b, for generic, and a
 
13
    variety of special matrices.
 
14
 
 
15
    See the file COPYING.FPC, included in this distribution,
 
16
    for details about the copyright.
 
17
 
 
18
    This program is distributed in the hope that it will be useful,
 
19
    but WITHOUT ANY WARRANTY; without even the implied warranty of
 
20
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
21
 
 
22
 **********************************************************************}
 
23
 
 
24
{Solve set of linear equations of the type Ax=b, for generic, and a variety of
 
25
special matrices.
 
26
One (generic) function for overdetermined sets of this kind : slegls
 
27
 
 
28
overdetermined are sets that look like this: (I don't know if I
 
29
translated "overdetermined" right)
 
30
 
 
31
    6   1  2  3     9
 
32
    3   9  3  4     2
 
33
   17  27 42 15    62
 
34
   17  27 42 15    61
 
35
 
 
36
The two bottom rows look much alike, which introduces a big uncertainty in the
 
37
result, therefore these matrices need special treatment.
 
38
 
 
39
All procedures have similar procedure with a "L" appended to the name. We
 
40
didn't receive docs for those procedures. If you know what the difference is,
 
41
please mail us }
 
42
 
 
43
Unit sle;
 
44
interface
 
45
{$I DIRECT.INC}
 
46
 
 
47
uses typ, omv;
 
48
 
 
49
{solve for special tridiagonal matrices}
 
50
Procedure sledtr(n: ArbInt; Var l, d, u, b, x: ArbFloat; Var term: ArbInt);
 
51
 
 
52
{solve for generic bandmatrices}
 
53
Procedure slegba(n, l, r: ArbInt;
 
54
                 Var a, b, x, ca: ArbFloat; Var term:ArbInt);
 
55
 
 
56
Procedure slegbal(n, l, r: ArbInt;
 
57
                  Var a1; Var b1, x1, ca: ArbFloat; Var term: ArbInt);
 
58
 
 
59
{generic solve for all matrices}
 
60
Procedure slegen(n, rwidth: ArbInt; Var a, b, x, ca: ArbFloat;
 
61
                 Var term: ArbInt);
 
62
 
 
63
Procedure slegenl(    n: ArbInt;
 
64
                  Var a1;
 
65
                  Var b1, x1, ca: ArbFloat;
 
66
                  Var term: ArbInt);
 
67
 
 
68
{solve for overdetermined matrices, see unit comments}
 
69
Procedure slegls(Var a: ArbFloat; m, n, rwidtha: ArbInt; Var b, x: ArbFloat;
 
70
                 Var term: ArbInt);
 
71
 
 
72
 
 
73
Procedure sleglsl(Var a1; m, n: ArbInt; Var b1, x1: ArbFloat;
 
74
                  Var term: ArbInt);
 
75
 
 
76
{Symmetrical positive definitive bandmatrices}
 
77
Procedure slegpb(n, l: ArbInt; Var a, b, x, ca: ArbFloat;
 
78
                 Var term: ArbInt);
 
79
 
 
80
Procedure slegpbl(n, l: ArbInt;
 
81
                  Var a1; Var b1, x1, ca: ArbFloat; Var term: ArbInt);
 
82
 
 
83
{Symmetrical positive definitive matrices}
 
84
Procedure slegpd(n, rwidth: ArbInt; Var a, b, x, ca: ArbFloat;
 
85
                 Var term: ArbInt);
 
86
 
 
87
Procedure slegpdl(n: ArbInt; Var a1; Var b1, x1, ca: ArbFloat;
 
88
                  Var term: ArbInt);
 
89
 
 
90
{Symmetrical matrices}
 
91
Procedure slegsy(n, rwidth: ArbInt; Var a, b, x, ca: ArbFloat;
 
92
                 Var term: ArbInt);
 
93
 
 
94
Procedure slegsyl(n: ArbInt; Var a1; Var b1, x1, ca: ArbFloat;
 
95
                  Var term: ArbInt);
 
96
 
 
97
{tridiagonal matrices}
 
98
Procedure slegtr(n:ArbInt; Var l, d, u, b, x, ca: ArbFloat;
 
99
                 Var term: ArbInt);
 
100
 
 
101
implementation
 
102
 
 
103
Uses DSL,MDT;
 
104
 
 
105
{Here originally stood an exact copy of mdtgtr from unit mdt}
 
106
{Here originally stood an exact copy of dslgtr from unit DSL}
 
107
 
 
108
Procedure decomp(Var qr: ArbFloat; m, n, rwidthq: ArbInt; Var alpha: ArbFloat;
 
109
                 Var pivot, term: ArbInt);
 
110
 
 
111
Var  i, j, jbar, k, ns, ii        : ArbInt;
 
112
     beta, sigma, alphak, qrkk, s : ArbFloat;
 
113
     pqr, pal, y, sum             : ^arfloat1;
 
114
     piv                          : ^arint1;
 
115
 
 
116
Begin
 
117
  term := 1;
 
118
  pqr := @qr;
 
119
  pal := @alpha;
 
120
  piv := @pivot;
 
121
  ns := n*sizeof(ArbFloat);
 
122
  getmem(y, ns);
 
123
  getmem(sum, ns);
 
124
  For j:=1 To n Do
 
125
    Begin
 
126
      s := 0;
 
127
      For i:=1 To m Do
 
128
        s := s+sqr(pqr^[(i-1)*rwidthq+j]);
 
129
      sum^[j] := s;
 
130
      piv^[j] := j
 
131
    End; {j}
 
132
  For k:=1 To n Do
 
133
    Begin
 
134
      sigma := sum^[k];
 
135
      jbar := k;
 
136
      For j:=k+1 To n Do
 
137
        If sigma < sum^[j] Then
 
138
          Begin
 
139
            sigma := sum^[j];
 
140
           jbar := j
 
141
          End;
 
142
      If jbar <> k
 
143
       Then
 
144
        Begin
 
145
          i := piv^[k];
 
146
          piv^[k] := piv^[jbar];
 
147
          piv^[jbar] := i;
 
148
          sum^[jbar] := sum^[k];
 
149
          sum^[k] := sigma;
 
150
          For i:=1 To m Do
 
151
            Begin
 
152
              ii := (i-1)*rwidthq;
 
153
              sigma := pqr^[ii+k];
 
154
              pqr^[ii+k] := pqr^[ii+jbar];
 
155
              pqr^[ii+jbar] := sigma
 
156
            End; {i}
 
157
        End; {column interchange}
 
158
      sigma := 0;
 
159
      For i:=k To m Do
 
160
        sigma := sigma+sqr(pqr^[(i-1)*rwidthq+k]);
 
161
      If sigma=0 Then
 
162
        Begin
 
163
          term := 2;
 
164
          freemem(y, ns);
 
165
          freemem(sum, ns);
 
166
          exit
 
167
        End;
 
168
      qrkk := pqr^[(k-1)*rwidthq+k];
 
169
      If qrkk < 0 Then
 
170
        alphak := sqrt(sigma)
 
171
      Else
 
172
        alphak := -sqrt(sigma);
 
173
      pal^[k] := alphak;
 
174
      beta := 1/(sigma-qrkk*alphak);
 
175
      pqr^[(k-1)*rwidthq+k] := qrkk-alphak;
 
176
      For j:=k+1 To n Do
 
177
        Begin
 
178
          s := 0;
 
179
          For i:=k To m Do
 
180
            Begin
 
181
              ii := (i-1)*rwidthq;
 
182
              s := s+pqr^[ii+k]*pqr^[ii+j]
 
183
            End; {i}
 
184
          y^[j] := beta*s
 
185
        End; {j}
 
186
      For j:=k+1 To n Do
 
187
        Begin
 
188
          For i:=k To m Do
 
189
            Begin
 
190
              ii := (i-1)*rwidthq;
 
191
              pqr^[ii+j] := pqr^[ii+j]-pqr^[ii+k]*y^[j]
 
192
            End; {i}
 
193
          sum^[j] := sum^[j]-sqr(pqr^[(k-1)*rwidthq+j])
 
194
        End {j}
 
195
    End; {k}
 
196
  freemem(y, ns);
 
197
 freemem(sum, ns);
 
198
End; {decomp}
 
199
 
 
200
Procedure decomp1(Var qr1; m, n: ArbInt; Var alpha1: ArbFloat;
 
201
                  Var pivot1, term: ArbInt);
 
202
 
 
203
Var             i, j, jbar, k, ns : ArbInt;
 
204
     beta, sigma, alphak, qrkk, s : ArbFloat;
 
205
     qr                           : ar2dr1 absolute qr1;
 
206
     alpha                        : arfloat1 absolute alpha1;
 
207
     pivot                        : arint1 absolute pivot1;
 
208
     y, sum                       : ^arfloat1;
 
209
Begin
 
210
  term := 1;
 
211
  ns := n*sizeof(ArbFloat);
 
212
  getmem(y, ns);
 
213
 getmem(sum, ns);
 
214
  For j:=1 To n Do
 
215
    Begin
 
216
      s := 0;
 
217
      For i:=1 To m Do
 
218
       s := s+sqr(qr[i]^[j]);
 
219
      sum^[j] := s;
 
220
     pivot[j] := j
 
221
    End; {j}
 
222
  For k:=1 To n Do
 
223
    Begin
 
224
      sigma := sum^[k];
 
225
     jbar := k;
 
226
      For j:=k+1 To n Do
 
227
        If sigma < sum^[j]
 
228
         Then
 
229
          Begin
 
230
            sigma := sum^[j];
 
231
           jbar := j
 
232
          End;
 
233
      If jbar <> k
 
234
       Then
 
235
        Begin
 
236
          i := pivot[k];
 
237
         pivot[k] := pivot[jbar];
 
238
         pivot[jbar] := i;
 
239
          sum^[jbar] := sum^[k];
 
240
         sum^[k] := sigma;
 
241
          For i:=1 To m Do
 
242
            Begin
 
243
              sigma := qr[i]^[k];
 
244
             qr[i]^[k] := qr[i]^[jbar];
 
245
              qr[i]^[jbar] := sigma
 
246
            End; {i}
 
247
        End; {column interchange}
 
248
      sigma := 0;
 
249
      For i:=k To m Do
 
250
       sigma := sigma+sqr(qr[i]^[k]);
 
251
      If sigma=0
 
252
       Then
 
253
        Begin
 
254
          term := 2;
 
255
         freemem(y, ns);
 
256
         freemem(sum, ns);
 
257
         exit
 
258
        End;
 
259
      qrkk := qr[k]^[k];
 
260
      If qrkk < 0 Then alphak := sqrt(sigma)
 
261
     Else alphak := -sqrt(sigma);
 
262
      alpha[k] := alphak;
 
263
      beta := 1/(sigma-qrkk*alphak);
 
264
      qr[k]^[k] := qrkk-alphak;
 
265
      For j:=k+1 To n Do
 
266
        Begin
 
267
          s := 0;
 
268
         For i:=k To m Do
 
269
          s := s+qr[i]^[k]*qr[i]^[j];
 
270
         y^[j] := beta*s
 
271
        End; {j}
 
272
      For j:=k+1 To n Do
 
273
        Begin
 
274
          For i:=k To m Do
 
275
           qr[i]^[j] := qr[i]^[j]-qr[i]^[k]*y^[j];
 
276
          sum^[j] := sum^[j]-sqr(qr[k]^[j])
 
277
        End {j}
 
278
    End; {k}
 
279
  freemem(y, ns);
 
280
 freemem(sum, ns);
 
281
End; {decomp1}
 
282
 
 
283
Procedure solve(Var qr: ArbFloat; m, n, rwidthq: ArbInt; Var alpha: ArbFloat;
 
284
                Var pivot: ArbInt; Var r, y: ArbFloat);
 
285
 
 
286
Var    i, j, ii            : ArbInt;
 
287
       gamma, s            : ArbFloat;
 
288
       pqr, pal, pr, py, z : ^arfloat1;
 
289
       piv                 : ^arint1;
 
290
Begin
 
291
  pqr := @qr;
 
292
  pal := @alpha;
 
293
  piv := @pivot;
 
294
  pr := @r;
 
295
  py := @y;
 
296
  getmem(z, n*sizeof(ArbFloat));
 
297
  For j:=1 To n Do
 
298
    Begin
 
299
      gamma := 0;
 
300
      For i:=j To m Do
 
301
        gamma := gamma+pqr^[(i-1)*rwidthq+j]*pr^[i];
 
302
      gamma := gamma/(pal^[j]*pqr^[(j-1)*rwidthq+j]);
 
303
      For i:=j To m Do
 
304
        pr^[i] := pr^[i]+gamma*pqr^[(i-1)*rwidthq+j]
 
305
    End; {j}
 
306
  z^[n] := pr^[n]/pal^[n];
 
307
  For i:=n-1 Downto 1 Do
 
308
    Begin
 
309
      s := pr^[i];
 
310
      ii := (i-1)*rwidthq;
 
311
      For j:=i+1 To n Do
 
312
        s := s-pqr^[ii+j]*z^[j];
 
313
      z^[i] := s/pal^[i]
 
314
    End; {i}
 
315
  For i:=1 To n Do
 
316
    py^[piv^[i]] := z^[i];
 
317
  freemem(z, n*sizeof(ArbFloat));
 
318
End; {solve}
 
319
 
 
320
Procedure solve1(Var qr1; m, n: ArbInt; Var alpha1: ArbFloat;
 
321
                 Var pivot1: ArbInt; Var r1, y1: ArbFloat);
 
322
 
 
323
Var    i, j                : ArbInt;
 
324
       gamma, s            : ArbFloat;
 
325
       qr                  : ar2dr1 absolute qr1;
 
326
       alpha               : arfloat1 absolute alpha1;
 
327
       r                   : arfloat1 absolute r1;
 
328
       y                   : arfloat1 absolute y1;
 
329
       pivot               : arint1 absolute pivot1;
 
330
       z                   : ^arfloat1;
 
331
Begin
 
332
  getmem(z, n*sizeof(ArbFloat));
 
333
  For j:=1 To n Do
 
334
    Begin
 
335
      gamma := 0;
 
336
      For i:=j To m Do
 
337
       gamma := gamma+qr[i]^[j]*r[i];
 
338
      gamma := gamma/(alpha[j]*qr[j]^[j]);
 
339
      For i:=j To m Do
 
340
       r[i] := r[i]+gamma*qr[i]^[j]
 
341
    End; {j}
 
342
  z^[n] := r[n]/alpha[n];
 
343
  For i:=n-1 Downto 1 Do
 
344
    Begin
 
345
      s := r[i];
 
346
      For j:=i+1 To n Do
 
347
       s := s-qr[i]^[j]*z^[j];
 
348
      z^[i] := s/alpha[i]
 
349
    End; {i}
 
350
  For i:=1 To n Do
 
351
   y[pivot[i]] := z^[i];
 
352
  freemem(z, n*sizeof(ArbFloat));
 
353
End; {solve1}
 
354
 
 
355
Procedure sledtr(n: ArbInt; Var l, d, u, b, x: ArbFloat; Var term: ArbInt);
 
356
 
 
357
Var               i, j, sr : ArbInt;
 
358
                    lj, di : ArbFloat;
 
359
        pd, pu, pb, px, dd : ^arfloat1;
 
360
                        pl : ^arfloat2;
 
361
Begin
 
362
  If n<1
 
363
   Then
 
364
    Begin
 
365
      term := 3;
 
366
     exit
 
367
    End; {wrong input}
 
368
  pl := @l;
 
369
 pd := @d;
 
370
 pu := @u;
 
371
 pb := @b;
 
372
 px := @x;
 
373
  sr := sizeof(ArbFloat);
 
374
  getmem(dd, n*sr);
 
375
  move(pb^, px^, n*sr);
 
376
  j := 1;
 
377
 di := pd^[j];
 
378
 dd^[j] := di;
 
379
  If di=0
 
380
   Then
 
381
    term := 2
 
382
  Else
 
383
    term := 1;
 
384
  while (term=1) and (j <> n) Do
 
385
    Begin
 
386
      i := j;
 
387
     j := j+1;
 
388
     lj := pl^[j]/di;
 
389
      di := pd^[j]-lj*pu^[i];
 
390
     dd^[j] := di;
 
391
      If di=0
 
392
       Then
 
393
        term := 2
 
394
      Else
 
395
        px^[j] := px^[j]-lj*px^[i]
 
396
    End; {j}
 
397
  If term=1
 
398
   Then
 
399
    Begin
 
400
      px^[n] := px^[n]/dd^[n];
 
401
      For i:=n-1 Downto 1 Do
 
402
        px^[i] := (px^[i]-pu^[i]*px^[i+1])/dd^[i]
 
403
    End; {term=1}
 
404
  freemem(dd, n*sr);
 
405
End; {sledtr}
 
406
 
 
407
Procedure slegba(n, l, r: ArbInt;
 
408
                 Var a, b, x, ca: ArbFloat; Var term:ArbInt);
 
409
 
 
410
Var
 
411
  sr, i, j, k, ipivot, lbj, lbi, ubi, ls,
 
412
         ii, jj, ll, ubj, rwidth       : ArbInt;
 
413
  normr, sumrowi, pivot, normt, maxim, h  : ArbFloat;
 
414
  pa, pb, px, au, sumrow, t, row              : ^arfloat1;
 
415
Begin
 
416
  If (n<1) Or (l<0) Or (r<0) Or (l>n-1) Or (r>n-1)
 
417
   Then
 
418
    Begin
 
419
      term := 3;
 
420
     exit
 
421
    End; {term=3}
 
422
  sr := sizeof(ArbFloat);
 
423
  pa := @a;
 
424
 pb := @b;
 
425
 px := @x;
 
426
  ll := l+r+1;
 
427
  ls := ll*sr;
 
428
  getmem(au, ls*n);
 
429
  getmem(sumrow, n*sr);
 
430
  getmem(t, n*sr);
 
431
  getmem(row, ls);
 
432
  move(pb^, px^, n*sr);
 
433
  jj := 1;
 
434
 ii := 1;
 
435
  For i:=1 To n Do
 
436
    Begin
 
437
      If i <= l+1 Then
 
438
        Begin
 
439
          If i <= n-r Then rwidth := r+i
 
440
         Else rwidth := n
 
441
        End
 
442
     Else
 
443
          If i <= n-r Then rwidth := ll
 
444
     Else rwidth := n-i+l+1;
 
445
      move(pa^[jj], au^[ii], rwidth*sr);
 
446
      fillchar(au^[ii+rwidth], (ll-rwidth)*sr, 0);
 
447
      jj := jj+rwidth;
 
448
     ii := ii+ll;
 
449
    End; {i}
 
450
  lbi := n-r+1;
 
451
 lbj := 0;
 
452
  normr := 0;
 
453
 term := 1;
 
454
  ii := 1;
 
455
  For i:=1 To n Do
 
456
    Begin
 
457
      sumrowi := omvn1v(au^[ii], ll);
 
458
      ii := ii+ll;
 
459
      sumrow^[i] := sumrowi;
 
460
      h := 2*random-1;
 
461
     t^[i] := sumrowi*h;
 
462
      h := abs(h);
 
463
     If normr<h Then normr := h;
 
464
      If sumrowi=0 Then term := 2
 
465
    End; {i}
 
466
  ubi := l;
 
467
 k := 0;
 
468
 jj := 1;
 
469
  while (k<n) and (term=1) Do
 
470
    Begin
 
471
      maxim := 0;
 
472
     k := k+1;
 
473
     ipivot := k;
 
474
     ii := jj;
 
475
      If ubi<n
 
476
       Then ubi := ubi+1;
 
477
      For i:=k To ubi Do
 
478
        Begin
 
479
          sumrowi := sumrow^[i];
 
480
          If sumrowi <> 0
 
481
           Then
 
482
            Begin
 
483
              h := abs(au^[ii])/sumrowi;
 
484
              ii := ii+ll;
 
485
              If maxim<h
 
486
               Then
 
487
                Begin
 
488
                  maxim := h;
 
489
                 ipivot := i
 
490
                End {maxim<h}
 
491
            End {sumrowi <> 0}
 
492
        End; {i}
 
493
      If maxim=0
 
494
       Then
 
495
        term := 2
 
496
      Else
 
497
        Begin
 
498
          If ipivot <> k
 
499
           Then
 
500
            Begin
 
501
              ii := (ipivot-1)*ll+1;
 
502
              move(au^[ii], row^, ls);
 
503
              move(au^[jj], au^[ii], ls);
 
504
              move(row^, au^[jj], ls);
 
505
              h := t^[ipivot];
 
506
             t^[ipivot] := t^[k];
 
507
             t^[k] := h;
 
508
              h := px^[ipivot];
 
509
             px^[ipivot] := px^[k];
 
510
             px^[k] := h;
 
511
              sumrow^[ipivot] := sumrow^[k]
 
512
            End; {ipivot <> k}
 
513
          pivot := au^[jj];
 
514
         ii := jj;
 
515
          For i:=k+1 To ubi Do
 
516
            Begin
 
517
              ii := ii+ll;
 
518
              h := au^[ii]/pivot;
 
519
              For j:=0 To ll-2 Do
 
520
                au^[ii+j] := au^[ii+j+1]-h*au^[jj+j+1];
 
521
              au^[ii+ll-1] := 0;
 
522
              t^[i] := t^[i]-h*t^[k];
 
523
              px^[i] := px^[i]-h*px^[k];
 
524
            End {i}
 
525
        End; {maxim <> 0}
 
526
        jj := jj+ll
 
527
    End; {k}
 
528
  If term=1
 
529
   Then
 
530
    Begin
 
531
      normt := 0;
 
532
     ubj := -l-1;
 
533
      jj := n*ll+1;
 
534
      For i:=n Downto 1 Do
 
535
        Begin
 
536
          jj := jj-ll;
 
537
          If ubj<r
 
538
           Then
 
539
            ubj := ubj+1;
 
540
          h := t^[i];
 
541
          For j:=1 To ubj+l Do
 
542
            h := h-au^[jj+j]*t^[i+j];
 
543
          t^[i] := h/au^[jj];
 
544
          h := px^[i];
 
545
          For j:=1 To ubj+l Do
 
546
            h := h-au^[jj+j]*px^[i+j];
 
547
          px^[i] := h/au^[jj];
 
548
          h := abs(t^[i]);
 
549
          If normt<h
 
550
           Then
 
551
            normt := h
 
552
        End; {i}
 
553
        ca := normt/normr
 
554
    End; {term=1}
 
555
  freemem(au, ls*n);
 
556
  freemem(sumrow, n*sr);
 
557
  freemem(t, n*sr);
 
558
  freemem(row, ls)
 
559
End; {slegba}
 
560
 
 
561
Procedure slegbal(n, l, r: ArbInt;
 
562
                  Var a1; Var b1, x1, ca: ArbFloat; Var term:ArbInt);
 
563
 
 
564
Var
 
565
  sr, i, j, k, ipivot, ubi, ls,
 
566
                 ll, ubj, rwidth       : ArbInt;
 
567
  normr, sumrowi, pivot, normt, maxim, h  : ArbFloat;
 
568
  a                                           : ar2dr1 absolute a1;
 
569
  b                                           : arfloat1 absolute b1;
 
570
  x                                           : arfloat1 absolute x1;
 
571
  au                                          : par2dr1;
 
572
  sumrow, t, row                              : ^arfloat1;
 
573
Begin
 
574
  If (n<1) Or (l<0) Or (r<0) Or (l>n-1) Or (r>n-1)
 
575
   Then
 
576
    Begin
 
577
      term := 3;
 
578
     exit
 
579
    End; {term=3}
 
580
  sr := sizeof(ArbFloat);
 
581
 ll := l+r+1;
 
582
 ls := ll*sr;
 
583
  AllocateAr2dr(n, ll, au);
 
584
  getmem(sumrow, n*sr);
 
585
 getmem(t, n*sr);
 
586
 getmem(row, ls);
 
587
  move(b[1], x[1], n*sr);
 
588
  For i:=1 To n Do
 
589
    Begin
 
590
      If i <= l+1 Then
 
591
        Begin
 
592
          If i <= n-r Then rwidth := r+i
 
593
         Else rwidth := n
 
594
        End
 
595
     Else
 
596
          If i <= n-r Then rwidth := ll
 
597
     Else rwidth := n-i+l+1;
 
598
      move(a[i]^, au^[i]^, rwidth*sr);
 
599
      fillchar(au^[i]^[rwidth+1], (ll-rwidth)*sr, 0);
 
600
    End; {i}
 
601
  normr := 0;
 
602
 term := 1;
 
603
  For i:=1 To n Do
 
604
    Begin
 
605
      sumrowi := omvn1v(au^[i]^[1], ll);
 
606
     sumrow^[i] := sumrowi;
 
607
      h := 2*random-1;
 
608
     t^[i] := sumrowi*h;
 
609
      h := abs(h);
 
610
     If normr<h Then normr := h;
 
611
      If sumrowi=0 Then term := 2
 
612
    End; {i}
 
613
  ubi := l;
 
614
 k := 0;
 
615
  while (k<n) and (term=1) Do
 
616
    Begin
 
617
      maxim := 0;
 
618
     k := k+1;
 
619
     ipivot := k;
 
620
      If ubi<n Then ubi := ubi+1;
 
621
      For i:=k To ubi Do
 
622
        Begin
 
623
          sumrowi := sumrow^[i];
 
624
          If sumrowi <> 0 Then
 
625
            Begin
 
626
              h := abs(au^[i]^[1])/sumrowi;
 
627
              If maxim<h Then
 
628
                Begin
 
629
                  maxim := h;
 
630
                 ipivot := i
 
631
                End {maxim<h}
 
632
            End {sumrowi <> 0}
 
633
        End; {i}
 
634
      If maxim=0 Then term := 2
 
635
     Else
 
636
        Begin
 
637
          If ipivot <> k Then
 
638
            Begin
 
639
              move(au^[ipivot]^, row^, ls);
 
640
              move(au^[k]^, au^[ipivot]^, ls);
 
641
              move(row^, au^[k]^, ls);
 
642
              h := t^[ipivot];
 
643
             t^[ipivot] := t^[k];
 
644
             t^[k] := h;
 
645
              h := x[ipivot];
 
646
             x[ipivot] := x[k];
 
647
             x[k] := h;
 
648
              sumrow^[ipivot] := sumrow^[k]
 
649
            End; {ipivot <> k}
 
650
          pivot := au^[k]^[1];
 
651
          For i:=k+1 To ubi Do
 
652
            Begin
 
653
              h := au^[i]^[1]/pivot;
 
654
              For j:=0 To ll-2 Do
 
655
                au^[i]^[j+1] := au^[i]^[j+2]-h*au^[k]^[j+2];
 
656
              au^[i]^[ll] := 0;
 
657
              t^[i] := t^[i]-h*t^[k];
 
658
              x[i] := x[i]-h*x[k];
 
659
            End {i}
 
660
        End; {maxim <> 0}
 
661
    End; {k}
 
662
  If term=1 Then
 
663
    Begin
 
664
      normt := 0;
 
665
     ubj := -l-1;
 
666
      For i:=n Downto 1 Do
 
667
        Begin
 
668
          If ubj<r Then ubj := ubj+1;
 
669
          h := t^[i];
 
670
          For j:=1 To ubj+l Do
 
671
           h := h-au^[i]^[j+1]*t^[i+j];
 
672
          t^[i] := h/au^[i]^[1];
 
673
          h := x[i];
 
674
          For j:=1 To ubj+l Do
 
675
           h := h-au^[i]^[j+1]*x[i+j];
 
676
          x[i] := h/au^[i]^[1];
 
677
          h := abs(t^[i]);
 
678
         If normt<h Then normt := h
 
679
        End; {i}
 
680
        ca := normt/normr
 
681
    End; {term=1}
 
682
  freemem(sumrow, n*sr);
 
683
 freemem(t, n*sr);
 
684
 freemem(row, ls);
 
685
  DeAllocateAr2dr(n, ll, au);
 
686
End; {slegbal}
 
687
 
 
688
Procedure slegen(n, rwidth: ArbInt; Var a, b, x, ca: ArbFloat;
 
689
                 Var term: ArbInt);
 
690
 
 
691
Var
 
692
          nsr, i, j, k, ipiv, ip, ik, i1n, k1n : ArbInt;
 
693
                                      singular : boolean;
 
694
           normr, pivot, l, normt, maxim, h, s : ArbFloat;
 
695
                pa, px, pb, au, sumrow, t, row : ^arfloat1;
 
696
 
 
697
Begin
 
698
  If (n<1) Or (rwidth<1)
 
699
   Then
 
700
    Begin
 
701
      term := 3;
 
702
     exit
 
703
    End; {wrong input}
 
704
  getmem(au, sqr(n)*sizeof(ArbFloat));
 
705
  nsr := n*sizeof(ArbFloat);
 
706
  getmem(t, nsr);
 
707
  getmem(row, nsr);
 
708
  getmem(sumrow, nsr);
 
709
  pa := @a;
 
710
 pb := @b;
 
711
 px := @x;
 
712
  For i:= 1 To n Do
 
713
    move(pa^[1+(i-1)*rwidth], au^[1+(i-1)*n], nsr);
 
714
  move(pb^[1], px^[1], nsr);
 
715
  normr := 0;
 
716
 singular := false ;
 
717
 i := 0;
 
718
 j := 0;
 
719
  while (i<n) and  (Not singular) Do
 
720
    Begin
 
721
      i := i+1;
 
722
     sumrow^[i] := omvn1v(au^[1+(i-1)*n], n);
 
723
      If sumrow^[i]=0
 
724
       Then
 
725
        singular := true
 
726
      Else
 
727
        Begin
 
728
          h := 2*random-1;
 
729
         t^[i] := sumrow^[i]*h;
 
730
         h := abs(h);
 
731
          If normr<h
 
732
           Then
 
733
            normr := h
 
734
        End
 
735
    End;
 
736
  k := 0;
 
737
  while (k<n) and  not singular Do
 
738
    Begin
 
739
      k := k+1;
 
740
     maxim := 0;
 
741
     ipiv := k;
 
742
      For i:=k To n Do
 
743
        Begin
 
744
          h := abs(au^[k+(i-1)*n])/sumrow^[i];
 
745
          If maxim<h
 
746
           Then
 
747
            Begin
 
748
              maxim := h;
 
749
             ipiv := i
 
750
            End
 
751
        End;
 
752
      If maxim=0
 
753
       Then
 
754
        singular := true
 
755
      Else
 
756
        Begin
 
757
          k1n := (k-1)*n;
 
758
          If ipiv <> k
 
759
           Then
 
760
            Begin
 
761
              ip := 1+(ipiv-1)*n;
 
762
             ik := 1+k1n;
 
763
              move(au^[ip], row^[1], nsr);
 
764
             move(au^[ik], au^[ip], nsr);
 
765
              move(row^[1], au^[ik], nsr);
 
766
              h := t^[ipiv];
 
767
             t^[ipiv] := t^[k];
 
768
             t^[k] := h;
 
769
              h := px^[ipiv];
 
770
             px^[ipiv] := px^[k];
 
771
             px^[k] := h;
 
772
              sumrow^[ipiv] := sumrow^[k]
 
773
            End;
 
774
          pivot := au^[k+k1n];
 
775
          For i:=k+1 To n Do
 
776
            Begin
 
777
              i1n := (i-1)*n;
 
778
             l := au^[k+i1n]/pivot;
 
779
              If l <> 0
 
780
               Then
 
781
                Begin
 
782
                  For j:=k+1 To n Do
 
783
                    au^[j+i1n] := au^[j+i1n]-l*au^[j+k1n];
 
784
                  t^[i] := t^[i]-l*t^[k];
 
785
                  px^[i] := px^[i]-l*px^[k]
 
786
                End
 
787
            End
 
788
        End
 
789
    End;
 
790
  If  Not singular
 
791
   Then
 
792
    Begin
 
793
      normt := 0;
 
794
      For i:=n Downto 1 Do
 
795
        Begin
 
796
          s := 0;
 
797
         i1n := (i-1)*n;
 
798
          For j:=i+1 To n Do
 
799
            s := s+t^[j]*au^[j+i1n];
 
800
          t^[i] := (t^[i]-s)/au^[i+i1n];
 
801
          s := 0;
 
802
          For j:=i+1 To n Do
 
803
            s := s+px^[j]*au^[j+i1n];
 
804
          px^[i] := (px^[i]-s)/au^[i+i1n];
 
805
          h := abs(t^[i]);
 
806
          If normt<h
 
807
           Then
 
808
            normt := h
 
809
        End;
 
810
      ca := normt/normr
 
811
    End;
 
812
   If singular
 
813
    Then
 
814
     term := 2
 
815
   Else
 
816
     term := 1;
 
817
  freemem(au, sqr(n)*sizeof(ArbFloat));
 
818
  freemem(t, nsr);
 
819
  freemem(row, nsr);
 
820
  freemem(sumrow, nsr);
 
821
End; {slegen}
 
822
 
 
823
Procedure slegenl(     n: ArbInt;
 
824
                  Var a1;
 
825
                  Var b1, x1, ca: ArbFloat;
 
826
                  Var term: ArbInt);
 
827
 
 
828
Var
 
829
     nsr, i, j, k, ipiv : ArbInt;
 
830
               singular : boolean;
 
831
     normr, pivot, l, normt, maxim, h, s : ArbFloat;
 
832
     a : ar2dr1 absolute a1;
 
833
     x : arfloat1 absolute x1;
 
834
     b : arfloat1 absolute b1;
 
835
     au: par2dr1;
 
836
     sumrow, t, row : ^arfloat1;
 
837
Begin
 
838
  If n<1 Then
 
839
    Begin
 
840
      term := 3;
 
841
     exit
 
842
    End; {wrong input}
 
843
  AllocateAr2dr(n, n, au);
 
844
  nsr := n*sizeof(ArbFloat);
 
845
  getmem(t, nsr);
 
846
  getmem(row, nsr);
 
847
  getmem(sumrow, nsr);
 
848
  For i:= 1 To n Do
 
849
   move(a[i]^, au^[i]^, nsr);
 
850
  move(b[1], x[1], nsr);
 
851
  normr := 0;
 
852
 singular := false ;
 
853
 i := 0;
 
854
 j := 0;
 
855
  while (i<n) and  (Not singular) Do
 
856
    Begin
 
857
      i := i+1;
 
858
     sumrow^[i] := omvn1v(au^[i]^[1], n);
 
859
      If sumrow^[i]=0
 
860
       Then
 
861
        singular := true
 
862
      Else
 
863
        Begin
 
864
          h := 2*random-1;
 
865
         t^[i] := sumrow^[i]*h;
 
866
         h := abs(h);
 
867
          If normr<h
 
868
           Then
 
869
            normr := h
 
870
        End
 
871
    End;
 
872
  k := 0;
 
873
  while (k<n) and  not singular Do
 
874
    Begin
 
875
      k := k+1;
 
876
     maxim := 0;
 
877
     ipiv := k;
 
878
      For i:=k To n Do
 
879
        Begin
 
880
          h := abs(au^[i]^[k])/sumrow^[i];
 
881
          If maxim<h
 
882
           Then
 
883
            Begin
 
884
              maxim := h;
 
885
             ipiv := i
 
886
            End
 
887
        End;
 
888
      If maxim=0
 
889
       Then
 
890
        singular := true
 
891
      Else
 
892
        Begin
 
893
          If ipiv <> k
 
894
           Then
 
895
            Begin
 
896
              move(au^[ipiv]^, row^, nsr);
 
897
              move(au^[k]^, au^[ipiv]^, nsr);
 
898
              move(row^, au^[k]^, nsr);
 
899
              h := t^[ipiv];
 
900
             t^[ipiv] := t^[k];
 
901
             t^[k] := h;
 
902
              h := x[ipiv];
 
903
             x[ipiv] := x[k];
 
904
             x[k] := h;
 
905
              sumrow^[ipiv] := sumrow^[k]
 
906
            End;
 
907
          pivot := au^[k]^[k];
 
908
          For i:=k+1 To n Do
 
909
            Begin
 
910
              l := au^[i]^[k]/pivot;
 
911
              If l <> 0
 
912
               Then
 
913
                Begin
 
914
                  For j:=k+1 To n Do
 
915
                    au^[i]^[j] := au^[i]^[j]-l*au^[k]^[j];
 
916
                  t^[i] := t^[i]-l*t^[k];
 
917
                  x[i] := x[i]-l*x[k]
 
918
                End
 
919
            End
 
920
        End
 
921
    End;
 
922
  If  Not singular
 
923
   Then
 
924
    Begin
 
925
      normt := 0;
 
926
      For i:=n Downto 1 Do
 
927
        Begin
 
928
          s := 0;
 
929
          For j:=i+1 To n Do
 
930
            s := s+t^[j]*au^[i]^[j];
 
931
          t^[i] := (t^[i]-s)/au^[i]^[i];
 
932
          s := 0;
 
933
          For j:=i+1 To n Do
 
934
            s := s+x[j]*au^[i]^[j];
 
935
          x[i] := (x[i]-s)/au^[i]^[i];
 
936
          h := abs(t^[i]);
 
937
          If normt<h
 
938
           Then
 
939
            normt := h
 
940
        End;
 
941
      ca := normt/normr
 
942
    End;
 
943
   If singular
 
944
    Then
 
945
     term := 2
 
946
   Else
 
947
     term := 1;
 
948
  freemem(t, nsr);
 
949
  freemem(row, nsr);
 
950
  freemem(sumrow, nsr);
 
951
  DeAllocateAr2dr(n, n, au);
 
952
End; {slegenl}
 
953
 
 
954
Procedure slegls(Var a: ArbFloat; m, n, rwidtha: ArbInt; Var b, x: ArbFloat;
 
955
                 Var term: ArbInt);
 
956
 
 
957
Var     i, j, ns, ms, ii                : ArbInt;
 
958
        normy0, norme1, s       : ArbFloat;
 
959
        pa, pb, px, qr, alpha, e, y, r  : ^arfloat1;
 
960
        pivot                           : ^arint1;
 
961
Begin
 
962
  If (n<1) Or (m<n)
 
963
   Then
 
964
    Begin
 
965
      term := 3;
 
966
     exit
 
967
    End;
 
968
  pa := @a;
 
969
 pb := @b;
 
970
 px := @x;
 
971
  ns := n*sizeof(ArbFloat);
 
972
 ms := m*sizeof(ArbFloat);
 
973
  getmem(qr, m*ns);
 
974
 getmem(alpha, ns);
 
975
 getmem(e, ns);
 
976
 getmem(y, ns);
 
977
  getmem(r, m*sizeof(ArbFloat));
 
978
 getmem(pivot, n*sizeof(ArbInt));
 
979
  For i:=1 To m Do
 
980
    move(pa^[(i-1)*rwidtha+1], qr^[(i-1)*n+1], ns);
 
981
  decomp(qr^[1], m, n, n, alpha^[1], pivot^[1], term);
 
982
  If term=2
 
983
   Then
 
984
    Begin
 
985
      freemem(qr, m*ns);
 
986
     freemem(alpha, ns);
 
987
     freemem(e, ns);
 
988
     freemem(y, ns);
 
989
      freemem(r, m*sizeof(ArbFloat));
 
990
     freemem(pivot, n*sizeof(ArbInt));
 
991
      exit
 
992
    End;
 
993
  move(pb^[1], r^[1], ms);
 
994
  solve(qr^[1], m, n, n, alpha^[1], pivot^[1], r^[1], y^[1]);
 
995
  For i:=1 To m Do
 
996
    Begin
 
997
      s := pb^[i];
 
998
     ii := (i-1)*rwidtha;
 
999
      For j:=1 To n Do
 
1000
        s := s-pa^[ii+j]*y^[j];
 
1001
      r^[i] := s
 
1002
    End; {i}
 
1003
  solve(qr^[1], m, n, n, alpha^[1], pivot^[1], r^[1], e^[1]);
 
1004
  normy0 := 0;
 
1005
 norme1 := 0;
 
1006
  For i:=1 To n Do
 
1007
    Begin
 
1008
      normy0 := normy0+sqr(y^[i]);
 
1009
     norme1 := norme1+sqr(e^[i])
 
1010
    End; {i}
 
1011
  If norme1 > 0.0625*normy0
 
1012
   Then
 
1013
    Begin
 
1014
      term := 2;
 
1015
      freemem(qr, m*ns);
 
1016
     freemem(alpha, ns);
 
1017
     freemem(e, ns);
 
1018
     freemem(y, ns);
 
1019
      freemem(r, m*sizeof(ArbFloat));
 
1020
     freemem(pivot, n*sizeof(ArbInt));
 
1021
      exit
 
1022
    End;
 
1023
  For i:=1 To n Do
 
1024
    px^[i] := y^[i];
 
1025
  freemem(qr, m*ns);
 
1026
 freemem(alpha, ns);
 
1027
 freemem(e, ns);
 
1028
 freemem(y, ns);
 
1029
  freemem(r, m*sizeof(ArbFloat));
 
1030
 freemem(pivot, n*sizeof(ArbInt));
 
1031
End; {slegls}
 
1032
 
 
1033
Procedure sleglsl(Var a1; m, n: ArbInt; Var b1, x1: ArbFloat;
 
1034
                  Var term: ArbInt);
 
1035
 
 
1036
Var     i, j, ns, ms                    : ArbInt;
 
1037
        normy0, norme1, s       : ArbFloat;
 
1038
        a                               : ar2dr1 absolute a1;
 
1039
        b                               : arfloat1 absolute b1;
 
1040
        x                               : arfloat1 absolute x1;
 
1041
        alpha, e, y, r                  : ^arfloat1;
 
1042
        qr                              : par2dr1;
 
1043
        pivot                           : ^arint1;
 
1044
Begin
 
1045
  If (n<1) Or (m<n)
 
1046
   Then
 
1047
    Begin
 
1048
      term := 3;
 
1049
     exit
 
1050
    End;
 
1051
  AllocateAr2dr(m, n, qr);
 
1052
  ns := n*sizeof(ArbFloat);
 
1053
 ms := m*sizeof(ArbFloat);
 
1054
  getmem(alpha, ns);
 
1055
 getmem(e, ns);
 
1056
 getmem(y, ns);
 
1057
  getmem(r, ms);
 
1058
 getmem(pivot, n*sizeof(ArbInt));
 
1059
  For i:=1 To m Do
 
1060
    move(a[i]^, qr^[i]^, ns);
 
1061
  decomp1(qr^[1], m, n, alpha^[1], pivot^[1], term);
 
1062
  If term=2
 
1063
   Then
 
1064
    Begin
 
1065
      freemem(qr, m*ns);
 
1066
     freemem(alpha, ns);
 
1067
     freemem(e, ns);
 
1068
     freemem(y, ns);
 
1069
      freemem(r, ms);
 
1070
     freemem(pivot, n*sizeof(ArbInt));
 
1071
      exit
 
1072
    End;
 
1073
  move(b[1], r^, ms);
 
1074
  solve1(qr^[1], m, n, alpha^[1], pivot^[1], r^[1], y^[1]);
 
1075
  For i:=1 To m Do
 
1076
    Begin
 
1077
      s := b[i];
 
1078
      For j:=1 To n Do
 
1079
       s := s-a[i]^[j]*y^[j];
 
1080
      r^[i] := s
 
1081
    End; {i}
 
1082
  solve1(qr^[1], m, n, alpha^[1], pivot^[1], r^[1], e^[1]);
 
1083
  normy0 := 0;
 
1084
 norme1 := 0;
 
1085
  For i:=1 To n Do
 
1086
    Begin
 
1087
      normy0 := normy0+sqr(y^[i]);
 
1088
     norme1 := norme1+sqr(e^[i])
 
1089
    End; {i}
 
1090
  If norme1 > 0.0625*normy0
 
1091
   Then
 
1092
    Begin
 
1093
      term := 2;
 
1094
      freemem(qr, m*ns);
 
1095
     freemem(alpha, ns);
 
1096
     freemem(e, ns);
 
1097
     freemem(y, ns);
 
1098
      freemem(r, m*sizeof(ArbFloat));
 
1099
     freemem(pivot, n*sizeof(ArbInt));
 
1100
      exit
 
1101
    End;
 
1102
  For i:=1 To n Do
 
1103
   x[i] := y^[i];
 
1104
  freemem(alpha, ns);
 
1105
 freemem(e, ns);
 
1106
 freemem(y, ns);
 
1107
  freemem(r, ms);
 
1108
 freemem(pivot, n*sizeof(ArbInt));
 
1109
  DeAllocateAr2dr(m, n, qr);
 
1110
End; {sleglsl}
 
1111
 
 
1112
Procedure slegpb(n, l: ArbInt; Var a, b, x, ca: ArbFloat;
 
1113
                 Var term: ArbInt);
 
1114
 
 
1115
Var
 
1116
    posdef                                                 : boolean;
 
1117
    i, j, k, r, p, q, jmin1, ii, jj, ri, ind,
 
1118
                                      ll, llm1, sr, rwidth : ArbInt;
 
1119
    h, normr, normt, sumrowi, hh, alim, norma              : ArbFloat;
 
1120
    pa, pb, px, al, t, v                                   : ^arfloat1;
 
1121
 
 
1122
    Procedure decomp(i, r: ArbInt);
 
1123
 
 
1124
    Var k: ArbInt;
 
1125
    Begin
 
1126
      ri := (r-1)*ll;
 
1127
      h := al^[ii+j];
 
1128
     q := ll-j+p;
 
1129
      For k:=p To jmin1 Do
 
1130
        Begin
 
1131
          h := h-al^[ii+k]*al^[ri+q];
 
1132
         q := q+1
 
1133
        End ;
 
1134
      If j<ll
 
1135
       Then
 
1136
        al^[ii+j] := h/al^[ri+ll];
 
1137
    End; {decomp}
 
1138
 
 
1139
Begin
 
1140
  If (n<1) Or (l<0) Or (l>n-1)
 
1141
   Then
 
1142
    Begin
 
1143
      term := 3;
 
1144
     exit
 
1145
    End; {wrong input}
 
1146
  sr := sizeof(ArbFloat);
 
1147
  pa := @a;
 
1148
  pb := @b;
 
1149
  px := @x;
 
1150
  ll := l+1;
 
1151
  getmem(al, ll*n*sr);
 
1152
  getmem(t, n*sr);
 
1153
  getmem(v, ll*sr);
 
1154
  move(pb^, px^, n*sr);
 
1155
  jj := 1;
 
1156
  ii := 1;
 
1157
  For i:=1 To n Do
 
1158
    Begin
 
1159
      If i>l Then rwidth := ll
 
1160
     Else rwidth := i;
 
1161
      move(pa^[jj], al^[ii+ll-rwidth], rwidth*sr);
 
1162
      jj := jj+rwidth;
 
1163
     ii := ii+ll
 
1164
    End; {i}
 
1165
  normr := 0;
 
1166
 p := ll+1;
 
1167
 norma := 0;
 
1168
  For i:=1 To n Do
 
1169
    Begin
 
1170
      If p>1
 
1171
       Then
 
1172
        p := p-1;
 
1173
      For j:=p To ll Do
 
1174
        v^[j] := al^[j+(i-1)*ll];
 
1175
      sumrowi := omvn1v(v^[p], ll-p+1);
 
1176
      r := i;
 
1177
     j := ll;
 
1178
      while (r<n) and (j>1) Do
 
1179
        Begin
 
1180
          r := r+1;
 
1181
         j := j-1;
 
1182
          sumrowi := sumrowi+abs(al^[j+(r-1)*ll])
 
1183
        End; {r,j}
 
1184
      If norma<sumrowi
 
1185
       Then
 
1186
        norma := sumrowi;
 
1187
      h := 2*random-1;
 
1188
     t^[i] := h;
 
1189
      h := abs(h);
 
1190
      If normr<h
 
1191
       Then
 
1192
        normr := h
 
1193
    End; {i}
 
1194
  llm1 := ll-1;
 
1195
 p := ll+1;
 
1196
 i := 0;
 
1197
 posdef := true ;
 
1198
  while (i<n) and posdef Do
 
1199
    Begin
 
1200
      i := i+1;
 
1201
     If p>1 Then p := p-1;
 
1202
     r := i-ll+p;
 
1203
     j := p-1;
 
1204
      ii := (i-1)*ll;
 
1205
      while j<llm1 Do
 
1206
        Begin
 
1207
          jmin1 := j;
 
1208
         j := j+1;
 
1209
          decomp(i, r);
 
1210
          r := r+1
 
1211
        End ; {j}
 
1212
      jmin1 := llm1;
 
1213
     j := ll;
 
1214
      decomp(i, i);
 
1215
      If h <= 0
 
1216
       Then
 
1217
        posdef := false
 
1218
      Else
 
1219
        Begin
 
1220
          alim := sqrt(h);
 
1221
         al^[ii+ll] := alim;
 
1222
          h := t^[i];
 
1223
         q := i;
 
1224
          For k:=llm1 Downto p Do
 
1225
            Begin
 
1226
              q := q-1;
 
1227
             h := h-al^[ii+k]*t^[q]
 
1228
            End ;
 
1229
          t^[i] := h/alim;
 
1230
          h := px^[i];
 
1231
         q := i;
 
1232
          For k:=llm1 Downto p Do
 
1233
            Begin
 
1234
              q := q-1;
 
1235
             h := h-al^[ii+k]*px^[q]
 
1236
            End; {k}
 
1237
          px^[i] := h/alim
 
1238
        End {posdef}
 
1239
    End; {i}
 
1240
    If posdef
 
1241
     Then
 
1242
      Begin
 
1243
        normt := 0;
 
1244
       p := ll+1;
 
1245
        For i:=n Downto 1 Do
 
1246
          Begin
 
1247
            If p>1
 
1248
             Then
 
1249
              p := p-1;
 
1250
            q := i;
 
1251
           h := t^[i];
 
1252
           hh := px^[i];
 
1253
            For k:=llm1 Downto p Do
 
1254
              Begin
 
1255
                q := q+1;
 
1256
                ind := (q-1)*ll+k;
 
1257
                h := h-al^[ind]*t^[q];
 
1258
               hh := hh-al^[ind]*px^[q]
 
1259
              End; {k}
 
1260
            ind := i*ll;
 
1261
            t^[i] := h/al^[ind];
 
1262
           px^[i] := hh/al^[ind];
 
1263
            h := abs(t^[i]);
 
1264
            If normt<h
 
1265
             Then
 
1266
              normt := h
 
1267
         End; {i}
 
1268
       ca := norma*normt/normr
 
1269
     End ; {posdef}
 
1270
  If posdef
 
1271
   Then
 
1272
    term := 1
 
1273
  Else
 
1274
    term := 2;
 
1275
  freemem(al, ll*n*sr);
 
1276
  freemem(t, n*sr);
 
1277
  freemem(v, ll*sr);
 
1278
End;  {slegpb}
 
1279
 
 
1280
Procedure slegpbl(n, l: ArbInt;
 
1281
                  Var a1; Var b1, x1, ca: ArbFloat; Var term: ArbInt);
 
1282
 
 
1283
Var
 
1284
    posdef                                    : boolean;
 
1285
    i, j, k, r, p, q, ll, sr, rwidth          : ArbInt;
 
1286
    h, normr, normt, sumrowi, hh, alim, norma : ArbFloat;
 
1287
    a                                         : ar2dr1 absolute a1;
 
1288
    b                                         : arfloat1 absolute b1;
 
1289
    x                                         : arfloat1 absolute x1;
 
1290
    al                                        : par2dr1;
 
1291
    t, v                                      : ^arfloat1;
 
1292
 
 
1293
    Procedure decomp(r: ArbInt);
 
1294
 
 
1295
    Var k: ArbInt;
 
1296
    Begin
 
1297
      h := al^[i]^[j];
 
1298
     q := ll-j+p;
 
1299
      For k:=p To j-1 Do
 
1300
        Begin
 
1301
          h := h-al^[i]^[k]*al^[r]^[q];
 
1302
         Inc(q)
 
1303
        End ;
 
1304
      If j<ll Then al^[i]^[j] := h/al^[r]^[ll];
 
1305
    End; {decomp}
 
1306
 
 
1307
Begin
 
1308
  If (n<1) Or (l<0) Or (l>n-1)
 
1309
   Then
 
1310
    Begin
 
1311
      term := 3;
 
1312
     exit
 
1313
    End; {wrong input}
 
1314
  sr := sizeof(ArbFloat);
 
1315
  ll := l+1;
 
1316
  AllocateAr2dr(n, ll, al);
 
1317
  getmem(t, n*sr);
 
1318
 getmem(v, ll*sr);
 
1319
  move(b[1], x[1], n*sr);
 
1320
  For i:=1 To n Do
 
1321
    Begin
 
1322
      If i>l Then rwidth := ll
 
1323
     Else rwidth := i;
 
1324
      move(a[i]^, al^[i]^[ll-rwidth+1], rwidth*sr);
 
1325
    End; {i}
 
1326
  normr := 0;
 
1327
 p := ll+1;
 
1328
 norma := 0;
 
1329
  For i:=1 To n Do
 
1330
    Begin
 
1331
      If p>1 Then Dec(p);
 
1332
      For j:=p To ll Do
 
1333
       v^[j] := al^[i]^[j];
 
1334
      sumrowi := omvn1v(v^[p], ll-p+1);
 
1335
      r := i;
 
1336
     j := ll;
 
1337
      while (r<n) and (j>1) Do
 
1338
        Begin
 
1339
          Inc(r);
 
1340
         Dec(j);
 
1341
          sumrowi := sumrowi+abs(al^[r]^[j])
 
1342
        End; {r,j}
 
1343
      If norma<sumrowi Then norma := sumrowi;
 
1344
      h := 2*random-1;
 
1345
     t^[i] := h;
 
1346
      h := abs(h);
 
1347
     If normr<h Then normr := h
 
1348
    End; {i}
 
1349
  p := ll+1;
 
1350
 i := 0;
 
1351
 posdef := true ;
 
1352
  while (i<n) and posdef Do
 
1353
    Begin
 
1354
      Inc(i);
 
1355
     If p>1 Then Dec(p);
 
1356
     r := i-ll+p;
 
1357
     j := p-1;
 
1358
      while j<ll-1 Do
 
1359
        Begin
 
1360
          Inc(j);
 
1361
         decomp(r);
 
1362
         Inc(r)
 
1363
        End ; {j}
 
1364
      j := ll;
 
1365
     decomp(i);
 
1366
      If h <= 0 Then posdef := false
 
1367
     Else
 
1368
        Begin
 
1369
          alim := sqrt(h);
 
1370
         al^[i]^[ll] := alim;
 
1371
          h := t^[i];
 
1372
         q := i;
 
1373
          For k:=ll-1 Downto p Do
 
1374
            Begin
 
1375
              q := q-1;
 
1376
             h := h-al^[i]^[k]*t^[q]
 
1377
            End ;
 
1378
          t^[i] := h/alim;
 
1379
          h := x[i];
 
1380
         q := i;
 
1381
          For k:=ll-1 Downto p Do
 
1382
            Begin
 
1383
              q := q-1;
 
1384
             h := h-al^[i]^[k]*x[q]
 
1385
            End; {k}
 
1386
          x[i] := h/alim
 
1387
        End {posdef}
 
1388
    End; {i}
 
1389
    If posdef
 
1390
     Then
 
1391
      Begin
 
1392
        normt := 0;
 
1393
       p := ll+1;
 
1394
        For i:=n Downto 1 Do
 
1395
          Begin
 
1396
            If p>1 Then Dec(p);
 
1397
            q := i;
 
1398
           h := t^[i];
 
1399
           hh := x[i];
 
1400
            For k:=ll-1 Downto p Do
 
1401
              Begin
 
1402
                Inc(q);
 
1403
                h := h-al^[q]^[k]*t^[q];
 
1404
               hh := hh-al^[q]^[k]*x[q]
 
1405
              End; {k}
 
1406
            t^[i] := h/al^[i]^[ll];
 
1407
           x[i] := hh/al^[i]^[ll];
 
1408
            h := abs(t^[i]);
 
1409
           If normt<h Then normt := h
 
1410
         End; {i}
 
1411
       ca := norma*normt/normr
 
1412
     End ; {posdef}
 
1413
  If posdef Then term := 1
 
1414
 Else term := 2;
 
1415
  freemem(t, n*sr);
 
1416
 freemem(v, ll*sr);
 
1417
  DeAllocateAr2dr(n, ll, al);
 
1418
End;  {slegpbl}
 
1419
 
 
1420
Procedure slegpd(n, rwidth: ArbInt; Var a, b, x, ca: ArbFloat;
 
1421
                 Var term: ArbInt);
 
1422
 
 
1423
Var
 
1424
    sr, i, j, k, kmin1, kk, k1n, i1n, ik, ii : ArbInt;
 
1425
                                          pd : boolean;
 
1426
        h, lkk, normr, normt, sumrowi, norma : ArbFloat;
 
1427
                           pa, pb, px, al, t : ^arfloat1;
 
1428
 
 
1429
Begin
 
1430
  If (n<1) Or (rwidth<1)
 
1431
   Then
 
1432
    Begin
 
1433
      term := 3;
 
1434
     exit
 
1435
    End;
 
1436
  sr := sizeof(ArbFloat);
 
1437
  getmem(al, sqr(n)*sr);
 
1438
 getmem(t, n*sr);
 
1439
  pa := @a;
 
1440
 pb := @b;
 
1441
 px := @x;
 
1442
  For i:=1 To n Do
 
1443
    move(pa^[1+(i-1)*rwidth], al^[1+(i-1)*n], i*sr);
 
1444
  move(pb^[1], px^[1], n*sr);
 
1445
  normr := 0;
 
1446
 pd := true ;
 
1447
 norma := 0;
 
1448
  For i:=1 To n Do
 
1449
    Begin
 
1450
      sumrowi := 0;
 
1451
      For j:=1 To i Do
 
1452
        sumrowi := sumrowi+abs(al^[j+(i-1)*n]);
 
1453
      For j:=i+1 To n Do
 
1454
        sumrowi := sumrowi+abs(al^[i+(j-1)*n]);
 
1455
      If norma<sumrowi
 
1456
       Then
 
1457
        norma := sumrowi;
 
1458
      t^[i] := 2*random-1;
 
1459
     h := abs(t^[i]);
 
1460
      If normr<h
 
1461
       Then
 
1462
        normr := h
 
1463
    End; {i}
 
1464
  k := 0;
 
1465
  while (k<n) and pd Do
 
1466
    Begin
 
1467
      kmin1 := k;
 
1468
     k := k+1;
 
1469
     k1n := (k-1)*n;
 
1470
     kk := k+k1n;
 
1471
     lkk := al^[kk];
 
1472
      For j:=1 To kmin1 Do
 
1473
        lkk := lkk-sqr(al^[j+k1n]);
 
1474
      If lkk<=0
 
1475
       Then
 
1476
        pd := false
 
1477
      Else
 
1478
        Begin
 
1479
          al^[kk] := sqrt(lkk);
 
1480
         lkk := al^[kk];
 
1481
          For i:=k+1 To n Do
 
1482
            Begin
 
1483
              i1n := (i-1)*n;
 
1484
             ik := k+i1n;
 
1485
             h := al^[ik];
 
1486
              For j:=1 To kmin1 Do
 
1487
                h := h-al^[j+k1n]*al^[j+i1n];
 
1488
              al^[ik] := h/lkk
 
1489
            End; {i}
 
1490
          h := t^[k];
 
1491
          For j:=1 To kmin1 Do
 
1492
            h := h-al^[j+k1n]*t^[j];
 
1493
          t^[k] := h/lkk;
 
1494
          h := px^[k];
 
1495
          For j:=1 To kmin1 Do
 
1496
            h := h-al^[j+k1n]*px^[j];
 
1497
          px^[k] := h/lkk
 
1498
        End {lkk > 0}
 
1499
    End; {k}
 
1500
    If pd
 
1501
     Then
 
1502
      Begin
 
1503
        normt := 0;
 
1504
        For i:=n Downto 1 Do
 
1505
          Begin
 
1506
            ii := i+(i-1)*n;
 
1507
           h := t^[i];
 
1508
            For j:=i+1 To n Do
 
1509
              h := h-al^[i+(j-1)*n]*t^[j];
 
1510
            t^[i] := h/al^[ii];
 
1511
            h := px^[i];
 
1512
            For j:=i+1 To n Do
 
1513
              h := h-al^[i+(j-1)*n]*px^[j];
 
1514
            px^[i] := h/al^[ii];
 
1515
            h := abs(t^[i]);
 
1516
            If normt<h
 
1517
              Then
 
1518
                normt := h
 
1519
          End; {i}
 
1520
        ca := norma*normt/normr
 
1521
      End; {pd}
 
1522
  If pd
 
1523
   Then
 
1524
    term := 1
 
1525
  Else
 
1526
    term := 2;
 
1527
  freemem(al, sqr(n)*sr);
 
1528
 freemem(t, n*sr);
 
1529
End; {slegpd}
 
1530
 
 
1531
Procedure slegpdl(n: ArbInt; Var a1; Var b1, x1, ca: ArbFloat;
 
1532
                  Var term: ArbInt);
 
1533
 
 
1534
Var                   sr, i, j, k, kmin1 : ArbInt;
 
1535
                                      pd : boolean;
 
1536
    h, lkk, normr, normt, sumrowi, norma : ArbFloat;
 
1537
                                       a : ar2dr1 absolute a1;
 
1538
                                       b : arfloat1 absolute b1;
 
1539
                                       x : arfloat1 absolute x1;
 
1540
                                      al : par2dr1;
 
1541
                                       t : ^arfloat1;
 
1542
 
 
1543
Begin
 
1544
  If n<1 Then
 
1545
    Begin
 
1546
      term := 3;
 
1547
     exit
 
1548
    End;
 
1549
  sr := sizeof(ArbFloat);
 
1550
  AllocateL2dr(n, al);
 
1551
 getmem(t, n*sr);
 
1552
  For i:=1 To n Do
 
1553
   move(a[i]^, al^[i]^, i*sr);
 
1554
  move(b[1], x[1], n*sr);
 
1555
  normr := 0;
 
1556
 pd := true ;
 
1557
 norma := 0;
 
1558
  For i:=1 To n Do
 
1559
    Begin
 
1560
      sumrowi := 0;
 
1561
      For j:=1 To i Do
 
1562
       sumrowi := sumrowi+abs(al^[i]^[j]);
 
1563
      For j:=i+1 To n Do
 
1564
       sumrowi := sumrowi+abs(al^[j]^[i]);
 
1565
      If norma<sumrowi Then norma := sumrowi;
 
1566
      t^[i] := 2*random-1;
 
1567
     h := abs(t^[i]);
 
1568
      If normr<h Then normr := h
 
1569
    End; {i}
 
1570
  k := 0;
 
1571
  while (k<n) and pd Do
 
1572
    Begin
 
1573
      kmin1 := k;
 
1574
     k := k+1;
 
1575
     lkk := al^[k]^[k];
 
1576
      For j:=1 To kmin1 Do
 
1577
       lkk := lkk-sqr(al^[k]^[j]);
 
1578
      If lkk<=0 Then pd := false
 
1579
     Else
 
1580
        Begin
 
1581
          al^[k]^[k] := sqrt(lkk);
 
1582
         lkk := al^[k]^[k];
 
1583
          For i:=k+1 To n Do
 
1584
            Begin
 
1585
              h := al^[i]^[k];
 
1586
              For j:=1 To kmin1 Do
 
1587
               h := h-al^[k]^[j]*al^[i]^[j];
 
1588
              al^[i]^[k] := h/lkk
 
1589
            End; {i}
 
1590
          h := t^[k];
 
1591
          For j:=1 To kmin1 Do
 
1592
           h := h-al^[k]^[j]*t^[j];
 
1593
          t^[k] := h/lkk;
 
1594
         h := x[k];
 
1595
          For j:=1 To kmin1 Do
 
1596
           h := h-al^[k]^[j]*x[j];
 
1597
          x[k] := h/lkk
 
1598
        End {lkk > 0}
 
1599
    End; {k}
 
1600
    If pd Then
 
1601
      Begin
 
1602
        normt := 0;
 
1603
        For i:=n Downto 1 Do
 
1604
          Begin
 
1605
            h := t^[i];
 
1606
            For j:=i+1 To n Do
 
1607
             h := h-al^[j]^[i]*t^[j];
 
1608
            t^[i] := h/al^[i]^[i];
 
1609
            h := x[i];
 
1610
            For j:=i+1 To n Do
 
1611
             h := h-al^[j]^[i]*x[j];
 
1612
            x[i] := h/al^[i]^[i];
 
1613
           h := abs(t^[i]);
 
1614
            If normt<h Then normt := h
 
1615
          End; {i}
 
1616
        ca := norma*normt/normr
 
1617
      End; {pd}
 
1618
  If pd Then term := 1
 
1619
 Else term := 2;
 
1620
  DeAllocateL2dr(n, al);
 
1621
 freemem(t, n*sr);
 
1622
End; {slegpdl}
 
1623
 
 
1624
Procedure slegsy(n, rwidth: ArbInt; Var a, b, x, ca: ArbFloat;
 
1625
                 Var term:ArbInt);
 
1626
 
 
1627
Var
 
1628
   i, j, kmin1, k, kplus1, kmin2, nsr, nsi, nsb,
 
1629
   imin1, jmin1, indexpivot, iplus1, indi, indj, indk, indp       : ArbInt;
 
1630
   h, absh, maxim, pivot, ct, norma, sumrowi, normt, normr, s : ArbFloat;
 
1631
              pa, pb, pb1, px, alt, l, d, t, u, v, l1, d1, u1, t1 : ^arfloat1;
 
1632
                                                                p : ^arint1;
 
1633
                                                                q : ^arbool1;
 
1634
Begin
 
1635
  If (n<1) Or (rwidth<1)
 
1636
   Then
 
1637
    Begin
 
1638
      term := 3;
 
1639
     exit
 
1640
    End; {if}
 
1641
  pa := @a;
 
1642
 pb := @b;
 
1643
 px := @x;
 
1644
  nsr := n*sizeof(ArbFloat);
 
1645
  nsi := n*sizeof(ArbInt);
 
1646
  nsb := n*sizeof(boolean);
 
1647
  getmem(alt, n*nsr);
 
1648
  getmem(l, nsr);
 
1649
  getmem(d, nsr);
 
1650
  getmem(t, nsr);
 
1651
  getmem(u, nsr);
 
1652
  getmem(v, nsr);
 
1653
  getmem(p, nsi);
 
1654
  getmem(q, nsb);
 
1655
  getmem(l1, nsr);
 
1656
  getmem(d1, nsr);
 
1657
  getmem(u1, nsr);
 
1658
  getmem(t1, nsr);
 
1659
  getmem(pb1, nsr);
 
1660
  move(pb^, pb1^, nsr);
 
1661
  For i:=1 To n Do
 
1662
    Begin
 
1663
      indi := (i-1)*n;
 
1664
      For j:=1 To i Do
 
1665
        alt^[indi+j] := pa^[(i-1)*rwidth+j];
 
1666
    End; {i}
 
1667
  norma := 0;
 
1668
  For i:=1 To n Do
 
1669
    Begin
 
1670
      indi := (i-1)*n;
 
1671
      p^[i] := i;
 
1672
     sumrowi := 0;
 
1673
      For j:=1 To i Do
 
1674
        sumrowi := sumrowi+abs(alt^[indi+j]);
 
1675
      For j:=i+1 To n Do
 
1676
        sumrowi := sumrowi+abs(alt^[(j-1)*n+i]);
 
1677
      If norma<sumrowi
 
1678
       Then
 
1679
        norma := sumrowi
 
1680
    End; {i}
 
1681
  kmin1 := -1;
 
1682
 k := 0;
 
1683
 kplus1 := 1;
 
1684
  while k<n Do
 
1685
    Begin
 
1686
      kmin2 := kmin1;
 
1687
     kmin1 := k;
 
1688
     k := kplus1;
 
1689
     kplus1 := kplus1+1;
 
1690
      indk := kmin1*n;
 
1691
      If k>3
 
1692
       Then
 
1693
        Begin
 
1694
          t^[2] := alt^[n+2]*alt^[indk+1]+alt^[2*n+2]*alt^[indk+2];
 
1695
          For i:=3 To kmin2 Do
 
1696
            Begin
 
1697
              indi := (i-1)*n;
 
1698
              t^[i] := alt^[indi+i-1]*alt^[indk+i-2]+alt^[indi+i]
 
1699
                       *alt^[indk+i-1]+alt^[indi+n+i]*alt^[indk+i]
 
1700
            End; {i}
 
1701
          t^[kmin1] := alt^[indk-n+kmin2]*alt^[indk+k-3]
 
1702
                       +alt^[indk-n+kmin1]*alt^[indk+kmin2]
 
1703
                       +alt^[indk+kmin1];
 
1704
          h := alt^[indk+k];
 
1705
          For j:=2 To kmin1 Do
 
1706
            h := h-t^[j]*alt^[indk+j-1];
 
1707
          t^[k] := h;
 
1708
          alt^[indk+k] := h-alt^[indk+kmin1]*alt^[indk+kmin2]
 
1709
        End {k>3}
 
1710
      Else
 
1711
       If k=3
 
1712
        Then
 
1713
        Begin
 
1714
          t^[2] := alt^[n+2]*alt^[2*n+1]+alt^[2*n+2];
 
1715
          h := alt^[2*n+3]-t^[2]*alt^[2*n+1];
 
1716
          t^[3] := h;
 
1717
          alt^[2*n+3] := h-alt^[2*n+2]*alt^[2*n+1]
 
1718
        End  {k=3}
 
1719
      Else
 
1720
       If k=2
 
1721
        Then
 
1722
        t^[2] := alt^[n+2];
 
1723
      maxim := 0;
 
1724
      For i:=kplus1 To n Do
 
1725
        Begin
 
1726
          indi := (i-1)*n;
 
1727
          h := alt^[indi+k];
 
1728
          For j:=2 To k Do
 
1729
            h := h-t^[j]*alt^[indi+j-1];
 
1730
          absh := abs(h);
 
1731
          If maxim<absh
 
1732
           Then
 
1733
            Begin
 
1734
              maxim := absh;
 
1735
             indexpivot := i
 
1736
            End; {if}
 
1737
          alt^[indi+k] := h
 
1738
        End; {i}
 
1739
      If maxim <> 0
 
1740
       Then
 
1741
        Begin
 
1742
          If indexpivot>kplus1
 
1743
           Then
 
1744
            Begin
 
1745
              indp := (indexpivot-1)*n;
 
1746
              indk := k*n;
 
1747
              p^[kplus1] := indexpivot;
 
1748
              For j:=1 To k Do
 
1749
                Begin
 
1750
                  h := alt^[indk+j];
 
1751
                  alt^[indk+j] := alt^[indp+j];
 
1752
                  alt^[indp+j] := h
 
1753
                End; {j}
 
1754
              For j:=indexpivot Downto kplus1 Do
 
1755
                Begin
 
1756
                  indj := (j-1)*n;
 
1757
                  h := alt^[indj+kplus1];
 
1758
                  alt^[indj+kplus1] := alt^[indp+j];
 
1759
                  alt^[indp+j] := h
 
1760
                End; {j}
 
1761
              For i:=indexpivot To n Do
 
1762
                Begin
 
1763
                  indi := (i-1)*n;
 
1764
                  h := alt^[indi+kplus1];
 
1765
                  alt^[indi+kplus1] := alt^[indi+indexpivot];
 
1766
                  alt^[indi+indexpivot] := h
 
1767
                End  {i}
 
1768
            End; {if}
 
1769
          pivot := alt^[k*n+k];
 
1770
          For i:=k+2 To n Do
 
1771
            alt^[(i-1)*n+k] := alt^[(i-1)*n+k]/pivot
 
1772
        End {maxim <> 0}
 
1773
    End; {k}
 
1774
  d^[1] := alt^[1];
 
1775
 i := 1;
 
1776
  while i<n Do
 
1777
    Begin
 
1778
      imin1 := i;
 
1779
     i := i+1;
 
1780
      u^[imin1] := alt^[(i-1)*n+imin1];
 
1781
      l^[imin1] := u^[imin1];
 
1782
     d^[i] := alt^[(i-1)*n+i]
 
1783
    End; {i}
 
1784
  mdtgtr(n, l^[1], d^[1], u^[1], l1^[1], d1^[1], u1^[1], v^[1],
 
1785
         q^[1], ct, term);
 
1786
  If term=1
 
1787
   Then
 
1788
    Begin
 
1789
      normr := 0;
 
1790
      For i:=1 To n Do
 
1791
        Begin
 
1792
          t^[i] := 2*random-1;
 
1793
         h := t^[i];
 
1794
          h := abs(h);
 
1795
          If normr<h
 
1796
           Then
 
1797
            normr := h
 
1798
        End; {i}
 
1799
      For i:=1 To n Do
 
1800
        Begin
 
1801
          indexpivot := p^[i];
 
1802
          If indexpivot <> i
 
1803
           Then
 
1804
            Begin
 
1805
              h := pb1^[i];
 
1806
             pb1^[i] := pb1^[indexpivot];
 
1807
              pb1^[indexpivot] := h
 
1808
            End {if}
 
1809
        End; {i}
 
1810
      i := 0;
 
1811
      while i<n Do
 
1812
        Begin
 
1813
          indi := i*n;
 
1814
          imin1 := i;
 
1815
         i := i+1;
 
1816
         j := 1;
 
1817
         h := t^[i];
 
1818
         s := pb1^[i];
 
1819
          while j<imin1 Do
 
1820
            Begin
 
1821
              jmin1 := j;
 
1822
             j := j+1;
 
1823
              s := s-alt^[indi+jmin1]*pb1^[j];
 
1824
              h := h-alt^[indi+jmin1]*t^[j]
 
1825
            End; {j}
 
1826
          t^[i] := h;
 
1827
         pb1^[i] := s
 
1828
        End; {i}
 
1829
      dslgtr(n, l1^[1], d1^[1], u1^[1], v^[1], q^[1], pb1^[1], px^[1], term);
 
1830
      dslgtr(n, l1^[1], d1^[1], u1^[1], v^[1], q^[1], t^[1], t1^[1], term);
 
1831
      i := n+1;
 
1832
     imin1 := n;
 
1833
     normt := 0;
 
1834
      while i>2 Do
 
1835
        Begin
 
1836
          iplus1 := i;
 
1837
         i := imin1;
 
1838
         imin1 := imin1-1;
 
1839
          h := t1^[i];
 
1840
         s := px^[i];
 
1841
          For j:=iplus1 To n Do
 
1842
            Begin
 
1843
              indj := (j-1)*n+imin1;
 
1844
              h := h-alt^[indj]*t1^[j];
 
1845
              s := s-alt^[indj]*px^[j]
 
1846
            End; {j}
 
1847
          px^[i] := s;
 
1848
          t1^[i] := h;
 
1849
         h := abs(h);
 
1850
          If normt<h
 
1851
           Then
 
1852
            normt := h
 
1853
        End; {i}
 
1854
      For i:=n Downto 1 Do
 
1855
        Begin
 
1856
          indexpivot := p^[i];
 
1857
          If indexpivot <> i
 
1858
           Then
 
1859
            Begin
 
1860
              h := px^[i];
 
1861
             px^[i] := px^[indexpivot];
 
1862
              px^[indexpivot] := h
 
1863
            End {if}
 
1864
        End; {i}
 
1865
      ca := norma*normt/normr
 
1866
    End {term=1}
 
1867
  Else
 
1868
    term := 2;
 
1869
  freemem(alt, n*nsr);
 
1870
  freemem(l, nsr);
 
1871
  freemem(d, nsr);
 
1872
  freemem(t, nsr);
 
1873
  freemem(u, nsr);
 
1874
  freemem(v, nsr);
 
1875
  freemem(p, nsi);
 
1876
  freemem(q, nsb);
 
1877
  freemem(l1, nsr);
 
1878
  freemem(d1, nsr);
 
1879
  freemem(u1, nsr);
 
1880
  freemem(t1, nsr);
 
1881
  freemem(pb1, nsr);
 
1882
End; {slegsy}
 
1883
 
 
1884
Procedure slegsyl(n: ArbInt; Var a1; Var b1, x1, ca: ArbFloat;
 
1885
                  Var term: ArbInt);
 
1886
 
 
1887
Var
 
1888
   i, j, k, nsr, nsi, nsb, indexpivot: ArbInt;
 
1889
   h, absh, maxim, pivot, ct, norma, sumrowi, normt, normr, s : ArbFloat;
 
1890
                                           a : ar2dr1 absolute a1;
 
1891
                                           b : arfloat1 absolute b1;
 
1892
                                           x : arfloat1 absolute x1;
 
1893
           b0, l, d, t, u, v, l1, d1, u1, t1 : ^arfloat1;
 
1894
                                         alt : par2dr1;
 
1895
                                           p : ^arint1;
 
1896
                                           q : ^arbool1;
 
1897
Begin
 
1898
  If n<1 Then
 
1899
    Begin
 
1900
      term := 3;
 
1901
     exit
 
1902
    End; {if}
 
1903
  nsr := n*sizeof(ArbFloat);
 
1904
 nsi := n*sizeof(ArbInt);
 
1905
 nsb := n*sizeof(boolean);
 
1906
  AllocateL2dr(n, alt);
 
1907
  getmem(l, nsr);
 
1908
 getmem(d, nsr);
 
1909
 getmem(t, nsr);
 
1910
  getmem(u, nsr);
 
1911
 getmem(v, nsr);
 
1912
 getmem(p, nsi);
 
1913
  getmem(q, nsb);
 
1914
 getmem(l1, nsr);
 
1915
 getmem(d1, nsr);
 
1916
  getmem(u1, nsr);
 
1917
 getmem(t1, nsr);
 
1918
 getmem(b0, nsr);
 
1919
  move(b[1], b0^, nsr);
 
1920
  For i:=1 To n Do
 
1921
   move(a[i]^, alt^[i]^, i*sizeof(ArbFloat));
 
1922
  norma := 0;
 
1923
  For i:=1 To n Do
 
1924
    Begin
 
1925
      p^[i] := i;
 
1926
     sumrowi := 0;
 
1927
      For j:=1 To i Do
 
1928
       sumrowi := sumrowi+abs(alt^[i]^[j]);
 
1929
      For j:=i+1 To n Do
 
1930
       sumrowi := sumrowi+abs(alt^[j]^[i]);
 
1931
      If norma<sumrowi Then norma := sumrowi
 
1932
    End; {i}
 
1933
  k := 0;
 
1934
  while k<n Do
 
1935
    Begin
 
1936
      Inc(k);
 
1937
      If k>3 Then
 
1938
        Begin
 
1939
          t^[2] := alt^[2]^[2]*alt^[k]^[1]+alt^[3]^[2]*alt^[k]^[2];
 
1940
          For i:=3 To k-2 Do
 
1941
            t^[i] := alt^[i]^[i-1]*alt^[k]^[i-2]+alt^[i]^[i]
 
1942
                     *alt^[k]^[i-1]+alt^[i+1]^[i]*alt^[k]^[i];
 
1943
          t^[k-1] := alt^[k-1]^[k-2]*alt^[k]^[k-3]
 
1944
                     +alt^[k-1]^[k-1]*alt^[k]^[k-2]+alt^[k]^[k-1];
 
1945
          h := alt^[k]^[k];
 
1946
          For j:=2 To k-1 Do
 
1947
           h := h-t^[j]*alt^[k]^[j-1];
 
1948
          t^[k] := h;
 
1949
          alt^[k]^[k] := h-alt^[k]^[k-1]*alt^[k]^[k-2]
 
1950
        End {k>3}
 
1951
      Else
 
1952
       If k=3
 
1953
        Then
 
1954
        Begin
 
1955
          t^[2] := alt^[2]^[2]*alt^[3]^[1]+alt^[3]^[2];
 
1956
          h := alt^[3]^[3]-t^[2]*alt^[3]^[1];
 
1957
          t^[3] := h;
 
1958
          alt^[3]^[3] := h-alt^[3]^[2]*alt^[3]^[1]
 
1959
        End  {k=3}
 
1960
      Else
 
1961
       If k=2 Then t^[2] := alt^[2]^[2];
 
1962
      maxim := 0;
 
1963
      For i:=k+1 To n Do
 
1964
        Begin
 
1965
          h := alt^[i]^[k];
 
1966
          For j:=2 To k Do
 
1967
           h := h-t^[j]*alt^[i]^[j-1];
 
1968
          absh := abs(h);
 
1969
          If maxim<absh Then
 
1970
            Begin
 
1971
              maxim := absh;
 
1972
             indexpivot := i
 
1973
            End; {if}
 
1974
          alt^[i]^[k] := h
 
1975
        End; {i}
 
1976
      If maxim <> 0
 
1977
       Then
 
1978
        Begin
 
1979
          If indexpivot>k+1 Then
 
1980
            Begin
 
1981
              p^[k+1] := indexpivot;
 
1982
              For j:=1 To k Do
 
1983
                Begin
 
1984
                  h := alt^[k+1]^[j];
 
1985
                  alt^[k+1]^[j] := alt^[indexpivot]^[j];
 
1986
                  alt^[indexpivot]^[j] := h
 
1987
                End; {j}
 
1988
              For j:=indexpivot Downto k+1 Do
 
1989
                Begin
 
1990
                  h := alt^[j]^[k+1];
 
1991
                  alt^[j]^[k+1] := alt^[indexpivot]^[j];
 
1992
                  alt^[indexpivot]^[j] := h
 
1993
                End; {j}
 
1994
              For i:=indexpivot To n Do
 
1995
                Begin
 
1996
                  h := alt^[i]^[k+1];
 
1997
                  alt^[i]^[k+1] := alt^[i]^[indexpivot];
 
1998
                  alt^[i]^[indexpivot] := h
 
1999
                End  {i}
 
2000
            End; {if}
 
2001
          pivot := alt^[k+1]^[k];
 
2002
          For i:=k+2 To n Do
 
2003
           alt^[i]^[k] := alt^[i]^[k]/pivot
 
2004
        End {maxim <> 0}
 
2005
    End; {k}
 
2006
  d^[1] := alt^[1]^[1];
 
2007
 i := 1;
 
2008
  while i<n Do
 
2009
    Begin
 
2010
      Inc(i);
 
2011
      u^[i-1] := alt^[i]^[i-1];
 
2012
      l^[i-1] := u^[i-1];
 
2013
     d^[i] := alt^[i]^[i]
 
2014
    End; {i}
 
2015
  mdtgtr(n, l^[1], d^[1], u^[1], l1^[1], d1^[1], u1^[1], v^[1],
 
2016
         q^[1], ct, term);
 
2017
  If term=1 Then
 
2018
    Begin
 
2019
      normr := 0;
 
2020
      For i:=1 To n Do
 
2021
        Begin
 
2022
          t^[i] := 2*random-1;
 
2023
         h := t^[i];
 
2024
          h := abs(h);
 
2025
          If normr<h Then normr := h
 
2026
        End; {i}
 
2027
      For i:=1 To n Do
 
2028
        Begin
 
2029
          indexpivot := p^[i];
 
2030
          If indexpivot <> i
 
2031
           Then
 
2032
            Begin
 
2033
              h := b0^[i];
 
2034
             b0^[i] := b0^[indexpivot];
 
2035
              b0^[indexpivot] := h
 
2036
            End {if}
 
2037
        End; {i}
 
2038
      i := 0;
 
2039
      while i<n Do
 
2040
        Begin
 
2041
          Inc(i);
 
2042
         j := 1;
 
2043
         h := t^[i];
 
2044
         s := b0^[i];
 
2045
          while j<i-1 Do
 
2046
            Begin
 
2047
              Inc(j);
 
2048
              s := s-alt^[i]^[j-1]*b0^[j];
 
2049
              h := h-alt^[i]^[j-1]*t^[j]
 
2050
            End; {j}
 
2051
          t^[i] := h;
 
2052
         b0^[i] := s
 
2053
        End; {i}
 
2054
      dslgtr(n, l1^[1], d1^[1], u1^[1], v^[1], q^[1], b0^[1], x[1], term);
 
2055
      dslgtr(n, l1^[1], d1^[1], u1^[1], v^[1], q^[1], t^[1], t1^[1], term);
 
2056
      i := n+1;
 
2057
     normt := 0;
 
2058
      while i>2 Do
 
2059
        Begin
 
2060
          Dec(i);
 
2061
          h := t1^[i];
 
2062
         s := x[i];
 
2063
          For j:=i+1 To n Do
 
2064
            Begin
 
2065
              h := h-alt^[j]^[i-1]*t1^[j];
 
2066
              s := s-alt^[j]^[i-1]*x[j]
 
2067
            End; {j}
 
2068
          x[i] := s;
 
2069
         t1^[i] := h;
 
2070
         h := abs(h);
 
2071
          If normt<h Then normt := h
 
2072
        End; {i}
 
2073
      For i:=n Downto 1 Do
 
2074
        Begin
 
2075
          indexpivot := p^[i];
 
2076
          If indexpivot <> i Then
 
2077
            Begin
 
2078
              h := x[i];
 
2079
             x[i] := x[indexpivot];
 
2080
             x[indexpivot] := h
 
2081
            End {if}
 
2082
        End; {i}
 
2083
      ca := norma*normt/normr
 
2084
    End {term=1}
 
2085
  Else
 
2086
    term := 2;
 
2087
  freemem(l, nsr);
 
2088
 freemem(d, nsr);
 
2089
 freemem(t, nsr);
 
2090
  freemem(u, nsr);
 
2091
 freemem(v, nsr);
 
2092
 freemem(p, nsi);
 
2093
  freemem(q, nsb);
 
2094
 freemem(l1, nsr);
 
2095
 freemem(d1, nsr);
 
2096
  freemem(u1, nsr);
 
2097
 freemem(t1, nsr);
 
2098
 freemem(b0, nsr);
 
2099
  DeAllocateL2dr(n, alt);
 
2100
End; {slegsyl}
 
2101
 
 
2102
Procedure slegtr(n:ArbInt; Var l, d, u, b, x, ca: ArbFloat;
 
2103
                 Var term: ArbInt);
 
2104
 
 
2105
Var                           singular, ch : boolean;
 
2106
               i, j, nm1, sr, n1s, ns, n2s : ArbInt;
 
2107
            normr, normt, h, lj, di, ui, m : ArbFloat;
 
2108
                                    pl, ll : ^arfloat2;
 
2109
    pd, pu, pb, px, dd, uu1, u2, t, sumrow : ^arfloat1;
 
2110
Begin
 
2111
  If n<1
 
2112
   Then
 
2113
    Begin
 
2114
      term := 3;
 
2115
     exit
 
2116
    End; {n<1}
 
2117
  sr := sizeof(ArbFloat);
 
2118
 n1s := (n-1)*sr;
 
2119
 ns := n1s+sr;
 
2120
 n2s := n1s;
 
2121
  getmem(ll, n1s);
 
2122
  getmem(uu1, n1s);
 
2123
  getmem(u2, n2s);
 
2124
  getmem(dd, ns);
 
2125
  getmem(t, ns);
 
2126
  getmem(sumrow, ns);
 
2127
 
 
2128
  pl := @l;
 
2129
 pd := @d;
 
2130
 pu := @u;
 
2131
 pb := @b;
 
2132
 px := @x;
 
2133
  move(pl^[2], ll^[2], n1s);
 
2134
  move(pd^[1], dd^[1], ns);
 
2135
  If n>1
 
2136
   Then
 
2137
    move(pu^[1], uu1^[1], n1s);
 
2138
  move(pb^[1], px^[1], ns);
 
2139
  normr := 0;
 
2140
 singular := false;
 
2141
  nm1 := n-1;
 
2142
 i := 0;
 
2143
  while (i<n) and not singular Do
 
2144
    Begin
 
2145
      i := i+1;
 
2146
      If i=1
 
2147
       Then
 
2148
        Begin
 
2149
          sumrow^[i] := abs(dd^[1]);
 
2150
          If n>1
 
2151
           Then
 
2152
            sumrow^[i] := sumrow^[i]+abs(uu1^[1])
 
2153
        End {i=1}
 
2154
      Else
 
2155
        If i=n
 
2156
         Then
 
2157
          sumrow^[i] := abs(ll^[n])+abs(dd^[n])
 
2158
        Else
 
2159
          sumrow^[i] := abs(ll^[i])+abs(dd^[i])+abs(uu1^[i]);
 
2160
      If sumrow^[i]=0
 
2161
       Then
 
2162
        singular := true
 
2163
      Else
 
2164
        Begin
 
2165
          h := 2*random-1;
 
2166
         t^[i] := sumrow^[i]*h;
 
2167
          h := abs(h);
 
2168
          If normr<h
 
2169
           Then
 
2170
            normr := h
 
2171
        End {sumrow^[i] <> 0}
 
2172
    End; {i}
 
2173
  j := 1;
 
2174
  while (j <> n) and  not singular Do
 
2175
    Begin
 
2176
      i := j;
 
2177
     j := j+1;
 
2178
     lj := ll^[j];
 
2179
      If lj <> 0
 
2180
       Then
 
2181
        Begin
 
2182
          di := dd^[i];
 
2183
          ch := abs(di/sumrow^[i])<abs(lj/sumrow^[j]);
 
2184
          If ch
 
2185
           Then
 
2186
            Begin
 
2187
              ui := uu1^[i];
 
2188
             dd^[i] := lj;
 
2189
             uu1^[i] := dd^[j];
 
2190
             m := di/lj;
 
2191
              dd^[j] := ui-m*dd^[j];
 
2192
              If i<nm1
 
2193
               Then
 
2194
                Begin
 
2195
                  u2^[i] := uu1^[j];
 
2196
                 uu1^[j] := -m*u2^[i]
 
2197
                End; {i<nm1}
 
2198
              sumrow^[j] := sumrow^[i];
 
2199
              h := t^[i];
 
2200
             t^[i] := t^[j];
 
2201
             t^[j] := h-m*t^[i];
 
2202
             h := px^[i];
 
2203
              px^[i] := px^[j];
 
2204
             px^[j] := h-m*px^[i]
 
2205
            End {ch}
 
2206
          Else
 
2207
            Begin
 
2208
              m := lj/di;
 
2209
             dd^[j] := dd^[j]-m*uu1^[i];
 
2210
              If i<nm1
 
2211
               Then
 
2212
                u2^[i] := 0;
 
2213
              t^[j] := t^[j]-m*t^[i];
 
2214
             px^[j] := px^[j]-m*px^[i]
 
2215
            End {not ch}
 
2216
        End {lj <> 0}
 
2217
      Else
 
2218
        Begin
 
2219
          If i < nm1
 
2220
            Then
 
2221
              u2^[i] := 0;
 
2222
          If dd^[i]=0
 
2223
           Then
 
2224
            singular := true
 
2225
        End {lj=0}
 
2226
    End; {j}
 
2227
  If dd^[n]=0
 
2228
   Then
 
2229
    singular := true;
 
2230
  If  Not singular
 
2231
    Then
 
2232
      Begin
 
2233
        normt := 0;
 
2234
       t^[n] := t^[n]/dd^[n];
 
2235
       px^[n] := px^[n]/dd^[n];
 
2236
       h := abs(t^[n]);
 
2237
        If normt<h
 
2238
         Then
 
2239
          normt := h;
 
2240
        If n>1
 
2241
         Then
 
2242
          Begin
 
2243
            t^[nm1] := (t^[nm1]-uu1^[nm1]*t^[n])/dd^[nm1];
 
2244
            px^[nm1] := (px^[nm1]-uu1^[nm1]*px^[n])/dd^[nm1];
 
2245
           h := abs(t^[nm1])
 
2246
          End; {n>1}
 
2247
        If normt<h
 
2248
         Then
 
2249
          normt := h;
 
2250
        For i:=n-2 Downto 1 Do
 
2251
          Begin
 
2252
            t^[i] := (t^[i]-uu1^[i]*t^[i+1]-u2^[i]*t^[i+2])/dd^[i];
 
2253
            px^[i] := (px^[i]-uu1^[i]*px^[i+1]-u2^[i]*px^[i+2])/dd^[i];
 
2254
            h := abs(t^[i]);
 
2255
            If normt<h
 
2256
             Then
 
2257
              normt := h
 
2258
          End; {i}
 
2259
        ca := normt/normr
 
2260
      End; {not singular}
 
2261
  If singular
 
2262
   Then
 
2263
    term := 2
 
2264
  Else
 
2265
    term := 1;
 
2266
  freemem(ll, n1s);
 
2267
  freemem(uu1, n1s);
 
2268
  freemem(u2, n2s);
 
2269
  freemem(dd, ns);
 
2270
  freemem(t, ns);
 
2271
  freemem(sumrow, ns);
 
2272
End; {slegtr}
 
2273
 
 
2274
End.