• #### A MATRIX type for Fortran

From Quadibloc@21:1/5 to All on Mon Oct 25 23:15:10 2021
I will have to admit that I'm not fully conversant
with the latest developments in Fortran.

But it is my understanding that starting around the
time of Fortran 90, constructs have been added to
the language so that programs can be written to
efficiently make use of the vector arithmetic features
of computers like the Cray I, and its successors and
imitators.

These days, except for the SX-Aurora TSUBASA from
NEC, most vector arithmetic on computers belongs to
a line of development that started from the AN/FSQ-7
for SAGE rather than the one that went through the
Cray I.

Be that as it may...

So, IIRC, one can write

A = B * C

in Fortran today and if A, B, and C are arrays having the
same dimensions (or if one of B or C is a scalar) it will
do an element-by-element multiplication, much as would
have been done by APL, without the need for a DO loop.

I remember that back in the 1970s I mentioned to someone
that I thought this sort of thing would be a nice feature to
have, even without vector hardware. The fellow replied that
I was being too slavish in following APL, matrix multiplication
would be better.

Well, I agreed that matrix multiplication was useful in
scientific and engineering computation.

to the construct of an array in Fortran.

Of course, one can add a built-in subroutine for matrix
multiplication.

But more recently, I thought that there might be a better
solution; and now on this newsgroup, I've thought to
present it more widely.

Fortran has a COMPLEX data type.

Why can't it have a MATRIX data type?

So there would be a MATRIX statement, having much the
same syntax as a DIMENSION statement. It could be used
to create one-dimensional vectors and matrices that were
not square.

If A, B, and C are declared as matrices (or one of B or C is
a scalar) then

A = B * C

would perform matrix multiplication as expected... if one of
B or C is a one-dimensional matrix, it would be treated as
either a row or column vector, whichever one would make
sense. (Unless that isn't practical; declaring 1 by n and n by 1
matrices instead to make it explicit is certainly an option as
well.)

Since a matrix is an a x b grid of floating-point numbers
combined into an object with specific mathematical
properties, while it would be possible to create an array
of matrices, it would *not* be possible to create a matrix
the elements of which are arrays.

For much the same reason that you can't make a
complex number the real and imaginary parts of which
are arrays of the same shape.

But, yes, it should be allowed to create a matrix the elements
of which are complex numbers. (Or quaternions, for that
matter!)

Oh, and by the way: yes, it should be possible to take a (4,4)
matrix the elements of which are REAL*16 and EQUIVALENCE
it to a (4,4) array the elements of which are REAL*16 with the
expected results of the same-numbered elements perfectly
corresponding. (This might even be a reason not to deprecate
EQUIVALENCE...)

Is this not a useful thing to include in Fortran?

Maybe someone else has even suggested it already?

John Savard

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From gah4@21:1/5 to Quadibloc on Tue Oct 26 03:47:27 2021
On Monday, October 25, 2021 at 11:15:12 PM UTC-7, Quadibloc wrote:
I will have to admit that I'm not fully conversant
with the latest developments in Fortran.

But it is my understanding that starting around the
time of Fortran 90, constructs have been added to
the language so that programs can be written to
efficiently make use of the vector arithmetic features
of computers like the Cray I, and its successors and
imitators.

Complicated question.

In the Cray days, there was much work on getting compilers to figure
out how to make good use of vector registers.

The Fortran added array expressions, which are, among others, supposed
to make it easier for vector processors. There is also FORALL, which was more specifically added to allow for vector processors.

Except that both don't work quite that well.

Both array expressions and FORALL are defined such that the whole
right hand side is evaluated before any changes are made to the left
side. Unless vector registers are big enough for the whole array, that
can mean temporary arrays.

Even worse, many array expressions, as written, are harder for compilers
to optimize well, than the DO form. That is especially true when one can
leave a computation early.

And then vector processors went out of style.

These days, except for the SX-Aurora TSUBASA from
NEC, most vector arithmetic on computers belongs to
a line of development that started from the AN/FSQ-7
for SAGE rather than the one that went through the
Cray I.

(snip)

I remember that back in the 1970s I mentioned to someone
that I thought this sort of thing would be a nice feature to
have, even without vector hardware. The fellow replied that
I was being too slavish in following APL, matrix multiplication
would be better.

Well, I agreed that matrix multiplication was useful in
scientific and engineering computation.

to the construct of an array in Fortran.

Of course, one can add a built-in subroutine for matrix
multiplication.

Fortran does have a built-in function for matrix
multiplication, MATMUL. Are expressions easier to read
with the * operator as a matrix multiply? I don't know.

But more recently, I thought that there might be a better
solution; and now on this newsgroup, I've thought to
present it more widely.

Fortran has a COMPLEX data type.

Why can't it have a MATRIX data type?

You can do it with defined type, and matrix multiply
through defined operations. (That is, define * to be
the matrix multiply operator for TYPE(MATRIX).)

So there would be a MATRIX statement, having much the
same syntax as a DIMENSION statement. It could be used
to create one-dimensional vectors and matrices that were
not square.

It seems a lot of work, just to get the * operator to change.

It would not be so hard to add an operator with different
syntax, that would be a matrix multiply operator.

This is completely separate from the question about vector,
or otherwise, processor. In either the current form, or one
with a specific operator, it is just as easy for compilers
to figure out.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From FortranFan@21:1/5 to Quadibloc on Tue Oct 26 06:41:55 2021
On Tuesday, October 26, 2021 at 2:15:12 AM UTC-4, Quadibloc wrote:
I will have to admit that I'm not fully conversant
with the latest developments in Fortran.
..

Why can't it have a MATRIX data type?
..

If A, B, and C are declared as matrices (or one of B or C is
a scalar) then

A = B * C

would perform matrix multiplication as expected... ..

Is this not a useful thing to include in Fortran?

Maybe someone else has even suggested it already?
..

@John Savard,

Given your first sentence, you may want to consider the following even if it is just an "academic" interest:
1) Grab a copy of "Modern Fortran Explained" that covers up to Fortran 2018, the current standard revision. With your experience it will be a quick read for you:
https://oxford.universitypressscholarship.com/view/10.1093/oso/9780198811893.001.0001/oso-9780198811893

2) Also post and follow modern Fortran discourse at this site: https://fortran-lang.discourse.group/

a) Fortran has long had native support for multidimensional arrays,

b) Toward ordinary operations assuming SHAPE conformance, the language supports the "vector" notation you show:
C = A + B

c) Starting with Fortran 90 introduced nearly 30 years ago, the language support for the outstanding "matrix" operation you show is

A = matmul( B, C )

