New intrinsic function proposed
Author Message
New intrinsic function proposed

One month in fortran 90 I have benifited greatly from its array notation and
intrinsic functions.  One common operation in statistics/maths, however, has
no easy solution.

real:: a(n),b(n),c(n,n)

do i=1,n
do j=1,n
c(i,j)=a(i)*b(j)
end do
end do

Is there a better way?  Or should we have a new intrinsic, say, out_product.
We can then write

c=out_product(a,b).

J.G. Liao

Sat, 30 Aug 1997 05:45:04 GMT
New intrinsic function proposed

Quote:
>One month in Fortran 90 I have benifited greatly from its array notation and
>intrinsic functions.  One common operation in statistics/maths, however, has
>no easy solution.

>real:: a(n),b(n),c(n,n)

>do i=1,n
>  do j=1,n
>    c(i,j)=a(i)*b(j)
>  end do
>end do
>Is there a better way?  Or should we have a new intrinsic, say, out_product.
>We can then write

>c=out_product(a,b).

Yes, such calculations are *very* common. In fact this is simply the
matrix multiplication:

T
C =  A * B    (if you consider 1-dim arrays as column vectors)

If F95 will allow the intrinsic routine MATMUL to let *both* arguments
be of rank 1 the syntax of such calculations could be:

C = MATMUL(A,B)

Has this obvious extention of MATMUL ever been considered?????????

This feature would be very important.

Best Regards,
Jens Bloch Helmers.

Sat, 30 Aug 1997 10:21:56 GMT
New intrinsic function proposed

Quote:

>One month in Fortran 90 I have benifited greatly from its array notation and
>intrinsic functions.  One common operation in statistics/maths, however, has
>no easy solution.

>real:: a(n),b(n),c(n,n)

>do i=1,n
>  do j=1,n
>    c(i,j)=a(i)*b(j)
>  end do
>end do

>Is there a better way?  Or should we have a new intrinsic, say, out_product.
>We can then write

>c=out_product(a,b).

>J.G. Liao

You can do an outer product using the intrinsic function spread and
element-wise matrix multiplication.  Your loop is replaced by the following:

spread(a,2,n) is a matrix of the form

a(1)    a(1)    a(1) . . .
a(2)    a(2)    a(2) . . .
.  .  .

and spread(b,1,n) is a matrix of the form

b(1)    b(2)    b(3) . . .
b(1)    b(2)    b(3) . . .
.  .  .

--
Ronald Sverdlove               Computational Science Research

Tel. 609-734-2517              CN 5300
FAX  609-734-2662              Princeton, NJ 08543-5300

Sun, 31 Aug 1997 02:41:34 GMT
New intrinsic function proposed

Quote:

> One month in Fortran 90 I have benifited greatly from its array notation and
> intrinsic functions.  One common operation in statistics/maths, however, has
> no easy solution.

> real:: a(n),b(n),c(n,n)

> do i=1,n
>   do j=1,n
>     c(i,j)=a(i)*b(j)
>   end do
> end do

> Is there a better way?  Or should we have a new intrinsic, say, out_product.
> We can then write

> c=out_product(a,b).

Even though FORTRAN 90 didn't directly implement this as an intrinsic
function, the second-level Basic Linear Algebra Subroutines (BLAS2) implement
a "rank-two update", which is a generalization of this kernel.

Any serious FORTRAN numerical programmer, whether using FORTRAN 90 or not,
netlib (at URL "http://netlib.att.com/netlib/blas/index.html").
--

Coordinated Science Laboratory     University of Illinois at Urbana-Champaign

Sun, 31 Aug 1997 04:26:22 GMT
New intrinsic function proposed

Quote:

> >One month in Fortran 90 I have benifited greatly from its array notation and
> >intrinsic functions.  One common operation in statistics/maths, however, has
> >no easy solution.

> >real:: a(n),b(n),c(n,n)

> >do i=1,n
> >  do j=1,n
> >    c(i,j)=a(i)*b(j)
> >  end do
> >end do

> >Is there a better way?  Or should we have a new intrinsic, say, out_product.
> >We can then write

> >c=out_product(a,b).

