← Back to branch summary

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

« back to all changes in this revision

Viewing changes to packages/extra/numlib/dsl.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:
packages/extra/numlib/dsl.pas
1
 
{
2
 
    $Id: dsl.pas,v 1.2 2002/09/07 15:43:01 peter Exp $
3
 
    This file is part of the Numlib package.
4
 
    Copyright (c) 1986-2000 by
5
 
     Kees van Ginneken, Wil Kortsmit and Loek van Reij of the
6
 
     Computational centre of the Eindhoven University of Technology
7
 
 
8
 
    FPC port Code          by Marco van de Voort (marco@freepascal.org)
9
 
             documentation by Michael van Canneyt (Michael@freepascal.org)
10
 
 
11
 
    Unknown unit. There doesn't exist any documentation for it, it isn't
12
 
    commented, and I don't recognize the algortism directly.
13
 
    I added some comments, since suffixes of the procedures seem to indicate
14
 
    some features of the matrixtype (from unit SLE)
15
 
    So probably Some pivot matrix?
16
 
 
17
 
    This code was probably internal in older libs, and only exported
18
 
    in later versions.
19
 
 
20
 
    See the file COPYING.FPC, included in this distribution,
21
 
    for details about the copyright.
22
 
 
23
 
    This program is distributed in the hope that it will be useful,
24
 
    but WITHOUT ANY WARRANTY; without even the implied warranty of
25
 
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
26
 
 
27
 
 **********************************************************************}
28
 
Unit dsl;
29
 
 
30
 
interface
31
 
{$I DIRECT.INC}
32
 
 
33
 
 
34
 
uses typ;
35
 
 
36
 
{Gen=generic, matrix without special or unknown ordering}
37
 
Procedure dslgen(n, rwidth: ArbInt; Var alu: ArbFloat; Var p: ArbInt;
38
 
                 Var b, x: ArbFloat; Var term: ArbInt);
39
 
 
40
 
{"tridiagonal matrix"}
41
 
Procedure dslgtr(n: ArbInt; Var l1, d1, u1, u2: ArbFloat;
42
 
                 Var p: boolean; Var b, x: ArbFloat; Var term: ArbInt);
43
 
 
44
 
{Symmetrical matrix}
45
 
Procedure dslgsy(n, rwidth: ArbInt; Var alt: ArbFloat; Var p: ArbInt;
46
 
                 Var q: boolean; Var b, x: ArbFloat; Var term: ArbInt);
47
 
 
48
 
{Symmetrical positive definitive matrix}
49
 
Procedure dslgpd(n, rwidth: ArbInt; Var al, b, x: ArbFloat;
50
 
                 Var term: ArbInt);
51
 
 
52
 
{Generic "band" matrix}
53
 
Procedure dslgba(n, lb, rb, rwa: ArbInt; Var au: ArbFloat; rwl: ArbInt;
54
 
                 Var l: ArbFloat; Var p: ArbInt; Var b, x: ArbFloat;
55
 
                 Var term: ArbInt);
56
 
 
57
 
{Positive definite bandmatrix}
58
 
Procedure dslgpb(n, lb, rwidth: ArbInt; Var al, b, x: ArbFloat;
59
 
                 Var term: ArbInt);
60
 
 
61
 
{Special tridiagonal matrix}
62
 
Procedure dsldtr(n:ArbInt; Var l, d, u, b, x: ArbFloat; Var term: ArbInt);
63
 
 
64
 
implementation
65
 
 
66
 
Procedure dslgen(n, rwidth: ArbInt; Var alu: ArbFloat; Var p: ArbInt;
67
 
                 Var b, x: ArbFloat; Var term: ArbInt);
68
 
 
69
 
Var
70
 
                          success : boolean;
71
 
    indk, j, k, indexpivot, kmin1 : ArbInt;
72
 
                      h, pivot, s : ArbFloat;
73
 
                               pp : ^arint1;
74
 
                     palu, pb, px : ^arfloat1;
75
 
 
76
 
Begin
77
 
  If (n<1) Or (rwidth<1) Then
78
 
    Begin
79
 
      term := 3;
80
 
     exit
81
 
    End; {wrong input}
82
 
  pp := @p;
83
 
 palu := @alu;
84
 
 pb := @b;
85
 
 px := @x;
86
 
  move(pb^, px^, n*sizeof(ArbFloat));