with the details of achieving efficient matrix multiplication left mostly to the processors.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From pehache@21:1/5 to All on Tue Oct 26 17:08:41 2021
Le 26/10/2021 à 12:47, gah4 a écrit :
On Monday, October 25, 2021 at 11:15:12 PM UTC-7, Quadibloc wrote:
I will have to admit that I'm not fully conversant
with the latest developments in Fortran.

But it is my understanding that starting around the
time of Fortran 90, constructs have been added to
the language so that programs can be written to
efficiently make use of the vector arithmetic features
of computers like the Cray I, and its successors and
imitators.

Complicated question.

In the Cray days, there was much work on getting compilers to figure
out how to make good use of vector registers.

The Fortran added array expressions, which are, among others, supposed
to make it easier for vector processors. There is also FORALL, which was more
specifically added to allow for vector processors.

Wasn't FORALL rather designed for general parallelism, allowing out of
order execution of the loop ?

Except that both don't work quite that well.

Both array expressions and FORALL are defined such that the whole
right hand side is evaluated before any changes are made to the left
side. Unless vector registers are big enough for the whole array, that
can mean temporary arrays.

Or unless the compiler can determine that this is safe to update the LHS
during the evaluation of the RHS. But yes, this means more complexity on
the compiler side.

Even worse, many array expressions, as written, are harder for compilers
to optimize well, than the DO form. That is especially true when one can leave a computation early.

And then vector processors went out of style.

Yes, and no. Most of modern CPUs have small vector registers.

--
"...sois ouvert aux idées des autres pour peu qu'elles aillent dans le
même sens que les tiennes.", ST sur fr.bio.medecine

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From gah4@21:1/5 to All on Tue Oct 26 11:42:41 2021
On Tuesday, October 26, 2021 at 8:08:45 AM UTC-7, pehache wrote:

(snip, I wrote)

Complicated question.

In the Cray days, there was much work on getting compilers to figure
out how to make good use of vector registers.

The Fortran added array expressions, which are, among others, supposed
to make it easier for vector processors. There is also FORALL, which was more
specifically added to allow for vector processors.

Wasn't FORALL rather designed for general parallelism, allowing out of
order execution of the loop ?

It was about the time of Cray, and I thought I knew that it was related
to vector processors. Now we have DO CONCURRENT that, as far as
I know, is related to out of order.

Except that both don't work quite that well.

Both array expressions and FORALL are defined such that the whole
right hand side is evaluated before any changes are made to the left
side. Unless vector registers are big enough for the whole array, that
can mean temporary arrays.

Or unless the compiler can determine that this is safe to update the LHS during the evaluation of the RHS. But yes, this means more complexity on
the compiler side.

Yes, I said "can". Maybe I should have made that more obvious.

But if you have infinite sized vector registers, then you wouldn't have
to consider that. But you never have those.

Even worse, many array expressions, as written, are harder for compilers
to optimize well, than the DO form. That is especially true when one can leave a computation early.

And then vector processors went out of style.

Yes, and no. Most of modern CPUs have small vector registers.

OK, but for it to make sense, vector registers have to keep getting
bigger, not smaller. Well, Cray registers were based on the processor pipeline, and that it could keep them full through a series of vector instructions.

Before Cray, there were a series of pipelined scalar processors,
such as the IBM 360/91, and the CDC 7600, but it is much easier
to keep them full with large vector registers.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Quadibloc@21:1/5 to All on Tue Oct 26 17:14:08 2021
On Tuesday, October 26, 2021 at 4:47:28 AM UTC-6, gah4 wrote:
It seems a lot of work, just to get the * operator to change.

element-by-element multiplication has the potential to cause
confusion, and, hence, errors in programming.

The reason I was doing it was not so much to have *access*
to matrix multiplication - as has been noted, Fortran now has a
function for that - but instead to have the language treat matrices consistently as mathematical objects which resemble numbers,
just as it does with complex numbers.

That would be a very natural style of writing certain types of
programs which made extensive use of matrices.

John Savard

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From gah4@21:1/5 to Quadibloc on Tue Oct 26 23:42:18 2021
On Tuesday, October 26, 2021 at 5:14:10 PM UTC-7, Quadibloc wrote:
On Tuesday, October 26, 2021 at 4:47:28 AM UTC-6, gah4 wrote:
It seems a lot of work, just to get the * operator to change.

using operator overloading for matrix multiplication versus element-by-element multiplication has the potential to cause
confusion, and, hence, errors in programming.

Octave and Matlab, two languages designed around matrix
operations, used * for the matrix multiply operator,
and .* for the element wise multiply operator.

That is slightly more obvious for an interpreted language,
and also with the matrix emphasis.

While matrix math is not rare in Fortran, it is also not the
primary use for the language.

Octave and Matlab, even more, have the / operator (multiply
by the inverse matrix), and ** (matrix exponential), again
with the ./ and .** element wise operators.

The reason I was doing it was not so much to have *access*
to matrix multiplication - as has been noted, Fortran now has a
function for that - but instead to have the language treat matrices consistently as mathematical objects which resemble numbers,
just as it does with complex numbers.

I don't know if there was any thought to that when array expressions
like such a bad idea, though.

I think I still remember, when I first learned C, how convenient it
was to have the % (modulo) operator, instead of the MOD function.
I suspect I do a lot more modulo than matrix multiplication, though.

That would be a very natural style of writing certain types of
programs which made extensive use of matrices.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Edmondo Giovannozzi@21:1/5 to All on Wed Oct 27 04:15:21 2021
Il giorno mercoledì 27 ottobre 2021 alle 08:42:20 UTC+2 gah4 ha scritto:
On Tuesday, October 26, 2021 at 5:14:10 PM UTC-7, Quadibloc wrote:
On Tuesday, October 26, 2021 at 4:47:28 AM UTC-6, gah4 wrote:
It seems a lot of work, just to get the * operator to change.

using operator overloading for matrix multiplication versus element-by-element multiplication has the potential to cause
confusion, and, hence, errors in programming.
Octave and Matlab, two languages designed around matrix
operations, used * for the matrix multiply operator,
and .* for the element wise multiply operator.

That is slightly more obvious for an interpreted language,
and also with the matrix emphasis.

While matrix math is not rare in Fortran, it is also not the
primary use for the language.

Octave and Matlab, even more, have the / operator (multiply
by the inverse matrix), and ** (matrix exponential), again
with the ./ and .** element wise operators.
The reason I was doing it was not so much to have *access*
to matrix multiplication - as has been noted, Fortran now has a
function for that - but instead to have the language treat matrices consistently as mathematical objects which resemble numbers,
just as it does with complex numbers.
I don't know if there was any thought to that when array expressions
like such a bad idea, though.

