C Library Functions - lala (3)
NAME
LALA(3f) - [M_matrix] initialize and/or pass commands to matrix laboratory interpreter LICENSE(MIT)
CONTENTS
Synopsis
Description
Options
Example
SYNOPSIS
subroutine lala(init,cmd)
integer,intent(in),optional :: init character(len=*),intent(in),optional :: cmd or character(len=*),intent(in),optional :: cmd(:)
DESCRIPTION
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
o as a stand-alone utility for working with lala() files and for basic computations. o embedded in a Fortran program, passing variables back and forth between the calling program and the utility. o to read configuration and data files that contain expressions and conditionally selected values. o for interactively inspecting data generated by the calling program. o for creating unit tests that allow for further interactive examination. The HELP command describes using the interpreter.
OPTIONS
INIT and CMD cannot be combined on a single call.
INIT indicate size of scratch space to allocate and (re)initialize LALA. CMD LALA command(s) to perform. May be CHARACTER scalar or vector 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
Example 1:
program demo_LALA use M_matrix, only : lalaExample 2: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
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 bigmatExample 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_userExample 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. > ]
| Nemo Release 3.1 | lala (3) | July 22, 2023 |