87
 
  For k:=1 To n Do
88
 
    Begin
89
 
      indexpivot := pp^[k];
90
 
      If indexpivot  <> k Then
91
 
        Begin
92
 
          h := px^[k];
93
 
         px^[k] := px^[indexpivot];
94
 
          px^[indexpivot] := h
95
 
        End {indexpivot <> k}
96
 
    End; {k}
97
 
  For k:=2 To n Do
98
 
    Begin
99
 
      s := px^[k];
100
 
     kmin1 := k-1;
101
 
      For j:=1 To kmin1 Do
102
 
        s := s-palu^[(k-1)*rwidth+j]*px^[j];
103
 
      px^[k] := s
104
 
    End; {k}
105
 
  success := true;
106
 
 k := n+1;
107
 
  while (k>1) and success Do
108
 
    Begin
109
 
      k := k-1;
110
 
     indk := (k-1)*rwidth;
111
 
      pivot := palu^[indk+k];
112
 
      If pivot=0 Then
113
 
        success := false
114
 
      Else
115
 
        Begin
116
 
          s := px^[k];
117
 
          For j:=k+1 To n Do
118
 
            s := s-palu^[indk+j]*px^[j];
119
 
          px^[k] := s/pivot
120
 
        End {pivot <> 0}
121
 
    End; {k}
122
 
  If success Then
123
 
    term := 1
124
 
  Else
125
 
    term := 2
126
 
End; {dslgen}
127
 
 
128
 
Procedure dslgtr(n: ArbInt; Var l1, d1, u1, u2: ArbFloat;
129
 
                 Var p: boolean; Var b, x: ArbFloat; Var term: ArbInt);
130
 
 
131
 
Var
132
 
                    i, j, nmin1 : ArbInt;
133
 
                          h, di : ArbFloat;
134
 
                        success : boolean;
135
 
          pd1, pu1, pu2, pb, px : ^arfloat1;
136
 
                            pl1 : ^arfloat2;
137
 
                             pp : ^arbool1;
138
 
Begin
139
 
  If n<1 Then
140
 
    Begin
141
 
      term := 3;
142
 
     exit
143
 
    End; {wrong input}
144
 
  pl1 := @l1;  pd1 := @d1;  pu1 := @u1; pu2 := @u2; pb := @b;  px := @x;
145
 
 pp := @p;
146
 
  move(pb^, px^, n*sizeof(ArbFloat));
147
 
  success := true;
148
 
 i := 0;
149
 
  while (i<>n) and success Do
150
 
    Begin
151
 
      i := i+1;
152
 
     success := pd1^[i]<>0
153
 
    End; {i}
154
 
  If success Then
155
 
    Begin
156
 
      nmin1 := n-1;
157
 
     j := 1;
158
 
      while j <> n Do
159
 
        Begin
160
 
          i := j;
161
 
         j := j+1;
162
 
          If pp^[i] Then
163
 
            Begin
164
 
              h := px^[i];
165
 
             px^[i] := px^[j];
166
 
             px^[j] := h-pl1^[j]*px^[i]
167
 
            End {pp^[i]}
168
 
          Else
169
 
            px^[j] := px^[j]-pl1^[j]*px^[i]
170
 
        End;  {j}
171
 
      di := pd1^[n];
172
 
      px^[n] := px^[n]/di;
173
 
      If n > 1 Then
174
 
        Begin
175
 
          di := pd1^[nmin1];
176
 
          px^[nmin1] := (px^[nmin1]-pu1^[nmin1]*px^[n])/di
177
 
        End; {n > 1}
178
 
      For i:=n-2 Downto 1 Do
179
 
        Begin
180
 
          di := pd1^[i];
181
 
          px^[i] := (px^[i]-pu1^[i]*px^[i+1]-pu2^[i]*px^[i+2])/di
182
 
        End {i}
183
 
    End; {success}
184
 
  If success Then
185
 
    term := 1
186
 
  Else
187
 
    term := 2
188
 
End; {dslgtr}
189
 
 
190
 
Procedure dslgsy(n, rwidth: ArbInt; Var alt: ArbFloat; Var p: ArbInt;
191
 
                 Var q: boolean; Var b, x: ArbFloat; Var term: ArbInt);
192
 
 
193
 
Var
194
 
    i, indexpivot, imin1, j, jmin1, iplus1, imin2, ns, ii  : ArbInt;