I think I still remember, when I first learned C, how convenient it
was to have the % (modulo) operator, instead of the MOD function.
I suspect I do a lot more modulo than matrix multiplication, though.
That would be a very natural style of writing certain types of
programs which made extensive use of matrices.

Python introduced the operator @ for matrix multiplication (in version 3.5 if I remember well), leaving * as an element by element multiplication.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From David Jones@21:1/5 to FortranFan on Wed Oct 27 12:30:16 2021
FortranFan wrote:

c) Starting with Fortran 90 introduced nearly 30 years ago, the
language support for the outstanding "matrix" operation you show is

A = matmul( B, C )

with the details of achieving efficient matrix multiplication
left mostly to the processors.

If there were to be a semi-serious attempt to make matrices more easily
used in Fortran, one might contemplate adding optional arguments to
MATMUL to allow optional transposition of the input arguments before
being used in the matrix-multiplication. This would then allow
avoidance of unnecessary construction of transposes.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Gary Scott@21:1/5 to Quadibloc on Wed Oct 27 08:36:17 2021
On 10/26/2021 7:14 PM, Quadibloc wrote:
On Tuesday, October 26, 2021 at 4:47:28 AM UTC-6, gah4 wrote:
It seems a lot of work, just to get the * operator to change.

using operator overloading for matrix multiplication versus element-by-element multiplication has the potential to cause
confusion, and, hence, errors in programming.

The reason I was doing it was not so much to have *access*
to matrix multiplication - as has been noted, Fortran now has a
function for that - but instead to have the language treat matrices consistently as mathematical objects which resemble numbers,
just as it does with complex numbers.

That would be a very natural style of writing certain types of
programs which made extensive use of matrices.

John Savard

I've wanted the same "natural" representation for boolean/logical

a .and. b .or. 5 .shift. 8

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Thomas Koenig@21:1/5 to David Jones on Wed Oct 27 17:07:56 2021
David Jones <dajhawkxx@nowherel.com> schrieb:
FortranFan wrote:

c) Starting with Fortran 90 introduced nearly 30 years ago, the
language support for the outstanding "matrix" operation you show is

A = matmul( B, C )

with the details of achieving efficient matrix multiplication
left mostly to the processors.

If there were to be a semi-serious attempt to make matrices more easily
used in Fortran, one might contemplate adding optional arguments to
MATMUL to allow optional transposition of the input arguments before
being used in the matrix-multiplication. This would then allow
avoidance of unnecessary construction of transposes.

Such arguments already exist, but you have to apply them to the
arguments.

If you write

A = matmul (transpose(b), c)

compilers can recognize this idiom and can

a) Use optimized loop nesting if the call to matmul is inlined

b) do something special in its own library

c) Call the relevant *GEMM routine (if the compiler does so)
and specify the TRANSA or TRANSB as needed

gfortran does a) and c) and, to a limited extent, b).

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Ron Shepard@21:1/5 to David Jones on Wed Oct 27 12:17:22 2021
On 10/27/21 7:30 AM, David Jones wrote:
FortranFan wrote:

c) Starting with Fortran 90 introduced nearly 30 years ago, the
language support for the outstanding "matrix" operation you show is

A = matmul( B, C )

with the details of achieving efficient matrix multiplication
left mostly to the processors.

If there were to be a semi-serious attempt to make matrices more easily
used in Fortran, one might contemplate adding optional arguments to
MATMUL to allow optional transposition of the input arguments before
being used in the matrix-multiplication. This would then allow
avoidance of unnecessary construction of transposes.

I remember having that same concern back when I started using f90. I
stumbled across the following. If you have a subprogram with an explicit interface with an intent(in) assumed shape matrix argument, then you can
use the loc() or c_loc() functions to look at the memory locations of
the actual and dummy arguments. What I found was that the transpose(A)
actual argument and the dummy argument memory locations are the same.
This is not required by the standard (any of them), but it was a nice
surprise to see that optimization being performed. I have looked at
several compilers over the years, and this seems to be a common
optimization.

What is happening, of course, is that the compiler is creating a new
data structure for the transposed array, and the row and column
addressing values for increments and offsets for the elements are being switched, but the actual memory is the same for the actual and the dummy arrays. Thus it is like a "shallow transpose" is being done, similar to
a shallow copy using pointer addresses rather than a deep copy that
actually references the memory. Thus, no matter how large the array
argument is, the calling sequence takes O(1) time. It is just the
metadata that is modified in the call, not the actual array data.

At the times I have looked at this, it did not also work for other declarations. Explicit shape dummy arguments don't work that way because
they must be contiguous, thus a deep transpose must be done, with copy-in/copy-out. Intent(inout) dummy arguments do not work this way
because the actual argument is considered to be an expression. However,
if you think about it, the same shallow transpose operation would also
work in this case too. For this reason, I have suggested that there
should be a new shallow transpose operator that is directly available to
the programmer that works the same way as this common optimization, but
one that would not be considered as an expression. It's only task is to manipulate the metadata for arrays. It would always be an O(1)
operation, no matter the size of the array. In addition to actual
arguments, it could also be used in expressions and on the left hand
side of expressions, the same way that the transpose operation can be
used in mathematical expressions. There is currently no way to do a
pointer assignemnt to the transpose of an array -- this would allow that
in a straightforward way. The programmer could do things like

B => shallow_transpose(A)
shallow_transpose(B) => A
B = expression...shallow_transpose(A)...expression
B = function(,..., shallow_transpose(A), ...)
call some_subroutine(..., shallow_transpose(A), ...)

and know that all of those transpose operations would be done with O(1)
effort. Currently, it is up to the compiler to decide which transpose operations are shallow (and cheap) and which ones are deep (and
expensive). This would give the programmer more control over the
low-level operations that are being performed.

BTW, I do not think there is any need to introduce an intrinsic MATRIX
type in fortran. That would introduce much complication to the language, requiring, for example, separate subprograms to be written for various combinations of array and MATRIX arguments, new conversion routines to
be defined, new semantics for assignments, pointer assignments, actual
and dummy arguments, and mixed-mode arithmetic, and so on. Consider for example, the nightmare of maintaining the LAPACK library that allows all combinations of array and MATRIX arguments. Instead, I think exposing
the array metadata to the programmer, even if it is limited to 2D arrays
in some cases, would be generally more productive.

\$.02 -Ron Shepard

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From pehache@21:1/5 to All on Fri Oct 29 00:08:15 2021
Le 26/10/2021 à 20:42, gah4 a écrit :

In the Cray days, there was much work on getting compilers to figure
out how to make good use of vector registers.

The Fortran added array expressions, which are, among others, supposed
to make it easier for vector processors. There is also FORALL, which was more
specifically added to allow for vector processors.

