Co-array error

[JohnNichols]

Steve:

I am using some structures and I keep getting this 8363 error and I cannot see why.

Code enclosed in zip file.

Looks legal to me.

JMN

[John_Campbell]

I always get annoyed when I see these discussions about the use of SELECTED_INT_KIND or SELECTED_REAL_KIND, as it is approaching the basic issue backwards.

The idea that you can dial up a numeric format based on this approach is like taking us back to the approaches of COBOL.

Who would have devised a numerical algorithm then decided they needed precision based on "selected_real_kind(15, 307)". Did they decide that selected_real_kind(14, 300) wasn't enough so we will go the extra yard and shoot for 307 ? Or is it just that these are the documented requirements to select real*8.

The reality is that there is  typically 4-byte or 8-byte precision available and there is a balance between accuracy and storage requirement. When you have to resort to higher precision, real*16 is available but that is used when the numerical algorithm does not manage it's precision; time for a better algorithm !

For most real calculations, it is difficult to provide input measurements to better than 4 or 5 significant figures, so providing input data for a real*16 calculation is questionable.

The better approach is to know your computer hardware, find out what number formats are supported and choose the one most suitable. "selected_real_kind(15, 307)" says it all as this declaration says, what values do I need to provide to get real*8. Why not just say real*8, which more clearly selects the numeric format to be used.