195
 
                                          success, regular : boolean;
196
 
                                                 h, ct, di : ArbFloat;
197
 
                               palt, pb, px, y, l, d, u, v : ^arfloat1;
198
 
                                                        pp : ^arint1;
199
 
                                                        pq : ^arbool1;
200
 
 
201
 
Begin
202
 
  If (n<1) Or (rwidth<1) Then
203
 
    Begin
204
 
      term := 3;
205
 
     exit
206
 
    End; {wrong input}
207
 
  palt := @alt;
208
 
 pp := @p;
209
 
 pq := @q;
210
 
 pb := @b;
211
 
 px := @x;
212
 
  ns := n*sizeof(ArbFloat);
213
 
  getmem(l, ns);
214
 
  getmem(d, ns);
215
 
  getmem(u, ns);
216
 
  getmem(v, ns);
217
 
  getmem(y, ns);
218
 
  move(pb^, y^, ns);
219
 
  success := true;
220
 
 i := 0;
221
 
 ii := 1;
222
 
  while (i<>n) and success Do
223
 
    Begin
224
 
      i := i+1;
225
 
     success := palt^[ii]<>0;
226
 
     ii := ii+rwidth+1
227
 
    End; {i}
228
 
  If success Then
229
 
    Begin
230
 
      For i:=1 To n Do
231
 
        Begin
232
 
          indexpivot := pp^[i];
233
 
          If indexpivot <> i Then
234
 
            Begin
235
 
              h := y^[i];
236
 
             y^[i] := y^[indexpivot];
237
 
              y^[indexpivot] := h
238
 
            End {indexpivot <> i}
239
 
        End; {i}
240
 
      i := 0;
241
 
      while i<n Do
242
 
        Begin
243
 
          imin1 := i;
244
 
         i := i+1;
245
 
         j := 1;
246
 
         h := y^[i];
247
 
          while j<imin1 Do
248
 
            Begin
249
 
              jmin1 := j;
250
 
             j := j+1;
251
 
              h := h-palt^[(i-1)*rwidth+jmin1]*y^[j]
252
 
            End; {j}
253
 
          y^[i] := h
254
 
        End; {i}
255
 
      d^[1] := palt^[1];
256
 
     di := d^[1];
257
 
      If n>1 Then
258
 
        Begin
259
 
          l^[1] := palt^[rwidth+1];
260
 
         d^[2] := palt^[rwidth+2];
261
 
          di := d^[2];
262
 
          u^[1] := palt^[2]
263
 
        End; {n>1}
264
 
      imin1 := 1;
265
 
     i := 2;
266
 
      while i<n Do
267
 
        Begin
268
 
          imin2 := imin1;
269
 
         imin1 := i;
270
 
         i := i+1;
271
 
          ii := (i-1)*rwidth;
272
 
          l^[imin1] := palt^[ii+imin1];
273
 
         d^[i] := palt^[ii+i];
274
 
         di := d^[i];
275
 
          u^[imin1] := palt^[ii-rwidth+i];
276
 
         v^[imin2] := palt^[ii-2*rwidth+i]
277
 
        End; {i}
278
 
      dslgtr(n, l^[1], d^[1], u^[1], v^[1], pq^[1], y^[1], px^[1], term);
279
 
      i := n+1;
280
 
     imin1 := n;
281
 
      while i>2 Do
282
 
        Begin
283
 
          iplus1 := i;
284
 
         i := imin1;
285
 
         imin1 := imin1-1;
286
 
         h := px^[i];
287
 
          For j:=iplus1 To n Do
288
 
            h := h-palt^[(j-1)*rwidth+imin1]*px^[j];
289
 
          px^[i] := h
290
 
        End; {i}
291
 
      For i:=n Downto 1 Do
292
 
        Begin
293
 
          indexpivot := pp^[i];
294
 
          If indexpivot <> i Then
295
 
            Begin
296
 
              h := px^[i];
297
 
             px^[i] := px^[indexpivot];
298
 
              px^[indexpivot] := h
299
 
            End {indexpivot <> i}
300
 
        End {i}
301
 
   End; {success}
302
 
  If success Then
303
 
    term := 1
304
 
  Else
305
 
    term := 2;
306
 
  freemem(l, ns);
307
 
  freemem(d, ns);
308
 
  freemem(u, ns);
309
 
  freemem(v, ns);
310
 
  freemem(y, ns)
311
 
