(j3.2006) MATMUL

Van Snyder van.snyder
Fri Nov 30 21:47:09 EST 2012

Bill Long wrote:
> Well, for 3 arguments, I'd use matmul(a1, matmul(a2,a3)).  For 4 
> arguments, I would use matmul(matmul(a1,a2),matmul(a3,a4)), etc. 
> Personally, I've never seem a case with more than three matrices.

Just because you've never seen it doesn't mean nobody ever needs it.  
I've recently seen a case with four, and another engineer says he needs 

As I mentioned in a recent message, which Bill might have missed, the 
amount of work for matrix multiplication can be vastly different between 
optimal and suboptimal parenthesizations.

For three factors,  matmul(a1, matmul(a2,a3)) and 
matmul(matmul(a1,a2),a3) are the only possibilities, so one could work 
out dynamically which one is better and then have an IF block.

For n factors, the number of parenthesizations is the Catalan number 
C(n-1), where C(n) = (2n)!/((n+1)!n!).  This has super-exponential 
growth, asymptotically O((2*n)**n).  Using 2-argument MATMUL, one could 
compute the optimum parenthesization, and then have a C(n-1)-way case 
selector to carry out the multiplications.  The values of C(n-1) for n = 
3..10 are 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786.  What are 
the chances of coding that correctly?  What are the chances of compiling it?

The point of a multiple-argument MATMUL is that the processor could 
determine the optimal parenthesization dynamically, based upon the 
extents of the arguments.  Optimal parenthesization is a simple 
dynamic-programming problem, for which an algorithm is described on 
pages 67--69 of "The Design and Analysis of Computer Algorithms" by Aho, 
Hopcroft and Ullman.

The desire for multi-argument MATMUL was not a trivial whim.

> Cheers,
> Bill
> On 11/30/12 3:43 PM, Van Snyder wrote:
>> On Wed, 2012-11-28 at 06:14 +0000, Whitlock, Stan wrote:
>>>     1.  MATMUL with more than 2 arguments - why put this in the
>>> standard except to force the vendors to do for free what you should
>>> buy from some math library company?  If enough customers want it,
>>> someone will sell it to you.
>> I took Stan's advice and started out to write a function with MATMUL's
>> functionality -- each argument can independently be any kind of any
>> numeric type --  with an arbitrary number of arguments.
>> I defined an abstract type with a deferred type-bound multiply
>> procedure.  I extended that with an "internal" type that has two class
>> pointers, and ten "leaf" types that have a rank-1 or rank-2 component of
>> default real, double precision, default complex, double complex, and
>> default integer.  I didn't use PDTs for the different kinds.  I wrote a
>> type-bound multiply procedure for each type.  The module is large.
>> I know how to compute the optimimum parenthesization, and build an
>> evaluation tree from it.
>> The one remaining problem is writing the function statement.
>> Fortran doesn't have a "..." dummy argument, and this wouldn't do the
>> right thing anyway if all it did was to specify an arbitrary number of
>> arguments with the same type, kind, and rank as the last named one.
>> I briefly contemplated a collection of generic routines, but the number
>> of them for three numeric types with two kinds each, plus logical, with
>> all arguments of rank 2 or one argument of rank 1, exploded as a
>> function of the number of arguments.  For 3..10 arguments, the numbers
>> of functions are 867, 6483, 46659, 326595, 2239491, 15116547, 100776963,
>> or 665127939, respectively.
>> I could allow, say, ten arguments of my base class, but that would
>> require users to start off with their matrices in one of my extension
>> types, or to put the target attribute on their matrices.  I wouldn't get
>> compile-time checking that at most one of the arguments was of rank one.
>> Any other ideas?
>> _______________________________________________
>> J3 mailing list
>> J3 at mailman.j3-fortran.org
>> http://mailman.j3-fortran.org/mailman/listinfo/j3

More information about the J3 mailing list