Wasn't FORALL rather designed for general parallelism, allowing out of
order execution of the loop ?

It was about the time of Cray, and I thought I knew that it was related
to vector processors.

Actually both parallelisation and vectorisation of loops require some indepency, so what was good for vectorisation was also good for
parallelism. FORALL is just a kind of generalized array syntax, ans
basically the array syntax allows the execution in any order (the price
to pay being the temporary arrays)

Now we have DO CONCURRENT that, as far as
I know, is related to out of order.

Yep. And is more general than array syntax, and without hidden temporary
arrays if used to replace an array syntax or a FORALL.

--
"...sois ouvert aux idées des autres pour peu qu'elles aillent dans le
même sens que les tiennes.", ST sur fr.bio.medecine

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Jos Bergervoet@21:1/5 to Ron Shepard on Fri Oct 29 19:47:56 2021
On 21/10/27 7:17 PM, Ron Shepard wrote:
On 10/27/21 7:30 AM, David Jones wrote:
FortranFan wrote:
...
...
... There is currently no way to do a
pointer assignemnt to the transpose of an array -- this would allow that
in a straightforward way. The programmer could do things like >
B => shallow_transpose(A)
shallow_transpose(B) => A
B = expression...shallow_transpose(A)...expression
B = function(,..., shallow_transpose(A), ...)
call some_subroutine(..., shallow_transpose(A), ...)

and know that all of those transpose operations would be done with O(1) effort.

This came up some time ago in the thread "Mathematics of arrays".
The minimum "hacking" solution (dangerous!) I found in gfortran
was the following: (just to play with it, not really useful..)
To be included in a module or "contains" of main program
!

function shallow_transpose(p1) result(p2)
real(wp), target, intent(in) :: p1(:,:)
real(wp), pointer :: p2(:,:)

type dummy
real(wp), pointer :: p(:,:)
end type
type(dummy) :: d
integer :: ar(22) ! to fetch descriptor.. yes, that big!

d%p => p1; ! look at pointer p1
ar = transfer(d, ar) ! fetch descriptor in ar(:)

ar(11:22) = [ar(17:22), ar(11:16)] ! swap index information

d = transfer(ar, d) ! put back in the descriptor
p2 => d%p ! copy pointer from wrapper
end

--
Jos

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Robin Vowels@21:1/5 to Quadibloc on Fri Oct 29 17:44:28 2021
On Tuesday, October 26, 2021 at 5:15:12 PM UTC+11, Quadibloc wrote:
I will have to admit that I'm not fully conversant
with the latest developments in Fortran.

But it is my understanding that starting around the
time of Fortran 90, constructs have been added to
the language so that programs can be written to
efficiently make use of the vector arithmetic features
of computers like the Cray I, and its successors and
imitators.

These days, except for the SX-Aurora TSUBASA from
NEC, most vector arithmetic on computers belongs to
a line of development that started from the AN/FSQ-7
for SAGE rather than the one that went through the
Cray I.

Be that as it may...

So, IIRC, one can write

A = B * C

in Fortran today and if A, B, and C are arrays having the
same dimensions (or if one of B or C is a scalar) it will
do an element-by-element multiplication, much as would
have been done by APL, without the need for a DO loop.
.
That has been available in PL/I since it was introduced in 1966.
It has been the mainstay of array operations on the Pilot ACE
from 1955 (and probably earlier) and on DEUCE from 1955,
as GIP. The specific program was LZ61B and its successor LZ61BM/1.
.
In IBM's PL/I-F (1966), the elements are accessed using a single loop (regardless of whether the operands are vectors or matrices),
without the need to manipulate subscripts, and thus allows
for highly efficient array operations.
.
I remember that back in the 1970s I mentioned to someone
that I thought this sort of thing would be a nice feature to
have, even without vector hardware. The fellow replied that
I was being too slavish in following APL, matrix multiplication
would be better.
.
for at least 4 years.
.
Well, I agreed that matrix multiplication was useful in
scientific and engineering computation.

to the construct of an array in Fortran.

Of course, one can add a built-in subroutine for matrix
multiplication.
.
That was done in 1990 with Fortran 90.
.
But more recently, I thought that there might be a better
solution; and now on this newsgroup, I've thought to
present it more widely.

Fortran has a COMPLEX data type.

Why can't it have a MATRIX data type?
.
Why?
.
So there would be a MATRIX statement, having much the
same syntax as a DIMENSION statement. It could be used
to create one-dimensional vectors and matrices that were
not square.
.
That can be done in FORTRAN (pre 1990) and from Fortran 90.
.
If A, B, and C are declared as matrices (or one of B or C is
a scalar) then

A = B * C

would perform matrix multiplication as expected...
.
Might that not be confusing?
An in any case, Fortran post 1990 already has a built-in subroutine
for matrix multiplication.
.
if one of
B or C is a one-dimensional matrix, it would be treated as
either a row or column vector, whichever one would make
sense. (Unless that isn't practical; declaring 1 by n and n by 1
matrices instead to make it explicit is certainly an option as
well.)

Since a matrix is an a x b grid of floating-point numbers
.
Matrices can be composed of integers.
.
combined into an object with specific mathematical
properties, while it would be possible to create an array
of matrices, it would *not* be possible to create a matrix
the elements of which are arrays.

For much the same reason that you can't make a
complex number the real and imaginary parts of which
are arrays of the same shape.

But, yes, it should be allowed to create a matrix the elements
of which are complex numbers.
.
That can be done in FORTRAN, and has been the case for
some 60+ years.
.
(Or quaternions, for that
matter!)

Oh, and by the way: yes, it should be possible to take a (4,4)
matrix the elements of which are REAL*16 and EQUIVALENCE
it to a (4,4) array the elements of which are REAL*16
.
Why on earth would you want to do that?
.
with the
expected results of the same-numbered elements perfectly
corresponding. (This might even be a reason not to deprecate
EQUIVALENCE...)

Is this not a useful thing to include in Fortran?
.
No.
.
Maybe someone else has even suggested it already?

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Ron Shepard@21:1/5 to Jos Bergervoet on Sat Oct 30 15:35:46 2021
On 10/29/21 12:47 PM, Jos Bergervoet wrote:
On 21/10/27 7:17 PM, Ron Shepard wrote:
On 10/27/21 7:30 AM, David Jones wrote:
FortranFan wrote:
...
...
... There is currently no way to do a pointer assignemnt to the
transpose of an array -- this would allow that in a straightforward
way. The programmer could do things like >
B => shallow_transpose(A)
shallow_transpose(B) => A
B = expression...shallow_transpose(A)...expression
B = function(,..., shallow_transpose(A), ...)
call some_subroutine(..., shallow_transpose(A), ...)

