← Back to branch summary

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

« back to all changes in this revision

Viewing changes to fpcsrc/packages/extra/numlib/eig.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/eig.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
 
 
11
    See the file COPYING.FPC, included in this distribution,
 
12
    for details about the copyright.
 
13
 
 
14
    This program is distributed in the hope that it will be useful,
 
15
    but WITHOUT ANY WARRANTY; without even the implied warranty of
 
16
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
17
 
 
18
 **********************************************************************}
 
19
 
 
20
unit eig;
 
21
{$I DIRECT.INC}
 
22
 
 
23
interface
 
24
 
 
25
uses typ;
 
26
 
 
27
const versie = 'augustus 1993';
 
28
 
 
29
procedure eiggs1(var a: ArbFloat; n, rwidth: ArbInt; var lam: ArbFloat;
 
30
                 var term: ArbInt);
 
31
 
 
32
procedure eiggs2(var a: ArbFloat; n, rwidth, k1, k2: ArbInt;
 
33
                 var lam: ArbFloat; var term: ArbInt);
 
34
 
 
35
procedure eiggs3(var a: ArbFloat; n, rwidtha: ArbInt; var lam, x: ArbFloat;
 
36
                 rwidthx: ArbInt; var term: ArbInt);
 
37
 
 
38
procedure eiggs4(var a: ArbFloat; n, rwidtha, k1, k2: ArbInt; var lam, x: ArbFloat;
 
39
                 rwidthx: ArbInt; var m2, term: ArbInt);
 
40
 
 
41
procedure eigts1(var d, cd: ArbFloat; n: ArbInt; var lam: ArbFloat;
 
42
                 var term: ArbInt);
 
43
 
 
44
procedure eigts2(var d, cd: ArbFloat; n, k1, k2: ArbInt; var lam: ArbFloat;
 
45
                 var term: ArbInt);
 
46
 
 
47
procedure eigts3(var d, cd: ArbFloat; n: ArbInt; var lam, x: ArbFloat;
 
48
                 rwidth: ArbInt; var term: ArbInt);
 
49
 
 
50
procedure eigts4(var d, cd: ArbFloat; n, k1, k2: ArbInt; var lam, x: ArbFloat;
 
51
                 rwidth: ArbInt; var m2, term: ArbInt);
 
52
 
 
53
procedure eigbs1(var a: ArbFloat; n, l: ArbInt; var lam: ArbFloat;
 
54
                 var term: ArbInt);
 
55
 
 
56
procedure eigbs2(var a: ArbFloat; n, l, k1, k2: ArbInt; var lam: ArbFloat;
 
57
                 var term: ArbInt);
 
58
 
 
59
procedure eigbs3(var a: ArbFloat; n, l: ArbInt; var lam, x: ArbFloat;
 
60
                 rwidthx: ArbInt; var term: ArbInt);
 
61
 
 
62
procedure eigbs4(var a: ArbFloat; n, l, k1, k2: ArbInt;
 
63
                 var lam, x: ArbFloat;  rwidthx: ArbInt;
 
64
                 var m2, term: ArbInt);
 
65
 
 
66
procedure eigge1(var a: ArbFloat; n, rwidth: ArbInt; var lam: complex;
 
67
                 var term: ArbInt);
 
68
 
 
69
procedure eigge3(var a: ArbFloat; n, rwidtha: ArbInt; var lam, x: complex;
 
70
                 rwidthx: ArbInt; var term: ArbInt);
 
71
 
 
72
procedure eiggg1(var a: ArbFloat; n, rwidtha: ArbInt; var b: ArbFloat;
 
73
                 rwidthb: ArbInt; var lam: ArbFloat; var term: ArbInt);
 
74
 
 
75
procedure eiggg2(var a: ArbFloat; n, rwidtha, k1, k2: ArbInt; var b: ArbFloat;
 
76
                 rwidthb: ArbInt; var lam: ArbFloat; var term: ArbInt);
 
77
 
 
78
procedure eiggg3(var a: ArbFloat; n, rwidtha: ArbInt; var b: ArbFloat;
 
79
                 rwidthb: ArbInt; var lam, x: ArbFloat; rwidthx: ArbInt;
 
80
                 var term: ArbInt);
 
81
 
 
82
procedure eiggg4(var a: ArbFloat; n, rwidtha, k1, k2: ArbInt; var b: ArbFloat;
 
83
                 rwidthb: ArbInt; var lam, x: ArbFloat; rwidthx: ArbInt;
 
84
                 var m2, term: ArbInt);
 
85
 
 
86
procedure eigsv1(var a: ArbFloat; m, n, rwidth: ArbInt; var sig: ArbFloat;
 
87
                 var term: ArbInt);
 
88
 
 
89
procedure eigsv3(var a: ArbFloat; m, n, rwidtha: ArbInt; var sig, u: ArbFloat;
 
90
                 rwidthu: ArbInt; var v: ArbFloat; rwidthv: ArbInt;
 
91
                 var term: ArbInt);
 
92
 
 
93
implementation
 
94
 
 
95
uses eigh1, eigh2;
 
96
 
 
97
procedure eiggs1(var a: ArbFloat; n, rwidth: ArbInt; var lam: ArbFloat;
 
98
                 var term: ArbInt);
 
99
var            i, sr, nsr : ArbInt;
 
100
    d, cd, dh, cdh, u, pa : ^arfloat1;
 
101
begin
 
102
  if n<1 then
 
103
    begin
 
104
      term:=3; exit
 
105
    end; {wrong input}
 
106
  pa:=@a;
 
