Overview
| Comment: | Initial revision |
|---|---|
| Downloads: | Tarball | ZIP archive | SQL archive |
| Timelines: | family | ancestors | descendants | both | origin/master | trunk |
| Files: | files | file ages | folders |
| SHA3-256: |
205a8be22a72d9fe022f89478f5b868c |
| User & Date: | gawthrop@users.sourceforge.net on 1998-08-13 14:58:35 |
| Other Links: | branch diff | manifest | tags |
Context
|
1998-08-13
| ||
| 15:10:47 | Change temp var name to mtt_t1 etc check-in: 8ea77c2604 user: gawthrop@users.sourceforge.net tags: origin/master, trunk | |
| 14:58:35 | Initial revision check-in: 205a8be22a user: gawthrop@users.sourceforge.net tags: origin/master, trunk | |
| 12:43:19 | Cahnged tmp variable to mtt_1 etc check-in: 648535a8c7 user: gawthrop@users.sourceforge.net tags: origin/master, trunk | |
Changes
Added mttroot/mtt/bin/trans/p/mtt_sparse.p version [5e780d3d5b].
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 |
PROCEDURE mtt_sparse( b : StateVector;
n : integer;
VAR x : StateVector;
VAR rsq : real);
{*
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.1 1998/08/13 08:40:40 peterg
## Initial revision
##
###############################################################
*}
(* Programs using routine SPARSE must define the type
TYPE
glnarray = ARRAY [1..n] OF real;
in the main routine. They must also provide two routines,
PROCEDURE asub(x: glnarray; VAR y: glnarray; n: integer);
and
PROCEDURE atsub(x: glnarray; VAR z: glnarray; n: integer);
which calculate A*x and (A transpose)*x *)
LABEL 1,99;
CONST
eps=0.000001;
VAR
j,iter,irst: integer;
rp,gg,gam,eps2,dgg,bsq,anum,aden: real;
g,h,xi,xj: glnarray;
BEGIN
eps2 := n*sqr(eps);
irst := 0;
1: irst := irst+1;
mtt_asub(x,xi,n);
rp := 0.0;
bsq := 0.0;
FOR j := 1 TO n DO BEGIN
bsq := bsq+sqr(b[j]);
xi[j] := xi[j]-b[j];
rp := rp+sqr(xi[j])
END;
mtt_atsub(xi,g,n);
FOR j := 1 TO n DO BEGIN
g[j] := -g[j];
h[j] := g[j]
END;
FOR iter := 1 TO 10*n DO BEGIN
mtt_asub(h,xi,n);
anum := 0.0;
aden := 0.0;
FOR j := 1 TO n DO BEGIN
anum := anum+g[j]*h[j];
aden := aden+sqr(xi[j])
END;
IF (aden = 0.0) THEN BEGIN
writeln('pause in routine SPARSE');
writeln('very singular matrix'); readln
END;
anum := anum/aden;
FOR j := 1 TO n DO BEGIN
xi[j] := x[j];
x[j] := x[j]+anum*h[j]
END;
mtt_asub(x,xj,n);
rsq := 0.0;
FOR j := 1 TO n DO BEGIN
xj[j] := xj[j]-b[j];
rsq := rsq+sqr(xj[j])
END;
IF ((rsq = rp) OR (rsq <= bsq*eps2)) THEN GOTO 99;
IF (rsq > rp) THEN BEGIN
FOR j := 1 TO n DO BEGIN
x[j] := xi[j]
END;
IF (irst >= 3) THEN GOTO 99;
GOTO 1
END;
rp := rsq;
mtt_atsub(xj,xi,n);
gg := 0.0;
dgg := 0.0;
FOR j := 1 TO n DO BEGIN
gg := gg+sqr(g[j]);
dgg := dgg+(xi[j]+g[j])*xi[j]
END;
IF (gg = 0.0) THEN GOTO 99;
gam := dgg/gg;
FOR j := 1 TO n DO BEGIN
g[j] := -xi[j];
h[j] := g[j]+gam*h[j]
END
END;
writeln('pause in routine SPARSE');
writeln('too many iterations'); readln;
99: END;
|