and know that all of those transpose operations would be done with
O(1) effort.

This came up some time ago in the thread "Mathematics of arrays".
The minimum "hacking" solution (dangerous!) I found in gfortran
was the following: (just to play with it, not really useful..)
To be included in a module or "contains" of main program
!

function shallow_transpose(p1)     result(p2)
real(wp), target, intent(in)  :: p1(:,:)
real(wp), pointer             :: p2(:,:)

type dummy
real(wp), pointer :: p(:,:)
end type
type(dummy) :: d
integer     :: ar(22)    ! to fetch descriptor.. yes, that big!

d%p => p1;                           ! look at pointer p1
ar = transfer(d, ar)                 ! fetch descriptor in ar(:)

ar(11:22) = [ar(17:22), ar(11:16)]   ! swap index information

d = transfer(ar, d)                  ! put back in the descriptor
p2 => d%p                            ! copy pointer from wrapper
end

Neat trick. This doesn't allow the function on the lhs of an expression,
but I think it is close for all of the other cases (apart from expression/intent(inout) mismatches). Unfortunately, it does not seem to
work in gfortran 10. That is not surprising because this would be highly nonportable code. But just from looking at what it does, that is exactly
what I have in mind for the shallow_transpose() function, but in a
standard and portable way.

\$.02 -Ron Shepard

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Steve Lionel@21:1/5 to Ron Shepard on Sun Oct 31 10:25:57 2021
On 10/30/2021 4:35 PM, Ron Shepard wrote:
But just from looking at what it does, that is exactly what I have in
mind for the shallow_transpose() function, but in a standard and
portable way.

Um, no. This is the absolute antithesis of "standard and portable". It
relies on a specific compiler's implementation of pointers, something
that can change even between versions of a compiler. It will not work
with other compilers.

--
Steve Lionel
ISO/IEC JTC1/SC22/WG5 (Fortran) Convenor
Retired Intel Fortran developer/support
Email: firstname at firstnamelastname dot com
Blog: https://stevelionel.com/drfortran
WG5: https://wg5-fortran.org

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Quadibloc@21:1/5 to Robin Vowels on Sun Oct 31 19:56:10 2021
On Friday, October 29, 2021 at 6:44:30 PM UTC-6, Robin Vowels wrote:

for at least 4 years.

I wasn't claiming originality. It was a good idea in PL/I,
and so it would be good for Fortran to have it too.

Just as later on, in Fortran 77, block-structured IF statements

John Savard

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From gah4@21:1/5 to Quadibloc on Mon Nov 1 09:10:01 2021
On Sunday, October 31, 2021 at 7:56:12 PM UTC-7, Quadibloc wrote:
On Friday, October 29, 2021 at 6:44:30 PM UTC-6, Robin Vowels wrote:

for at least 4 years.

I wasn't claiming originality. It was a good idea in PL/I,
and so it would be good for Fortran to have it too.

Just as later on, in Fortran 77, block-structured IF statements

One interesting difference, though.

Fortran array assignments require the "as if" the whole right hand side
is evaluated before the left side is changed. If the compiler can't be sure,
it has to use a temporary array.

PL/I array assignment is done element by element. If an array element
changes, that change is used in later assignments. That was before
vector processors, maybe before people thought about building them.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Robin Vowels@21:1/5 to All on Mon Nov 1 17:51:03 2021
On Tuesday, November 2, 2021 at 3:10:03 AM UTC+11, gah4 wrote:
On Sunday, October 31, 2021 at 7:56:12 PM UTC-7, Quadibloc wrote:
On Friday, October 29, 2021 at 6:44:30 PM UTC-6, Robin Vowels wrote:

for at least 4 years.

I wasn't claiming originality. It was a good idea in PL/I,
and so it would be good for Fortran to have it too.

Just as later on, in Fortran 77, block-structured IF statements
.
Algol also had block-structured IF statements.
.
One interesting difference, though.

Fortran array assignments require the "as if" the whole right hand side
is evaluated before the left side is changed. If the compiler can't be sure, it has to use a temporary array.

PL/I array assignment is done element by element.

As is done in Fortran. However, in PL/I

If an array element changes,
.
No. If an array element on the LHS of the assignment
is also used on the RHS,
.
that change is used in later assignments. That was before
vector processors, maybe before people thought about building them.
.
Vector processors were in operation long before PL/I (1966).
The Pilot ACE (1951) did vector arithmetic and vector moves.
So did the DEUCE (1955) which had a similar architecture to Pilot ACE.
I suspect (but do not know for sure) that ACE (c. 1955) also had these
array features.

The rule in PL/I is that in an array assignment,
it is treated as if the loop were written out as one (or more nested
DO loops, as required), with the subscripts being varied in the
corresponding row-major order.
.
However, that rule* doesn't stop a PL/I compiler from reducing
nested loops to a SINGLE loop (as was done in IBM's PL/I (F)
compiler and the DEUCE PL/I compiler) to produce a highly-
optimised loop.
_________
* It is necessary only for the programmer's understanding of the
order in which the computations are carried out.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Jos Bergervoet@21:1/5 to Steve Lionel on Tue Nov 2 09:04:46 2021
On 21/10/31 3:25 PM, Steve Lionel wrote:
On 10/30/2021 4:35 PM, Ron Shepard wrote:
But just from looking at what it does, that is exactly what I have in
mind for the shallow_transpose() function, but in a standard and
portable way.

Um, no. This is the absolute antithesis of "standard and portable".

perhaps have used equivalence, instead of the transfer function. But
I admit this was close to what is conceivably possible (in terms of approximating the antithesis, I mean..)

It
relies on a specific compiler's implementation of pointers, something
that can change even between versions of a compiler.

That is true, but it did not rely on possible alignment problems, which equivalence could perhaps have introduced (with some care..) That could
also easily have provided out-of-bounds addressing, which was ow not
used at all!

It will not work
with other compilers.

OK, it did at least satisfy that basic requirement. :^)

[But seriously, I think Ron is right that it could be portable if each
version of a compiler would do this in the way appropriate for, well..
that version of the compiler! After all that is how all code generation
is kept portable, even though "a specific compiler's implementation"
that "can change between versions" is used for almost everything.]

--
Jos

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Ron Shepard@21:1/5 to Jos Bergervoet on Tue Nov 2 10:35:27 2021
On 11/2/21 3:04 AM, Jos Bergervoet wrote:
On 21/10/31 3:25 PM, Steve Lionel wrote:
On 10/30/2021 4:35 PM, Ron Shepard wrote:
But just from looking at what it does, that is exactly what I have in
mind for the shallow_transpose() function, but in a standard and
portable way.