End; {dslgsy}
312
 
 
313
 
Procedure dslgpd(n, rwidth: ArbInt; Var al, b, x: ArbFloat;
314
 
                 Var term: ArbInt);
315
 
 
316
 
Var
317
 
       ii, imin1, i, j : ArbInt;
318
 
                h, lii : ArbFloat;
319
 
               success : boolean;
320
 
           pal, pb, px : ^arfloat1;
321
 
Begin
322
 
  If (n<1) Or (rwidth<1) Then
323
 
    Begin
324
 
      term := 3;
325
 
     exit
326
 
    End; {wrong input}
327
 
  pal := @al;
328
 
 pb := @b;
329
 
 px := @x;
330
 
  move(pb^, px^, n*sizeof(ArbFloat));
331
 
  success := true;
332
 
 i := 0;
333
 
 ii := 1;
334
 
  while (i<>n) and success Do
335
 
    Begin
336
 
      i := i+1;
337
 
     success := pal^[ii]<>0;
338
 
     ii := ii+rwidth+1
339
 
    End; {i}
340
 
  If success Then
341
 
    Begin
342
 
      For i:=1 To n Do
343
 
        Begin
344
 
          ii := (i-1)*rwidth;
345
 
          h := px^[i];
346
 
         imin1 := i-1;
347
 
          For j:=1 To imin1 Do
348
 
            h := h-pal^[ii+j]*px^[j];
349
 
          lii := pal^[ii+i];
350
 
          px^[i] := h/lii
351
 
        End; {i}
352
 
      For i:=n Downto 1 Do
353
 
        Begin
354
 
          h := px^[i];
355
 
          For j:=i+1 To n Do
356
 
            h := h-pal^[(j-1)*rwidth+i]*px^[j];
357
 
          px^[i] := h/pal^[(i-1)*rwidth+i]
358
 
        End {i}
359
 
    End; {success}
360
 
  If success Then
361
 
    term := 1
362
 
  Else
363
 
    term := 2
364
 
End;  {dslgpd}
365
 
 
366
 
Procedure dslgba(n, lb, rb, rwa: ArbInt; Var au: ArbFloat; rwl: ArbInt;
367
 
                 Var l: ArbFloat; Var p: ArbInt; Var b, x: ArbFloat;
368
 
                 Var term: ArbInt);
369
 
 
370
 
Var
371
 
   i, j, k, ipivot, ubi, ubj : ArbInt;
372
 
   h, pivot                  : ArbFloat;
373
 
   pau, pl, px, pb           : ^arfloat1;
374
 
   pp                        : ^arint1;
375
 
 
376
 
Begin
377
 
  If (n<1) Or (lb<0) Or (rb<0) Or (lb>n-1)
378
 
        Or (rb>n-1) Or (rwa<1) Or (rwl<0) Then
379
 
    Begin
380
 
      term := 3;
381
 
     exit
382
 
    End; {term=3}
383
 
  pau := @au;
384
 
 pl := @l;
385
 
 pb := @b;
386
 
 px := @x;
387
 
 pp := @p;
388
 
  move(pb^, px^, n*sizeof(ArbFloat));
389
 
  ubi := lb;
390
 
  For k:=1 To n Do
391
 
    Begin
392
 
      ipivot := pp^[k];
393
 
      If ipivot <> k Then
394
 
        Begin
395
 
          h := px^[k];
396
 
         px^[k] := px^[ipivot];
397
 
          px^[ipivot] := h
398
 
        End; {ipivot <> k}
399
 
      If ubi<n Then
400
 
        ubi := ubi+1;
401
 
      For i:=k+1 To ubi Do
402
 
        px^[i] := px^[i]-px^[k]*pl^[(k-1)*rwl+i-k]
403
 
    End; {k}
404
 
  ubj := 0;
405
 
 i := n;
406
 
 term := 1;
407
 
  while (i >= 1) and (term=1) Do
408
 
    Begin
409
 
      If ubj<rb+lb+1 Then
410
 
        ubj := ubj+1;
411
 
      h := px^[i];
412
 
      For j:=2 To ubj Do
413
 
        h := h-pau^[(i-1)*rwa+j]*px^[i+j-1];
414
 
      pivot := pau^[(i-1)*rwa+1];
415
 
      If pivot=0 Then
416
 
        term := 2
417
 
      Else
418
 
        px^[i] := h/pivot;
