LALA(3f) - [M_matrix] initialize and/or pass commands to matrix
laboratory interpreter
LICENSE(MIT)
subroutine lala(init,cmd)
integer,intent(in),optional :: init
character(len=*),intent(in),optional :: cmd
or
character(len=*),intent(in),optional :: cmd(:)
LALA(3f) is modeled on MATLAB(3f) (MATrix LABoratory), a FORTRAN
package developed by Argonne National Laboratories for in-house use.
It provides comprehensive vector and tensor operations in a package
which may be programmed, either through a macro language or through
execution of script files.
LALA(3f) Functions supported include (but are not by any means limited
to) sin, cos, tan, arcfunctions, upper triangular, lower triangular,
determinants, matrix multiplication, identity, Hilbert matrices,
eigenvalues and eigenvectors, matrix roots and products, inversion
and so on and so forth.
LALA() can be used
+ as a stand-alone utility for working with lala() files and
for basic computations.
+ embedded in a Fortran program, passing variables back and forth
between the calling program and the utility.
+ to read configuration and data files that contain expressions
and conditionally selected values.
+ for interactively inspecting data generated by the calling program.
+ for creating unit tests that allow for further interactive examination.
The HELP command describes using the interpreter.
INIT indicate size of scratch space to allocate and (re)initialize
LALA.
CMD LALA command(s) to perform. May be CHARACTER scalar or vector
INIT and CMD cannot be combined on a single call.
The first call may be an initialization declaring the number of
doubleprecision complex values to allocate for the combined scratch
and variable storage area. This form may be repeated and reinitializes
the utility at each call. A size of zero will deallocate any allocated
storage (after which the routine cannot be called with commands until
reallocated by another call to lala()).
If no parameters are supplied interactive mode is entered.
If a CMD is passed and no previous initialization call was made the
scratch space will be allocated to 200000.
Example 1:
program demo_LALA
use M_matrix, only : lala
write(*,'(a)')'optionally initialize scratch area size'
call LALA(20000)
write(*,'(a)')'do some commands'
call LALA([character(len=80) :: &
& 'semi; ',&
& 'a=magic(4),b=-a ',&
& 'a+b;a;b ',&
& "display('That is all Folks!') "])
write(*,'(a)')'do a single command'
call LALA('who')
write(*,'(a)')'enter interactive mode'
call LALA()
write(*,'(a)')'ending program'
end program demo_LALA
Example 2:
program bigmat
use M_matrix, only : lala
! pass strings to LALA but do not enter interactive mode
call lala(20000) ! initialize silently
call lala( 'a=[1 2 3 4; 5 6 7 8]')
call lala( [character(len=80) :: &
& 'semi;lines(999999) ',&
& '// create a magic square and add 100 to all the values',&
& 'A=magic(4),<X,Y>=shape(A) ',&
& 'B=A+ones(X,Y)*100 ',&
& '// save all current values to a file ',&
& "save('sample.laf') ",&
& '// clear all user values ',&
& 'clear ',&
& '// show variable names, load values from file ',&
& '// and show again to show the variables are restored ',&
& "who;load('sample.laf');who "])
end program bigmat
Example 3: Sample program with custom user function
program custom_user
use M_matrix
implicit none
call set_usersub(lala_user)
call lala()
contains
!-------------------------------------------------------------
subroutine lala_user(a,m,n,s,t) ! sample user routine
! Allows personal Fortran subroutines to be linked into
! LALA. The subroutine should have the heading
!
! subroutine name(a,m,n,s,t)
! integer :: m,n
! doubleprecision a(:),s,t
!
! The LALA statement Y = USER(X,s,t) results in a call to
! the subroutine with a copy of the matrix X stored in the
! argument A, its column and row dimensions in M and N,
! and the scalar parameters S and T stored in S and T.
! If S and T are omitted, they are set to 0.0. After
! the return, A is stored in Y. The dimensions M and
! N may be reset within the subroutine. The statement Y =
! USER(K) results in a call with M = 1, N = 1 and A(1,1) =
! FLOAT(K). After the subroutine has been written, it must
! be compiled and linked to the LALA object code within the
! local programming environment.
!
implicit none
integer :: m,n
doubleprecision :: a(:)
doubleprecision :: s,t
integer :: i, j, k
write(*,*)'MY ROUTINE'
write(*,*)'M=',m
write(*,*)'N=',n
write(*,*)'S=',s
write(*,*)'T=',t
k=0
do i = 1, m
do j = 1, n
k=k+1
write(*,*)i,j,a(k)
enddo
enddo
k=0
if(s.eq.0)s=1
do i = 1, m
do j = 1, n
k=k+1
a(k)=a(k)*s+t
enddo
enddo
end subroutine lala_user
end program custom_user
Example inputs
>:avg:
>for i = 2:2:n, for j = 2:2:n, t = (a(i-1,j-1)+a(i-1,j)+a(i,j-1)+a(i,j))/4; ...
>a(i-1,j-1) = t; a(i,j-1) = t; a(i-1,j) = t; a(i,j) = t;
>:cdiv:
>// ======================================================
>// cdiv
>a=sqrt(random(8))
>ar = real(a); ai = imag(a); br = real(b); bi = imag(b);
>p = bi/br;
>t = (ai - p*ar)/(br + p*bi);
>cr = p*t + ar/br;
>ci = t;
>p2 = br/bi;
>t2 = (ai + p2*ar)/(bi + p2*br);
>ci2 = p2*t2 - ar/bi;
>cr2 = t2;
>s = abs(br) + abs(bi);
>ars = ar/s;
>ais = ai/s;
>brs = br/s;
>bis = bi/s;
>s = brs**2 + bis**2;
>cr3 = (ars*brs + ais*bis)/s;
>ci3 = (ais*brs - ars*bis)/s;
>[cr ci; cr2 ci2; cr3 ci3]
>// ======================================================
>:exp:
>t = 0*x + eye; s = 0*eye(x); n = 1;
>while abs(s+t-s) > 0, s = s+t, t = x*t/n, n = n + 1
>:four:
> n
> pi = 4*atan(1);
> i = sqrt(-1);
> w = exp(2*pi*i/n);
> F = [];
> for k = 1:n, for j = 1:n, F(k,j) = w**((j-1)*(k-1));
> F = F/sqrt(n);
> alpha = r*pi;
> rho = exp(i*alpha);
> S = log(rho*F)/i - alpha*eye;
> serr = norm(imag(S),1);
> S = real(S);
> serr = serr + norm(S-S',1)
> S = (S + S')/2;
> ferr = norm(F-exp(i*S),1)
> :gs:
> for k = 1:n, for j = 1:k-1, d = x(k,:)*x(j,:)'; x(k,:) = x(k,:) - d*x(j,:); ...
> end, s = norm(x(k,:)), x(k,:) = x(k,:)/s;
> :jacobi:
> [n, n] = shape(A);
> X = eye(n);
> anorm = norm(A,'fro');
> cnt = 1;
> while cnt > 0, ...
> cnt = 0; ...
> for p = 1:n-1, ...
> for q = p+1:n, ...
> if anorm + abs(a(p,q)) > anorm, ...
> cnt = cnt + 1; ...
> exec('jacstep'); ...
> end, ...
> end, ...
> end, ...
> display(rat(A)), ...
> end
> :jacstep:
> d = (a(q,q)-a(p,p))*0.5/a(p,q);
> t = 1/(abs(d)+sqrt(d*d+1));
> if d < 0, t = -t; end;
> c = 1/sqrt(1+t*t); s = t*c;
> R = eye(n); r(p,p)=c; r(q,q)=c; r(p,q)=s; r(q,p)=-s;
> X = X*R;
> A = R'*A*R;
> :kron:
> // C = Kronecker product of A and B
> [m, n] = shape(A);
> for i = 1:m, ...
> ci = a(i,1)*B; ...
> for j = 2:n, ci = [ci a(i,j)*B]; end ...
> if i = 1, C = ci; else, C = [C; ci];
> :lanczos:
> [n,n] = shape(A);
> q1 = rand(n,1);
> ort
> alpha = []; beta = [];
> q = q1/norm(q1); r = A*q(:,1);
> for j = 1:n, exec('lanstep',0);
> :lanstep:
> alpha(j) = q(:,j)'*r;
> r = r - alpha(j)*q(:,j);
> if ort <> 0, for k = 1:j-1, r = r - r'*q(:,k)*q(:,k);
> beta(j) = norm(r);
> q(:,j+1) = r/beta(j);
> r = A*q(:,j+1) - beta(j)*q(:,j);
> if j > 1, T = diag(beta(1:j-1),1); T = diag(alpha) + T + T'; eig(T)
> :mgs:
> for k = 1:n, s = norm(x(k,:)), x(k,:) = x(k,:)/s; ...
> for j = k+1:n, d = x(j,:)*x(k,:)'; x(j,:) = x(j,:) - d*x(k,:);
> :net:
> C = [
> 1 2 15 . . .
> 2 1 3 . . .
> 3 2 4 11 . .
> 4 3 5 . . .
> 5 4 6 7 . .
> 6 5 8 . . .
> 7 5 9 30 . .
> 8 6 9 10 11 .
> 9 7 8 30 . .
> 10 8 12 30 31 34
> 11 3 8 12 13 .
> 12 10 11 34 36 .
> 13 11 14 . . .
> 14 13 15 16 38 .
> 15 1 14 . . .
> 16 14 17 20 35 37
> 17 16 18 . . .
> 18 17 19 . . .
> 19 18 20 . . .
> 20 16 19 21 . .
> 21 20 22 . . .
> 22 21 23 . . .
> 23 22 24 35 . .
> 24 23 25 39 . .
> 25 24 . . . .
> 26 27 33 39 . .
> 27 26 32 . . .
> 28 29 32 . . .
> 29 28 30 . . .
> 30 7 9 10 29 .
> 31 10 32 . . .
> 32 27 28 31 34 .
> 33 26 34 . . .
> 34 10 12 32 33 35
> 35 16 23 34 36 .
> 36 12 35 38 . .
> 37 16 38 . . .
> 38 14 36 37 . .
> 39 24 26 . . .
> ];
> [n, m] = shape(C);
> A = 0*ones(n,n);
> for i=1:n, for j=2:m, k=c(i,j); if k>0, a(i,k)=1;
> check = norm(A-A',1), if check > 0, quit
> [X,D] = eig(A+eye);
> D = diag(D); D = D(n:-1:1)
> X = X(:,n:-1:1);
> [x(:,1)/sum(x(:,1)) x(:,2) x(:,3) x(:,19)]
> :pascal:
> //Generate next Pascal matrix
> [k,k] = shape(L);
> k = k + 1;
> L(k,1:k) = [L(k-1,:) 0] + [0 L(k-1,:)];
> :pdq:
> alpha = []; beta = 0; q = []; p = p(:,1)/norm(p(:,1));
> t = A'*p(:,1);
> alpha(1) = norm(t);
> q(:,1) = t/alpha(1);
> X = p(:,1)*(alpha(1)*q(:,1))'
> e(1) = norm(A-X,1)
> for j = 2:r, exec('pdqstep',ip); ...
> X = X + p(:,j)*(alpha(j)*q(:,j)+beta(j)*q(:,j-1))', ...
> e(j) = norm(A-X,1)
> :pdqstep:
> t = A*q(:,j-1) - alpha(j-1)*p(:,j-1);
> if ort>0, for i = 1:j-1, t = t - t'*p(:,i)*p(:,i);
> beta(j) = norm(t);
> p(:,j) = t/beta(j);
> t = A'*p(:,j) - beta(j)*q(:,j-1);
> if ort>0, for i = 1:j-1, t = t - t'*q(:,i)*q(:,i);
> alpha(j) = norm(t);
> q(:,j) = t/alpha(j);
> :pop:
> y = [ 75.995 91.972 105.711 123.203 ...
> 131.669 150.697 179.323 203.212]'
> t = [ 1900:10:1970 ]'
> t = (t - 1940*ones(t))/40; [t y]
> n = 8; A(:,1) = ones(t); for j = 2:n, A(:,j) = t .* A(:,j-1);
> A
> c = A\y
> :qr:
> scale = s(m);
> sm = s(m)/scale; smm1 = s(m-1)/scale; emm1 = e(m-1)/scale;
> sl = s(l)/scale; el = e(l)/scale;
> b = ((smm1 + sm)*(smm1 - sm) + emm1**2)/2;
> c = (sm*emm1)**2;
> shift = sqrt(b**2+c); if b < 0, shift = -shift;
> shift = c/(b + shift)
> f = (sl + sm)*(sl-sm) - shift
> g = sl*el
> for k = l: m-1, exec('qrstep',ip)
> e(m-1) = f
> :qrstep:
> exec('rot');
> if k <> l, e(k-1) = f
> f = cs*s(k) + sn*e(k)
> e(k) = cs*e(k) - sn*s(k)
> g = sn*s(k+1)
> s(k+1) = cs*s(k+1)
> exec('rot');
> s(k) = f
> f = cs*e(k) + sn*s(k+1)
> s(k+1) = -sn*e(k) + cs*s(k+1)
> g = sn*e(k+1)
> e(k+1) = cs*e(k+1)
> :rho:
> //Conductivity example.
> //Parameters ---
> rho //radius of cylindrical inclusion
> n //number of terms in solution
> m //number of boundary points
> //initialize operation counter
> flop = [0 0];
> //initialize variables
> m1 = round(m/3); //number of points on each straight edge
> m2 = m - m1; //number of points with Dirichlet conditions
> pi = 4*atan(1);
> //generate points in Cartesian coordinates
> //right hand edge
> for i = 1:m1, x(i) = 1; y(i) = (1-rho)*(i-1)/(m1-1);
> //top edge
> for i = m2+1:m, x(i) = (1-rho)*(m-i)/(m-m2-1); y(i) = 1;
> //circular edge
> for i = m1+1:m2, t = pi/2*(i-m1)/(m2-m1+1); ...
> x(i) = 1-rho*sin(t); y(i) = 1-rho*cos(t);
> //convert to polar coordinates
> for i = 1:m-1, th(i) = atan(y(i)/x(i)); ...
> r(i) = sqrt(x(i)**2+y(i)**2);
> th(m) = pi/2; r(m) = 1;
> //generate matrix
> //Dirichlet conditions
> for i = 1:m2, for j = 1:n, k = 2*j-1; ...
> a(i,j) = r(i)**k*cos(k*th(i));
> //Neumann conditions
> for i = m2+1:m, for j = 1:n, k = 2*j-1; ...
> a(i,j) = k*r(i)**(k-1)*sin((k-1)*th(i));
> //generate right hand side
> for i = 1:m2, b(i) = 1;
> for i = m2+1:m, b(i) = 0;
> //solve for coefficients
> c = A\b
> //compute effective conductivity
> c(2:2:n) = -c(2:2:n)
> sigma = sum(c)
> //output total operation count
> ops = flop(2)
> :rogers.exec:
> exec('d.boug'); // reads data
> [g,k] = shape(p); // p is matrix of gene frequencies
> wv = ncen/sum(ncen); // ncen contains population sizes
> pbar = wv*p; // weighted average of p
> p = p - ones(g,1)*pbar; // deviations from mean
> p = sqrt(diag(wv)) * p; // weight rows of p by sqrt of pop size
> h = diag(pbar); h = h*(eye-h); // diagonal contains binomial variance: p*(1-p)
> r = p*inv(h)*p'/k; // normalized covariance matrix
> eig(r)'
> :rosser:
> A = [
> 611. 196. -192. 407. -8. -52. -49. 29.
> 196. 899. 113. -192. -71. -43. -8. -44.
> -192. 113. 899. 196. 61. 49. 8. 52.
> 407. -192. 196. 611. 8. 44. 59. -23.
> -8. -71. 61. 8. 411. -599. 208. 208.
> -52. -43. 49. 44. -599. 411. 208. 208.
> -49. -8. 8. 59. 208. 208. 99. -911.
> 29. -44. 52. -23. 208. 208. -911. 99. ];
> :rot:
> // subexec rot(f,g,cs,sn)
> rho = g; if abs(f) > abs(g), rho = f;
> cs = 1.0; sn = 0.0; z = 1.0;
> r = norm([f g]); if rho < 0, r = -r; r
> if r <> 0.0, cs = f/r
> if r <> 0.0, sn = g/r
> if abs(f) > abs(g), z = sn;
> if abs(g) >= abs(f), if cs <> 0, z = 1/cs;
> f = r;
> g = z;
> :rqi:
> rho = (x'*A*x)
> x = (A-rho*eye)\x;
> x = x/norm(x)
> :setup:
> diary('xxx')
> !tail -f xxx > /dev/tty1 &
> !tail -f xxx > /dev/tty2 &
> :sigma:
> RHO = .5 M = 20 N = 10 SIGMA = 1.488934271883534
> RHO = .5 M = 40 N = 20 SIGMA = 1.488920312974229
> RHO = .5 M = 60 N = 30 SIGMA = 1.488920697912116
> :strut.laf:
> // Structure problem, Forsythe, Malcolm and Moler, p. 62
> s = sqrt(2)/2;
> A = [
> -s . . 1 s . . . . . . . . . . . .
> -s . -1 . -s . . . . . . . . . . . .
> . -1 . . . 1 . . . . . . . . . . .
> . . 1 . . . . . . . . . . . . . .
> . . . -1 . . . 1 . . . . . . . . .
> . . . . . . -1 . . . . . . . . . .
> . . . . -s -1 . . s 1 . . . . . . .
> . . . . s . 1 . s . . . . . . . .
> . . . . . . . -1 -s . . 1 s . . . .
> . . . . . . . . -s . -1 . -s . . . .
> . . . . . . . . . -1 . . . 1 . . .
> . . . . . . . . . . 1 . . . . . .
> . . . . . . . . . . . -1 . . . s .
> . . . . . . . . . . . . . . -1 -s .
> . . . . . . . . . . . . -s -1 . . 1
> . . . . . . . . . . . . s . 1 . .
> . . . . . . . . . . . . . . . -s -1];
> b = [
> . . . 10 . . . 15 . . . . . . . 10 .]';
> :test1:
> // -----------------------------------------------------------------
> // start a new log file
> sh rm -fv log.txt
> diary('log.txt')
> // -----------------------------------------------------------------
> titles=['GNP deflator'
> 'GNP '
> 'Unemployment'
> 'Armed Force '
> 'Population '
> 'Year '
> 'Employment '];
> data = ...
> [ 83.0 234.289 235.6 159.0 107.608 1947 60.323
> 88.5 259.426 232.5 145.6 108.632 1948 61.122
> 88.2 258.054 368.2 161.6 109.773 1949 60.171
> 89.5 284.599 335.1 165.0 110.929 1950 61.187
> 96.2 328.975 209.9 309.9 112.075 1951 63.221
> 98.1 346.999 193.2 359.4 113.270 1952 63.639
> 99.0 365.385 187.0 354.7 115.094 1953 64.989
> 100.0 363.112 357.8 335.0 116.219 1954 63.761
> 101.2 397.469 290.4 304.8 117.388 1955 66.019
> 104.6 419.180 282.2 285.7 118.734 1956 67.857
> 108.4 442.769 293.6 279.8 120.445 1957 68.169
> 110.8 444.546 468.1 263.7 121.950 1958 66.513
> 112.6 482.704 381.3 255.2 123.366 1959 68.655
> 114.2 502.601 393.1 251.4 125.368 1960 69.564
> 115.7 518.173 480.6 257.2 127.852 1961 69.331
> 116.9 554.894 400.7 282.7 130.081 1962 70.551];
> short
> X = data;
> [n,p] = shape(X)
> mu = ones(1,n)*X/n
> X = X - ones(n,1)*mu; X = X/diag(sqrt(diag(X'*X)))
> corr = X'*X
> y = data(:,p); X = [ones(y) data(:,1:p-1)];
> long e
> beta = X\y
> expected = [ ...
> -3.482258634594421D+03
> 1.506187227124484D-02
> -3.581917929257409D-02
> -2.020229803816908D-02
> -1.033226867173703D-02
> -5.110410565317738D-02
> 1.829151464612817D+00
> ]
> display('EXPE and BETA should be the same')
> :tryall:
> diary('log.txt')
> a=magic(8)
> n=3
> exec('avg')
> b=random(8,8)
> exec('cdiv')
> exec('exp')
> exec('four')
> exec('gs')
> exec('jacobi')
> // jacstep
> exec('kron')
> exec('lanczos')
> // lanstep
> exec('longley')
> exec('mgs')
> exec('net')
> exec('pascal')
> exec('pdq')
> // pdqstep
> exec('pop')
> exec('qr')
> // qrstep
> exec('rho')
> exec('rosser')
> // rot
> exec('rqi')
> exec('setup')
> exec('sigma')
> exec('strut.laf')
> exec('w5')
> exec('rogers.exec
> exec('rogers.load
> :w5:
> w5 = [
> 1. 1. 0. 0. 0.
> -10. 1. 1. 0. 0.
> 40. 0. 1. 1. 0.
> 205. 0. 0. 1. 1.
> 024. 0. 0. 0. -4.
> ]
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | init | |||
logical, | intent(in), | optional | :: | echo |
Type | Visibility | Attributes | Name | Initial | |||
---|---|---|---|---|---|---|---|
integer, | public, | parameter | :: | EPS(GG_MAX_NAME_LENGTH) | = | [iachar(['e', 'p', 's', ' ', ' ']), GG_PAD(6:)] | |
integer, | public, | parameter | :: | EYE(GG_MAX_NAME_LENGTH) | = | [iachar(['e', 'y', 'e', ' ', ' ']), GG_PAD(6:)] | |
integer, | public, | parameter | :: | FLOPS(GG_MAX_NAME_LENGTH) | = | [iachar(['f', 'l', 'o', 'p', 's']), GG_PAD(6:)] | |
integer, | public, | parameter | :: | RAND(GG_MAX_NAME_LENGTH) | = | [iachar(['r', 'a', 'n', 'd', ' ']), GG_PAD(6:)] | |
doubleprecision, | public | :: | s | ||||
doubleprecision, | public | :: | t |
subroutine LALA_init(init,echo)
! ident_1="@(#) M_matrix lala(3f) initialize and/or pass commands to matrix laboratory interpreter"
integer,intent(in) :: init
logical,intent(in),optional :: echo
doubleprecision :: s,t
integer,parameter :: EPS(GG_MAX_NAME_LENGTH)= [iachar(['e','p','s',' ',' ']),GG_PAD(6:)]
integer,parameter :: FLOPS(GG_MAX_NAME_LENGTH)= [iachar(['f','l','o','p','s']),GG_PAD(6:)]
integer,parameter :: EYE(GG_MAX_NAME_LENGTH)= [iachar(['e','y','e',' ',' ']),GG_PAD(6:)]
integer,parameter :: RAND(GG_MAX_NAME_LENGTH)= [iachar(['r','a','n','d',' ']),GG_PAD(6:)]
if(present(echo)) G_ECHO=echo
G_PROMPT=.true.
G_ERR=0
if(allocated(G_PSEUDO_FILE))deallocate(G_PSEUDO_FILE)
allocate(G_PSEUDO_FILE(0))
G_LIN=blank
G_VAR_IDS=blank
GM_BIGMEM=INIT
if(GM_BIGMEM.lt.0)GM_BIGMEM=200000
if(allocated(GM_REALS) )deallocate(GM_REALS)
if(allocated(GM_IMAGS) )deallocate(GM_IMAGS)
allocate(GM_REALS(GM_BIGMEM),GM_IMAGS(GM_BIGMEM)) ! set to size of GM_BIGMEM
G_INPUT_LUN = STDIN ! unit number for terminal input
call mat_files(G_INPUT_LUN,G_BUF)
G_RIO = G_INPUT_LUN ! current file to read commands from
call mat_files(STDOUT,G_BUF)
call mat_help_text() ! initialize help text
G_CURRENT_RANDOM_SEED = 0 ! random number seed
G_CURRENT_RANDOM_TYPE = 0 ! set the type of random numbers to compute
G_LINECOUNT(2) = 23 ! initial line limit for paging output
G_TOP_OF_SAVED = GG_MAX_NUMBER_OF_NAMES-3 ! move up to allow room for the built-in values eps, flops, eye, rand
call mat_wset(5,0.0D0,0.0d0,GM_REALS(GM_BIGMEM-4),GM_IMAGS(GM_BIGMEM-4),1)
call update('eps',1,1,GM_BIGMEM-4)
!=============================================================
call mat_copyid(G_VAR_IDS(1,GG_MAX_NUMBER_OF_NAMES-3),EPS)
G_VAR_DATALOC(GG_MAX_NUMBER_OF_NAMES-3) = GM_BIGMEM-4
G_VAR_ROWS(GG_MAX_NUMBER_OF_NAMES-3) = 1
G_VAR_COLS(GG_MAX_NUMBER_OF_NAMES-3) = 1
!=============================================================
! interesting way to calculate the epsilon value of a machine
s = 1.0d0
SET_ST: do
s = s/2.0D0
t = 1.0d0 + s
if (t .LE. 1.0d0) exit
enddo SET_ST
GM_REALS(GM_BIGMEM-4) = 2.0d0*s
call update('flops',1,2,GM_BIGMEM-3)
!=============================================================
call mat_copyid(G_VAR_IDS(1,GG_MAX_NUMBER_OF_NAMES-2),flops)
G_VAR_DATALOC(GG_MAX_NUMBER_OF_NAMES-2) = GM_BIGMEM-3
G_VAR_ROWS(GG_MAX_NUMBER_OF_NAMES-2) = 1
G_VAR_COLS(GG_MAX_NUMBER_OF_NAMES-2) = 2
!=============================================================
call update('eye',-1,-1,GM_BIGMEM-1)
!=============================================================
call mat_copyid(G_VAR_IDS(1,GG_MAX_NUMBER_OF_NAMES-1), eye)
G_VAR_DATALOC(GG_MAX_NUMBER_OF_NAMES-1) = GM_BIGMEM-1
G_VAR_ROWS(GG_MAX_NUMBER_OF_NAMES-1) = -1
G_VAR_COLS(GG_MAX_NUMBER_OF_NAMES-1) = -1
!=============================================================
GM_REALS(GM_BIGMEM-1) = 1.0D0
call update('rand',1,1,GM_BIGMEM)
!=============================================================
call mat_copyid(G_VAR_IDS(1,GG_MAX_NUMBER_OF_NAMES), rand)
G_VAR_DATALOC(GG_MAX_NUMBER_OF_NAMES) = GM_BIGMEM
G_VAR_ROWS(GG_MAX_NUMBER_OF_NAMES) = 1
G_VAR_COLS(GG_MAX_NUMBER_OF_NAMES) = 1
!=============================================================
G_FMT = 1
G_FLOP_COUNTER(1) = 0
G_FLOP_COUNTER(2) = 0
G_DEBUG_LEVEL = 0
G_PTZ = 0
G_PT = G_PTZ
G_FORTRAN_TEXT=help_intrinsics('manual',m_help=.true.) ! load Fortran documentation
end subroutine LALA_init