Um, no. This is the absolute antithesis of "standard and portable".

perhaps have used equivalence, instead of the transfer function. But
I admit this was close to what is conceivably possible (in terms of approximating the antithesis, I mean..)

It relies on a specific compiler's implementation of pointers,
something that can change even between versions of a compiler.

That is true, but it did not rely on possible alignment problems, which equivalence could perhaps have introduced (with some care..) That could
also easily have provided out-of-bounds addressing, which was ow not
used at all!

It will not work with other compilers.

OK, it did at least satisfy that basic requirement. :^)

[But seriously, I think Ron is right that it could be portable if each version of a compiler would do this in the way appropriate for, well..
that version of the compiler! After all that is how all code generation
is kept portable, even though "a specific compiler's implementation"
that "can change between versions" is used for almost everything.]

I guess my wording in the previous post must have been confusing. I
certainly was not saying that that implementation of the
shallow_transpose() function was portable. In fact, I noted that it did
not work with an older version of gfortran, so it was not even portable
between different versions of the same compiler. But that implementation
does demonstrate that such an intrinsic function could provide that functionality in a standard way that would be portable to the
programmer. The low-level implementation of that intrinsic function
would of course depend on the underlying metadata that the compiler uses
to describe arrays.

\$.02 -Ron Shepard

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From gah4@21:1/5 to Ron Shepard on Tue Nov 2 17:42:54 2021
On Tuesday, November 2, 2021 at 8:35:31 AM UTC-7, Ron Shepard wrote:

(snip)
I guess my wording in the previous post must have been confusing. I
certainly was not saying that that implementation of the
shallow_transpose() function was portable. In fact, I noted that it did
not work with an older version of gfortran, so it was not even portable between different versions of the same compiler. But that implementation
does demonstrate that such an intrinsic function could provide that functionality in a standard way that would be portable to the
programmer. The low-level implementation of that intrinsic function
would of course depend on the underlying metadata that the compiler uses
to describe arrays.

I wasn't especially trying to follow it, but I thought I remembered
something about standardizing the descriptors, such that one could
do things like this.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Steve Lionel@21:1/5 to All on Thu Nov 4 12:14:17 2021
On 11/2/2021 8:42 PM, gah4 wrote:
I wasn't especially trying to follow it, but I thought I remembered
something about standardizing the descriptors, such that one could
do things like this.

"C descriptors" are somewhat standardized, but they're not visible to
Fortran code and the standard permits layout differences. C code can
access the descriptors through macros and structure declarations in ISO_Fortran_binding.h. Some compilers might use C descriptors to
represent pointers in Fortran code, but there's no requirement that this
be done. Even the specific layout of a C descriptor's bounds section is implementation-dependent.

--
Steve Lionel
ISO/IEC JTC1/SC22/WG5 (Fortran) Convenor
Retired Intel Fortran developer/support
Email: firstname at firstnamelastname dot com
Blog: https://stevelionel.com/drfortran
WG5: https://wg5-fortran.org

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Lynn McGuire@21:1/5 to Quadibloc on Thu Nov 4 13:38:56 2021
On 10/26/2021 1:15 AM, Quadibloc wrote:
I will have to admit that I'm not fully conversant
with the latest developments in Fortran.

But it is my understanding that starting around the
time of Fortran 90, constructs have been added to
the language so that programs can be written to
efficiently make use of the vector arithmetic features
of computers like the Cray I, and its successors and
imitators.

These days, except for the SX-Aurora TSUBASA from
NEC, most vector arithmetic on computers belongs to
a line of development that started from the AN/FSQ-7
for SAGE rather than the one that went through the
Cray I.

Be that as it may...

So, IIRC, one can write

A = B * C

in Fortran today and if A, B, and C are arrays having the
same dimensions (or if one of B or C is a scalar) it will
do an element-by-element multiplication, much as would
have been done by APL, without the need for a DO loop.

I remember that back in the 1970s I mentioned to someone
that I thought this sort of thing would be a nice feature to
have, even without vector hardware. The fellow replied that
I was being too slavish in following APL, matrix multiplication
would be better.

Well, I agreed that matrix multiplication was useful in
scientific and engineering computation.

to the construct of an array in Fortran.

Of course, one can add a built-in subroutine for matrix
multiplication.

But more recently, I thought that there might be a better
solution; and now on this newsgroup, I've thought to
present it more widely.

Fortran has a COMPLEX data type.

Why can't it have a MATRIX data type?

So there would be a MATRIX statement, having much the
same syntax as a DIMENSION statement. It could be used
to create one-dimensional vectors and matrices that were
not square.

If A, B, and C are declared as matrices (or one of B or C is
a scalar) then

A = B * C

would perform matrix multiplication as expected... if one of
B or C is a one-dimensional matrix, it would be treated as
either a row or column vector, whichever one would make
sense. (Unless that isn't practical; declaring 1 by n and n by 1
matrices instead to make it explicit is certainly an option as
well.)

Since a matrix is an a x b grid of floating-point numbers
combined into an object with specific mathematical
properties, while it would be possible to create an array
of matrices, it would *not* be possible to create a matrix
the elements of which are arrays.

For much the same reason that you can't make a
complex number the real and imaginary parts of which
are arrays of the same shape.

But, yes, it should be allowed to create a matrix the elements
of which are complex numbers. (Or quaternions, for that
matter!)

Oh, and by the way: yes, it should be possible to take a (4,4)
matrix the elements of which are REAL*16 and EQUIVALENCE
it to a (4,4) array the elements of which are REAL*16 with the
expected results of the same-numbered elements perfectly
corresponding. (This might even be a reason not to deprecate
EQUIVALENCE...)

Is this not a useful thing to include in Fortran?

Maybe someone else has even suggested it already?

John Savard

Hi John,

When was the last time you wrote anything in Fortran ?

I write software in Fortran almost every day. Or C++. Or if I am
lucky, both. I could never use a matrix object package in our software
since the computations are so long and usually dependent upon the
previous computations. We do use some of the Harwell routines, heavily modified.

Lynn

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Quadibloc@21:1/5 to Lynn McGuire on Thu Nov 4 22:15:50 2021
On Thursday, November 4, 2021 at 12:39:03 PM UTC-6, Lynn McGuire wrote:

When was the last time you wrote anything in Fortran ?

I must admit, quite a while ago.

I write software in Fortran almost every day. Or C++. Or if I am
lucky, both. I could never use a matrix object package in our software
since the computations are so long and usually dependent upon the
previous computations. We do use some of the Harwell routines, heavily modified.

I know that Harwell is in Britain, and has something to do with atomic energy...

http://www.chilton-computing.org.uk/acl/society/impact/hsl.htm

