program f90test code98: compare ranf() and random_number pseudo random number generators code97: update by removing old comments to cmfortran code96: retest=f90test.f redone on borg = convex spp1200/xa-16 integer, parameter :: m = 6 integer, parameter :: n = 4 integer :: i,j integer, dimension(2) :: s2, ctr1, ctr2, ctr3, b2 integer, dimension(3) :: s3 ,at ,ar1 ,ar2 ,br1 ,br2 integer, dimension(4) :: as(4) integer, dimension(2,2) :: c ,bi integer, dimension(2,3) :: b, a integer, dimension(3,2) :: ct integer, dimension(3,4) :: cs integer, dimension(4,3) :: cst logical, dimension(2,3) :: test logical, dimension(64,64) :: inmask real, parameter :: tol = 0.5e-5 integer, parameter :: niter = 5000 real :: diffav real, dimension(8,8) :: us real, dimension(64,64) :: u , du real :: ranf, xran real, dimension(m,n) :: uniranf, uniran real, dimension(n,m) :: truniranf, truniran intrinsic sum,maxval,minval,product & ,dot_product,matmul,transpose & ,cshift,eoshift,spread data b/1,2,3,4,5,6/ !replace constructors initialization data as/2,3,4,5/ data at/2,3,4/ c --------------------Array Constructors: b(1,1:3) = (/1, 3, 5/) ! initialize first row, along dimension 2. b(2,1:3) = (/2, 4, 6/) ! initialize second row, along dimension 2. print*,'Note: constructors like "(/1,2/)" allowed in fc9.5' br1 = b(1,:) br2 = b(2,:) print60,br1,br2 60 format(' b(2,3)'/(3i3)) c --------------------Sum Function sum: isum = sum(b) ! => isum = 21; i.e., Front-End scalar. print61,' isum=sum(b)=',isum 61 format(1x,a36,i4) isum = sum(b(:,1:3:2)) ! => isum = 14; sole ':' means all values '1:2'. print61,' isum = sum("b(:,1:3:2)")=',isum bi=b(:,1:3:2) isum=sum(bi) print61,' isum = sum("b(:,1:3:2)")=',isum print*,'CAUTION: "dim=", etc., markers= NOT allowed in intrinsics' s2 = sum(b,2) ! redeclared with the correct array section shape. print62,' s2 = sum(b,2)=',s2 ! => s2 = (/9,12/), row sums 62 format(1x,a32,2i3) s3 = sum(b,1) ! => s3 = (/3,7,11/); column sums. print63,' s3 = sum(b,1)=',s3 63 format(1x,a32,3i3) print*,'CAUTION: "mask=" marker= STILL not allowed either.' s3 = sum(b,1,b.gt.3) ! => s3 = (/0,4,11/); i.e., conditional col sum print63,' s3 = sum(b,1,"b.gt.3") =',s3 test=b.gt.3 s3 = sum(b,1,test) ! => s3 = (/0,4,11/); i.e., conditional col sum print63,' s3 = sum(b,1,"b.gt.3") =',s3 s2 = sum(b,2,test) ! => s2 = (/5,10/); i.e., conditional row sum print62,' s2 = sum(b,2,b.gt.3) =',s2 cf8er:isum = sum(b,0,test) ! => isum = 18; i.e., add only elements cf8er:print61,' isum = sum(b,0,b.gt.3) =',isum ! that are greater than three. print*,' CAUTION: If "sum(array[dim[,mask]])", CANT use zero (0)' & ,' for [dim] for whole array when there is a mask.' c --------------------Maximum Value Function maxval: imax = maxval(b) ! => imax = 6; array maximum value. print61,' imax = maxval(b)=',imax s3 = maxval(b,1) ! => s3 = (/2,4,6/); column maximums. print63,' s3 = maxval(b,1)=',s3 s2 = maxval(b,2) ! => s2 = (/5,6/); row maximums. print62,' s2 = maxval(b,2)=',s2 c --------------------Minimum Value Function minval: imin = minval(b) ! => imin = 1; array minimum value. print61,' imin = minval(b)=',imin c --------------------Product Function product: s2 = product(b,2) ! => s2 = (/15,48/); products of column elements. print62,' s2 = product(b,2)=',s2 c --------------------Dot Product Function dot_product: idot = dot_product(br1,br2) ! => idot = 44; dot product of row print61,' idot = dot_product(b(1,:),b(2,:))=',idot ! vectors of b. print*,' CAUTION: Array syntax not allowed in actual arguments.' c --------------------Matrix Multiplication Function matmul: ! assuming array b of the previous section. ![Ans] = matmul([Array_1],[Array_2]) ! computes matrix multiplication ! of two rank two matrices. c = matmul(b(:,1:2),b(:,2:3)) ! => c(1,:)=(/15,23/);c(2,:)=(/22,34/). c=transpose(c) print623,'c=matmul(b(:,1:2),b(:,2:3))=',c 623 format(1x,a36/(2i3)) ![Ans] = transpose([Array]) ! transforms an array to its transpose. ct = transpose(b) ! => ct(1,:)=(/1,2/);ct(2,:)=(/3,4/);ct(3,:)=(/5,6/). ctr1 = ct(1,:) ctr2 = ct(2,:) ctr3 = ct(3,:) print623,'ct = transpose(b)=',ctr1,ctr2,ctr3 c --------------------Circular Shift Function cshift: ! assume b is again initialized as ! b = 1 3 5 ! 2 4 6 a = cshift(b,1,2) ! => a = 3 5 1 ! 4 6 2 cshift EG1: ar1 = a(1,:) ar2 = a(2,:) print633,'a = cshift(a,1,2)=',ar1,ar2 633 format(1x,a36/(3i3)) ! i.e., b(i,(j+shift) "mod" n) -> a(i,j) for j=1:2, etc.; ! nonstandard modulus fn: 0 "mod" n = n; 1 "mod" n = 1; ...; n "mod" n = n ! i.e., the result is computed from shifting subscript in specified ! dimension of the source array by the specified shift. a = cshift(b,-1,2) ! => a = 5 1 3 ! 6 2 4 cshift EG2: ar1 = a(1,:) ar2 = a(2,:) print633,'a = cshift(b,-1,2)=',ar1,ar2 ! i.e., b(i,(j+shift) "mod" n) -> a(i,j) for j=2:3, etc. cshift EG3: s2(1) = 1 s2(2) = 2 a = cshift(b,s2,2) ! a = 3 5 1 ! 6 2 4 ! i.e., an array-valued shift, or shift per row. ar1 = a(1,:) ar2 = a(2,:) print633,'a = cshift(b,(/1,2/),2)=',ar1,ar2 cshift Laplace Example: ! Jacobi Iteration for a 5-star discretization of ! 2D Laplace's equation: u = 0 u(1,:)=2 u(64,:)=2 u(:,1)=2 u(:,64)=1 inmask = .FALSE. inmask(2:63,2:63) = .TRUE. diffav = 1 iter=0 do while (diffav.gt.tol.and.iter.lt.niter) iter=iter+1 du = 0 where(inmask) du = 0.25*(cshift(u,1,1)+cshift(u,-1,1)+cshift(u,1,2) & +cshift(u,-1,2)) - u u = u + du end where du = du*du diffav = sqrt(sum(du)/(62*62)) end do ! which is the main program fragment of laplace.fcm. print*,'CAUTION: array sections not allowed in print' us = u(1:64:9,1:64:9) us=transpose(us) print66,'u = laplace-shift(u)= ; iter=',iter,'; av-diff =' & ,diffav,us 66 format(1x,a36,i5,a11,e10.3/(8f8.4)) c --------------------End Off Shift Function eoshift: a = eoshift(b,-1,0,1) ! a = 0 0 0 note default boundary value is 0. ! 1 3 5 ar1 = a(1,:) ar2 = a(2,:) print633,'a = eoshift(b,-1,0,1)=',ar1,ar2 s2=(/-1,0/) b2=(/7,8/) a = eoshift(b,s2,b2,2) ! => a = 7 1 3 ! 2 4 6 ar1 = a(1,:) ar2 = a(2,:) print633,'a = eoshift(b,(/-1,0/),(/7,8/),2)=',ar1,ar2 a = eoshift(b,2,0,2) ! => a = 5 0 0 ! => 6 0 0 ar1 = a(1,:) ar2 = a(2,:) print623,'a = eoshift(b,2,2)=',ar1,ar2 c --------------------Spread Function spread: cs = spread(as,1,3) ! contents of cs: ! 2 3 4 5 ! 2 3 4 5 ! 2 3 4 5 cst = transpose(cs) print64,'as =',as 64 format(1x,a32,4i3) print643,'cs = spread(as,1,3)=',cst 643 format(1x,a36/(4i3)) c -------------------- cs = spread(at,2,4) ! contents of c: ! 2 2 2 2 ! 3 3 3 3 ! 4 4 4 4 cst = transpose(cs) print63,'at =',at print643,'cs = spread(at,2,4)=',cst c --------------------------------------------------------------------------- ! i.e., b=spread(a,d,c) => ! a(n_1,n_2,...,n_(d-1),n_d,...,n_r) -> b(n_1,n_2,...,n_(d-1),c,n_d,...,n_r) ! where r is the rank of source array a and n_i is the size of dimension i; ! noting that a new dimension of size c is added before dimension d. c --------------------------------------------------------------------------- ! Initialize scalar xran with a pseudo random number call random_number(harvest=xran) call random_number(uniran) ! xran and uniran contain uniformly distributed random numbers truniran = transpose(uniran) write(6,65) xran, truniran 65 format(' f90 uniform random_number(): xran =',f14.10/ & ' and f90 subroutine random_number() uniform random array:' & /(4f14.10)) ! standard UNICOS random number generator ranf: do i = 1, m do j = 1, n uniranf(i,j) = ranf() enddo enddo truniranf = transpose(uniranf) write(6,651) truniranf 651 format(' UNICOS function ranf() uniform random array:'/(4f14.10)) stop end