107
  sr:=sizeof(ArbFloat); nsr:=n*sr;
 
108
  getmem(d, nsr); getmem(cd, nsr); getmem(dh, nsr); getmem(cdh, nsr);
 
109
  getmem(u, n*nsr);
 
110
  for i:=1 to n do move(pa^[(i-1)*rwidth+1], u^[(i-1)*n+1], i*sr);
 
111
  tred1(u^[1], n, n, d^[1], cd^[1], term);
 
112
  move(d^[1], dh^[1], nsr); move(cd^[1], cdh^[1], nsr);
 
113
  tql1(d^[1], cd^[1], n, lam, term);
 
114
  if term=2 then bisect(dh^[1], cdh^[1], n, 1, n, 0, lam, term);
 
115
  freemem(d, nsr); freemem(cd, nsr); freemem(dh, nsr); freemem(cdh, nsr);
 
116
  freemem(u, n*nsr);
 
117
end; {eiggs1}
 
118
 
 
119
procedure eiggs2(var a: ArbFloat; n, rwidth, k1, k2: ArbInt;
 
120
                 var lam: ArbFloat; var term: ArbInt);
 
121
var          i, sr, nsr : ArbInt;
 
122
           d, cd, u, pa : ^arfloat1;
 
123
begin
 
124
  if (n<1) or (k1<1) or (k2<k1) or (k2>n) then
 
125
    begin
 
126
      term:=3; exit
 
127
    end; {wrong input}
 
128
  pa:=@a;
 
129
  sr:=sizeof(ArbFloat); nsr:=n*sr;
 
130
  getmem(d, nsr); getmem(cd, nsr); getmem(u, n*nsr);
 
131
  for i:=1 to n do move(pa^[(i-1)*rwidth+1], u^[(i-1)*n+1], i*sr);
 
132
  tred1(u^[1], n, n, d^[1], cd^[1], term);
 
133
  bisect(d^[1], cd^[1], n, k1, k2, 0, lam, term);
 
134
  freemem(d, nsr); freemem(cd, nsr); freemem(u, n*nsr);
 
135
end; {eiggs2}
 
136
 
 
137
procedure eiggs3(var a: ArbFloat; n, rwidtha: ArbInt; var lam, x: ArbFloat;
 
138
                 rwidthx: ArbInt; var term: ArbInt);
 
139
var   nsr : ArbInt;
 
140
    d, cd : ^arfloat1;
 
141
begin
 
142
  if n<1 then
 
143
    begin
 
144
      term:=3; exit
 
145
    end;
 
146
  nsr:=n*sizeof(ArbFloat);
 
147
  getmem(d, nsr); getmem(cd, nsr);
 
148
  tred2(a, n, rwidtha, d^[1], cd^[1], x, rwidthx, term);
 
149
  tql2(d^[1], cd^[1], n, lam, x, rwidthx, term);
 
150
  freemem(d, nsr); freemem(cd, nsr)
 
151
end;  {eiggs3}
 
152
 
 
153
procedure eiggs4(var a: ArbFloat; n, rwidtha, k1, k2: ArbInt; var lam, x: ArbFloat;
 
154
                 rwidthx: ArbInt; var m2, term: ArbInt);
 
155
var      i, sr, nsr : ArbInt;
 
156
       pa, d, cd, u : ^arfloat1;
 
157
begin
 
158
  if (n<1) or (k1<1) or (k2<k1) or (k2>n) then
 
159
    begin
 
160
      term:=3; exit
 
161
    end; {wrong input}
 
162
  pa:=@a;
 
163
  sr:=sizeof(ArbFloat); nsr:=n*sr;
 
164
  getmem(d, nsr); getmem(cd, nsr); getmem(u, n*nsr);
 
165
  for i:=1 to n do move(pa^[(i-1)*rwidtha+1], u^[(i-1)*n+1], i*sr);
 
166
  tred1(u^[1], n, n, d^[1], cd^[1], term);
 
167
  trsturm1(d^[1], cd^[1], n, k1, k2, lam, x, rwidthx, m2, term);
 
168
  trbak1(u^[1], n, n, cd^[1], k1, k2, x, rwidthx);
 
169
  freemem(d, nsr); freemem(cd, nsr); freemem(u, n*nsr) { toegevoegd 3 apr 92 }
 
170
end; {eiggs4}
 
171
 
 
172
procedure eigts1(var d, cd: ArbFloat; n: ArbInt; var lam: ArbFloat;
 
173
                 var term: ArbInt);
 
174
var               sr, nsr : ArbInt;
 
175
         pd, pcd, dh, cdh : ^arfloat1;
 
176
begin
 
177
  if n<1 then
 
178
    begin
 
179
      term:=3; exit
 
180
    end; {wrong input}
 
181
  sr:=sizeof(ArbFloat); nsr:=n*sr;
 
182
  pd:=@d; pcd:=@cd; getmem(dh, nsr); getmem(cdh, nsr);
 
183
  move(pd^[1], dh^[1], nsr); move(pcd^[1], cdh^[2], (n-1)*sr);
 
184
  tql1(dh^[1], cdh^[1], n, lam, term);
 
185
  if term=2 then
 
186
    begin
 
187
      move(pd^[1], dh^[1], nsr); move(pcd^[1], cdh^[2], (n-1)*sr);
 
188
      bisect(dh^[1], cdh^[1], n, 1, n, 0, lam, term)
 
189
    end;
 
190
  freemem(dh, nsr); freemem(cdh, nsr);
 
191
end; {eigts1}
 
192
 
 
193
procedure eigts2(var d, cd: ArbFloat; n, k1, k2: ArbInt; var lam: ArbFloat;
 
194
                 var term: ArbInt);
 