on a site I went to for nice pictures of the front panel of the IBM 370/195 explains what the subroutine library is.

And it's the Atomic Energy Research Establishment that was known as Harwell, from the place it was located; it has now been renamed the Harwell Science and Innovation Campus.

John Savard

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From gah4@21:1/5 to Steve Lionel on Thu Nov 4 23:56:03 2021
On Thursday, November 4, 2021 at 9:14:21 AM UTC-7, Steve Lionel wrote:
On 11/2/2021 8:42 PM, gah4 wrote:
I wasn't especially trying to follow it, but I thought I remembered something about standardizing the descriptors, such that one could
do things like this.

"C descriptors" are somewhat standardized, but they're not visible to
Fortran code and the standard permits layout differences. C code can
access the descriptors through macros and structure declarations in ISO_Fortran_binding.h. Some compilers might use C descriptors to
represent pointers in Fortran code, but there's no requirement that this
be done. Even the specific layout of a C descriptor's bounds section is implementation-dependent.

So with some strange mix of Fortran and C, it should be possible to
do what was wanted. Maybe not so convenient, though.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Edmondo Giovannozzi@21:1/5 to All on Fri Nov 5 03:51:40 2021
Il giorno venerdì 5 novembre 2021 alle 07:56:05 UTC+1 gah4 ha scritto:
On Thursday, November 4, 2021 at 9:14:21 AM UTC-7, Steve Lionel wrote:
On 11/2/2021 8:42 PM, gah4 wrote:
I wasn't especially trying to follow it, but I thought I remembered something about standardizing the descriptors, such that one could
do things like this.

"C descriptors" are somewhat standardized, but they're not visible to Fortran code and the standard permits layout differences. C code can access the descriptors through macros and structure declarations in ISO_Fortran_binding.h. Some compilers might use C descriptors to
represent pointers in Fortran code, but there's no requirement that this be done. Even the specific layout of a C descriptor's bounds section is implementation-dependent.
So with some strange mix of Fortran and C, it should be possible to
do what was wanted. Maybe not so convenient, though.

I don't think that a manipulation of the descriptor should be allowed except when using pointers.
If I write something like:

B = shallowTranspose(A)

Now B and A are pointing to the same memory and now one have aliasing without the processor really knowing. While one should have something like:
..., target :: A
..., pointer :: B

B => shallowTranspose(A)

Or whatever manipulation of the descriptor one want to allow. Now the standard Fortran rules apply and the processor know that the two arrays may share the same memory.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Steve Lionel@21:1/5 to All on Fri Nov 5 10:40:26 2021
On 11/5/2021 2:56 AM, gah4 wrote:
"C descriptors" are somewhat standardized, but they're not visible to
Fortran code and the standard permits layout differences. C code can
access the descriptors through macros and structure declarations in
ISO_Fortran_binding.h. Some compilers might use C descriptors to
represent pointers in Fortran code, but there's no requirement that this
be done. Even the specific layout of a C descriptor's bounds section is
implementation-dependent.
So with some strange mix of Fortran and C, it should be possible to
do what was wanted. Maybe not so convenient, though.

You can't use the CFI_ functions to remap a Fortran pointer, but you
could create a new pointer with dimensions and bounds as desired and
pass that back to Fortran. It would have to be done as a sort of
callback, because the descriptor structure would have to be created in
the C code. It would be awkward.

--
Steve Lionel
ISO/IEC JTC1/SC22/WG5 (Fortran) Convenor
Retired Intel Fortran developer/support
Email: firstname at firstnamelastname dot com
Blog: https://stevelionel.com/drfortran
WG5: https://wg5-fortran.org

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Ron Shepard@21:1/5 to Edmondo Giovannozzi on Fri Nov 5 19:19:49 2021
On 11/5/21 5:51 AM, Edmondo Giovannozzi wrote:
Il giorno venerdì 5 novembre 2021 alle 07:56:05 UTC+1 gah4 ha scritto:
On Thursday, November 4, 2021 at 9:14:21 AM UTC-7, Steve Lionel wrote:
On 11/2/2021 8:42 PM, gah4 wrote:
I wasn't especially trying to follow it, but I thought I remembered
something about standardizing the descriptors, such that one could
do things like this.

"C descriptors" are somewhat standardized, but they're not visible to
Fortran code and the standard permits layout differences. C code can
access the descriptors through macros and structure declarations in
ISO_Fortran_binding.h. Some compilers might use C descriptors to
represent pointers in Fortran code, but there's no requirement that this >>> be done. Even the specific layout of a C descriptor's bounds section is
implementation-dependent.
So with some strange mix of Fortran and C, it should be possible to
do what was wanted. Maybe not so convenient, though.

I don't think that a manipulation of the descriptor should be allowed except when using pointers.
If I write something like:

B = shallowTranspose(A)

Now B and A are pointing to the same memory and now one have aliasing without the processor really knowing. While one should have something like:
..., target :: A
..., pointer :: B

B => shallowTranspose(A)

Or whatever manipulation of the descriptor one want to allow. Now the standard Fortran rules apply and the processor know that the two arrays may share the same memory.

Yes, pointer assignment was exactly what was suggested in the earlier
post, along with ability to use shallow_transpose() in expressions and
as an actual argument to a subprogram. In that latter case, it would be
the programmer's responsibility to avoid aliases with other arguments --
I don't think the compiler could know and check for those in all cases.

I also suggested the possibility of using the function on the left hand
side of expressions.

shallow_transpose(B) = <some expression>

That would be beyond what is currently allowed in fortran, but it would
match what is allowed in formal mathematical expressions. In terms of
low-level operations, that would require the allocation of memory to
evaluate the expression, association of that memory to the pointer, and adjusting the pointer metadata to account for the transpose. That seems complicated, knowing when that memory can be deallocated by the compiler
and what other things can be done with it (e.g. subsequent pointer
assignments, etc.). So maybe that isn't really practical, despite its
apparent usefulness.

In any case, in the end the shallow_transpose() would be an intrinsic
function defined as part of the fortran standard. The programmer would
not need to fiddle around with any descriptors using C interop or
anything like that.

\$.02 -Ron Shepard

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Jos Bergervoet@21:1/5 to Ron Shepard on Mon Nov 8 22:58:47 2021
On 21/11/06 1:19 AM, Ron Shepard wrote:
On 11/5/21 5:51 AM, Edmondo Giovannozzi wrote:
Il giorno venerdì 5 novembre 2021 alle 07:56:05 UTC+1 gah4 ha scritto:
On Thursday, November 4, 2021 at 9:14:21 AM UTC-7, Steve Lionel wrote:
On 11/2/2021 8:42 PM, gah4 wrote:
I wasn't especially trying to follow it, but I thought I remembered
something about standardizing the descriptors, such that one could
do things like this.

