Hello,

I was reading peric CFD book and in forming the matrix A for the momentum eqn, the order of entries start from the southwest corner of the domain, proceeding northwards 1st along the grid lines, then eastwards across the domain.

however, i thought that in fortran, the movement should be

do j=1,nj do i=1,ni ie eastwards 1st then northwards.

which will result in faster operation

So which one will be better? I'm getting confused.

Thank you!

Peric stores the 2D array in a single 1d array so he doesn't need a nested do loop over (i,j).

Ya I know. What I mean is that is there a difference if the matrix is computed column 1st then row, or row 1st then column.

In fortran, doing row then column ie

do j=1,y do i=1,x is much faster than the reverse. Hence, I'm wondering if what Peric did influence the speed.

No it doesn't - he hasn't stored anything in a 2D array so there is no j loop in what he has done. What you are talking about is the fastest way to access data in a 2D array - in a 1D array this is not an issue.

the best way of programming is the one that others understand and is readable, especially if you are working in a big company. Nowadays, It is the compiler's and optimizer's job to take care of the speed. If you are a programmer and do a code that you are the only one to understand, you can put a chair and a laptop in your backgarden and stay there.

Those kind of programmers are mainly a waste of money for both research and industry.

Paul

The userfriendliness of a program is second to its structure when aiming at the optimization of computing time...

the structure being itself second to its maintainability and readability...

non readable program nobody wants to look at: hours spent fixing a bug = infinite optimized CPU run time = code is not being fixed total = infinite * 2

readable program: hours spent fixing a bug = 2.0 non optimized CPU time after bugfix = 2.0 total = 4.0

4.0 < infinite * 2

...so which program do you use at the end?

You are right if life is so simple that one can always use numbers to qualify things.. I am talking abt such people who are confident enough in his program before trying to optimize the computational time!

... instead of presenting irrelevant and unprompted opinions on the meaning of life?

Matrices in Fortran are "column major" as opposed to the "row major" structure of arrays in C. Internally all contiguous arrays are handled as 1D chunks of memory, anyway. The question is just how the 2D elements are ordered in this 1D array. For example, in a Fortran array x(i,j), the elements x(i,j) and x(i+1,j) are neighbors in the internal 1D representation, while elements x(i,j) and x(i,j+1) are separated by the size of the i-dimension. That's why running the first index in the inner-most loop is generally more efficient for large Fortran arrays, because of a higher chance of cache hits. (And yes, a good compiler may take care of that for you, although it's really not that difficult to write "do j; do i;" instead of "do i; do j;") The situation is reversed in C.

"... instead of presenting irrelevant and unprompted opinions on the meaning of life? "

Actually I did answer the question - my point is all this 2D (i,j) indexing is not what Peric does - the array x(M,N) is stored as y(M*N) and you are free to choose how x(i,j) relates to y(k); i.e. you could store x(i,j) in y(j*M+i) (the way FORTAN would internally store x) or, as Peric does, y(i*N+j).

The important point is that there is only the 1D y(k) array in Peric's code so that there is only 1 loop variable and not 2! So the question of which is faster for the 2D array case is irrelevant to What Peric is discussing in the book.

If your code is for academic puporse you can have fun with thisgeeky approach.

If you want a code that is aimed to last long or commercial, don't do this.

You can bet after reading Mani's answer, at least one nerd will play in inverting the indices of the vector to earn fantastic CPU speed-up, then won't put any comments or documentation, and you can say bye-bye to the maintenance and the parallelizatoion.

In theory it is fine. If you are realistic, don't do this y(i*N+j) thing, it is a waste of time.

You're Missing the point of the original posters question. The original question was concerned with the discussion in the book by Ferziger and Peric and the ordering of the 2D arrays. If you download the sample code from the book webpage (and read the relevant section of the book) you'll find that the ordering I've explained is exactly what Peric has done.

I never actually said it was a good idea to do things this way.

[Ford Prefect]

Hey,

After some digging through the archives I found this thread.

Does anyone know if the 1d array ordering of matrices is generally more difficult for the compiler to optimize, with regard to speed? I would assume that it all boils down to memory chunks and how they fit in the L1 cache, or?

Cheers!