c=========================================================== c fdemo1: Program which demonstrates many of the c essential features of Fortran 77. Some 'safe' language c extensions are used; all extensions are valid c Fortran 90. c=========================================================== c=========================================================== c Source code formatting rules: c c Columns Use c c 1-5 numeric statement label c 6 continuation character: '&' recommended c 7-72 statement c c BE EXTREMELY CAREFUL NOT TO TYPE BEYOND COLUMN 72! c=========================================================== C COMMENT LINES: Use 'c' 'C' or '*' IN FIRST COLUMN *=========================================================== c----------------------------------------------------------- c The 'program' statement names a Fortran main routine. c Optional, but recommended and note that there can c only be one 'program' (main routine) per executable. c----------------------------------------------------------- program fdemo1 c=========================================================== c BEGINNING OF DECLARATION STATEMENTS c c Declarations (or specification statements) must c ALWAYS appear before ANY executable statements. c=========================================================== c----------------------------------------------------------- c The 'implicit none' statement is an extension which c forces us to explicitly declare all variables and c functions (apart from Fortran built in functions). c HIGHLY RECOMMENDED. c----------------------------------------------------------- implicit none c----------------------------------------------------------- c PARAMETERS c----------------------------------------------------------- c The parameter declaration effectively assigns a c CONSTANT value to a name. Note that each c parameter statement must be accompanied by an c appropriate declaration of the type of the c parameter. Also note that, except in strings, c blanks (spaces) are ignored in Fortran---you can c use this fact to make code more readable. c----------------------------------------------------------- integer zero parameter ( zero = 0 ) c----------------------------------------------------------- c Always specify floating point constants using c scientific notation. Use 'd' (instead of 'e') for c real*8 constants. c----------------------------------------------------------- real*8 pi parameter ( pi = 3.141 5926 5358 9793 d0 ) real*8 tiny parameter ( tiny = 1.0 d-50 ) c----------------------------------------------------------- c VARIABLES c----------------------------------------------------------- c The main data types we will be using are c c integer, real*8, logical, c character*1, character*2, ... etc., character*(*) c c but note that Fortran has support for complex c arithmetic. Note that complex*16 means real*8 c values are used for both the real and imaginary c parts of the variable. c----------------------------------------------------------- c----------------------------------------------------------- c (a) SCALARS c----------------------------------------------------------- real*8 a, b, c real*8 res1, res2, res3, res4 integer i, j, k, n integer ires1, ires2, ires3, ires4 logical switch logical lres1, lres2, lres3 complex*16 ca, cb c----------------------------------------------------------- c (b) ARRAYS c----------------------------------------------------------- integer n1, n2, n3 parameter ( n1 = 4, n2 = 3, n3 = 2) c----------------------------------------------------------- c (b.1) 1-D ARRAYS: Note, in a main program, all c dimension bounds must be integer parameters or c integer constants. c----------------------------------------------------------- real*8 r1a(n1), r1b(n2) integer i1i(n1) c----------------------------------------------------------- c (b.2) 2-D ARRAYS: c----------------------------------------------------------- real*8 r2a(n1,n2) c----------------------------------------------------------- c (b.3) 3-D ARRAYS: c----------------------------------------------------------- real*8 r3a(n1,n2,n3) c=========================================================== c END OF DECLARATION STATEMENTS c=========================================================== c=========================================================== c BEGINNING OF EXECUTABLE STATEMENTS c=========================================================== c*********************************************************** c Assignment statements and simple arithmetic c expressions c*********************************************************** c----------------------------------------------------------- c Assignment to scalar variables ... again, note c the use of scientific notation (d0) to specify c a real*8 constant. c c The only valid logical constants are .true. and c .false. (don't forget to include the .'s) c----------------------------------------------------------- a = 0.025d0 b = -1.234d-16 c = 1.0d0 i = 3000 switch = .true. c----------------------------------------------------------- c Note the use of the continuation character in c column 6 to continue a statement on a second line. c----------------------------------------------------------- write(*,*) 'a = ', a,' b = ', b write(*,*) ' c = ', c, ' i = ', i, & ' switch = ', switch call prompt('Through scalar assignment') c----------------------------------------------------------- c Arithmetic expressions. Fortran has standard c operator precedences except that the exponentiation c operator '**' associates RIGHT to LEFT: e.g. c c i ** j ** k is equivalent to i ** (j ** k) c c Parentheses force evaluation of subexpressions. c----------------------------------------------------------- a = 2.0d0 b = 3.0d0 c = 3.0d0 res1 = a + b res2 = a**2 + b**2 res3 = (a**2 + b**2)**(0.5d0) write(*,*) 'res1 = ', res1,' res2 = ', res2 write(*,*) ' res3 = ', res3 call prompt('Through real*8 arithmetic expressions') c----------------------------------------------------------- c Notice the integer truncation which occurs when c dividing the integer 2 by the integer 3. c----------------------------------------------------------- i = 2 j = 3 k = 2 ires1 = 2 + 3 ires2 = 2 / 3 ires3 = i ** j ** k ires4 = (i ** j) ** k write(*,*) 'ires1 = ', ires1, ' ires2 = ', ires2 write(*,*) 'ires3 = ', ires3, ' ires4 = ', ires4 call prompt('Through integer arithmetic expressions') c----------------------------------------------------------- c "Mixed-mode" computations c----------------------------------------------------------- c----------------------------------------------------------- c i + j is computed using integer arithmetic and c the result is converted to a real*8 value before being c assigned to res2. c----------------------------------------------------------- res1 = i + j c----------------------------------------------------------- c 3 / 4 is evaluated using integer arithmetic (yielding c 0) and then the value is converted to real*8. c----------------------------------------------------------- res2 = 3 / 4 c----------------------------------------------------------- c The appearance of a double precision constant c forces the division to be computed using real*8 c arithmetic c----------------------------------------------------------- res3 = 3 / 4.0d0 write(*,*) 'res1 = ', res1, ' res2 = ', res2 write(*,*) ' res3 = ', res3 call prompt('Through mixed-mode arithmetic') c*********************************************************** c CONTROL STATEMENTS c*********************************************************** c*********************************************************** c DO LOOPS c c Note that 'end do' is not Fortran 77, but a safe c extension (it is legal Fortran 90). c*********************************************************** do i = 1 , 3 write(*,*) 'Loop 1: i = ', i end do call prompt('Through loop 1') c----------------------------------------------------------- c The same do loop with the optional loop increment c specified explicitly c----------------------------------------------------------- do i = 1 , 3 , 1 write(*,*) 'Loop 2: i = ', i end do call prompt('Through loop 2') c----------------------------------------------------------- c Another do-loop with a non-default loop increment ... c----------------------------------------------------------- do i = 1 , 7 , 2 write(*,*) 'Loop 3: i = ', i end do call prompt('Through loop 3') c----------------------------------------------------------- c ... and one with a negative increment c----------------------------------------------------------- do i = 3 , 1 , -1 write(*,*) 'Loop 4: i = ', i end do call prompt('Through loop 4') c----------------------------------------------------------- c Nested do-loops. c----------------------------------------------------------- do i = 1 , 3 do j = 1 , 2 write(*,*) 'Loop 5: i, j = ', i, j end do end do call prompt('Through loop 5') c----------------------------------------------------------- c Any of the do-loop parameters can be variables, c expressions or parameters: safest to ALWAYS use c integer values. c----------------------------------------------------------- n = 6 do i = 2 , n , n / 3 write(*,*) 'Loop 6: i = ', i end do call prompt('Through loop 6') c*********************************************************** c LOGICAL EXPRESSIONS c c Note that the Fortran comparison and logical c operators all have the form: .operator. c c Comparison: .eq. .ne. .gt. .lt. c .ge. .le. c Logical: .not. (unary) c .and. .or. c*********************************************************** a = 25.0d0 b = 12.0d0 lres1 = a .gt. b lres2 = (a .lt. b) .or. (b .ge. 0.0d0) lres3 = a .eq. b write(*,*) 'lres1 = ', lres1, ' lres2 = ', lres2, & ' lres3 = ', lres3 call prompt('Through basic conditionals') c*********************************************************** c IF-THEN-ELSE STATEMENTS. c*********************************************************** if( a .gt. b ) then write(*,*) a, ' > ', b end if call prompt('Through if 1') if( b .gt. a ) then write(*,*) b, ' > ', a else write(*,*) a, ' > ', b end if call prompt('Through if 2') c----------------------------------------------------------- c Nested IF statement. c----------------------------------------------------------- if( a .gt. b ) then if( a .gt. 2 * b ) then write(*,*) a, ' > ', 2 * b else write(*,*) a, ' <= ', 2 * b end if else write(*,*) a, ' <= ', b end if call prompt('Through nested if') c----------------------------------------------------------- c IF ... ELSE IF .. IF construct can be used in lieu c of 'CASE' statement. c----------------------------------------------------------- do i = 1 , 4 if( i .eq. 1 ) then write(*,*) 'Case 1' else if( i .eq. 2 ) then write(*,*) 'Case 2' else if( i .eq. 3 ) then write(*,*) 'Case 3' else write(*,*) 'Default case' end if end do call prompt('Through case via if') c*********************************************************** c WHILE LOOPS c c The do while( ... ) ... end do construct is valid c Fortran 90, and a safe Fortran 77 extension. c*********************************************************** a = 0.1d0 b = 0.0d0 do while ( b .le. 1.0d0 ) write(*,*) 'Do while loop: b = ', b b = b + a end do call prompt('Through while loop') c*********************************************************** c USING BUILT-IN (INTRINSIC) FUNCTIONS c*********************************************************** res1 = sin(0.3d0 * Pi) res2 = cos(0.3d0 * Pi) res3 = res1**2 + res2**2 res4 = sqrt(res3) write(*,*) 'res1 = ', res1, ' res2 = ', res2 write(*,*) 'res3 = ', res3, ' res4 = ', res4 call prompt('Through built-in fcn 1') c----------------------------------------------------------- c atan, acos, asin, etc. return arctangent, arcosine, c arcsine etc. in RADIANS c----------------------------------------------------------- res1 = atan(1.0d0) write(*,*) 'res1 = ', res1 call prompt('Through built-in fcn 2') c----------------------------------------------------------- c min and max will return the minimum and maximum c respectively of an arbitrary number of arguments c of any UNIQUE data type. Do NOT mix types in c a single statement as in c c write(*,*) min(1,2.0d0) c----------------------------------------------------------- write(*,*) 'min(3.0d0,2.0d0) = ', min(3.0d0,2.0d0) write(*,*) 'min(1,-3,5,0) = ', min(1,-3,5,0) call prompt('Through built-in fcn 3') c----------------------------------------------------------- c mod is particularly useful for calculating when one c integer divides another evenly c----------------------------------------------------------- do i = 0 , 1000 if( mod(i,100) .eq. 0 ) then write(*,*) 'i = ', i end if end do call prompt('Through built-in fcn 4') c----------------------------------------------------------- c Stop program execution c----------------------------------------------------------- call prompt('Through fdemo1') stop c=========================================================== c END OF EXECUTABLE STATEMENTS c=========================================================== c----------------------------------------------------------- c End of program unit (fdemo1) c----------------------------------------------------------- end c=========================================================== c Prints a message on stdout and then waits for input c from stdin. c=========================================================== c----------------------------------------------------------- c A new program unit (prompt). c----------------------------------------------------------- subroutine prompt(pstring) implicit none character*(*) pstring integer rc character*1 resp write(*,*) pstring write(*,*) 'Enter any non-blank character & '// & 'enter to continue' read(*,*,iostat=rc,end=900) resp c----------------------------------------------------------- c Return to calling program. c----------------------------------------------------------- return 900 continue c----------------------------------------------------------- c Stop program execution. This section of code is c the "end-of-file" handler for standard input c (via the end=900 clause of the read statement). c In this case, it is perfectly good style simply c to quit. c----------------------------------------------------------- stop c----------------------------------------------------------- c End of program unit (prompt). c----------------------------------------------------------- end c=========================================================== c fdemo2: Program which demonstrates basic usage c of character variables in Fortran 77. c=========================================================== program fdemo2 implicit none c----------------------------------------------------------- c See below for definition of integer function c 'indlnb'. Note that this and other useful routines c are available in the 'p410f' library. c----------------------------------------------------------- integer indlnb c----------------------------------------------------------- c Define some character variables of various lengths c c Note that c c character*1 foo c c and c c character foo c c are synonymous, i.e. if an explicit length c specification is not given, the variable will c be a single character long. c----------------------------------------------------------- character*1 c1 character*2 c2 character*4 c4 character*26 lcalph character cc1*1, cc2*2, cc4*4 character*60 buffer c----------------------------------------------------------- c Assignment of constant strings to char. variables. c If length of character expression being assigned c is less than length of character variable, variable c is 'right-padded' with blanks. c----------------------------------------------------------- c1 = 'a' c2 = 'bc' c4 = 'defg' lcalph = 'abcdefghijklmnopqrstuvwxyz' write(*,*) 'c1 = ', c1 write(*,*) 'c2 = ', c2 write(*,*) 'c4 = ', c4 write(*,*) 'lcalph = ', lcalph call prompt('Through constant assignment') c----------------------------------------------------------- c // is the string concatentation operator c----------------------------------------------------------- write(*,*) 'c1 // c2 // c4 = ', c1 // c2 // c4 call prompt('Through concatenation') c----------------------------------------------------------- c The integer intrinsic (built-in) function 'len' c returns the length of its string argument c----------------------------------------------------------- write(*,*) 'len(c1) = ', len(c1) write(*,*) 'len(buffer) = ', len(buffer) call prompt('Through string length') c----------------------------------------------------------- c Substring extraction c----------------------------------------------------------- write(*,*) 'lcalph(1:13) = ', lcalph(1:13) write(*,*) 'lcalph(18:18) = ', lcalph(18:18) call prompt('Through substring extraction') c----------------------------------------------------------- c Substring assignment c----------------------------------------------------------- c4(4:4) = 'Z' write(*,*) 'c4 = ', c4 call prompt('Through substring assignment') c----------------------------------------------------------- c Use of 'indlnb' c----------------------------------------------------------- buffer = 'somefilename' write(*,*) '<' // buffer // '>' write(*,*) '<' // buffer(1:indlnb(buffer)) // '>' buffer = 'Some multi-word message' write(*,*) '<' // buffer // '>' write(*,*) '<' // buffer(1:indlnb(buffer)) // '>' buffer = ' ' write(*,*) 'indlnb(buffer) = ', indlnb(buffer) call prompt('Through indlnb usage') call prompt('Through fdemo2') stop end c----------------------------------------------------------- c Prints a message on stdout and then waits for input c from stdin. c----------------------------------------------------------- subroutine prompt(pstring) implicit none character*(*) pstring integer rc character*1 resp write(*,*) pstring write(*,*) 'Enter any non-blank character & '// & 'enter to continue' read(*,*,iostat=rc,end=900) resp return 900 continue stop end c----------------------------------------------------------- c Returns index of last non-blank character in 's', c or 0 if the string is completely blank. c----------------------------------------------------------- integer function indlnb(s) character*(*) s do indlnb = len(s) , 1 , -1 if( s(indlnb:indlnb) .ne. ' ' ) return end do indlnb = 0 return end write(*,*) 'Hello World!' stop end program greeting call sayhello() stop end c=========================================================== c Computes and reports estimate of machine epsilon. c c Recall: machine epsilon is smallest positive 'eps' c such that c c (1.0d0 + eps ) .ne. (1.0d0) c c Program accepts optional argument which specifies c division factor: values close to 1.0 will result c in more accurate estimate of machine epsilon. c=========================================================== program meps implicit none c----------------------------------------------------------- c Note use of 'r8arg', available in 'libp410f.a' which c works exactly like 'i4arg' except that it returns c a real*8 value parsed from the specified command-line c argument c----------------------------------------------------------- real*8 r8arg real*8 default_fac parameter ( default_fac = 2.0d0 ) real*8 eps, neweps, fac fac = r8arg(1,default_fac) write(0,*) 'meps: using division factor: ', fac eps = 1.0d0 neweps = 1.0d0 do while( 1.0d0 .ne. (1.0d0 + neweps) ) eps = neweps neweps = neweps / fac end do write(*,*) eps stop end c=========================================================== c mysum: reads numbers one per line from stdin c and writes sum on stdout. Ignores invalid inputs c but counts number encountered and reports on stderr. c=========================================================== program mysum implicit none c----------------------------------------------------------- c vi: Current number read from stdin c sum: Current sum of numbers read c rc: For storing return status from READ c nbad: Count of number of bad inputs c----------------------------------------------------------- real*8 vi, sum integer rc, nbad c----------------------------------------------------------- c Initialize ... c----------------------------------------------------------- nbad = 0 sum = 0.0d0 c----------------------------------------------------------- c The following construct is roughly equivalent to c a while loop, execution keeps returning to the c top of the loop until end of file is detected on c stdin. c----------------------------------------------------------- 100 continue read(*,*,iostat=rc,end=200) vi if( rc .eq. 0 ) then c----------------------------------------------------------- c Read a bona fide real*8 value, update sum. c----------------------------------------------------------- sum = sum + vi else c----------------------------------------------------------- c Input was invalid. c----------------------------------------------------------- nbad = nbad + 1 end if go to 100 200 continue c----------------------------------------------------------- c Write sum on standard output. c----------------------------------------------------------- write(*,*) sum c----------------------------------------------------------- c Report # of invalid inputs only if there were some. c----------------------------------------------------------- if( nbad .gt. 0 ) then c----------------------------------------------------------- c Unit 0 is stderr (standard error) on most Unix c systems: if you redirect stdin using '>' and this c message is tripped, it will still appear on the c terminal. c----------------------------------------------------------- write(0,*) nbad, ' invalid inputs' end if stop end c=========================================================== c go-to-less version of mysum. No more or less correct c than go-to-full version, of course, but some would assert c that this version is "cleaner"/"safer" in some sense, c or at the very least, better stylistically. c=========================================================== program mysum implicit none real*8 vi, sum integer rc, rceof, nbad nbad = 0 sum = 0.0d0 c----------------------------------------------------------- c 'rc' will be set to a NEGATIVE value when end of file c is encountered. Initialize 'rc' to 0 to ensure that c we get into the 'do while' loop. c----------------------------------------------------------- rc = 0 do while ( rc .ge. 0 ) read(*,*,iostat=rc) vi if( rc .eq. 0 ) then sum = sum + vi c----------------------------------------------------------- c New logic to ensure that EOF doesn't get counted c as a bad input. c----------------------------------------------------------- else if ( rc .gt. 0 ) then nbad = nbad + 1 end if end do write(*,*) sum if( nbad .gt. 0 ) then write(0,*) nbad, ' invalid inputs' end if stop end c=========================================================== c Less-commented (i.e. more reasonable level of c comments) version of mysum. c=========================================================== c mysum_s: reads numbers one per line from stdin c and writes sum on stdout. Ignores invalid inputs c but counts number encountered and reports on stderr. c=========================================================== program mysum implicit none real*8 vi, sum integer rc, nbad nbad = 0 sum = 0.0d0 100 continue read(*,*,iostat=rc,end=200) vi if( rc .eq. 0 ) then sum = sum + vi else nbad = nbad + 1 end if go to 100 200 continue write(*,*) sum if( nbad .gt. 0 ) then write(0,*) nbad, ' invalid inputs' end if stop end subroutine sayhello() write(*,*) 'Hello World!' return end c=========================================================== c Program to test understanding of Fortran 77 array c storage scheme. c c EXERCISE: Predict program output. c=========================================================== program t2d implicit none integer n, m parameter ( n = 3, m = 4 ) real*8 a(n,m), aij integer i, j do i = 1 , n do j = 1 , m a(i,j) = 100.0d0 * i + m end do end do write(0,*) 'main: a(2,2) = ', a(2,2) call sub(a,2,2) stop end subroutine sub(a,n,m) implicit none integer n, m real*8 a(n,m) write(0,*) 'sub: a(2,2) = ', a(2,2) return end c=========================================================== c Demonstration main program, subroutines and functions c to illustrate RECOMMENDED use of common blocks c using 'include' statement. Safe Fortran 77 c extension. c=========================================================== program tcommon1 implicit none c----------------------------------------------------------- c By convention, I use the extension '.inc' for c Fortran source files which are to be included. c----------------------------------------------------------- include 'coma.inc' string = 'foo' v(1) = 1.0d0 v(2) = 2.0d0 v(3) = 3.0d0 x = 10.0d0 y = 20.0d0 z = 30.0d0 i = 314 call subcom() stop end c=========================================================== c This subroutine dumps information passed to it in c a common block. c=========================================================== subroutine subcom() include 'coma.inc' write(0,*) 'In subcom:' write(0,*) 'string = ', string write(0,*) 'v = ', v write(0,*) 'x = ', x, ' y = ', y, ' z = ', z write(0,*) 'i = ', i return end c=========================================================== c Demonstration main program and subroutine c to illustrate use of COMMON blocks for creating c 'global' storage. Common blocks should always c be labelled (named) and should be used sparingly. c=========================================================== program tcommon implicit none c----------------------------------------------------------- c Declare variables to be placed in common block c----------------------------------------------------------- character*16 string real*8 v(3), & x, y, z integer i c----------------------------------------------------------- c Variables are stored in a common block in the c order in which they are specified in the 'common' c statement. ALWAYS order variables from longest to c shortest to avoid "alignment problems". Don't c try to put a variable in more than one common block c and note that entire arrays (such as 'v') are placed c in the common block by simply specifying the name c of the array. Finally, note that variables in a c common block CAN NOT be initialized with a 'data' c statement. c----------------------------------------------------------- common / coma / & string, & v, & x, y, z, & i string = 'foo' v(1) = 1.0d0 v(2) = 2.0d0 v(3) = 3.0d0 x = 10.0d0 y = 20.0d0 z = 30.0d0 i = 314 call subcom() stop end c=========================================================== c This subroutine dumps information passed to it in c a common block. c=========================================================== subroutine subcom() c----------------------------------------------------------- c Overall layout of common block should be identical c in all program units which use the common block. c----------------------------------------------------------- character*16 string real*8 v(3), & x, y, z integer i common / coma / & string, & v, & x, y, z, & i write(0,*) 'In subcom:' write(0,*) 'string = ', string write(0,*) 'v = ', v write(0,*) 'x = ', x, ' y = ', y, ' z = ', z write(0,*) 'i = ', i return end c------------------------------------------------- c Simple illustration of use of libraries c c dmach is a function defined in the LINPACK c library c c /usr/local/PGI/lib/liblinpack.a c c It computes and returns information c concerning the invoking machine's floating c point model. c------------------------------------------------- program tdmach implicit none real*8 dmach write(0,*) dmach(1) stop end c=========================================================== c Test program for subroutine 'dmfrom', 'dmto' and c 'dmmmult' (see 'dmroutines.f') c c Program expects one argument, the name of a file which c contains a real*8 square matrix written as descibed c in the documentation for 'dmfrom' in 'dmroutines.f' c Use '-' to read from stdin. Program then computes c square of matrix and outputs result to stdout. c=========================================================== program tdm implicit none integer iargc character*256 fname c----------------------------------------------------------- c Maximum size for input and output arrays (matrices). c----------------------------------------------------------- integer maxsize parameter ( maxsize = 100 000 ) real*8 a(maxsize), asq(maxsize) integer d1a, d2a if( iargc() .ne. 1 ) go to 900 call getarg(1,fname) c----------------------------------------------------------- c Read matrix ... c----------------------------------------------------------- call dmfrom(fname,a,d1a,d2a,maxsize) if( d1a .gt. 0 .and. d2a .gt. 0 ) then if( d1a .eq. d2a ) then c----------------------------------------------------------- c Compute square ... c----------------------------------------------------------- call dmmmult(a,a,asq,d1a,d1a) c----------------------------------------------------------- c ... and output. c----------------------------------------------------------- call dmto('-',asq,d1a,d1a) else write(0,*) 'tdm: Input array not square' end if else write(0,*) 'tdm: dmfrom() failed' end if stop 900 continue write(0,*) 'usage: tdm <file name>' write(0,*) write(0,*) ' Use ''tdm -'' to read ', & 'from standard input' stop end c=========================================================== c Demonstrates use of the real*8 (pseudo-)random number c generator, 'drand48' available on the lnx machines via c the 'p410f' utility library. c c usage: tdrand48 <nrand> [<seed>] c c where <nrand> is the number of random numbers on the c interval (0..1) to be generated, and <seed> is the c optional, integer-valued "seed" for the random number c generator. c c The program outputs the random numbers generated to c standard output (one per line), and their average to c standard error. In the limit of very large <nrand>, c this average should approach 0.5 since 'drand48' c generates numbers uniformly distributed on the unit c interval. c c Note carefully that for fixed seed, 'drand48' (and c most other pseudo-random number generators) will c ALWAYS RETURN THE SAME SEQUENCE OF NUMBERS! In c fact, the purpose of seeding the generator is precisely c to produce variation ("randomness") in the sequence c of iterates produced. If 'drand48' is not explicitly c seeded via 'srand48' a specific default seed value is c used. c=========================================================== program tdrand48 implicit none integer iargc, i4arg real*8 drand48 real*8 ranval, sum integer i, nrand if( iargc() .lt. 1 ) go to 900 nrand = i4arg(1,-1) if( nrand .le. 0 ) go to 900 if( iargc() .gt. 1 ) then call srand48(i4arg(2,0)) end if sum = 0.0d0 do i = 1 , nrand c----------------------------------------------------------- c Generate a random number c----------------------------------------------------------- ranval = drand48() sum = sum + ranval end do write(0,*) write(0,*) 'Average: ', sum / nrand stop 900 continue write(0,*) 'usage: tdrand48 <nrand> [<seed>]' stop end c=========================================================== c Returns a double precision vector (one-dimensional c array) read from file 'fname'. If 'fname' is the c string '-', the vector is read from standard input. c c The file should contain one number per line; invalid c input is ignored. c c This routine illustrates a general technique for c reading data from a FORMATTED (ASCII) file. In c Fortran, one associates a "logical unit number" c (an integer) with a file via the OPEN statement. c The unit number can then be used as the first c "argument" of the READ and WRITE statements to c perform input and output on the file. c c Fortran reserves the following unit numbers: c c 5 terminal input (stdin) c 6 terminal output (stdout) c 0 error output on Unix systems (stderr) c=========================================================== subroutine dvfrom(fname,v,n,maxn) c----------------------------------------------------------- c Arguments: c c fname: (I) File name c v: (O) Return vector c n: (O) Length of v (# read) c maxn: (I) Maximum number to read c----------------------------------------------------------- implicit none c----------------------------------------------------------- c The integer functions 'indlnb' and 'getu' are c defined in the 'p410f' library. c----------------------------------------------------------- integer indlnb, getu c----------------------------------------------------------- c Declaration of routine arguments: note c "adjustable dimensioning" of v; any array which c is declared with adjustable dimensions must be c a subroutine argument; any adjustable dimensions c must also be subroutine arguments. c----------------------------------------------------------- character*(*) fname integer n, maxn real*8 v(maxn) c----------------------------------------------------------- c Programming style: Use parameter (ustdin) rather c than constant value (5) for stdin logical unit # c----------------------------------------------------------- integer ustdin parameter ( ustdin = 5 ) c----------------------------------------------------------- c Local variables: c c vn: Current number read from input c ufrom: Logical unit number for READ c rc: For storing return status from READ c----------------------------------------------------------- real*8 vn integer ufrom, rc c----------------------------------------------------------- c Intialize c----------------------------------------------------------- n = 0 c----------------------------------------------------------- c Read from stdin? c----------------------------------------------------------- if( fname .eq. '-' ) then c----------------------------------------------------------- c Set unit number to stdin default c----------------------------------------------------------- ufrom = ustdin else c----------------------------------------------------------- c Get an available unit number c----------------------------------------------------------- ufrom = getu() c----------------------------------------------------------- c Open the file for formatted I/O c----------------------------------------------------------- open(ufrom,file=fname(1:indlnb(fname)), & form='formatted',status='old',iostat=rc) if( rc .ne. 0 ) then c----------------------------------------------------------- c Couldn't open the file, print error message c and return. c----------------------------------------------------------- write(0,*) 'dvfrom: Error opening ', & fname(1:indlnb(fname)) return end if end if c----------------------------------------------------------- c Input numbers into vector (one per line) until c EOF or maximum allowable number read c----------------------------------------------------------- 100 continue read(ufrom,*,iostat=rc,end=200) vn if( rc .eq. 0 ) then n = n + 1 if( n .gt. maxn ) then write(0,*) 'dvfrom: Read maximum of ', & maxn, ' from ', & fname(1:indlnb(fname)) n = maxn go to 200 end if v(n) = vn end if go to 100 200 continue c----------------------------------------------------------- c If we are reading from a file, close the file. c This releases the unit number for subsequent use. c----------------------------------------------------------- if( ufrom .ne. ustdin ) then close(ufrom) end if return end c=========================================================== c Test program for subroutine 'dvfrom'. c c Program expects one argument which is the filename c to be passed to 'dvfrom' c=========================================================== program tdvfrom implicit none c----------------------------------------------------------- c The integer function 'iargc' returns the number of c arguments supplied to the program. It is c automatically available to all Fortran programs on c most Unix systems, as is 'getarg' (see below). c----------------------------------------------------------- integer iargc, indlnb integer maxn parameter ( maxn = 100 000 ) real*8 v(maxn) integer n character*256 fname c----------------------------------------------------------- c Unless exactly one argument is supplied, print usage c message and exit. c----------------------------------------------------------- if( iargc() .ne. 1 ) then write(0,*) 'usage: tdvfrom <file name>' write(0,*) write(0,*) ' Use ''tdvfrom -'' to read ', & 'from standard input' stop end if c----------------------------------------------------------- c The subroutine 'getarg' (Unix) takes 2 arguments. c The first is an integer input argument specifying c which argument is to be fetched, the second is c a character output argument which, on return, c contains the fetched argument. c c Get the filename. c----------------------------------------------------------- call getarg(1,fname) c----------------------------------------------------------- c Call the routine ... c----------------------------------------------------------- call dvfrom(fname,v,n,maxn) c----------------------------------------------------------- c ... and report how many numbers were read. c----------------------------------------------------------- write(0,*) 'tdvfrom: ', n, ' read from '// & fname(1:indlnb(fname)) stop end c=========================================================== c Writes a double precision vector to file 'fname'. c If fname is the string '-' then the vector is written c to standard output. c=========================================================== subroutine dvto(fname,v,n) c----------------------------------------------------------- c Arguments: c c fname: (I) File name c v: (I) Vector to be written c n: (I) Length of vector c----------------------------------------------------------- implicit none integer getu, indlnb character*(*) fname integer n real*8 v(n) integer ustdout parameter ( ustdout = 6 ) integer i, uto, rc if( fname .eq. '-' ) then uto = ustdout else uto = getu() open(uto,file=fname(1:indlnb(fname)), & form='formatted',iostat=rc) if( rc .ne. 0 ) then write(0,*) 'dvto: Error opening ', & fname(1:indlnb(fname)) return end if end if do i = 1 , n write(uto,*) v(i) end do if( uto .ne. ustdout ) then close(uto) end if return end c=========================================================== c Test program for subroutine 'dvto'. c c Program expects two arguments, the name of a file c for output ('-' for stdout) and the length of the c test vector to be written. c=========================================================== program tdvto implicit none c----------------------------------------------------------- c The integer function 'i4arg' is defined in the c 'p410f' library. It takes two arguments, the first c is an integer specifying which program argument is c to be parsed as an integer, and the second is a c default value which will be returned if the argument c was not supplied or could not be converted to an c integer. c----------------------------------------------------------- integer iargc, i4arg integer maxn parameter ( maxn = 100 000 ) real*8 v(maxn) integer n integer i character*256 fname c----------------------------------------------------------- c Unless exactly two arguments are supplied, print usage c message and exit. c c Note the use of the "logical-if" statement (no then) c----------------------------------------------------------- if( iargc() .ne. 2 ) go to 900 call getarg(1,fname) n = i4arg(2,-1) if( n .eq. -1 ) go to 900 c----------------------------------------------------------- c Limit the value of n c----------------------------------------------------------- n = min(n,maxn) c----------------------------------------------------------- c Define test vector c----------------------------------------------------------- do i = 1 , n v(i) = i end do c----------------------------------------------------------- c Call the routine .. c----------------------------------------------------------- call dvto(fname,v,n) c----------------------------------------------------------- c Normal exit c----------------------------------------------------------- stop c----------------------------------------------------------- c Usage exit c----------------------------------------------------------- 900 continue write(0,*) 'usage: tdvto <file name> <n>' write(0,*) write(0,*) ' Use ''tdvto -'' to write ', & 'to standard output' stop end c=========================================================== c Demonstration main program and subprograms c illustrating the 'EXTERNAL' statement and how c subprograms may be passed as ARGUMENTS to other c subprograms. This technique is often used to c pass "user-defined" functions to routines which c can do generic things with such functions (such c as integrating or differentiating them, for example). c=========================================================== program texternal c----------------------------------------------------------- c The 'external' statement tells the compiler that the c specified names are names of externally-defined c subprograms (i.e. subroutines or functions) c----------------------------------------------------------- real*8 r8fcn external r8fcn, r8sub2 c----------------------------------------------------------- c Call 'r8fcncaller' which then invokes 'r8fcn' c----------------------------------------------------------- call r8fcncaller(r8fcn) c----------------------------------------------------------- c Call 'r8subcaller' which then invokes 'r8sub2' c----------------------------------------------------------- call subcaller(r8sub2) stop end c=========================================================== c Input 'fcn' is the name of an externally defined c real*8 function. This routine invokes that function c with argument 10.0d0 and writes the result on c standard error c=========================================================== subroutine r8fcncaller(fcn) implicit none real*8 fcn external fcn real*8 fcnval fcnval = fcn(10.0d0) write(0,*) 'r8caller: ', fcnval return end c=========================================================== c Input 'sub' is the name of an externally defined c subroutine. This routine invokes that subroutine c with arguments 10.0d0 and 20.0d0. c=========================================================== subroutine subcaller(sub) implicit none external sub call sub(10.0d0,20.0d0) return end c----------------------------------------------------------- c Demonstration real*8 function c----------------------------------------------------------- real*8 function r8fcn(x) implicit none real*8 x r8fcn = x**2 return end c=========================================================== c Demonstration subroutine c=========================================================== subroutine r8sub2(x,y) implicit none real*8 x, y write(0,*) 'r8sub: x = ', x, ' y = ', y return end c=========================================================== c Demonstration main program and subroutine to c illustrate use of SAVE and DATA statements. c=========================================================== program tsavedata implicit none integer i do i = 1 , 10 call sub1() end do stop end c=========================================================== c Subprogram 'sub1': writes a message to standard c error the FIRST time it is called, and writes c the number of times it has been called so far to c standard output EVERY time it is called. c=========================================================== subroutine sub1() implicit none logical first integer ncall c----------------------------------------------------------- c Strict f77 statement ordering demands that c ANY DATA statements appear after ALL variable c declarations. Note the use of '/' to delimit the c initialization value. c----------------------------------------------------------- data first / .true. / c----------------------------------------------------------- c This 'save' statement guarantees that ALL local c storage is preserved between calls. c----------------------------------------------------------- save if( first ) then ncall = 1 write(0,*) 'First call to sub1' first = .false. end if write(*,*) 'sub1: Call ', ncall ncall = ncall + 1 return end c=========================================================== c Demonstration main program, subroutines and functions c to illustrate argument passing (call by address) in c Fortran. c=========================================================== program tsub real*8 r8side integer n parameter ( n = 6 ) real*8 v1(n) real*8 a, b a = -1.0d0 b = 1.0d0 write(*,*) 'Pre r8swap: a = ', a, ' b = ', b call r8swap(a,b) write(*,*) 'Post r8swap: a = ', a, ' b = ', b call prompt('Through r8swap') a = 10.0d0 b = r8side(a) write(*,*) 'Post r8side: a = ', a, ' b = ', b call prompt('Through r8side') c----------------------------------------------------------- c Load 'v1' with 0.0d0 c----------------------------------------------------------- call dvloadsc(v1,n,0.0d0) call dvstderr('v1 loaded with 0.0',v1,n) call prompt('Through dvloadsc') c----------------------------------------------------------- c 'v1' and 'v1(1)' have the SAME ADDRESS and thus c this call to 'dvloadsc' has precisely the same effect c as the previous one. c----------------------------------------------------------- call dvloadsc(v1(1),n,0.0d0) call dvstderr('v1 loaded with 0.0',v1,n) call prompt('Through dvloadsc (second time)') c----------------------------------------------------------- c Load v(2:n-1) with 1.0d0, values 'v(1)' and 'v(n)' c are unchanged c----------------------------------------------------------- call dvloadsc(v1(2),n-2,1.0d0) call dvstderr('v1 loaded with 0.0 and 1.0',v1,n) call prompt('Through dvloadsc (third time)') c----------------------------------------------------------- c It is actually a violation of strict F77 to pass c the same address more than once to a subroutine c or argument, but in many cases, such as this one c it is perfectly safe. This sequence uses the c routine 'dvaddsc' to increment each value of 'v1' c by 2.0d0. c----------------------------------------------------------- call dvaddsc(v1,v1,n,2.0d0) call dvstderr('v1 incremented by 2.0',v1,n) call prompt('Through dvaddsc') call prompt('Through tsub') stop end c=========================================================== c This routine swaps its two real*8 arguments c=========================================================== subroutine r8swap(val1,val2) implicit none real*8 val1, val2 real*8 temp temp = val1 val1 = val2 val2 = temp return end c=========================================================== c Real*8 function 'r8side' which has the 'side effect' c of overwriting its argument with 0.0d0. As a general c matter of style, Fortran FUNCTION subprograms should c act like real functions (i.e. NO side-effects) where c possible. c c Also note that the name of a Fortran c function is treated as a local variable in the c subprogram source code and MUST be assigned a value c before any 'return' statements are encountered. c=========================================================== real*8 function r8side(x) implicit none real*8 x r8side = x * x * x x = 0.0d0 return end c=========================================================== c Loads output real*8 vector 'v' with input scalar c value 'sc'. c=========================================================== subroutine dvloadsc(v,n,sc) implicit none integer n real*8 v(n) real*8 sc integer i do i = 1 , n v(i) = sc end do return end c=========================================================== c Adds real*8 scalar to input real*8 vector 'v1', c and returns results in output real*8 vector 'v2' c=========================================================== subroutine dvaddsc(v1,v2,n,sc) implicit none integer n real*8 v1(n), v2(n) real*8 sc integer i do i = 1 , n v2(i) = v1(i) + sc end do return end c=========================================================== c Dumps 'string' and the real*8 vector 'v' to stderr. c=========================================================== subroutine dvstderr(string,v,n) implicit none character*(*) string integer n real*8 v(n) integer i write(0,*) string do i = 1 , n write(0,*) v(i) end do return end c=========================================================== c Prints a message on stdout and then waits for input c from stdin. c=========================================================== subroutine prompt(pstring) implicit none character*(*) pstring integer rc character*1 resp write(*,*) pstring write(*,*) 'Enter anything & <CR> to continue' read(*,*,iostat=rc,end=900) resp return 900 continue stop end !=========================================================== ! arraydemo90.f: Fortran 90 version of arraydemo.f ! ! Fortran 90 implements run-time allocation of arrays ! and other data objects. ! ! Array names are identified as dynamically allocated ! via the 'allocatable' attribute in the array ! declaration. Storage is allocated via the 'allocate' ! statement, and freed with 'deallocate'. !=========================================================== program arraydemo90 implicit none integer iargc, i4arg !----------------------------------------------------------- ! Identify a1, a2 and a3 as rank-2 allocatable arrays. ! ! Alternate equivalent declaration: ! ! real*8, allocatable:: a1(:,:), a2(:,:), a3(:,:) !----------------------------------------------------------- real*8, allocatable, dimension( : , : ) :: & a1, a2, a3 !----------------------------------------------------------- ! Run-time array bounds. !----------------------------------------------------------- integer n1, n2 !----------------------------------------------------------- ! Get the desired array bounds from the command-line ! and perform superficial check for validity. !----------------------------------------------------------- if( iargc() .ne. 2 ) go to 900 n1 = i4arg(1,-1) n2 = i4arg(2,-1) if( n1 .le. 0 .or. n2 .le. 0 ) go to 900 !----------------------------------------------------------- ! Allocate the arrays !----------------------------------------------------------- allocate( a1(n1,n2) ) allocate( a2(n1,n2) ) allocate( a3(n1,n2) ) !----------------------------------------------------------- ! Define and manipulate the 2-d arrays using various ! subroutines. !----------------------------------------------------------- call load2d( a1, n1, n2, 1.0d0 ) call load2d( a2, n1, n2, -1.0d0 ) call add2d( a1, a2, a3, n1, n2 ) !----------------------------------------------------------- ! Dump the 3 arrays to standard error. !----------------------------------------------------------- call dump2d( a1, n1, n2, 'a1' ) call dump2d( a2, n1, n2, 'a2' ) call dump2d( a3, n1, n2, 'a1 + a2' ) !----------------------------------------------------------- ! Deallocate the arrays !----------------------------------------------------------- deallocate(a1) deallocate(a2) deallocate(a3) stop 900 continue write(0,*) 'usage: arraydemo90 <n1> <n2>' stop end !----------------------------------------------------------- ! Loads a 2-D array with the values: ! ! a(i,j) = sc * (100 * j + i) !----------------------------------------------------------- subroutine load2d(a,d1,d2,sc) implicit none integer d1, d2 real*8 a(d1,d2) real*8 sc integer i, j do j = 1 , d2 do i = 1 , d1 a(i,j) = sc * (100.0d0 * j + i) end do end do return end !----------------------------------------------------------- ! Adds 2-D arrays 'a1' and 'a2' element-wise and returns ! result in 'a3' !----------------------------------------------------------- subroutine add2d(a1,a2,a3,d1,d2) implicit none integer d1, d2 real*8 a1(d1,d2), a2(d1,d2), a3(d1,d2) integer i, j do j = 1 , d2 do i = 1 , d1 a3(i,j) = a1(i,j) + a2(i,j) end do end do return end !----------------------------------------------------------- ! Dumps 2-d array labelled with 'label' on stderr !----------------------------------------------------------- subroutine dump2d(a,d1,d2,label) implicit none integer d1, d2 real*8 a(d1,d2) character*(*) label integer i, j, st if( d1 .gt. 0 .and. d2 .gt. 0 ) then write(0,100) label 100 format( /' <<< ',A,' >>>'/) do j = 1 , d2 st = 1 110 continue write(0,120) ( a(i,j) , i = st , min(st+7,d1)) 120 format(' ',8F9.3) st = st + 8 if( st .le. d1 ) go to 110 if( j .lt. d2 ) write(0,*) end do end if return end c=========================================================== c arraydemo.f: Program which demonstrates manipulation c of 'run-time' dimensioned arrays in Fortran. c c The program accepts two integer arguments which c specify the bounds for the two-dimensional arrays c which are to be defined and manipulated. c c The basic guidelines are as follows: c c (1) To deal with run-time defined dimensions, c perform all array manipulation (including c input and output) in SUBPROGRAMS rather c that the main program. c c (2) Always pass ALL bounds of an array, along c with the array itself, to subprograms which c are to manipulate the array. c c (3) Declare sufficient storage in the main routine c to deal with the largest array(s) you c anticipate dealing with, but make sure that c you always check that the size of the storage c is sufficient c c (4) An address of a location in a ONE dimensional c array can be passed to a subprogram expecting c a multi-dimensional array. c=========================================================== program arraydemo implicit none integer iargc, i4arg c----------------------------------------------------------- c Single-dimensioned array which can be used to provide c storage for the multi-dimensional array manipulation. c ("Poor-man's memory allocation") c----------------------------------------------------------- integer maxq parameter ( maxq = 100 000 ) real*8 q(maxq) c----------------------------------------------------------- c 'Pointer' to next available location in 'q' c----------------------------------------------------------- integer qnext c----------------------------------------------------------- c 'Pointers' for three 2-D arrays ('a1', 'a2', and 'a3') c----------------------------------------------------------- integer narray parameter ( narray = 3 ) integer a1, a2, a3 c----------------------------------------------------------- c Array bounds which are to be defined at run time c----------------------------------------------------------- integer n1, n2 c----------------------------------------------------------- c Get the desired array bounds from the command-line c and check that there is sufficient 'main-storage'. c----------------------------------------------------------- if( iargc() .ne. 2 ) go to 900 n1 = i4arg(1,-1) n2 = i4arg(2,-1) if( n1 .le. 0 .or. n2 .le. 0 ) go to 900 if( narray * n1 * n2 .gt. maxq ) then write(0,*) 'arraydemo: Insufficient main storage' stop end if c----------------------------------------------------------- c Initialize the main storage pointer ... c----------------------------------------------------------- qnext = 1 c----------------------------------------------------------- c ... and set up the 'pointers' for the two arrays c with bounds (n1,n2). c----------------------------------------------------------- a1 = qnext qnext = qnext + n1 * n2 a2 = qnext qnext = qnext + n1 * n2 a3 = qnext c----------------------------------------------------------- c Define and manipulate the 2-d arrays using various c subroutines. c----------------------------------------------------------- call load2d( q(a1), n1, n2, 1.0d0 ) call load2d( q(a2), n1, n2, -1.0d0 ) call add2d( q(a1), q(a2), q(a3), n1, n2 ) c----------------------------------------------------------- c Dump the 3 arrays to standard error. c----------------------------------------------------------- call dump2d( q(a1), n1, n2, 'a1' ) call dump2d( q(a2), n1, n2, 'a2' ) call dump2d( q(a3), n1, n2, 'a1 + a2' ) stop 900 continue write(0,*) 'usage: arraydemo <n1> <n2>' stop end c----------------------------------------------------------- c Loads a 2-D array with the values: c c a(i,j) = sc * (100 * j + i) c----------------------------------------------------------- subroutine load2d(a,d1,d2,sc) implicit none integer d1, d2 real*8 a(d1,d2) real*8 sc integer i, j do j = 1 , d2 do i = 1 , d1 a(i,j) = sc * (100.0d0 * j + i) end do end do return end c----------------------------------------------------------- c Adds 2-D arrays 'a1' and 'a2' element-wise and returns c result in 'a3' c----------------------------------------------------------- subroutine add2d(a1,a2,a3,d1,d2) implicit none integer d1, d2 real*8 a1(d1,d2), a2(d1,d2), a3(d1,d2) integer i, j do j = 1 , d2 do i = 1 , d1 a3(i,j) = a1(i,j) + a2(i,j) end do end do return end c----------------------------------------------------------- c Dumps 2-d array labelled with 'label' on stderr c----------------------------------------------------------- subroutine dump2d(a,d1,d2,label) implicit none integer d1, d2 real*8 a(d1,d2) character*(*) label integer i, j, st if( d1 .gt. 0 .and. d2 .gt. 0 ) then write(0,100) label 100 format( /' <<< ',A,' >>>'/) do j = 1 , d2 st = 1 110 continue write(0,120) ( a(i,j) , i = st , min(st+7,d1)) 120 format(' ',8F9.3) st = st + 8 if( st .le. d1 ) go to 110 if( j .lt. d2 ) write(0,*) end do end if return end c=========================================================== c Program illustrating "catastrophic" loss of precision c resulting from the subtraction of two nearly equal c floating point values. c=========================================================== program catprec implicit none real*8 x parameter ( x = 0.2d0 ) integer i real*8 h, dsinx write(*,*) ' h d(sin) approx '// & 'd(sin) exact d(sin) err' write(*,*) h = 0.5d0 do i = 1 , 16 c----------------------------------------------------------- c Algebraically, in the limit h -> 0, dsinx should c approach cos(x), but sin(x+h) -> sin(x) so c catastrophic loss of precision occurs. c----------------------------------------------------------- dsinx = (sin(x+h) - sin(x)) / h write(*,1000) h, dsinx, cos(x), dsinx - cos(x) 1000 format(1P,E12.3,2E16.8,E12.3) h = 0.125d0 * h end do stop end c=========================================================== c Implements matrix-matrix multiply c c c = a b c c where a, b and c are n x n (square) real*8 matrices. c=========================================================== subroutine dmmmult(a,b,c,n) implicit none integer n real*8 a(n,n), b(n,n), c(n,n) integer i, j, k do j = 1 , n do i = 1 , n c(i,j) = 0.0d0 do k = 1 , n c(i,j) = c(i,j) + a(i,k) * b(k,j) end do end do end do return end c=========================================================== c Writes a double precision matrix (two dimensional c array) to file 'fname'. If 'fname' is the c string '-', the matrix is written to standard input. c c This routine is modelled on 'dvto' previously c discussed in class: see ~phys410/f77/ex3/dvto.f c=========================================================== subroutine dmto(fname,a,d1,d2) c----------------------------------------------------------- c Arguments: c c fname: (I) File name c a: (I) Input matrix c d1: (I) First dimension of a c d2: (I) Second dimension of a c----------------------------------------------------------- implicit none integer indlnb, getu character*(*) fname integer d1, d2 real*8 a(d1,d2) integer ustdout parameter ( ustdout = 6 ) integer uto, rc c----------------------------------------------------------- c Parse fname: either "attach" 'uto' to stdout or c get a unit number using 'getu', and open the c file 'fname' for formatted I/O via 'uto' c----------------------------------------------------------- if( fname .eq. '-' ) then uto = ustdout else uto = getu() open(uto,file=fname(1:indlnb(fname)), & form='formatted',iostat=rc) if( rc .ne. 0 ) then write(0,*) 'dmto: Error opening ', & fname(1:indlnb(fname)) return end if end if c----------------------------------------------------------- c Write dimensions, then array elements c----------------------------------------------------------- write(uto,*,iostat=rc) d1, d2 if( rc .ne. 0 ) then write(0,*) 'dmto: Error writing dimensions' go to 500 end if write(uto,*,iostat=rc) a if( rc .ne. 0 ) then write(0,*) 'dmto: Error reading matrix' end if c----------------------------------------------------------- c Exit: Close file and return c----------------------------------------------------------- 500 continue if( uto .ne. ustdout ) then close(uto) end if return end c=========================================================== c Returns a double precision matrix (two dimensional c array) read from file 'fname'. If 'fname' is the c string '-', the matrix is read from standard input. c c The dimensions of the matrix must precede the matrix c elements themselves in the file. Specifically, the c file should have been created using the following c list-directed, formatted READ statement c (or equivalent): c c integer d1, d2 c real*8 a(d1,d2) c integer uout c c write(uout,*) d1, d2 c write(uout,*) a c c This routine is modelled on 'dvfrom' previously c discussed in class: see ~phys410/f77/ex3/dvfrom.f c c Note the use of helper routine 'dmfrom1' which c reads actual array values once bounds have been c extracted from file. c=========================================================== subroutine dmfrom(fname,a,d1,d2,asize) c----------------------------------------------------------- c Arguments: c c fname: (I) File name c a: (O) Return matrix c d1: (O) First dimension of a c d2: (O) Second dimension of a c asize: (I) Maximum size (d1 * d2) of a c----------------------------------------------------------- implicit none integer indlnb, getu character*(*) fname integer d1, d2, asize real*8 a(d1,d2) integer ustdin parameter ( ustdin = 5 ) integer ufrom, rc c----------------------------------------------------------- c Parse fname: either "attach" 'ufrom' to stdin or c get a unit number using 'getu', and open the c file 'fname' for formatted I/O via 'ufrom' c----------------------------------------------------------- if( fname .eq. '-' ) then ufrom = ustdin else ufrom = getu() open(ufrom,file=fname(1:indlnb(fname)), & form='formatted',iostat=rc,status='old') if( rc .ne. 0 ) then write(0,*) 'dmfrom: Error opening ', & fname(1:indlnb(fname)) return end if end if c----------------------------------------------------------- c Read dimensions and abort if there is insufficient c storage for the entire matrix. Note the 'go to' c to the 'exit block' since we've opened a file now c and should close it, even if there's an error. c Also, we set the dimensions to 0 for all error c conditions as a way of communicating failure to c the calling routine. c----------------------------------------------------------- read(ufrom,*,iostat=rc) d1, d2 if( rc .ne. 0 ) then write(0,*) 'dmfrom: Error reading dimensions' d1 = 0 d2 = 0 go to 500 end if if( (d1 * d2) .gt. asize ) then write(0,*) 'dmfrom: Insufficient storage' d1 = 0 d2 = 0 go to 500 end if c----------------------------------------------------------- c Now that dimensions have been determined call c helper routine to read values c----------------------------------------------------------- call dmfrom1(ufrom,a,d1,d2,rc) if( rc .ne. 0 ) then write(0,*) 'dmfrom: Error reading matrix' d1 = 0 d2 = 0 end if c----------------------------------------------------------- c Exit: Close file and return c----------------------------------------------------------- 500 continue if( ufrom .ne. ustdin ) then close(ufrom) end if return end c=========================================================== c Helper routine for dmfrom: Reads array values, returns c I/O status to calling routine via 'rc' c=========================================================== subroutine dmfrom1(ufrom,a,d1,d2,rc) implicit none integer d1, d2, ufrom, rc real*8 a(d1,d2) read(ufrom,*,iostat=rc) a return end
last update: Wed May 11, 2016