"C descriptors" are somewhat standardized, but they're not visible to
Fortran code and the standard permits layout differences. C code can
access the descriptors through macros and structure declarations in
ISO_Fortran_binding.h. Some compilers might use C descriptors to
represent pointers in Fortran code, but there's no requirement that
this
be done. Even the specific layout of a C descriptor's bounds section is >>>> implementation-dependent.
So with some strange mix of Fortran and C, it should be possible to
do what was wanted. Maybe not so convenient, though.

I don't think that a manipulation of the descriptor should be allowed
except when using pointers.
If I write something like:

B = shallowTranspose(A)

Now B and A are pointing to the same memory and now one have aliasing
without the processor really knowing. While one should have something
like:
..., target :: A
..., pointer :: B

B => shallowTranspose(A)

Or whatever manipulation of the descriptor one want to allow. Now the
standard Fortran rules apply and the processor know that the two
arrays may share the same memory.

Yes, pointer assignment was exactly what was suggested in the earlier
post, along with ability to use shallow_transpose() in expressions and
as an actual argument to a subprogram. In that latter case, it would be
the programmer's responsibility to avoid aliases with other arguments --
I don't think the compiler could know and check for those in all cases.

I also suggested the possibility of using the function on the left hand
side of expressions.

shallow_transpose(B) = <some expression>

That would be beyond what is currently allowed in fortran,

Why? The lhs is allowed to be a pointer (assuming that it
has been allocated the correct array size!) And if B has
been allocated the correct transposed array size then the
result of shallow_transpose(B) has the correct size, so
the assignement is no problem (and does indeed work nicely

but it would
match what is allowed in formal mathematical expressions. In terms of low-level operations, that would require the allocation of memory to
evaluate the expression,

That is needed regardless of what the lhs is..

association of that memory to the pointer, and

That's the programmers responsibiliy, not the compilers
(we are talking about pointers, not allocatable arrays).

That's what the function shallow_transpose does, which (at
this point of the thread) we assume to exist of course..

That seems
complicated, knowing when that memory can be deallocated by the compiler
and what other things can be done with it (e.g. subsequent pointer assignments, etc.).

I see no difference between the situation after:
B = <some expression>
and
shallow_transpose(B) = <some expression>

At least not for what can be done with the memory region
involved. To me it would seem to be more tricky to have
B(i:j, n:m) = <some expression>
where only a slice of the total allocated memory of B is
being used.

But why would anything have to be "deallocated by the
compiler" in these examples? Again: we aren't talking
about allocatable arrays, but just pointers, so the
burden is on the programmer.

So maybe that isn't really practical, despite its
apparent usefulness.

In any case, in the end the shallow_transpose() would be an intrinsic function defined as part of the fortran standard.

I do think it might be useful, but the name is too
long (it's even worse than dot_product!)

The programmer would
not need to fiddle around with any descriptors using C interop or
anything like that.

Of course. But can it be generalized? For a two-
dimensional array we could have something like:

B(:,:, iorder=[2,1]) = <some expression>

And for a higher dimensional case for instance:

B(:,:,2:3,:, iorder=[2,1,4,3]) = <some expression>

where "iorder" just specifies a change in indexing
order! (Admittedly it would be error prone because of
the added complexity to what index slices now mean..)

--
Jos

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From gah4@21:1/5 to Jos Bergervoet on Mon Nov 8 19:44:02 2021
On Monday, November 8, 2021 at 2:00:04 PM UTC-8, Jos Bergervoet wrote:

(snip)

I see no difference between the situation after:
B = <some expression>
and
shallow_transpose(B) = <some expression>

Other than it isn't usual to use functions on the left side of an assignment.

PL/I has pseudo-variables, which do look like functions, and can appear
on the left side. In the case of a function that is its own inverse, though
I don't see the need.

PL/I has the REAL and IMAG functions to extract the real and imaginary
parts of a complex value, and COMPLEX to make a complex value from
real and imaginary parts.

Then there are the REAL and IMAG pseudo-variables to change the real
or imaginary part of a complex value, and COMPLEX to extract both.

X=REAL(Z);
Y=IMAG(Z);

can be written as:

COMPLEX(X,Y) = Z;

IMAG(Z) = Y;

to change only the imaginary part. Much more obvious than:

Z = COMPLEX(REAL(Z), Y);

and even:

DO IMAG(Z) = 1 TO 100;

where it is both a function and pseudo-variable.

There is also SUBSTR(), as a function to extract a substring, and pseudo-variable to change a substring in a string.

In any case, interpreted languages commonly have the shallow
transpose, where all of memory allocation is not visible to the
programmer.

At least not for what can be done with the memory region
involved. To me it would seem to be more tricky to have
B(i:j, n:m) = <some expression>
where only a slice of the total allocated memory of B is
being used.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From Robin Vowels@21:1/5 to All on Mon Nov 8 22:39:31 2021
On Tuesday, November 9, 2021 at 2:44:03 PM UTC+11, gah4 wrote:
On Monday, November 8, 2021 at 2:00:04 PM UTC-8, Jos Bergervoet wrote:

(snip)
I see no difference between the situation after:
B = <some expression>
and
shallow_transpose(B) = <some expression>
Other than it isn't usual to use functions on the left side of an assignment.

PL/I has pseudo-variables, which do look like functions,
.
Most of them do not.
Most are simple names, such as ONSOURCE, ONCHAR, etc.
One pseudo-variable that is likely to be seen is SUBSTR,
which can replace zero or more characters in an existing string. E.g., SUBSTR(S, 4, 7) = 'abcdefg';
which replaces 7 characters in string S, commencing with the fourth
character.
Fortran has a comparable facility.
.
and can appear
on the left side. In the case of a function that is its own inverse, though
I don't see the need.

PL/I has the REAL and IMAG functions to extract the real and imaginary
parts of a complex value, and COMPLEX to make a complex value from
real and imaginary parts.
.
COMPLEX as a pseudo-variable has long been dropped from the language.
.
Then there are the REAL and IMAG pseudo-variables to change the real
or imaginary part of a complex value, and COMPLEX to extract both.

X=REAL(Z);
Y=IMAG(Z);

can be written as:

COMPLEX(X,Y) = Z;
.
Not any more.
.
IMAG(Z) = Y;

to change only the imaginary part. Much more obvious than:

Z = COMPLEX(REAL(Z), Y);

and even:

DO IMAG(Z) = 1 TO 100;

where it is both a function and pseudo-variable.

There is also SUBSTR(), as a function to extract a substring, and pseudo-variable to change a substring in a string.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)