> Yes, such calculations are *very* common. In fact this is simply the
> matrix multiplication:

>                 T
>       C =  A * B    (if you consider 1-dim arrays as column vectors)

> If F95 will allow the intrinsic routine MATMUL to let *both* arguments
> be of rank 1 the syntax of such calculations could be:

>       C = MATMUL(A,B)

Of course, this is also the syntax of a dot product:
T
D = A * B    (if you consider 1-dim arrays as column vectors)

How would you tell the difference?

Quote:
> Has this obvious extention of MATMUL ever been considered?????????

> This feature would be very important.

It is less important than you might think when you consider the SPREAD
intrinsic.  Consider the following example:

\$ cat outprod.f90
program outprod
real a(4),b(3),c(4,3)
a = (/ (i, i = 1, 4) /)
b = (/ (i, i = 1, 3) /)
print *, c
end program outprod
\$
\$
\$ f90 outprod.f90
\$ a.out
1.,  2.,  3.,  4.,  2.,  4.,  6.,  8.,  3.,  6.,  9.,  12.
\$#^^^^^col 1^^^^^^   ^^^^^^col 2^^^^^^   ^^^^^^col 3^^^^^^^

Of course, this is *really* ugly, so the original suggestion of an
outer product intrinsic might be useful after all.

----------------------------------------------- Roger Glover
XXXX\ \ / \ /XXX  \ / \ X   \ /\\\  X///X /\\\  Cray Research, Inc.
/ \ / \/ /\ / \ / \X /\ X  \  / \   X\  \ X     DISCLAIMER HAIKU:
//X/  X\\\X //X/  X \ X X\\   / \   X/X \ X \\\   C R I may not
/ \   X///X / \/  X//XX X  \  / \   X  \\ X   \   Share these opinions with me
/ \   X   X /\\\/ X   X X///X /XXX/ X///X /XXX/   This is not my fault

Tue, 02 Sep 1997 01:49:09 GMT
New intrinsic function proposed
Quote:
>> Yes, such calculations are *very* common. In fact this is simply the
>> matrix multiplication:

>>                 T
>>       C =  A * B    (if you consider 1-dim arrays as column vectors)

>> If F95 will allow the intrinsic routine MATMUL to let *both* arguments
>> be of rank 1 the syntax of such calculations could be:

>>       C = MATMUL(A,B)

>Of course, this is also the syntax of a dot product:
>         T
>        D = A * B    (if you consider 1-dim arrays as column vectors)

>How would you tell the difference?

T
F95 could simply *define* MATMUL(A,B) to be  A * B  if A *and* B has rank 1.
There is no conflict of interest because

T
A * B can (and should) be calculated as DOT_PRODUCT(A,B)

Quote:
>It is less important than you might think when you consider the SPREAD
>intrinsic. .......

real a(m), b(n), c(m,n)

c = matmul( reshape(a,(/m,1/)), reshape(b,(/1,n/)) )

Anyway we both use the MATMUL function, so why don't extend the functionality
of this function instead of inventing a complete new intrinsic?

Quote:
>Of course, this is *really* ugly, so the original suggestion of an
>outer product intrinsic might be useful after all.

Usually I don't like the idea of adding a whole lot of new intrinsics,
but the outer product realy belongs among the other native basic array
features of this wonderful language.
(IMHO)

Jens B. Helmers.

Tue, 02 Sep 1997 12:46:30 GMT
New intrinsic function proposed
It should be noted that I believe that FORALL, an HPF construct that is
part of the F9x draft partially addresses this problem while providing a
more general capability.  Correct me if I am wrong, but I believe the
construct
real:: a(n),b(n),c(n,n)

do i=1,n
do j=1,n
c(i,j)=a(i)*b(j)
end do
end do

is written as

real:: a(n),b(n),c(n,n)

FORALL (i=1:n, j=1:n)  c(i,j) = a(i)*b(j)

While not as succinct as the outer product, it is more general and more
succinct than the nested do constructs.

Thu, 04 Sep 1997 06:31:56 GMT
New intrinsic function proposed
What you need is the FORALL construct:

