Path: utzoo!utgpu!news-server.csri.toronto.edu!rpi!usc!elroy.jpl.nasa.gov!jato!csi.jpl.nasa.gov!mwette
From: mwette@csi.jpl.nasa.gov (Matt Wette)
Newsgroups: comp.lang.fortran
Subject: Re: matrix multiplication
Keywords: parameters, pass-by-value, pass-by-reference
Message-ID: <1991May6.142020.29744@jato.jpl.nasa.gov>
Date: 6 May 91 14:20:20 GMT
References: <1991May04.170203.22222@ariel.unm.edu>
Sender: news@jato.jpl.nasa.gov
Reply-To: mwette@csi.jpl.nasa.gov (Matt Wette)
Organization: Jet Propulsion Laboratory, Pasadena, CA
Lines: 64
Nntp-Posting-Host: mr-ed.jpl.nasa.gov

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
 -----------------------------------------------------------------