home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / TOEPLZ.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1991-04-29  |  1.7 KB  |  63 lines
  1. PROCEDURE toeplz(r: gltwon; VAR x: glnarray; y: glnarray; n: integer);
  2. (* Programs using routine TOEPLZ must define the types
  3.    glnarray = ARRAY [1..n] OF real;
  4.    gltwon = ARRAY [1..2*n] OF real;
  5. in the main routine. *)
  6. LABEL 98,99;
  7. VAR
  8.    m2,m1,m,k,j: integer;
  9.    sxn,shn,sgn,sgd,sd,qt2,qt1,qq,pt2,pt1,pp: real;
  10.    g,h: glnarray;
  11. BEGIN
  12.    IF (r[n] = 0.0) THEN  GOTO 99;
  13.    x[1] := y[1]/r[n];
  14.    IF (n = 1) THEN GOTO 99;
  15.    g[1] := r[n-1]/r[n];
  16.    h[1] := r[n+1]/r[n];
  17.    FOR m := 1 TO n DO BEGIN
  18.       m1 := m+1;
  19.       sxn := -y[m1];
  20.       sd := -r[n];
  21.       FOR j := 1 TO m DO BEGIN
  22.          sxn := sxn+r[n+m1-j]*x[j];
  23.          sd := sd+r[n+m1-j]*g[m-j+1]
  24.       END;
  25.       IF (sd = 0.0) THEN GOTO 98;
  26.       x[m1] := sxn/sd;
  27.       FOR j := 1 TO m DO BEGIN
  28.          x[j] := x[j]-x[m1]*g[m-j+1]
  29.       END;
  30.       IF (m1 = n) THEN GOTO 99;
  31.       sgn := -r[n-m1];
  32.       shn := -r[n+m1];
  33.       sgd := -r[n];
  34.       FOR j := 1 TO m DO BEGIN
  35.          sgn := sgn+r[n+j-m1]*g[j];
  36.          shn := shn+r[n+m1-j]*h[j];
  37.          sgd := sgd+r[n+j-m1]*h[m-j+1]
  38.       END;
  39.       IF ((sd = 0.0) OR (sgd = 0.0)) THEN GOTO 98;
  40.       g[m1] := sgn/sgd;
  41.       h[m1] := shn/sd;
  42.       k := m;
  43.       m2 := (m+1) DIV 2;
  44.       pp := g[m1];
  45.       qq := h[m1];
  46.       FOR j := 1 TO m2 DO BEGIN
  47.          pt1 := g[j];
  48.          pt2 := g[k];
  49.          qt1 := h[j];
  50.          qt2 := h[k];
  51.          g[j] := pt1-pp*qt2;
  52.          g[k] := pt2-pp*qt1;
  53.          h[j] := qt1-qq*pt2;
  54.          h[k] := qt2-qq*pt1;
  55.          k := k-1
  56.       END
  57.    END;
  58.    writeln('pause in TOEPLZ - should not arrive here!'); readln;
  59.    GOTO 99;
  60. 98:   writeln('pause in TOEPLZ - Levinson method fails');
  61.    writeln('matrix has a singular principal minor'); readln;
  62. 99:   END;
  63.