I have always made my precision decisions based on the number of bytes, and when converting from alternative precision hardware (and I've seen most of them) the most important consideration was how many bytes is equivalent to how the number was previously stored. Typically there is no choice, as the answer is nearly always real*8, as the two other alternatives are not suitable.

Integer formats are also interesting, as the main influence is probably the compiler or operating system, being either 16-bit, 32-bit or 64-bit, with the move to 64-bit requiring more use of integer*8. This requires all the SELECTED_INT_KIND arguments to be changed from 9 to 10, when what we really want is 8 bytes.

John

[FortranFan]

John Campbell wrote:

I always get annoyed ...

The idea that you can dial up a numeric format based on this approach is like taking us back to the approaches of COBOL.

Who would have devised a numerical algorithm..

John

What is worth a consideration is the thought that these comments in Quote #22 may actually be misleading to many innocent Fortranner readers of this forum and they, in fact, may be more annoying to way more than the one and only person who feels annoyed by Quotes such as #12 and #20.  There are now several books and publications which suggest the same "best practice" based on SELECTED_REAL_KIND and SELECTED_INT_KIND as suggested by Steve above.  Even on this forum, Steve, Craig Dedo, and others have offered nice explanations many a times that provide the basis for said practice.

I think John Campbell should take his argument elsewhere, on a totally different platform such as peer-reviewed publications in computing journals or papers in ACM Fortran Forum, or write book(s) that argue his point - this is not the right avenue.   And if John (and anyone else) feels so strongly about the supposed absurdity of SELECTED_REAL_KIND, etc., they should look deeply into ISO_FORTRAN_ENV intrinsic module and what it offers.  ISO_FORTRAN_ENV module has been part of standard Fortran for quite a while now and I think it makes all of John's points moot: so he should offer a critical analysis in a scientific manner on a peer-reviewed platform with regard to any deficiencies in this intrinsic module relative to his continued insistence of REAL*4, ..., INTEGER*8 type of declarations which have long been communicated (since the Fortran 90 standard?) as being non-standard and for damn good reasons too.

Until then, if not SELECTED_REAL_KIND and SELECTED_INT_KIND, then I will find it far more preferable that fellow programmers make use of constructs such as those shown below than having to read/work with code that has REAL*8 foo, INTEGER*4 I, ... type of declarations:

   ..
   USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY : INT32, REAL64
   ..
   INTEGER(KIND=INT32) :: I
   REAL(KIND=REAL64)   :: foo
   ..

[JohnNichols]

John:

I do not disagree with you. This all started as I have been developing the code for a Water Supply Analysis program using the MKL library. One can follow the tortuous path on that forum. The de-facto standard for Water Supply distribution software is EPANET, this is coded in Pascal and C. This code is based the Gradient Method, which uses the Newton Ralphson method ( takes me back to 1980 and Conte and de Boore).  

I have done a lot of water supply modelling and have some large models, with solvers written using the Hardy Cross method and Optimization Techniques (developed by Kurczera at the University of Newcastle). The models are saved as AutoCAD drawings and it is quite painful to write programs to write the model data into the different forms used by different programs. I had presented a paper on this topic in 1988 to a conference on the need for a standard for exchanging data. There are sound economic reasons why a Water Supply Modelling program developer is not going to want to make it easy to transfer data to other models, some of these modellers are free and undercut the market.

I had always wanted to get the Gradient Method running on the large data files, but I was not interested in taking the data from the model into EPANET. It makes maintaining the models hard, and whilst the EPA GUI has a number of interesting features it is in PASCAL and the developer has retired. I stumbled across the original C code from 1991 for the EPA program and it was only 2000 lines. As mecej4 has shown it compiles nicely on the command line IFORT compiler, but it does not play nice with the VS compiler. Also I wanted to get the program to access the data directly from the AutoCAD drawing and return the output as a visual representation to the AutoCAD drawing. There are very  many sound engineering design reasons for this procedure.

So I coded the program in C#, and got it running (actually mecej4 solved the PARDISO problem for the C# code, I can not claim that step) . I also wrote a simple program in C# to tear the DXF apart and put it into null format. Rossman who had developed the 91 code had an input file in that code that is essentially a single line for each entity in the model. This is a great idea and probably the best thing he could have done, but for the 98 version he moved to multiple lines for the same data, so he had the node data on one line and several hundred lines later for a large model he would code the co-ordinates for the node. Makes it very hard to read the models and check for topological errors. So I used the 91 format and really it should be adopted as the default null standard, but that is beyond my pay grade.

The cycle time to fix the errors in the largest model is down to a few minutes now that I have this code. mecej4 asked for the large models as samples, at the moment I have only published a small model of about 14 nodes. In checking the large Tamworth model I found the key error in the EPANET procedure is the handling of the demand. Rossman had created a Qbase, but in reality one wants to know the underlying cause of the Qbase, so I took a step backwards and amended Rossman's format to read in the different demand types in the data file. It meant that the C# dxf converter needs to be fixed for this code change.

mecej4 commented that the Fortran version of the program would run faster, so in a fit of madness, I started to code the program in Fortran. The main issue is getting the data into the format required by the PARDISO solver, not an straightforward problem I can assure you. It seems simple, but it is painful. I stuck with the null format for the data file and so reading this data into Fortran is interesting, mecej4 suggested a cool method I am now using, Reid suggested an alternative that is also interesting. The problem is the 4 alpha group at the start of the line to say what each line is, this needs to be recovered so the right code is read for the input data. Rossman's 98 change was to put in a two alpha group at the start of each section of data. Everyone else has different formats for this type of data. Rossman's 91 method is like the original DXF format compared to his 98 method, which is a hodgepodge.

So over the weekend I had in most of the code- with quite a few logic mistakes, taking the square root of a negative number for instance.  I asked a few questions on this forum and Steve made the comment that I should use the SELECTED_REAL_KIND. At the time I was trying to get the PARDISO solver to accept the input, so stupidly I changed the code for the SELECTED_REAL_KIND, and I copied the suggested numbers from Reid's book. I created many more problems than I solved, the biggest problem was that the reals were all now 8's instead of 4's. mecej4 has made quite a few comments on this point in relation to the PARDISO solver.  At the moment I am using the 64 bit PARDISO solver (again thanks to mecej4) using the 5.0.0 version from the developer. I compile with VS.  I got it running last night.

I also noted as I was debugging the logic and coding errors in the Newton Ralphson model that the VS listed the output as REAL*8. Steve originally commented on my use of the double precision.

So I have been thinking about this issue, the selected real kind can pick only a few numbers of interest so why go to all the trouble, particularly as it impacts on the library files.  Can we put in a statement to check the real type and then link with the correct library call?  The SELECTED_REAL_KIND can give the error code -1, so what happens when that happens.

mecej4 comments strongly suggest hard coding the REAL*4 or REAL*8 into the code and forget the SELECTED_REAL_KIND issue, and return to the PARDISO solver not the 64 bit routine. Like a lot of coding stuff, you try and then see what breaks. Reid's comments on SELECTED_REAL_KIND do not give it a lot of support - I can see the use, but really you are back to "KNOW YOUR MACHINE".

In the end this whole exercise is about publishing a paper on the speeds of the different solution methods and the use of the null format. So it is  not about a correct answer, as we try and teach our students, there is always a correct answer, but the best answer for the problem at hand.

My preferred answer is to find the fastest solution, which at the moment for the small problem is 122 ms using Fortran and the PARDISO_64, as mecej4 said it would be. I was surprised the code speed went from 177-200 for the C# calling the PARDISO to 122 ms.

The problem set is of course constrained by physical reality:

a single pipe might supply a million people, 300,000 homes, or 45,000 litres per second, at 3 metres per second that is 15 square metres of pipe which has a diameter less than 4.5 metres or thereabout.  So the diameter is going to be less than 4500, the highest point in the world is 27000 feet as measured by the British who did not believe the results and made it a bit higher - or 9000 metres, and the longest model pipe is really not going to be over 10 km, the friction factors will be in the range of 0.02 to 200 depending on the model and g is 9.806 or 32. Assuming that there is 1 node for each 10 people as there is in the large models I have built so far, and we serve 300,000 homes that is 1 million people, or 100,000 nodes. So number of pipes is about equal to the number of nodes so we have no number bigger than 100,000.  The array size if it was dense is a problem but it is not dense.

A real 4 is more than adequate as is the default integer for the INTEL Fortran compiler.

As you rightly pointed out the problem is knowing your problem. The question then reminds one of the age old rivalary between an applied mathematician who develops a a solution and a pure mathematician who says "It is not elegant". Both are correct, which is the problem both are correct, but the PM is free to seek an elegant solution.

So you are the applied mathematician and Steve is the PM, both of you are correct. 

I love computers

JMN

[mecej4]

John Campbell wrote:
For most real calculations, it is difficult to provide input measurements to better than 4 or 5 significant figures,...

"Real calculations" by what criterion/definition?

Perhaps this is a bit off-topic, but here is a counterexample: Current uncorrected GPS uses six (or more) significant figures, which corresponds to a resolution of about 11 m. Map applications report positions in degrees to six places after the decimal point (after applying differential correction), which gives a resolution of 0.11 m. A cat trying to catch a mouse using GPS would need more precision than REAL*4 would provide.

[Steven_L_Intel1]

You don't have to check for a -1 value as the compiler will give an error if you use that for a kind. That's the design.

Not all implementations have the same integer or real kinds or sizes.

[FortranFan]

John Nichols wrote:

..

mecej4 comments strongly suggest hard coding the REAL*4 or REAL*8 into the code and forget the SELECTED_REAL_KIND issue, and return to the PARDISO solver not the 64 bit routine. ..

..

John Nichols,

I think you're totally mistaken: I highly doubt mecej4 ever suggested to you to hard-code REAL*4 or REAL*8 in your code.  He may have suggested to you to use 4-byte integers, instead of 8-byte integers, so you can avoid having to use PARDISO_64.  But that is a totally separate matter from Campbell's remarks in Quote #22.

Regardless, I'll be utterly astonished if mecej4 were to disavow the use of either SELECTED_INT_KIND or ISO_FORTRAN_ENV to define 4-byte integers and instead follow the integer*4 style touted by Campbell.

[JohnNichols]

I stand corrected, apologies to mecej4.

JMN

[mecej4]

Me too, I tried to apologize to myself, but the message came back marked "return to sender". Reminds me of a chapter heading decorated with a motto in the MacPascal manual: "RECURSIVE: see recursive", and of Groucho Marx's famous statement, "I wouldn't join any club that would have me for a member!"

[JohnNichols]

This problem would be solved if we all LISPED using a quantum computer that had a cycle time in femto seconds

In Oz we would call that a bl..dy beauty

JMN

[John_Campbell]

Fortran Fan,

My point is that when selecting a number precision, I think in terms of 4 bytes or 8 bytes ( 32 bits or 64 bits ) ; the commonly supported real number formats. I certainly don't think of an exponent of 307.

Your example code also demonstrates another common problem with the use of SELECTED_REAL_KIND, where this Kind value is defined in a module and not listed with the code. I prefer the numerical precision to be more accessible and clearly identified.

I don't agree that my point about clarity of coding should be excluded from this forum, as clarity and simplicity are important for minimising mistakes when developing new code. There is a significant issue here about how best to develop robust error free code, which has been a strength of Fortran for many years.

I do agree with you that this discussion has been addressed in the past and I should better control my annoyance.

John

[mecej4]

John Campbell has some valid concerns, and our thoughts about precision are definitely colored by our past experience. I do agree with what he wrote in the second paragraph of  #31, because I have often found myself miffed at having to search other files than the current source file, written by someone else, to find the definition 'WP' or some such vague mnemonic for 'working precision'. The program author's need for easy changing of precision (to see if double precision is really needed, for example) may lead to a coding style that is probably not good for a reader of the code who does not share the same needs as the author.

The kind numbers used in an implementation do not tell us much about precision, beyond "kind=3 probably provides more precision than kind=2". Likewise, REAL*4 and REAL*8 do not tell us exactly how much precision is available (24 and 53 bits in IEEE arithmetic). 

We may have found a magnitude range of 10^-38 to 10⁺³⁸ to be more than enough for most of our calculations, but there exist many examples where that range is insufficient, and having n-digits precision in the input does not gives n (or even one) correct digits in the output. One famous example is Wilkinson's Perfidious Polynomial, where changing the least significant bit in one coefficient (i.e, 1+2^-23 in place of 1) changes one root from 20 to 20.8 (see http://en.wikipedia.org/wiki/Wilkinson's_polynomial ). Another instance of the need for care in controlling precision arises in stiff differential equation IVP-s.