195
var               sr, nsr : ArbInt;
 
196
                 pcd, cdh : ^arfloat1;
 
197
begin
 
198
  if (n<1) or (k1<1) or (k2<k1) or (k2>n) then
 
199
    begin
 
200
      term:=3; exit
 
201
    end; {wrong input}
 
202
  pcd:=@cd;
 
203
  term:=1; sr:=sizeof(ArbFloat); nsr:=n*sr; getmem(cdh, nsr);
 
204
  move(pcd^[1], cdh^[2], (n-1)*sr);
 
205
  bisect(d, cdh^[1], n, k1, k2, 0, lam, term);
 
206
  freemem(cdh, nsr)
 
207
end; {eigts2}
 
208
 
 
209
procedure eigts3(var d, cd: ArbFloat; n: ArbInt; var lam, x: ArbFloat;
 
210
                 rwidth: ArbInt; var term: ArbInt);
 
211
var             i, sr, nsr : ArbInt;
 
212
              px, pcd, cdh : ^arfloat1;
 
213
begin
 
214
  if n<1 then
 
215
    begin
 
216
      term:=3; exit
 
217
    end;
 
218
  px:=@x; pcd:=@cd;
 
219
  sr:=sizeof(ArbFloat); nsr:=n*sr;
 
220
  getmem(cdh, nsr);
 
221
  move(pcd^[1], cdh^[2], (n-1)*sr);
 
222
  for i:=1 to n do fillchar(px^[(i-1)*rwidth+1], nsr, 0);
 
223
  for i:=1 to n do px^[(i-1)*rwidth+i]:=1;
 
224
  tql2(d, cdh^[1], n, lam, px^[1], rwidth, term);
 
225
  freemem(cdh, nsr);
 
226
end;  {eigts3}
 
227
 
 
228
procedure eigts4(var d, cd: ArbFloat; n, k1, k2: ArbInt; var lam, x: ArbFloat;
 
229
                 rwidth: ArbInt; var m2, term: ArbInt);
 
230
var                    sr : ArbInt;
 
231
                 pcd, cdh : ^arfloat1;
 
232
begin
 
233
  if (n<1) or (k1<1) or (k2<k1) or (k2>n) then
 
234
    begin
 
235
      term:=3; exit
 
236
    end; {wrong input}
 
237
  term:=1;
 
238
  pcd:=@cd; sr:=sizeof(ArbFloat);
 
239
  getmem(cdh, n*sr);
 
240
  move(pcd^[1], cdh^[2], (n-1)*sr);
 
241
  trsturm1(d, cdh^[1], n, k1, k2, lam, x, rwidth, m2, term);
 
242
  freemem(cdh, n*sr)
 
243
end; {eigts4}
 
244
 
 
245
procedure eigbs1(var a: ArbFloat; n, l: ArbInt; var lam: ArbFloat;
 
246
                 var term: ArbInt);
 
247
var             u, d, cd : ^arfloat1;
 
248
      uwidth, i, sr, nsr : ArbInt;
 
249
begin
 
250
  if (n<1) or (l<0) or (l>n-1) then
 
251
    begin
 
252
      term:=3; exit
 
253
    end; {wrong input}
 
254
  sr:=sizeof(ArbFloat); nsr:=n*sr; uwidth:=l+1;
 
255
  getmem(u, uwidth*nsr); getmem(d, nsr); getmem(cd, nsr);
 
256
  transf(a, n, l, u^[1], uwidth);
 
257
  bandrd1(u^[1], n, l, uwidth, d^[1], cd^[1]);
 
258
  eigts1(d^[1], cd^[2], n, lam, term);
 
259
  freemem(u, uwidth*nsr); freemem(d, nsr); freemem(cd, nsr);
 
260
end; {eigbs1}
 
261
 
 
262
procedure eigbs2(var a: ArbFloat; n, l, k1, k2: ArbInt; var lam: ArbFloat;
 
263
                 var term: ArbInt);
 
264
var                  u, d, cd : ^arfloat1;
 
265
           i, sr, nsr, uwidth : ArbInt;
 
266
begin
 
267
  if (n<1) or (k1<1) or (k2<k1) or (k2>n) or (l<0) or (l>n-1) then
 
268
    begin
 
269
      term:=3; exit
 
270
    end; {wrong input}
 
271
  sr:=sizeof(ArbFloat); nsr:=n*sr; uwidth:=l+1;
 
272
  getmem(u, uwidth*nsr); getmem(d, nsr); getmem(cd, nsr);
 
273
  transf(a, n, l, u^[1], uwidth);
 
274
  bandrd1(u^[1], n, l, uwidth, d^[1], cd^[1]);
 
275
  eigts2(d^[1], cd^[2], n, k1, k2, lam, term);
 
276
  freemem(u, uwidth*nsr); freemem(d, nsr); freemem(cd, nsr)
 
277
end; {eigbs2}
 
278
 
 
279
procedure eigbs3(var a: ArbFloat; n, l: ArbInt; var lam, x: ArbFloat;
 
280
                 rwidthx: ArbInt; var term: ArbInt);
 
281
var                  u, d, cd : ^arfloat1;
 
282
           i, sr, nsr, uwidth : ArbInt;
 
283
begin
 
284
  if (n<1) or (l<0) or (l>n-1) then
 
285
    begin
 
286
      term:=3; exit
 
287
    end; {wrong input}
 
288
  sr:=sizeof(ArbFloat); nsr:=n*sr; uwidth:=l+1;
 
