home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / LUBKSB.DEM < prev    next >
Encoding:
Text File  |  1991-04-29  |  1.8 KB  |  73 lines
  1. PROGRAM d2r3(input,output,dfile);
  2. (* driver for routine LUBKSB *)
  3. LABEL 10,99;
  4. CONST
  5.    np=20;
  6. TYPE
  7.    glnpbynp=ARRAY [1..np,1..np] OF real;
  8.    glnarray=ARRAY [1..np] OF real;
  9.    glindx=ARRAY [1..np] OF integer;
  10. VAR
  11.    j,k,l,m,n : integer;
  12.    p : real;
  13.    a,b,c : glnpbynp;
  14.    indx : glindx;
  15.    x : glnarray;
  16.    dfile : text;
  17.  
  18. (*$I MODFILE.PAS *)
  19. (*$I LUDCMP.PAS *)
  20.  
  21. (*$I LUBKSB.PAS *)
  22.  
  23. BEGIN
  24.    glopen(dfile,'matrx1.dat');
  25. 10:   readln(dfile);
  26.    readln(dfile);
  27.    readln(dfile,n,m);
  28.    readln(dfile);
  29.    FOR k := 1 to n DO BEGIN
  30.       FOR l := 1 to n-1 DO read(dfile,a[k,l]);
  31.       readln(dfile,a[k,n])
  32.    END;
  33.    readln(dfile);
  34.    FOR l := 1 to m DO BEGIN
  35.       FOR k := 1 to n-1 DO read(dfile,b[k,l]);
  36.       readln(dfile,b[n,l])
  37.    END;
  38. (* save matrix a for later testing *)
  39.    FOR l := 1 to n DO BEGIN
  40.       FOR k := 1 to n DO BEGIN
  41.          c[k,l] := a[k,l]
  42.       END
  43.    END;
  44. (* do lu decomposition *)
  45.    ludcmp(c,n,np,indx,p);
  46. (* solve equations for each right-hand vector *)
  47.    FOR k := 1 to m DO BEGIN
  48.       FOR l := 1 to n DO BEGIN
  49.          x[l] := b[l,k]
  50.       END;
  51.       lubksb(c,n,np,indx,x);
  52. (* test results with original matrix *)
  53.       writeln('right-hand side vector:');
  54.       FOR l := 1 to n-1 DO write(b[l,k]:12:6);
  55.       writeln(b[n,k]:12:6); 
  56.       writeln ('result of matrix applied to sol''n vector');
  57.       FOR l := 1 to n DO BEGIN
  58.          b[l,k] := 0.0;
  59.          FOR j := 1 to n DO BEGIN
  60.             b[l,k] := b[l,k]+a[l,j]*x[j]
  61.          END
  62.       END;
  63.       FOR l := 1 to n-1 DO write(b[l,k]:12:6);
  64.       writeln(b[n,k]:12:6);
  65.       writeln('***********************************')
  66.    END;
  67.    IF eof(dfile) THEN GOTO 99;
  68.    writeln('press RETURN for next problem:');
  69.    readln;
  70.    GOTO 10;
  71. 99:   close(dfile)
  72. END.
  73.