scavo@cie.uoregon.edu (Tom Scavo) (05/05/91)
As you all know, a classic programming exercise is the matrix
multiplication problem
A = B C .
mxn mxp pxn
A robust algorithm will also handle the special case
X = Y X (1)
mxn mxm mxn
where the result is to replace one of the matrices being
multiplied. Such a calculation will break a naive algorithm,
however. One solution is to pass B and C by *value* (with
loss of efficiency in both time and space) but of course Fortran's
insistence that all parameters be passed by reference precludes
such an approach. On the other hand, pass-by-value can be
simulated using a local array (for the lhs) declared within the
matrix multiplication routine. Unfortunately, the dimensions of
this local array put an operational limit on the size of the input
matrices.
I already have an inadequate routine that uses such a local array.
How would I write a *general* matrix multiplication routine that
handles tricky computations such as those in (1), and accomodates
matrices of arbitrary size? Also, does anyone know of a Fortran
compiler that implements pass-by-value-result, a kind of delayed
pass-by-reference where the actual parameters are updated at the
point of return? Or is this non-standard?
Tom Scavo
scavo@cie.uoregon.edupmontgom@euphemia.math.ucla.edu (Peter Montgomery) (05/05/91)
In article <1991May04.170203.22222@ariel.unm.edu> scavo@cie.uoregon.edu (Tom Scavo) writes: >I already have an inadequate routine that uses such a local array. >How would I write a *general* matrix multiplication routine that >handles tricky computations such as those in (1), and accommodates >matrices of arbitrary size? Convince your vendor to support certain Fortran 90 features. Matrix multiplication is an F90 intrinsic function. My nominee for the more important F90 feature is automatic arrays: local arrays whose dimensions can depend upon the values of parameters passed to the subroutine (and upon common variables). [Also high on my list are the bit manipulation functions and IMPLICIT NONE.] With automatic arrays you could allocate a product array of the desired size upon entry, build the product there, and copy it before exiting. -- Peter L. Montgomery pmontgom@MATH.UCLA.EDU Department of Mathematics, UCLA, Los Angeles, CA 90024-1555 If I spent as much time on my dissertation as I do reading news, I'd graduate.
dik@cwi.nl (Dik T. Winter) (05/05/91)
In article <1991May04.170203.22222@ariel.unm.edu> scavo@cie.uoregon.edu (Tom Scavo) writes: > A robust algorithm will also handle the special case > X = Y X (1) > mxn mxm mxn > where the result is to replace one of the matrices being > multiplied. Not in Fortran. > On the other hand, pass-by-value can be > simulated using a local array (for the lhs) declared within the > matrix multiplication routine. Unfortunately, the dimensions of > this local array put an operational limit on the size of the input > matrices. Still your routine will not be conforming to the Fortran standard. > > How would I write a *general* matrix multiplication routine that > handles tricky computations such as those in (1), and accomodates > matrices of arbitrary size? If you are thinking about portrable Fortran you must forget it. It is not permitted. To give an example of the problem, suppose the following (simple) routine: SUBROUTINE ADD(A, B, C) A = B + C RETURN END The following calls are illegal: CALL ADD(P, P, Q) CALL ADD(P, Q, P) The reason is that *if* one of the parameters changes value within the routine that parameter may not be aliased to another on the call side. Some compilers allow this, but some choke on a number of cases. > Also, does anyone know of a Fortran > compiler that implements pass-by-value-result, a kind of delayed > pass-by-reference where the actual parameters are updated at the > point of return? Or is this non-standard? This is certainly non-standard in general. For simple variables it is implemented on a large number of compilers (IBM amongst them). -- dik t. winter, cwi, amsterdam, nederland dik@cwi.nl
mwette@csi.jpl.nasa.gov (Matt Wette) (05/06/91)
In article <1991May04.170203.22222@ariel.unm.edu>, scavo@cie.uoregon.edu (Tom Scavo) writes: |> As you all know, a classic programming exercise is the matrix |> multiplication problem |> |> A = B C . |> mxn mxp pxn |> |> A robust algorithm will also handle the special case |> |> X = Y X (1) |> mxn mxm mxn |> |> where the result is to replace one of the matrices being |> multiplied. Computing "X <- Y * X" is actually quite simple to do as a Fortran routine if you're willing to pass a "work vector". My version is included below. In fact, this case is quite efficient. A nastier one (in Fortran) is "X <- X * Y". A "general" routine could branch to code which could handle your particular case. Note that you have to take care of transposes also (i.e., "X <- Y' * X"). I think the convenience you want from pass-by-value is not worth the large loss in efficiecy you'll get from copying (possibly) large arrays. Matt SUBROUTINE MULBRR(L, M, N, A, NA, B, NB, WKV) C INTEGER L, M, N, NA, NB DOUBLE PRECISION A(NA,1), B(NB,1), WKV(1) C C B = A * B A IS (L X M), B(RHS) IS (M X N), B(LHS) IS (L X N). C WKV IS A WORK VECTOR OF SIZE L. C C $Id: mulbrr.f,v 1.1 1991/01/24 00:50:12 mwette Exp $ C INTEGER I,J,K C IF (L .LE. 0 .OR. M .LE. 0 .OR. N .LE. 0) RETURN DO 50 J = 1,N DO 10 I = 1,L WKV(I) = 0.0D0 10 CONTINUE DO 30 K = 1,M DO 20 I = 1,L WKV(I) = WKV(I) + A(I,K)*B(K,J) 20 CONTINUE 30 CONTINUE DO 40 I = 1,L B(I,J) = WKV(I) 40 CONTINUE 50 CONTINUE RETURN C C --- LAST LINE OF MULBRR --- END -- _________________________________________________________________ Matthew R. Wette | Jet Propulsion Laboratory, 198-326 mwette@csi.jpl.nasa.gov | 4800 Oak Grove Dr, Pasadena,CA 91109 -----------------------------------------------------------------
jlg@cochiti.lanl.gov (Jim Giles) (05/06/91)
In article <1991May04.170203.22222@ariel.unm.edu>, scavo@cie.uoregon.edu (Tom Scavo) writes: |> [...] Also, does anyone know of a Fortran |> compiler that implements pass-by-value-result, a kind of delayed |> pass-by-reference where the actual parameters are updated at the |> point of return? Or is this non-standard? It is quite within the standard to pass procedure arguments by value/result. It is also standard to pass by reference. The reason for the Fortran prohibition of aliasing parameters through procedure calls is to make the semantics of both forms of procedure call identical. That's why you can't do MATMLT(X,Y,X,M,N,M) - if the procedure assigns to X. If you were allowed to do so, the semantics of the two parameter passing methods would be different. Yes, you can implement Fortran with value/result parameters. The bad news is that it still doesn't do what you want since aliasing IN and INOUT (or OUT) parameters is illegal. J. Giles
scavo@cie.uoregon.edu (Tom Scavo) (05/07/91)
In an earlier article, I naively ask: > > Also, does anyone know of a Fortran > > compiler that implements pass-by-value-result, a kind of delayed > > pass-by-reference where the actual parameters are updated at the > > point of return? Or is this non-standard? In article <3444@charon.cwi.nl> dik@cwi.nl (Dik T. Winter) writes: >This is certainly non-standard in general. For simple variables it >is implemented on a large number of compilers (IBM amongst them). In article <23261@lanl.gov> jlg@cochiti.lanl.gov (Jim Giles) writes: >It is quite within the standard to pass procedure arguments by value/result. >It is also standard to pass by reference. The reason for the Fortran >prohibition of aliasing parameters through procedure calls is to make >the semantics of both forms of procedure call identical. That's why >you can't do MATMLT(X,Y,X,M,N,M) - if the procedure assigns to X. If >you were allowed to do so, the semantics of the two parameter passing >methods would be different. So what is the concensus on this issue? Is it within the standard to implement pass-by-value-result? If so, what compilers have gone this route? As an aside, how did it happen that two competing parameter passing mechanisms got into the Fortran 77 standard? (This question assumes, of course, that this is indeed the case.) Was it an error of omission on the part of the standards committee, or were the insecurities of pass-by-reference simply not known at the time? (I'm just interested in a proper historical perspective.) Looking to the future, how does Fortran 90 deal with the parameter passing problem? Tom Scavo scavo@cie.uoregon.edu
dik@cwi.nl (Dik T. Winter) (05/08/91)
In article <1991May07.155218.13741@ariel.unm.edu> scavo@cie.uoregon.edu (Tom Scavo) writes: > So what is the concensus on this issue? Is it within the > standard to implement pass-by-value-result? If so, what > compilers have gone this route? I fail to see any difference between what Jim Giles said and what I said. Pass-by-value-result is allowed, as is pass-by-reference. The standard allows both. As I said, IBM compilers use pass-by-value-result, but CDC compilers use pass-by-reference. > > As an aside, how did it happen that two competing parameter > passing mechanisms got into the Fortran 77 standard? (This > question assumes, of course, that this is indeed the case.) This is indeed *not* the case. The standard simply leaves the parameter passing mechanisms unspecified. It only specifies the result for non-aliased parameters. Any program assuming one way or the other depends on improper aliasing, and that makes the program illegal. If you want a detour to Ada: in Ada it *is* specified that simple parameters are passed-by-value-return. Non-simple parameters are to be passed by either value-return or reference. The reason is that if the subroutine you are calling resides on a different processor you really *want* value-return in most cases. -- dik t. winter, cwi, amsterdam, nederland dik@cwi.nl
khb@chiba.Eng.Sun.COM (Keith Bierman fpgroup) (05/09/91)
In article <1991May07.155218.13741@ariel.unm.edu> scavo@cie.uoregon.edu (Tom Scavo) writes:
As an aside, how did it happen that two competing parameter
passing mechanisms got into the Fortran 77 standard? (This
The standard was carefully crafted to allow compilers to do the
highest performing thing on their platform. If you code according to
the standard, you can't tell how parameter passing occurs.
On traditional "scalar" uniprocessors, passing addresses is often the
best way. On machines like the TMC and MASSPAR, copying is often
cheaper.
Looking to the future, how does Fortran 90 deal with the
parameter passing problem?
The same way; the standard doesn't mandate how you do it; just what
the effects are. With interface blocks, the coder can communicate with
the compiler a bit better (parameters passed just as IN, or OUT may
merit special treatment, for example).
--
----------------------------------------------------------------
Keith H. Bierman keith.bierman@Sun.COM| khb@chiba.Eng.Sun.COM
SMI 2550 Garcia 12-33 | (415 336 2648)
Mountain View, CA 94043