289
  getmem(u, uwidth*nsr); getmem(d, nsr); getmem(cd, nsr);
 
290
  transf(a, n, l, u^[1], uwidth);
 
291
  bandrd2(u^[1], n, l, uwidth, d^[1], cd^[1], x, rwidthx);
 
292
  tql2(d^[1], cd^[1], n, lam, x, rwidthx, term);
 
293
  freemem(u, uwidth*nsr); freemem(d, nsr); freemem(cd, nsr)
 
294
end; {eigbs3}
 
295
 
 
296
procedure eigbs4(var a: ArbFloat; n, l, k1, k2: ArbInt;
 
297
                 var lam, x: ArbFloat;  rwidthx: ArbInt;
 
298
                 var m2, term: ArbInt);
 
299
var  i, j, k, m, uwidth : ArbInt;
 
300
     plam, px, pa, v, u : ^arfloat1;
 
301
                s, norm : ArbFloat;
 
302
begin
 
303
  if (n<1) or (k1<1) or (k2<k1) or (k2>n) or (l<0) or (l>n-1) then
 
304
    begin
 
305
      term:=3; exit
 
306
    end; {wrong input}
 
307
  plam:=@lam; px:=@x; pa:=@a; getmem(v, n*sizeof(ArbFloat));
 
308
  uwidth:=l+1; getmem(u, n*uwidth*sizeof(ArbFloat));
 
309
  eigbs2(a, n, l, k1, k2, plam^[1], term);
 
310
  { kijk of norm(A-lambda.I)=0 }
 
311
  { zo ja, lever dan de eenheidsvectoren e(k1) t/m e(k2) af }
 
312
  norm:=0; j:=1;
 
313
  for i:=1 to n do
 
314
  begin
 
315
      if i<=l then m:=i else m:=l+1; s:=0;
 
316
      for k:=j to j+m-1 do
 
317
      if k=j+m-1 then s:=s+abs(pa^[k]-plam^[1]) else s:=s+abs(pa^[k]);
 
318
      if s>norm then norm:=s;
 
319
      j:=j+m
 
320
  end;
 
321
  if norm=0 then
 
322
  begin
 
323
      for i:=k1 to k2 do for j:=1 to n do
 
324
      if j=i then px^[(j-1)*rwidthx+i-k1+1]:=1
 
325
      else px^[(j-1)*rwidthx+i-k1+1]:=0;
 
326
      freemem(v, n*sizeof(ArbFloat)); freemem(u, n*uwidth*sizeof(ArbFloat));
 
327
      m2:=k2; term:=1; exit
 
328
  end;
 
329
  transf(a, n, l, u^[1], uwidth);
 
330
  i:=k1; m2:=k1-1;
 
331
  while (i <= k2) and (term=1) do
 
332
    begin
 
333
      bandev(u^[1], n, l, uwidth, plam^[i-k1+1], v^[1], term);
 
334
      if term=1 then
 
335
        begin
 
336
          m2:=i; for j:=1 to n do px^[(j-1)*rwidthx+i-k1+1]:=v^[j]
 
337
        end;
 
338
      i:=i+1
 
339
    end; {i}
 
340
  freemem(v, n*sizeof(ArbFloat));
 
341
  freemem(u, n*uwidth*sizeof(ArbFloat));
 
342
end; {eigbs4}
 
343
 
 
344
procedure eigge1(var a: ArbFloat; n, rwidth: ArbInt; var lam: complex;
 
345
                 var term: ArbInt);
 
346
var pa, h, dummy : ^arfloat1;
 
347
           i, ns : ArbInt;
 
348
begin
 
349
  if n<1 then
 
350
    begin
 
351
      term:=3; exit
 
352
    end;
 
353
  ns:=n*sizeof(ArbFloat); pa:=@a;
 
354
  getmem(dummy, ns); getmem(h, n*ns);
 
355
  for i:=1 to n do move(pa^[(i-1)*rwidth+1], h^[(i-1)*n+1], ns);
 
356
  orthes(h^[1], n, n, dummy^[1]);
 
357
  hessva(h^[1], n, n, lam, term);
 
358
  freemem(dummy, ns); freemem(h, n*ns);
 
359
end;  {eigge1}
 
360
 
 
361
procedure eigge3(var a: ArbFloat; n, rwidtha: ArbInt; var lam, x: complex;
 
362
                 rwidthx: ArbInt; var term: ArbInt);
 
363
var     pa, pd, u, v: ^arfloat1;
 
364
    m1, m2, i, j, ns: ArbInt;
 
365
begin
 
366
  if n<1 then
 
367
    begin
 
368
      term:=3; exit
 
369
    end;
 
370
  ns:=n*sizeof(ArbFloat); getmem(pd, ns); getmem(u, n*ns); getmem(v, n*ns);
 
371
  pa:=@a; for i:=1 to n do move(pa^[(i-1)*rwidtha+1], u^[(i-1)*n+1], ns);
 
372
  balance(u^[1], n, n, m1, m2, pd^[1]);
 
373
  orttrans(u^[1], n, n, v^[1], n);
 
374
  hessvec(u^[1], n, n, lam, v^[1], n, term);
 
375
  if term=1 then
 
376
    begin
 
377
      balback(pd^[1], n, m1, m2, 1, n, v^[1], n);
 
378
      normeer(lam, n, v^[1], n);
 
379
      transx(v^[1], n, n, lam, x, rwidthx)
 
380
    end;
 
381
  freemem(pd, ns); freemem(u, n*ns); freemem(v, n*ns);
 
382
end;  {eigge3}
 
