Question: **Is it considered a good practice to use automatic arrays and intrinsic array operations** (`matmul`

, `transpose`

, etc)?

The reasonable answer, I guess and I hope, is **yes, it is a good practice**.

For automatic arrays, let us focus on their usage as local variables. I do not intend to discuss what kind of arrays should be used as dummy variables. You may see Fortran Best Practices for such discussions.

Now let me explain why I have this question.

For me, **one of the core strengths of Fortran is its native support for matrices and matrix/vector operations**. By â€śnative support for matricesâ€ť, I mean particularly automatic arrays â€” to get a matrix, you do not need to do anything more than declare it. As for matrix/vector operations, I mean intrinsic array operations such as `matmul`

, `transpose`

, `dot_product`

, and broadcasting.

Thanks to such support, we can write scientific computing code in a way that is very similar to how the algorithms are formulated in mathematics. Without such support, a language should not be called a â€śFormula Translatorâ€ť.

Take matrix/vector operations as an example. Without intrinsic support for them, we have to either re-invent these wheels ourselves, or code them as loops whenever such operations appear in an algorithm. However, as a mathematician, I can tell that **mathematicians rarely consider matrix/vector operations as loops** â€” loops are only how we calculate them, but not how we define/understand them. Loops are used only when we describe truly iterative algorithms, e.g., CG, steepest descent, â€¦

Not surprisingly, most other popular languages in scientific computing are also strong at supporting matrices and matrix/vector operations. My most familiar language is MATLAB â€” if you call it a language. It supports almost all matrix/vector operations you can name. A well-knowledged MATLAB programmer can code a sophisticated algorithm without using any loop but only using matrix/vector operations.

Indeed, as a common principle in MATLAB programming, you should **never** code a matrix/vector operation in loops. **It is a bad practice to re-invent any tool that the language intrinsically provides.** You may argue that this is because MATLAB is slow at loops as a scripting language. But I understand it from a more philosophical perspective: **why does the language provide the functionality if we are not encouraged to use it?**

So, it seems that the only reasonable answer to my question is â€śyes, it is a good practice; do use automatic arrays and intrinsic array operations as long as you canâ€ť.

However, is it really the case?

Take automatic arrays as an example. These arrays are often allocated on the stack. In that case, a double-precision matrix of modest size will lead you to a stack overflow (I know that â€śstackâ€ť is not a specification of the language standard but a choice of the compiler implementations, but should we expect all users to know this? Anyway, **such knowledge does not relieve the hurt of a stack overflow**). On macOS, this size is about 256*256. Of course, we can tell the compiler not to save automatic arrays at the stack, or do so only for small arrays, or simply use allocatable arrays instead of automatic ones. However, let us acknowledge that **it is an inconvenience, which essentially does not exist in other popular languages that have native support for matrices**. Imagine that you are using a piece of sophisticated code that contains automatic arrays, then you will encounter disasters if you do not specify the correct compiler options. Surely, the author of the code should document the compiler options clearly, and the user should read the documentation carefully. But **why donâ€™t users of MATLAB/Python/Julia/R need to pay such attention when using intrinsic functionalities of the languages?**

Now take `matmul`

as another example. Consider the code

```
A = B + matmul(C, D)
```

A couple of problems may occur if you code Fortran in that way. First, a temporary array will likely be created for `matmul(C, D)`

; in addition to the performance penalty, this temporary array may again be allocated on the stack, which may again lead to a mysterious segmentation fault if you are unlucky even when the matrices are not so big (see my reply below for a simple example of such a segmentation fault). Second, according to my reading and testing, the performance of `matmul`

~~is sometimes inferior to a naive triple loop~~ (This is wrong, thank @JohnCampbell for pointing it out, and my apologies for misinformation), let alone optimized subroutines. Surely, if you are knowledged enough, you can choose the correct compiler option and link to the correct library so that `matmul`

is translated to, e.g., `gemm`

. But, again, **why donâ€™t users of MATLAB/Python/Julia/R need to pay such attention when using intrinsic functionalities of the languages?**