The key point is that errors in the output from a computation have more than one component, one part being attributable to errors in the input, and another to "round-off". When derivatives have been approximated using difference ratios, we have yet another source of error. The precision selected has to be sufficient to contend with all these causes of loss of precision.

[JohnNichols]

We are all trying to help people with our code, we must not lose sight of that fact, the tersest reply can often be a good learning lesson in - read the book dummy.

Although sometimes the book is quite confusing, I cannot say that I enjoy the 2003 Fortran explained book, Winston and Horne's Lisp is a much better read on a much harder topic, but one continues.

I serve on an IEEE committee, for a variety of reasons, but these IEEE people at the top of the game (and I am just a flunky) have much simpler issues to deal with, such as someone saying 16 inches in science terms is equivalent to 406.4 mm when we are talking about comparing two measurements. These people who wrote this are "scientists". At least this group is way beyond that level of physics.

JMN

[JohnNichols]

Steve:

When I try and upload the capture file showing the output files, I get the second capture file shown - I cannot proceed to actually use your picture icon to upload files, it always crashes now. The HTML offending code is a jscript with the name of the capture file embedded but not working button to close the window.

JMN

[JohnNichols]

Steve:

Is there a reason why I get wild swings in the run time for the code attached, it is the water supply program with a final name - a Norse God.  All my programs are ultimately named for the GODS.