383
 
 
384
procedure eiggg1(var a: ArbFloat; n, rwidtha: ArbInt; var b: ArbFloat;
 
385
                 rwidthb: ArbInt; var lam: ArbFloat; var term: ArbInt);
 
386
var u, v, pa, pb : ^arfloat1;
 
387
        i, j, ns : ArbInt;
 
388
begin
 
389
  if n<1 then
 
390
    begin
 
391
      term:=3; exit
 
392
    end;
 
393
  pa:=@a; pb:=@b; ns:=n*sizeof(ArbFloat); getmem(u, n*ns); getmem(v, n*ns);
 
394
  for i:=1 to n do move(pa^[(i-1)*rwidtha+1], u^[(i-1)*n+1], ns);
 
395
  for i:=1 to n do move(pb^[(i-1)*rwidthb+1], v^[(i-1)*n+1], ns);
 
396
  reduc1(u^[1], n, n, v^[1], n, term);
 
397
  if term=1 then eiggs1(u^[1], n, n, lam, term);
 
398
  freemem(u, n*ns); freemem(v, n*ns);
 
399
end; {eiggg1}
 
400
 
 
401
procedure eiggg2(var a: ArbFloat; n, rwidtha, k1, k2: ArbInt; var b: ArbFloat;
 
402
                 rwidthb: ArbInt; var lam: ArbFloat; var term: ArbInt);
 
403
var u, v, pa, pb : ^arfloat1;
 
404
        i, j, ns : ArbInt;
 
405
begin
 
406
  if (n<1) or (k1<1) or (k2<k1) or (k2>n) then
 
407
    begin
 
408
      term:=3; exit
 
409
    end;
 
410
  pa:=@a; pb:=@b; ns:=n*sizeof(ArbFloat); getmem(u, n*ns); getmem(v, n*ns);
 
411
  for i:=1 to n do move(pa^[(i-1)*rwidtha+1], u^[(i-1)*n+1], ns);
 
412
  for i:=1 to n do move(pb^[(i-1)*rwidthb+1], v^[(i-1)*n+1], ns);
 
413
  reduc1(u^[1], n, n, v^[1], n, term);
 
414
  if term=1 then eiggs2(u^[1], n, n, k1, k2, lam, term);
 
415
  freemem(u, n*ns); freemem(v, n*ns)
 
416
end; {eiggg2}
 
417
 
 
418
procedure eiggg3(var a: ArbFloat; n, rwidtha: ArbInt; var b: ArbFloat;
 
419
                 rwidthb: ArbInt; var lam, x: ArbFloat; rwidthx: ArbInt;
 
420
                 var term: ArbInt);
 
421
var u, v, pa, pb : ^arfloat1;
 
422
        i, j, ns : ArbInt;
 
423
begin
 
424
  if n<1 then
 
425
    begin
 
426
      term:=3; exit
 
427
    end;
 
428
  pa:=@a; pb:=@b;
 
429
  ns:=n*sizeof(ArbFloat);
 
430
  getmem(u, n*ns); getmem(v, n*ns);
 
431
  for i:=1 to n do move(pa^[(i-1)*rwidtha+1], u^[(i-1)*n+1], ns);
 
432
  for i:=1 to n do move(pb^[(i-1)*rwidthb+1], v^[(i-1)*n+1], ns);
 
433
  reduc1(u^[1], n, n, v^[1], n, term);
 
434
  if term=1 then
 
435
    begin
 
436
      eiggs3(u^[1], n, n, lam, x, rwidthx, term);
 
437
      if term=1 then rebaka(v^[1], n, n, 1, n, x, rwidthx, term)
 
438
    end;
 
439
  freemem(u, n*ns); freemem(v, n*ns)
 
440
end; {eiggg3}
 
441
 
 
442
procedure eiggg4(var a: ArbFloat; n, rwidtha, k1, k2: ArbInt; var b: ArbFloat;
 
443
                 rwidthb: ArbInt; var lam, x: ArbFloat; rwidthx: ArbInt;
 
444
                 var m2, term: ArbInt);
 
445
 
 
446
var u, v, pa, pb : ^arfloat1;
 
447
     i, j, ns, t : ArbInt;
 
448
begin
 
449
  if (n<1) or (k1<1) or (k2<k1) or (k2>n) then
 
450
    begin
 
451
      term:=3; exit
 
452
    end;
 
453
  pa:=@a; pb:=@b; ns:=n*sizeof(ArbFloat); getmem(u, n*ns); getmem(v, n*ns);
 
454
  for i:=1 to n do move(pa^[(i-1)*rwidtha+1], u^[(i-1)*n+1], ns);
 
455
  for i:=1 to n do move(pb^[(i-1)*rwidthb+1], v^[(i-1)*n+1], ns);
 
456
  reduc1(u^[1], n, n, v^[1], n, term);
 
457
  if term=1 then
 
458
    begin
 
459
      eiggs4(u^[1], n, n, k1, k2, lam, x, rwidthx, m2, term);
 
460
      if m2 < k2 then term:=4;
 
461
      if m2 > k1-1 then
 
462
        begin
 
463
          rebaka(v^[1], n, n, k1, m2, x, rwidthx, t);
 
464
          if t=2 then
 
465
            begin
 
466
              term:=4; m2:=k1-1
 
467
            end
 
468
        end
 
469
    end;
 
470
  freemem(u, n*ns); freemem(v, n*ns)
 
471
end; {eiggg4}
 
472
 
 
473
procedure eigsv1(var a: ArbFloat; m, n, rwidth: ArbInt; var sig: ArbFloat;
 
474
                 var term: ArbInt);
 