419
 
      i := i-1
420
 
    End {i}
421
 
End; {dslgba}
422
 
 
423
 
Procedure dslgpb(n, lb, rwidth: ArbInt; Var al, b, x: ArbFloat;
424
 
                 Var term: ArbInt);
425
 
 
426
 
Var
427
 
   ll, ii, llmin1, p, i, q, k : ArbInt;
428
 
            h, hh, alim       : ArbFloat;
429
 
                  pal, pb, px : ^arfloat1;
430
 
Begin
431
 
  If (lb<0) Or (lb>n-1) Or (n<1) Or (rwidth<1) Then
432
 
    Begin
433
 
      term := 3;
434
 
     exit
435
 
    End; {wrong input}
436
 
  pal := @al;
437
 
 pb := @b;
438
 
 px := @x;
439
 
  move(pb^, px^, n*sizeof(ArbFloat));
440
 
  ll := lb+1;
441
 
  llmin1 := ll-1;
442
 
 p := ll+1;
443
 
 term := 1;
444
 
 i := 1;
445
 
  while (i <= n) and (term=1) Do
446
 
    Begin
447
 
      ii := (i-1)*rwidth;
448
 
      If p>1 Then
449
 
        p := p-1;
450
 
      h := px^[i];
451
 
     q := i;
452
 
      For k:=llmin1 Downto p Do
453
 
        Begin
454
 
          q := q-1;
455
 
         h := h-pal^[ii+k]*px^[q]
456
 
        End; {k}
457
 
      alim := pal^[ii+ll];
458
 
      If alim=0 Then
459
 
        term := 2
460
 
      Else
461
 
        px^[i] := h/alim;
462
 
      i := i+1
463
 
    End; {i}
464
 
  If term=1 Then
465
 
    Begin
466
 
      p := ll+1;
467
 
      For i:=n Downto 1 Do
468
 
        Begin
469
 
          If p>1 Then
470
 
            p := p-1;
471
 
          q := i;
472
 
         h := px^[i];
473
 
          For k:=llmin1 Downto p Do
474
 
            Begin
475
 
              q := q+1;
476
 
             h := h-pal^[(q-1)*rwidth+k]*px^[q]
477
 
            End; {k}
478
 
          px^[i] := h/pal^[(i-1)*rwidth+ll]
479
 
        End {i}
480
 
    End {term=1}
481
 
End; {dslgpb}
482
 
 
483
 
Procedure dsldtr(n:ArbInt; Var l, d, u, b, x: ArbFloat; Var term: ArbInt);
484
 
 
485
 
Var
486
 
                   i, j : ArbInt;
487
 
                     di : ArbFloat;
488
 
         pd, pu, pb, px : ^arfloat1;
489
 
                     pl : ^arfloat2;
490
 
Begin
491
 
  If n<1 Then
492
 
    Begin
493
 
      term := 3;
494
 
     exit
495
 
    End; {wrong input}
496
 
  pl := @l;
497
 
 pd := @d;
498
 
 pu := @u;
499
 
 pb := @b;
500
 
 px := @x;
501
 
  move(pb^, px^, n*sizeof(ArbFloat));
502
 
  j := 1;
503
 
  while j <> n Do
504
 
    Begin
505
 
      i := j;
506
 
     j := j+1;
507
 
     px^[j] := px^[j]-pl^[j]*px^[i]
508
 
    End;
509
 
  di := pd^[n];
510
 
  If di=0 Then
511
 
    term := 2
512
 
  Else
513
 
    term := 1;
514
 
  If term=1 Then
515
 
    px^[n] := px^[n]/di;
516
 
  i := n-1;
517
 
  while (i >= 1) and (term=1) Do
518
 
    Begin
519
 
      di := pd^[i];
520
 
      If di=0 Then
521
 
        term := 2
522
 
      Else
523
 
        px^[i] := (px^[i]-pu^[i]*px^[i+1])/di;
524
 
      i := i-1
525
 
    End; {i}
526
 
End; {dsldtr}
527
 
 
528
 
End.
529
 
{
530
 
  $Log: dsl.pas,v $
531
 
  Revision 1.2  2002/09/07 15:43:01  peter
532
 
    * old logs removed and tabs fixed
533
 
 
534
 
  Revision 1.1  2002/01/29 17:55:18  peter
535
 
    * splitted to base and extra
536
 
 
537
 
}