FORALL(i=1:n, j=1:n) c(i,j)=a(i)*b(j)

This isn't the exact form of anyone's FORALL proposal, but it captures the
concept.  FORALL is a subject of debate (not whether to include it, but what
syntax to use).  It will almost certainly show up in the next standard.

Personally, I prefer to have a class of iteration variables with limited
scope which would be self-declared (the form of their identifier denotes
their type).  The above would then be written:

C(%i,%j) = a(%i)*b(%j)

The iteration variables would automatically vary over the range of the index
they were used in place of (unless a triplet notation is also used).  The
use of each iteration variable would have to be conformable with all other
uses in the same statement.  The scope of the variables is one statement
only - so the same iteration variables can be used on other statements with
different limits.  The above written with explicit triplet notation is:

C(%i:1:n,%j:1:n) = a(%i:1:n)*b(%j:1:n)

You could also stride.  The only restriction is that each variable (%i for
example) must iterate over the same number of elements in each place it's
used.  Reduction indices are also possible.  The sum reduction is:

sum = a(#i)

Here, #i iterates over all the values of a and sums over them.  The precedence
of the summation operation is below multiplication and above addition.  So,
common matrix multiply is:

A(%i,%k) = B(%i,#j) * C(#j,%k)

Like I said, I would prefer this kind of operation, but don't expect the
committee to ever adopt anything this simple.

Cheers.

Thu, 04 Sep 1997 07:01:33 GMT
New intrinsic function proposed
: Yes, such calculations are *very* common. In fact this is simply the
: matrix multiplication:

:                 T
:       C =  A * B    (if you consider 1-dim arrays as column vectors)

: If F95 will allow the intrinsic routine MATMUL to let *both* arguments
: be of rank 1 the syntax of such calculations could be:

:       C = MATMUL(A,B)

: Has this obvious extention of MATMUL ever been considered?????????

: This feature would be very important.

: Best Regards,
: Jens Bloch Helmers.

if such extension is allowed, there must be a way to distingush it from
dot_product(a,b).  Maybe a new intrinsic function is a better way.

J.G. Liao

Fri, 05 Sep 1997 01:17:15 GMT
New intrinsic function proposed

Helmers> matrix multiplication:
Helmers>   C =  A * B    (if you consider 1-dim arrays as column vectors)

Helmers> If F95 will allow the intrinsic routine MATMUL to let *both* arguments
Helmers> be of rank 1 the syntax of such calculations could be:
Helmers>       C = MATMUL(A,B)
Helmers> Has this obvious extention of MATMUL ever been considered?????????

Helmers> This feature would be very important.

Liao>   if such extension is allowed, there must be a way to distingush it from
Liao>   dot_product(a,b).  Maybe a new intrinsic function is a better way.

I'd have to agree that a new intrinsic would a less confusing
approach.  The name outer_product comes to mind.  Otherwise, it
is not clear which product is intended.

However, there is no reason why this *has* to be an intrinsic.
It could just as easily be a user-written function - a rather
trivial one at that.  And the user could even define it as an
operator, allowing you to write things like

C = A .outer. B

if that syntax strikes your fancy.

For dot product, it makes special sense for the function to be
intrinsic because some vendors may be able to do "smart" things
with it - like use extra precision for storage of the intermediate
result, or even use special hardware instructions.  For the outer
product, I don't know off-hand of anything that the vendor could
probably do much better than the obvious (which I freely admit
could just reflect my ignorance - or perhaps I could get by with

It might make sense to have outer_product as an intrinsic just
to provide symmetry, but I don't see it as a very big issue given
how easy it is to write this as a user function.  Let's see,
off the top of my head, without even running it through a compiler
for syntax checking - so caveat emptor - something like
function outer_product(a,b) result(c)
real, intent(in) :: a(:), b(:)
real :: c(size(a),size(b))
integer :: i,j
do j = 1 , size(b)
do i = 1 , size(a)
c(i,j) = a(i)*b(j)
end do
end do
end function outer_product

--
--
Richard Maine

Sat, 06 Sep 1997 01:11:07 GMT

 Page 1 of 1 [ 10 post ]

Relevant Pages