475
 
 
476
var                     pa, pq, u, e : ^arfloat1;
 
477
          i, j, k, l, ns, ii, jj, kk : ArbInt;
 
478
 c, f, g, h, p, s, x, y, z, eps, tol : ArbFloat;
 
479
                  conv, goon, cancel : boolean;
 
480
begin
 
481
  if (n<1) or (m<n) then
 
482
    begin
 
483
      term:=3; exit
 
484
    end;
 
485
  pa:=@a; pq:=@sig; term:=1;
 
486
  ns:=n*sizeof(ArbFloat); getmem(e, ns); getmem(u, m*ns);
 
487
  for i:=1 to m do move(pa^[(i-1)*rwidth+1], u^[(i-1)*n+1], ns);
 
488
  g:=0; x:=0; tol:=midget/macheps;
 
489
  for i:=1 to n do
 
490
    begin
 
491
      ii:=(i-1)*n; e^[i]:=g;
 
492
      s:=0; for j:=i to m do s:=s+sqr(u^[(j-1)*n+i]);
 
493
      if s<tol then g:=0 else
 
494
        begin
 
495
          f:=u^[ii+i]; if f<0 then g:=sqrt(s) else g:=-sqrt(s);
 
496
          h:=f*g-s; u^[ii+i]:=f-g;
 
497
          for j:=i+1 to n do
 
498
            begin
 
499
              s:=0;
 
500
              for k:=i to m do
 
501
                begin
 
502
                  kk:=(k-1)*n; s:=s+u^[kk+i]*u^[kk+j]
 
503
                end; {k}
 
504
              f:=s/h;
 
505
              for k:=i to m do
 
506
                begin
 
507
                  kk:=(k-1)*n; u^[kk+j]:=u^[kk+j]+f*u^[kk+i]
 
508
                end {k}
 
509
            end {j}
 
510
        end; {s}
 
511
      pq^[i]:=g; s:=0;
 
512
      for j:=i+1 to n do s:=s+sqr(u^[ii+j]);
 
513
      if s < tol then g:=0 else
 
514
        begin
 
515
          f:=u^[ii+i+1]; if f < 0 then g:=sqrt(s) else g:=-sqrt(s);
 
516
          h:=f*g-s; u^[ii+i+1]:=f-g;
 
517
          for j:=i+1 to n do e^[j]:=u^[ii+j]/h;
 
518
          for j:=i+1 to m do
 
519
            begin
 
520
              s:=0; jj:=(j-1)*n;
 
521
              for k:=i+1 to n do s:=s+u^[jj+k]*u^[ii+k];
 
522
              for k:=i+1 to n do u^[jj+k]:=u^[jj+k]+s*e^[k]
 
523
            end {j}
 
524
        end; {s}
 
525
      y:=abs(pq^[i])+abs(e^[i]); if y > x then x:=y
 
526
    end; {i}
 
527
  eps:=macheps*x;
 
528
  for k:=n downto 1 do
 
529
    begin
 
530
      conv:=false;
 
531
      repeat
 
532
        l:=k; goon:=true;
 
533
        while goon do
 
534
          begin
 
535
            if abs(e^[l]) <= eps then
 
536
              begin
 
537
                goon:=false; cancel:=false
 
538
              end else
 
539
            if abs(pq^[l-1]) <= eps then
 
540
              begin
 
541
                goon:=false; cancel:=true
 
542
              end
 
543
            else l:=l-1
 
544
          end; {goon}
 
545
        if cancel then
 
546
          begin
 
547
            c:=0; s:=1;
 
548
            i:=l; goon:=true;
 
549
            while goon do
 
550
              begin
 
551
                f:=s*e^[i]; e^[i]:=c*e^[i]; goon:=abs(f) > eps;
 
552
                if goon then
 
553
                  begin
 
554
                    g:=pq^[i]; h:=sqrt(f*f+g*g); pq^[i]:=h;
 
555
                    c:=g/h; s:=-f/h;
 
556
                    i:=i+1; goon:=i <= k
 
557
                  end {goon}
 
558
              end {while goon}
 
559
          end; {cancel}
 
560
        z:=pq^[k];
 
561
        if k=l then conv:=true else
 
562
          begin
 
563
            x:=pq^[l]; y:=pq^[k-1]; g:=e^[k-1]; h:=e^[k];
 
564
            f:=((y-z)*(y+z)+(g-h)*(g+h))/(2*h*y); g:=sqrt(f*f+1);
 
565
            if f < 0 then s:=f-g else s:=f+g;
 
566
            f:=((x-z)*(x+z)+h*(y/s-h))/x;
 
567
            c:=1; s:=1;
 
568
            for i:=l+1 to k do
 
569
              begin
 
570
                g:=e^[i]; y:=pq^[i]; h:=s*g; g:=c*g;
 
571
                z:=sqrt(f*f+h*h); e^[i-1]:=z; c:=f/z; s:=h/z;
 
572
                f:=x*c+g*s; g:=-x*s+g*c; h:=y*s; y:=y*c;
 
573
                z:=sqrt(f*f+h*h); pq^[i-1]:=z; c:=f/z; s:=h/z;
 
574
                f:=c*g+s*y; x:=-s*g+c*y
 
575
              end; {i}
 
576
            e^[l]:=0; e^[k]:=f; pq^[k]:=x
 
577
          end {k <> l}
 
578
      until conv;
 
579
      if z < 0 then pq^[k]:=-z
 
580
    end; {k}
 
581
  for i:=1 to n do
 
582
    begin
 
583
      k:=i; p:=pq^[i];
 
584
      for j:=i+1 to n do
 
