program cmtest code: test if cm Fort90 array functions from cm-guide work using cmf include'/usr/include/cm/timer-fort.h' integer, array(2,3):: a,b,cu integer, array(2):: s2,ctr1,ctr2,ctr3,cr1,cr2,b2 integer, array(3):: s3,at,ar1,ar2,br1,br2 integer c(2,2),ct(3,2),cs(3,4),cst(4,3),as(4),d(2,9) logical inmask(64,64) real, array(64,64):: u,du real us(8,8),diffav Caution: Fortran90 array intrinsics need to be listed in INTRINSIC stmt cft77: intrinsic sum,maxval,minval,product,dot_product,matmul,transpose cft77: & ,cshift,eoshift,spread data as/2,3,4,5/ data at/2,3,4/ call cm_timer_clear(0) call cm_timer_start(0) print* print60 60 format(1x/1x,'CM TEST Program for Array Intrinsic Functions') c --------------------CONSTRUCTOR 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. CAUTION: no constructors like "[1,2]" allowed in cf77v6.0 br1 = b(1,:) br2 = b(2,:) print63,'constructor defined b(2,3) =',br1,br2 63 format(1x,a36/(12x,3i3)) c --------------------SUM isum = sum(b) ! => isum = 21; i.e., Front-End scalar. print61,'isum = sum(b) =',isum 61 format(1x,a36/(12x,i3)) isum = sum(b(:,1:3:2)) ! => isum = 14; sole ':' means all values '1:2'. print61,'isum = sum(b(:,1:3:2)) =',isum s2 = sum(b,dim=2) ! declared with the correct array section shape. print62,'s2 = sum(b,2) =',s2 ! => s2 = [9,12], row sums 62 format(1x,a36/(12x,2i3)) s3 = sum(b,dim=1) ! => s3 = [3,7,11]; column sums. print63,'s3 = sum(b,dim=1) =',s3 isum = sum(b,mask=b.gt.3) ! =>isum = 15; i.e., add only elements print61,'isum = sum(b,mask=b.gt.3) =',isum s3 = sum(b,dim=1,mask=b.gt.3) !=> s3 = [0,4,11], conditional col sum print63,'s3 = sum(b,dim=1,mask=b.gt.3) =',s3 s2 = sum(b,dim=2,mask=b.gt.3) ! => s2 = [5,10], conditional row sum print62,'s2 = sum(b,dim=2,mask=b.gt.3) =',s2 c --------------------MAXVAL imax = maxval(b) ! => imax = 6; array maximum value. print61,'imax = maxval(b) =',imax s3 = maxval(b,dim=1) ! => s3 = [2,4,6]; column maximums. print63,'s3 = maxval(b,dim=1) =',s3 s2 = maxval(b,dim=2) ! => s2 = [5,6]; row maximums. print62,'s2 = maxval(b,dim=2) =',s2 c --------------------MINVAL imin = minval(b) ! => imin = 1; array minimum value. print61,'imin = minval(b) =',imin c --------------------PRODUCT s2 = product(b,dim=2) ! => s2 = [15,48]; products of column elements. print62,'s2 = product(b,dim=2) =',s2 c --------------------DOTPRODUCT idot = dotproduct(b(1,:),b(2,:)) ! => idot = 44; dot product of row print61,'idot = dotproduct(b(1,:),b(2,:)) =',idot ! vectors of b. c --------------------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]. cr1=c(1,:) cr2=c(2,:) print62,'c = matmul(b(:,1:2),b(:,2:3)) =',cr1,cr2 c --------------------TRANSPOSE ![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,:) print62,'ct = transpose(b) =',ctr1,ctr2,ctr3 c -------------------- CCCCCCCCCCC FORALL egs? c -------------------- CCCCCCCCCCC WHERE egs? c --------------------CSHIFT ! assume b is again initialized as ! b = 1 3 5 ! 2 4 6 a = cshift(b,shift=1,dim=2) ! => a = 3 5 1 ! 4 6 2 print63,'a = cshift(a,shift=1,dim=2) =',transpose(a) ! i.e., b(i,j+shift) -> a(i,j) for j=1:2, etc.; ! i.e., the result is computed from shifting subscript in specified ! dimension of the source array by the specified shift. a = cshift(b,dim=2,shift=-1) ! => a = 5 1 3 ! 6 2 4 print63,'a = cshift(b,dim=2,shift=-1) =',transpose(a) ! i.e., b(i,j+shift) -> a(i,j) for j=2:3, etc. a = cshift(b,2,[1,2]) ! a = 3 5 1 ! 6 2 4 ! i.e., an array-valued shift, or shift per row. CAUTION: Keywords like "dim=" are optional, but if used, must be used ! for all but array arguments. print63,'a = cshift(b,2,[1,2]) =',transpose(a) c --------------------LAPLACE'S EQUATION ! 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.5.e-3.and.iter.lt.100) iter=iter+1 du = 0 where(inmask) du = 0.25*(cshift(u,1,1)+cshift(u,1,-1)+cshift(u,2,1) & +cshift(u,2,-1)) - 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. cf6er:print66,'u = laplace-shift(u) =',u(1:64:16,1:64:16) us=transpose(u(1:64:9,1:64:9)) print66,'u = laplace-shift(u)= ; iter=',iter,'; av-diff =' & ,diffav,us 66 format(1x,a36,i3,a7,e10.3/(8f7.3)) c --------------------EOSHIFT a = eoshift(b,1,-1) ! a = 0 0 0 note default boundary value is 0. ! 1 3 5 print63,'a = eoshift(b,1,-1) =',transpose(a) a = eoshift(b,dim=2,shift=[-1,0],boundary=[7,8]) ! => a = 7 1 3 ! 2 4 6 print63,'a = eoshift(b,2,[-1,0],[7,8]) =',transpose(a) a = eoshift(b,2,2) ! => a = 5 0 0 ! => 6 0 0 print63,'a = eoshift(b,2,2) =',transpose(a) c --------------------REPLICATE cst = replicate(b,dim=1,ncopies=2) ! contents of cst: ! 1 3 5 ! 2 4 6 ! 1 3 5 ! 2 4 6 print63,'cst=replicate(b,dim=1,ncopies=2) =',transpose(cst) c --------------------------------------------------------------------------- d = replicate(b,DIM=2,NCOPIES=3) ! contents of d: ! 1 3 5 1 3 5 1 3 5 ! 2 4 6 2 4 6 2 4 6 print69,'d = replicate(b,DIM=2,NCOPIES=3) =',transpose(d) 69 format(1x,a36/(12x,9i3)) d(:,2:4) = replicate(b,2,1) ! contents of d(:,2:4) is b print63,'d(:,2:4) = replicate(b,2,1) =',transpose(d(:,2:4)) c --------------------SPREAD as = [2,3,4,5] cs = spread(as,dim=1,ncopies=3) ! contents of cs: ! 2 3 4 5 ! 2 3 4 5 ! 2 3 4 5 print64,'as = [2,3,4,5] = ',as 64 format(1x,a36/(12x,4i3)) print64,'cs = spread(as,dim=1,ncopies=3) =',transpose(cs) c --------------------------------------------------------------------------- at = [2,3,4] cs = spread(at,dim=2,ncopies=4) ! contents of c: ! 2 2 2 2 ! 3 3 3 3 ! 4 4 4 4 print63,'at = [2,3,4] = ',at print64,'cs = spread(at,dim=2,ncopies=4) =',transpose(cs) 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 --------------------------------------------------------------- CCCCCCCCC LAYOUT egs? c --------------------------------------------------------------- CCCCCCCCC ALIGN egs? c --------------------------------------------------------------- call cm_timer_stop(0) print*,'total CM5 time' call cm_timer_print(0) stop end