Indeed, in MATLAB, `A = B + C*D`

is almost the only reasonable way to code the matrix calculation illustrated above. Any code involving loops is essentially wrong unless these matrices have special structures. In addition, MATLAB users normally (but not always) do not need to worry about stack overflow or which library to link to. They simply rely on the intrinsic functionality of the language.

With these questions in mind, I checked the code of `stdlib`

. Very interestingly, it seems that the contributors have diverse opinions.

Example 1: stdlib_stats_corr.fypp, Lines 119â€“226:

```
do i = 1, size(x, 2)
center(:, i) = x(:, i) - mean_
end do
#:if t1[0] == 'r'
res = matmul( center, transpose(center))
#:else
res = matmul( center, transpose(conjg(center)))
#:endif
```

We see `matmul`

and `transpose`

, which I like very much. However, we also see a loop that can be replaced with broadcasting. There must be some careful consideration behind the code, but it is simply interesting to see code of two different â€śstylesâ€ť side by side.

Example 2: stdlib_linalg.fypp, `diag`

versus `eye`

:

```
! Part of the code for `diag` (vector to matrix):
#:for k1, t1 in RCI_KINDS_TYPES
module function diag_${t1[0]}$${k1}$(v) result(res)
${t1}$, intent(in) :: v(:)
${t1}$ :: res(size(v),size(v))
end function diag_${t1[0]}$${k1}$
#:endfor
```

```
! Code for `eye`
pure function eye(dim1, dim2) result(result)
integer, intent(in) :: dim1
integer, intent(in), optional :: dim2
integer(int8), allocatable :: result(:, :)
integer :: dim2_
integer :: i
dim2_ = optval(dim2, dim1)
allocate(result(dim1, dim2_))
result = 0_int8
do i = 1, min(dim1, dim2_)
result(i, i) = 1_int8
end do
end function eye
```

We see that `diag`

declares its result as an automatic array, while `eye`

declares it as an allocatable array. Why are they different?

Before the end, let me stress that **the native support for matrices and matrix/vector operations is a core strength of Fortran**. It would be extremely wired if we are indifferent about its performance and even its occasional failure under default settings.

Sorry for this long post. What is your opinion on my question (see the very top line)? Thank you very much for your attention.

Related posts (to the authors of the posts, **I apologize** for bugging you with an @, but I just thought you might be interested in this topic since you discussed similar things before):

- Reply to â€śIs the Fortran Standard crucial for the future of the Fortran language?â€ť (see the comment on â€śDecent array supportâ€ť) by @Pap
- Mapping matrix & vector arithmetic to BLAS calls by @nshaffer
- Automatic vs. allocatable arrays for function results by @milancurcic and replies by @Beliavsky @certik @zoziha @themos
- Reply to â€śBest practice of allocating memory in Fortranâ€ť by @JohnCampbell
- When to use allocatable arrays? by @fortran4r
- Why a function returns an array is much slower than a subroutine returns an array? (real MWE included) by @CRquantum
- Reply to â€śWhy abandon Fortran for Linear Algebra?â€ť by @Beliavsky
- Julia: Fast as Fortran, Beautiful as Python by @certik
- Fortran: An ideal language for writing â€śTemplatesâ€ť of numerical algorithms (?) by zaikunzhang

Back in the good old days â€“ the â€śGolden Eraâ€ť of computers, it was easy to separate the men from the boys (sometimes called â€śReal Menâ€ť and â€śQuiche Eatersâ€ť in the literature). During this period, the Real Men were the ones that understood computer programming, and the Quiche Eaters were the ones that didnâ€™t. A real computer programmer said things like

â€śDO 10 I=1,10â€ťandâ€śABENDâ€ť(they actually talked in capital letters, you understand) â€¦

In the 1950â€™s von Neumann â€¦ was confronted with the FORTRAN concept; John Backus remembered von Neumann being unimpressed and that he asked â€śwhy would you want more than machine language?â€ť