it ranges from 109 to 280 ms using etime? is this normal?

Bugging the hell of me.

JMN

[FortranFan]

John Campbell wrote:

..

Your example code also demonstrates another common problem with the use of SELECTED_REAL_KIND, where this Kind value is defined in a module and not listed with the code. I prefer the numerical precision to be more accessible and clearly identified.

..

mecej4 wrote:

... Likewise, REAL*4 and REAL*8 do not tell us exactly how much precision is available ...

1. My example code was not with SELECTED_REAL_KIND, but with ISO_FORTRAN_ENV
2. REAL64 from ISO_FORTRAN_ENV is as identifiable as (or as vague as!) REAL*8, but  the latter style is now NON_STANDARD, so unless someone can make it part of the standard again, it is better not to bring REAL*8, INTEGER*4, etc. into discussions other than for legacy purposes.
3. But, when it comes to "clear" from a bits-and-bytes perspective, ISO_FORTRAN_ENV is intended to offer good foundation for writing portable code.  So I think it deserves a closer look by the numerically inclined.
4. Now when it comes to accessibility, I think the debate moves on to a separate terrain, one of personal preference.  Having to reference ISO_FORTRAN_ENV in my code via a USE statement (or WP/DP from some kinds module) is not something I consider a big deal; in fact, I think it offers benefits that are worth the extra lines of code I've to write.  And I only need to "look up" what REAL64 or INT32, etc. mean once or a few times; after that, "muscle memory" takes over.  But others may see it differently.  But again, the argument here is not about numerical precision, so it is a separate discussion.
5. Now to muddy up waters, the original thread is concerned about making use of external math libraries for sparse matrix solvers, etc. and since many such libraries are now written in C/C++, one should also think about ISO_C_BINDING module and making use of defined kinds therein such as C_DOUBLE, C_INT, etc. rather than assuming an intrinsic or defined kind in Fortran will always interoperate with the double, int, etc. from the companion "C" processor used by the external library.

[JohnNichols]

Steve:

In message 34, I include a screen shot, one shows the C# routine (lhs) output for the Water Supply Algorithm using 32 bit code and the MKL library PARDISO SOLVER and the rhs shows the Intel Fortran 64 bit and the PARDISO 5.0.0 solver.

The interesting observation is that the first inversion gives exactly the same relative error to 7 significant figures, it diverges in the second and subsequent iterations. As shown.  This means that there are several hypothesis to explain this behaviour?

1.  C# and Fortran give you slightly different answers - interestingly the final answers agree to 5 or 6 significant figures?

2. 64 and 32 bit differences in the code make a difference in accuracy of the algorithms?

3. The PARDISO solver from the MKL version to the latest one has made subtle changes?