585
        if pq^[j] < p then
 
586
          begin
 
587
            k:=j; p:=pq^[j]
 
588
          end;
 
589
        if k <> i then
 
590
          begin
 
591
            pq^[k]:=pq^[i]; pq^[i]:=p
 
592
          end
 
593
    end; {i}
 
594
  freemem(e, ns); freemem(u, m*ns)
 
595
end; {eigsv1}
 
596
 
 
597
procedure eigsv3(var a: ArbFloat; m, n, rwidtha: ArbInt; var sig, u: ArbFloat;
 
598
                 rwidthu: ArbInt; var v: ArbFloat; rwidthv: ArbInt;
 
599
                 var term: ArbInt);
 
600
 
 
601
var                pa, pu, pq, pv, e : ^arfloat1;
 
602
          i, j, k, l, ns, ii, jj, kk : ArbInt;
 
603
 c, f, g, h, p, s, x, y, z, eps, tol : ArbFloat;
 
604
                  conv, goon, cancel : boolean;
 
605
begin
 
606
  if (n<1) or (m<n)
 
607
  then
 
608
    begin
 
609
      term:=3; exit
 
610
    end;
 
611
  pa:=@a; pu:=@u; pq:=@sig; pv:=@v; term:=1;
 
612
  ns:=n*sizeof(ArbFloat); getmem(e, ns);
 
613
  for i:=1 to m do move(pa^[(i-1)*rwidtha+1], pu^[(i-1)*rwidthu+1], ns);
 
614
  g:=0; x:=0; tol:=midget/macheps;
 
615
  for i:=1 to n do
 
616
    begin
 
617
      ii:=(i-1)*rwidthu;
 
618
      e^[i]:=g; s:=0;
 
619
      for j:=i to m do s:=s+sqr(pu^[(j-1)*rwidthu+i]);
 
620
      if s<tol then g:=0 else
 
621
        begin
 
622
          f:=pu^[ii+i]; if f<0 then g:=sqrt(s) else g:=-sqrt(s);
 
623
          h:=f*g-s; pu^[ii+i]:=f-g;
 
624
          for j:=i+1 to n do
 
625
            begin
 
626
              s:=0;
 
627
              for k:=i to m do
 
628
                begin
 
629
                  kk:=(k-1)*rwidthu; s:=s+pu^[kk+i]*pu^[kk+j]
 
630
                end; {k}
 
631
              f:=s/h;
 
632
              for k:=i to m do
 
633
                begin
 
634
                  kk:=(k-1)*rwidthu; pu^[kk+j]:=pu^[kk+j]+f*pu^[kk+i]
 
635
                end {k}
 
636
            end {j}
 
637
        end; {s}
 
638
      pq^[i]:=g; s:=0; for j:=i+1 to n do s:=s+sqr(pu^[ii+j]);
 
639
      if s < tol then g:=0 else
 
640
        begin
 
641
          f:=pu^[ii+i+1];
 
642
          if f < 0 then g:=sqrt(s) else g:=-sqrt(s);
 
643
          h:=f*g-s; pu^[ii+i+1]:=f-g;
 
644
          for j:=i+1 to n do e^[j]:=pu^[ii+j]/h;
 
645
          for j:=i+1 to m do
 
646
            begin
 
647
              s:=0; jj:=(j-1)*rwidthu;
 
648
              for k:=i+1 to n do s:=s+pu^[jj+k]*pu^[ii+k];
 
649
              for k:=i+1 to n do pu^[jj+k]:=pu^[jj+k]+s*e^[k]
 
650
            end {j}
 
651
        end; {s}
 
652
      y:=abs(pq^[i])+abs(e^[i]); if y > x then x:=y
 
653
    end; {i}
 
654
  for i:=n downto 1 do
 
655
    begin
 
656
      ii:=(i-1)*rwidthu;
 
657
      if g <> 0 then
 
658
        begin
 
659
          h:=pu^[ii+i+1]*g;
 
660
          for j:=i+1 to n do pv^[(j-1)*rwidthv+i]:=pu^[ii+j]/h;
 
661
          for j:=i+1 to n do
 
662
            begin
 
663
              s:=0;
 
664
              for k:=i+1 to n do s:=s+pu^[ii+k]*pv^[(k-1)*rwidthv+j];
 
665
              for k:=i+1 to n do
 
666
                begin
 
667
                  kk:=(k-1)*rwidthv; pv^[kk+j]:=pv^[kk+j]+s*pv^[kk+i]
 
668
                end {k}
 
669
            end {j}
 
670
        end; {g}
 
671
      ii:=(i-1)*rwidthv;
 
672
      for j:=i+1 to n do
 
673
        begin
 
674
          pv^[ii+j]:=0; pv^[(j-1)*rwidthv+i]:=0
 
675
        end; {j}
 
676
      pv^[ii+i]:=1; g:=e^[i]
 
677
    end; {i}
 
678
  for i:=n downto 1 do
 
679
    begin
 
680
      g:=pq^[i]; ii:=(i-1)*rwidthu;
 
681
      for j:=i+1 to n do pu^[ii+j]:=0;
 
682
      if g <> 0 then
 
683
        begin
 
684
          h:=pu^[ii+i]*g;
 
685
          for j:=i+1 to n do
 
686
            begin
 
687
              s:=0;
 
688
              for k:=i+1 to m do
 
689
                begin
 
690
                  kk:=(k-1)*rwidthu; s:=s+pu^[kk+i]*pu^[kk+j]
 
691
                end; {k}
 
692
              f:=s/h;
 
693
              for k:=i to m do
 
694
                begin
 