4. Rounding differences in the use of REAL*8 etc. against the C# code precision? Sample of the C# below.

 }
                    double[] F = new double[Data.maxN];
                    double[,] J = new double[Data.maxN, Data.maxN];
                    int[] ipvt = new int[Data.maxN];
                    for (int j = 1; j <= 1; j++)

This is fairly simple algorithms that have been well studied for 50 years.

Any ideas?  It would be an interesting topic for an undergraduate class to explore.

regards

JMN

[Steven_L_Intel1]

I can't help regarding PARDISO. Feel free to ask in the MKL forum.

[JohnNichols]

But you should be able to comment on 1 or 2 - will they make a difference or it is all PARDISO.  

The people on MKL also generally read this forum.

So can I use the iso_fortran_env module in Intel Fortran and can you give us a short sample of how to declare a character, real and an integer that matches the new standard.  It ahs certainly generated some angst and the word from Dr Fortran should put this baby to bed.

The Fortran explained book appeared to suggest that if I had  Character*80 variable I could do variaible(j) or am I misunderstanding their point

JMN

[mecej4]

JMN: There is a bug in the program that you posted in #35 (and, possibly, earlier Fortran versions of your program). Here is what I see in magni.f90 (line numbers added):

330                  spipes(Npipes)%K1 = K1
382                  Read(iline,*, ERR = 160, END = 120)AF,dumPipeID, area,Hmin,Hmax,H0,K1
389                  stanks(Ntanks)%K1 = k1

The problem is that the component K1 of type SPIPE is integer, but the component K1 of type STANK is double precision. The convenience variable K1 used in the read statement is, however, is of type integer. Therefore, as Steve explains in https://software.intel.com/en-us/forums/topic/532551, reading the string "0.1" into an integer variable using list-directed input causes k1 to be 0. Thus, the value assigned to  stanks(Ntanks)%K1 is 0 instead of 0.1. Your C# program, however, would not be affected by this issue, and that may partly explain the discrepancies in the outputs, if the value assigned affects the output at all.

[FortranFan]

John Nichols wrote:

..

So can I use the iso_fortran_env module in Intel Fortran and can you give us a short sample of how to declare a character, real and an integer that matches the new standard.  It ahs certainly generated some angst and the word from Dr Fortran should put this baby to bed.

..

John Nichols,

Look here - https://software.intel.com/en-us/node/511041 - the online documentation for Intel Fortran that explains ISO_FORTRAN_ENV.  So yes, this intrinsic module is available in Intel Fortran.

Since it was I who brought this intrinsic module into the discussion, let me reiterate: most experts (and I dare say Steve included per his Quotes #12 and #20 above) still recommend user-defined kinds based on SELECTED_REAL_KIND and SELECTED_INT_KIND.  However, John Campbell (and perhaps some others) find this unappealing, especially the range part of the specification and prefer a more explicit invocation in terms of bytes such as REAL*4. It was in response to this that I wrote the Fortran 2003 standard includes the ISO_FORTRAN_ENV intrinsic that can allow Campbell and co. to continue somewhat similarly, in bits rather than bytes, and not have to "hard-wire" the KIND, a practice that Steve recommends we should avoid.  So this results in something like REAL(KIND=REAL32) instead of REAL*4 or REAL(4) and allows portability benefits to other compilers and platforms where ISO_FORTRAN_ENV is supported.  See Quote #23 for the example.  Now note one can also "rename" the intrinsic module parameters to something else, say a better mnemonic, and use it that way instead:

   ..
   USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY : R8B => REAL64
   ..
   REAL(R8B) :: foo

There are a lot of possible "good" practices here: the message is to avoid non-standard usage (/stand compiler flag can catch these) and not hard-wire the kinds.  I think Dr Fortran will be happy with this and leave it at that.  I don't think further dogma about using ISO_FORTRAN_ENV versus SELECTED_REAL_KIND,.. is warranted.

Also, as I mentioned, don't forget ISO_C_BINDING module for C interoperability.  I use it all the time with my C#-Fortran "mixed language, managed+unmanaged memory" solutions.

[Steven_L_Intel1]

I think there is some misunderstanding regarding SELECTED_INT_KIND and SELECTED_REAL_KIND. The idea is that you know what range and precision your application needs and you specify those in the intrinsics. The compiler then selects the "smallest" kind that meets the requirements. Many people will use 15,307 when they want IEEE double or something that acts like IEEE double. There's nothing magic about this - use (10,100) if you know that's all your application needs.

The point is to get away from putting implementation-specific details in your source code.