695
                  kk:=(k-1)*rwidthu;
 
696
                  pu^[kk+j]:=pu^[kk+j]+f*pu^[kk+i]
 
697
                end {k}
 
698
            end; {j}
 
699
          for j:=i to m do
 
700
            begin
 
701
              jj:=(j-1)*rwidthu+i; pu^[jj]:=pu^[jj]/g
 
702
            end {j}
 
703
        end {g}
 
704
      else
 
705
        for j:=i to m do pu^[(j-1)*rwidthu+i]:=0;
 
706
      pu^[ii+i]:=pu^[ii+i]+1
 
707
    end; {i}
 
708
  eps:=macheps*x;
 
709
  for k:=n downto 1 do
 
710
    begin
 
711
      conv:=false;
 
712
      repeat
 
713
        l:=k; goon:=true;
 
714
        while goon do
 
715
          begin
 
716
            if abs(e^[l]) <= eps then
 
717
              begin
 
718
                goon:=false; cancel:=false
 
719
              end else
 
720
            if abs(pq^[l-1]) <= eps then
 
721
              begin
 
722
                goon:=false; cancel:=true
 
723
              end else l:=l-1
 
724
          end; {goon}
 
725
        if cancel then
 
726
          begin
 
727
            c:=0; s:=1; i:=l; goon:=true;
 
728
            while goon do
 
729
              begin
 
730
                f:=s*e^[i]; e^[i]:=c*e^[i]; goon:=abs(f) > eps;
 
731
                if goon then
 
732
                  begin
 
733
                    g:=pq^[i]; h:=sqrt(f*f+g*g); pq^[i]:=h;
 
734
                    c:=g/h; s:=-f/h;
 
735
                    for j:=1 to m do
 
736
                      begin
 
737
                        jj:=(j-1)*rwidthu; y:=pu^[jj+l-1]; z:=pu^[jj+i];
 
738
                        pu^[jj+l-1]:=y*c+z*s; pu^[jj+i]:=-y*s+z*c
 
739
                      end; {j}
 
740
                    i:=i+1; goon:=i <= k
 
741
                  end {goon}
 
742
              end {while goon}
 
743
          end; {cancel}
 
744
        z:=pq^[k]; if k=l then conv:=true else
 
745
          begin
 
746
            x:=pq^[l]; y:=pq^[k-1]; g:=e^[k-1]; h:=e^[k];
 
747
            f:=((y-z)*(y+z)+(g-h)*(g+h))/(2*h*y); g:=sqrt(f*f+1);
 
748
            if f < 0 then s:=f-g else s:=f+g;
 
749
            f:=((x-z)*(x+z)+h*(y/s-h))/x;
 
750
            c:=1; s:=1;
 
751
            for i:=l+1 to k do
 
752
              begin
 
753
                g:=e^[i]; y:=pq^[i]; h:=s*g; g:=c*g;
 
754
                z:=sqrt(f*f+h*h); e^[i-1]:=z; c:=f/z; s:=h/z;
 
755
                f:=x*c+g*s; g:=-x*s+g*c; h:=y*s; y:=y*c;
 
756
                for j:=1 to n do
 
757
                  begin
 
758
                    jj:=(j-1)*rwidthv;
 
759
                    x:=pv^[jj+i-1]; z:=pv^[jj+i];
 
760
                    pv^[jj+i-1]:=x*c+z*s; pv^[jj+i]:=-x*s+z*c
 
761
                  end; {j}
 
762
                z:=sqrt(f*f+h*h); pq^[i-1]:=z; c:=f/z; s:=h/z;
 
763
                f:=c*g+s*y; x:=-s*g+c*y;
 
764
                for j:=1 to m do
 
765
                  begin
 
766
                    jj:=(j-1)*rwidthu;
 
767
                    y:=pu^[jj+i-1]; z:=pu^[jj+i];
 
768
                    pu^[jj+i-1]:=y*c+z*s; pu^[jj+i]:=-y*s+z*c
 
769
                  end {j}
 
770
              end; {i}
 
771
            e^[l]:=0; e^[k]:=f; pq^[k]:=x
 
772
          end {k <> l}
 
773
      until conv;
 
774
      if z < 0 then
 
775
        begin
 
776
          pq^[k]:=-z;
 
777
          for j:=1 to n do
 
778
            begin
 
779
              jj:=(j-1)*rwidthv+k; pv^[jj]:=-pv^[jj]
 
780
            end {j}
 
781
        end {z}
 
782
    end; {k}
 
783
  for i:=1 to n do
 
784
    begin
 
785
      k:=i; p:=pq^[i];
 
786
      for j:=i+1 to n do
 
787
        if pq^[j] < p then
 
788
          begin
 
789
            k:=j; p:=pq^[j]
 
790
          end;
 
791
        if k <> i then
 
792
          begin
 
793
            pq^[k]:=pq^[i]; pq^[i]:=p;
 
794
            for j:=1 to m do
 
795
              begin
 
796
                jj:=(j-1)*rwidthu;
 
797
                p:=pu^[jj+i]; pu^[jj+i]:=pu^[jj+k]; pu^[jj+k]:=p;
 
798
              end;
 
799
            for j:=1 to n do
 
800
              begin
 
801
                jj:=(j-1)*rwidthv;
 
802
                p:=pv^[jj+i]; pv^[jj+i]:=pv^[jj+k]; pv^[jj+k]:=p
 
803
              end { interchange in u and v column i with comlumn k }
 
804
          end
 
805
    end; {i}
 
806
  freemem(e, ns)
 
807
end; {eigsv3}
 
808
end.