Computer Languages for Engineers Book of Examples 2015 University of Duisburg-Essen Faculty of Engineering Department of Civil Engineering Structural Analysis and Construction Dr. E. Baeck 9.6.2015 Contents I FORTRAN 1 1 Development Tools 3 1.1 The Development Toolchain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Some History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.2 FORTRAN’s History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 The Fixed FORTRAN Format and it’s Roots . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Some free FORTRAN Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.5 The Open Watcom Development Suite . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.6 The MinGW Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.6.1 Running Compilers of the MinGW Package . . . . . . . . . . . . . . . . . . . . 10 1.6.2 Installing the MinGW Package . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.7 The G95 Compiler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.8 The Code::Blocks IDE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2 FORTRAN Basics 17 2.1 Structure of a FORTRAN Program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Format and Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3 Character Set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.4 Available Data Formats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4.1 Negative Numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4.2 Endianness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Data Types, Variables and Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.5.1 Data Types of FORTRAN 77 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.5.2 Data Types, Constants and KIND in FORTRAN 90 . . . . . . . . . . . . . . . . 23 2.5.3 Representation of a Float Number . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.5.4 Data Ranges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.5.5 Declaration of Variables and Constants in FORTRAN 77 . . . . . . . . . . . . . 26 2.5.6 Declaration of Variables and Constants in FORTRAN 90 . . . . . . . . . . . . . 27 2.5.7 Complex Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.6.1 Unary Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.6.2 Arithmetic Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.6.3 Comparison Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 File IO, Screen and Keyboard . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.7.1 32 2.5 2.6 2.7 Open a File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Page iv . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 34 35 35 36 36 37 38 39 40 42 42 43 43 44 44 46 48 48 49 50 50 50 51 52 53 56 58 58 58 59 61 63 . . . . . . 65 65 65 66 67 69 78 Linear Algebra, Vectors and Matrices 4.1 Helper Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Outlines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 Reset and List . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 85 85 85 2.8 2.9 2.10 2.11 2.12 2.13 2.14 3 4 Computer Languages for Engineering - SS 15 2.7.2 Writing Texts, write Statement . . . . . . . . . . . . . . . . . . . . . 2.7.3 Formatting, FORMAT Statement . . . . . . . . . . . . . . . . . . . 2.7.4 Read from Keyboard or File . . . . . . . . . . . . . . . . . . . . . . 2.7.5 Close a File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Loops . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8.1 Explicit Loops with Counter in 66, 77 and 90+ . . . . . . . . . . . . 2.8.2 Simple Nested Loop Example . . . . . . . . . . . . . . . . . . . . . 2.8.3 Quit a Cycle or a Loop . . . . . . . . . . . . . . . . . . . . . . . . . 2.8.4 Implicit, General Loop without a Control Structure in 90+ . . . . . . 2.8.5 Factorial in FORTRAN 90++ . . . . . . . . . . . . . . . . . . . . . Branching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9.1 if Statement, Fortran 66 like . . . . . . . . . . . . . . . . . . . . . . 2.9.2 Implementation of a Quadratic Equation Solver . . . . . . . . . . . . 2.9.2.1 Some Theory . . . . . . . . . . . . . . . . . . . . . . . . 2.9.2.2 A Flow-Chart of the QuadSolver . . . . . . . . . . . . . . 2.9.2.3 Quadratic Equation, Solver Implementation Fortran 66 like 2.9.2.4 Quadratic Equation, Solver Implementation Fortran 90 like Subroutines and Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10.1 Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10.2 Subroutines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10.3 Functions as Parameters . . . . . . . . . . . . . . . . . . . . . . . . Arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11.1 Static Array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11.2 Dynamical Array . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11.3 Automatic Array . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11.4 A little Array Example . . . . . . . . . . . . . . . . . . . . . . . . . 2.11.5 Pseudo Dynamic Arrays in FORTRAN77 . . . . . . . . . . . . . . . Global Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12.1 Classical Fortran and Common . . . . . . . . . . . . . . . . . . . . . 2.12.2 Some Aspects of the Module Concept of FORTRAN90 . . . . . . . . 2.12.3 Using global Data . . . . . . . . . . . . . . . . . . . . . . . . . . . Memory Mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Commandline Arguments . . . . . . . . . . . . . . . . . . . . . . . . . . . . Some Examples 3.1 Hello World . . . . . . . . . . . . . . . . . . 3.2 Simple Sum . . . . . . . . . . . . . . . . . . 3.3 Calculation of real*4/8 Precision . . . . . . . 3.4 Relative Precision with Functions . . . . . . 3.5 Newton’s Algorithm to calculate a Root . . . 3.6 Matrix Product with 77-Main and 90-Library E. Baeck . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CONTENTS 4.2 II 5 4.1.3 LU-Extract, MatMult and DiffMat . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.1.4 Text Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.1.5 Memory Manager . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.1.6 Multi Matrix Allocation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.1.7 Tracing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 Gauss-LU-Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.2.1 Gauss Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.2.2 Memory Manager . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 4.2.3 LU-Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.2.4 Tracing and Memory Pointers . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.2.4.1 Tracing gaussLU 4.2.4.2 Preparing Main . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.2.5 Gauss FB-Substitution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 4.2.6 Linare Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 C/C++ 123 Development Tools 5.1 6 Page v 127 The Code::Blocks IDE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 Basics of C/C++ 131 6.1 The Preprocessor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 6.2 Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 6.3 Data Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 6.4 Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.4.1 Assignment Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.4.2 Arthmetic Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.4.3 Compound Arithmetic Assignment . . . . . . . . . . . . . . . . . . . . . . . . 133 6.4.4 Increment - Decrement Operators . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.4.5 Relational and Equality Operators . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.4.6 Logical Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.4.7 Bitwise Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.4.8 Compound Bitwise Assignment . . . . . . . . . . . . . . . . . . . . . . . . . . 135 6.4.9 Explicit Type Casting Operator . . . . . . . . . . . . . . . . . . . . . . . . . . 135 6.4.10 sizeof Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 6.4.11 Address and Value Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 6.4.12 C++ operator synonyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 6.5 Taking about the Hello . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 6.6 Line Output with printf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 6.7 A For Loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 6.8 Static Arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 6.9 Branching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 6.9.1 if Branching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 6.9.2 switch Branching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 6.10 Exceptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 9.6.2015 Page vi 6.11 OOP with Classes . . . . . . . . . . . . 6.11.1 Some UML Diagrams . . . . . 6.11.2 C++ Class . . . . . . . . . . . . 6.11.2.1 Declaration . . . . . 6.11.2.2 Implementation . . . 6.11.2.3 Structures and Classes 7 III Computer Languages for Engineering - SS 15 . . . . . . . . . . . . Profile Example 7.1 Class Concept for the Tin-Walled Approach 7.2 Implementation . . . . . . . . . . . . . . . 7.2.1 Base, the Base Class of all Classes . 7.2.2 Node Class for Model Nodes . . . . 7.2.3 Checking the Node Class . . . . . . 7.2.4 Element Class for Model Elements . 7.2.5 Checking the Element Class . . . . 7.2.6 Profile Class for Model Profiles . . 7.2.7 Checking the Profile Class . . . . . 7.2.8 H-Profile Class for Model Profiles . 7.2.9 Checking the HProfile Class . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 148 149 149 151 152 . . . . . . . . . . . 157 157 157 159 161 163 165 169 171 180 183 186 Appendix A The Console’s Friends A.1 Directory Commands . . . . . . . . . . . . A.1.1 Selecting a Drive . . . . . . . . . . A.1.2 Listing the Content of a Directory . A.1.3 Creating and Removing a Directory A.1.4 Browsing through Directories . . . A.2 File Commands . . . . . . . . . . . . . . . A.3 Environment Commands . . . . . . . . . . 189 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 196 196 196 196 197 197 198 B Code::Blocks’s first Project 201 C Some Theory C.1 Section Properties . . . . . . . . . . . . . . . . . . . . . . . . C.1.1 The Area of a Profile Section . . . . . . . . . . . . . . C.1.2 First Moments of an Area . . . . . . . . . . . . . . . C.1.3 Second Moments of an Area . . . . . . . . . . . . . . C.1.4 Center of Mass . . . . . . . . . . . . . . . . . . . . . C.1.5 Moments of Inertia with Respect to the Center of Mass C.1.6 Main Axis Transformation . . . . . . . . . . . . . . . 209 209 209 209 210 211 211 212 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D Conventions 213 D.1 The Java Code Conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213 E. Baeck Part I FORTRAN 1 1 Development Tools 1.1 The Development Toolchain To develop a computer program there are three possible program types. • Interpreted Programs An interpreted program file is compiled, i.e. translated, into machine code at run time. That means, that repeated statements are compiled not only ones. To execute a interpreted program1 the source code is interpreted, i.e. compiled, line by line by an interpreter. A very simple interpreted program is a batch statement in the console window. The Basic language originally was designed as an interpreted language. The MS-Dos of the first years was shipping a Basic interpretor as development tool. Today we know the Basic language as Visual Basic for Application mostly from the MS software packages like MS Office. • Compiled Programs To create a compiled program, also called executable, there are some steps necessary. After having created the source files with an editor each of this source files has to be compiled by the so-called compiler into an object module. The object consists of binary native code, code which can be understood by the processor of the destination platform. Within a second step all object modules and the the used libraries are linked to an executable by the so-called linker. The libraries, especially the system libraries, are used to access the system resources like keyboard, screen and discs. This code is part of the develop system and can be used in an own application by linking. If you are developing software for a windows system (MS-Window, X-Windows, etc.) in an additional last step system resources are linked to the executable by an resource linkder. Program languages which are used to build an native executable are for example FORTRAN, C, C++2 . • Programs running on virtual machines Languages which are designed for bytecode are compiled in an neutral format. This format is not directly running an the processor. The compiled module can be executed on a real processor using a virtual machine, which translates the bytecode just in time before executing it. So bytecode may 1 An interpreted program often is called script. The languages C and C++ are very important, because the operating systems Linux and Windows are written in C, new parts of Windows are written in C++ and C# too 2 3 Page 4 Computer Languages for Engineering - SS 15 often be either directly executed on a virtual machine (i.e. interpreter), or it may be further compiled into machine code for better performance. Languages which are designed for the execution of bytecode are Java, Smalltalk, Python, Forth, TCL and C#. 1.2 1.2.1 Some History Motivation To understand the importance of the development of the FORTRAN language, the code for the calculation of the nth fibunacci number is given in machine and assembler code. The fibunacci number are defined as follows. fn = fn−1 + fn−2 for n > 2, with f0 = 1 and f1 = 1 (1.1) The implementation of the calculation of the nth fibunacci numberis given below in machine code 1 2 3 4 8B542408 FA027706 B9010000 C84AEBF1 83FA0077 06B80000 0000C383 B8010000 00C353BB 01000000 008D0419 83FA0376 078BD98B 5BC3 The implementation of the calculation of the nth fibunacci numberis given below in x86 assembler. 1 2 3 4 5 6 fib: mov edx, [esp+8] cmp edx, 0 ja @f mov eax, 0 ret 7 8 9 10 11 12 @@: cmp edx, 2 ja @f mov eax, 1 ret 13 14 15 16 17 @@: push ebx mov ebx, 1 mov ecx, 1 18 19 @@: 20 lea cmp jbe mov mov dec jmp @b 21 22 23 24 25 26 27 28 29 30 E. Baeck @@: pop ebx ret eax, edx, @f ebx, ecx, edx [ebx+ecx] 3 ecx eax 1.2. SOME HISTORY 1.2.2 Page 5 FORTRAN’s History The first compiler was written by Grace Hopper, in 1952, for the A-0 System language, which today has no relevance any more. The computer language FORTRAN (= FORmular TRANslator) was developed form IBM for the computer Type 704 in 1954-1957 (see figure 1.13 ) and was the first computer language which was able to handle mathematical formulas in nearly mathematical notation. In the time before FORTRAN only machine code or assembler was available. Therefore FORTRAN is a very important milestone in programming. The computer Type 704 was able to perform 4000 integer multiplications per second. A modern computer is able to perform some 100 millions of integer multiplications. The 1st FORTRAN version was followed in 1958 by FORTRAN II. FORTRAN IV was published in 1962 and became ANSI-standard as FORTRAN 66. The next standardized version is called FORTRAN 77. Today FORTRAN 77 is often used as computer language in engineering, especially in ad don software packages like ISML or NAG. Figure 1.1: IBM Type 704 [1] FORTRAN 90 was the next revision after FORTRAN 77. The concept of object orientated programming (OOP) was introduced in the FORTRAN language as well as the usual more flexible free format. Furthermore dynamical memory management, build in matrix arithmetic and the possibility of recursive functions were implemented. With FORTRAN 95 the next revision was published. Obsolete constructs are removed and besides automatic deallocation of arrays which go out of scope some new key words are introduced. With FORTRAN 2003 a better interoperability with the C programming language was introduced. A better integration into the host operating system is given (access to command line arguments and environment variables). The last revision is the revision FORTRAN 2008 which introduces Co-Arrays, a parallel processing model and the data type BIT. Since FORTRAN has been in use for more than fifty years, there is a vast body of FORTRAN in daily use throughout the scientific and engineering communities. It is the primary language for some of the most intensive supercomputing tasks, such as weather and climate modeling, computational fluid dynamics, computational chemistry, computational economics, plant breeding and computational physics. Even today, half a century later, many of the floating-point benchmarks to gauge the performance of new computer processors are still written in FORTRAN (e.g., CFP2006, the floating-point component of the SPEC CPU2006 benchmarks). 3 Note: The programmer sitting in front of his operator panel is not the author in his young days. 9.6.2015 Page 6 1.3 Computer Languages for Engineering - SS 15 The Fixed FORTRAN Format and it’s Roots The fixed FORTRAN format is a child of the punching cards time. Each line of code was punched onto one card. You had to punch n cards for a source code with n lines of code. To handle the cards in a batch, i.e. to avoid chaos in the order, the cards were numbered in the last 8 columns (see figure 1.2). Therefor in fixed FORTRAN format Figure 1.2: IBM FORTRAN Punching Card the last 8 columns of 80 available columns are columns of comment. Their information is not considered in the compile step. The leading 5 columns are columns for label numbers and must not used for statements. The column 6 is the column for the continuation line flag. A card punch (see figure 1.3) such as the IBM 3525 (not to be confused with keypunch), is an electronically mechanized output device used to punch data into punched cards. Sometimes combined with card readers and other functions to form multifunction machines, such as the IBM 2540 card reader-punch, such devices were attached to a computer. If you look at the properties of a console window on WindowsXP, you will find a buffer size per line for 80 characters as Figure 1.3: IBM Card Punch standard. The 80 character line length results from the format of a punch card with it’s 80 columns. E. Baeck 1.4. SOME FREE FORTRAN TOOLS 1.4 Page 7 Some free FORTRAN Tools Commercial FORTRAN compilers are often expensive and there for we try to find an open source development package for our lecture. Fortunately there are many FORTRAN development tools available for free in the Internet some of them are discussed below. If you use one of this packages it is recommended to update it from the original project site. • Watcom Fortran 77 Open Watcom is a project of the open source community to maintain and enhance the Watcom C, C++, and Fortran cross compilers and tools. An Open Source license from Sybase allows free commercial and non-commercial use of Open Watcom. More information about Open Watcom’s features, details about Open Watcom’s history and the supported platforms are given on the project site http://www.openwatcom.org/index.php. • MinGW Package MinGW 4 provides a complete Open Source programming tool set which is suitable for the development of native MS-Windows applications, and which do not depend on any 3rd-party C-Runtime DLLs (only the Microsoft C runtime, MSVCRT). More informations about MinGW are available on the project page http://www.mingw.org/. • G95 Package G95 is based on the G77 of the MinGW package. G95 provides a FORTRAN compiler for the versions 77, 90, 95 and 2003. The extension of source file selects the FORTRAN version. More informations about G95 are available on the project page http://g95.org/. • Code::Blocks Code::Blocks provides a cross-platform IDE for Linux, Mac-OS and Windows supporting a width range of available compilers and debuggers. More informations about Code::Blocks are available on the project site http://www.codeblocks.org. 1.5 The Open Watcom Development Suite The first FORTRAN 77 compiler of the Watcom Suite was published in 1985 for the IBM PC. Since 2003 the FORTRAN 77 and the C/C++ compilers are available as an open source project under Open Watcom. The Open Watcom IDE5 consists of the following three main modules. • IDE Project Manager Within the Project Manager application projects are initialized. The application type is set. Source files are added to the project. Tools are started out of the manager application. • Editor The Editor can be used to write source files for a project. 4 5 Minimalist GNU for Windows package is a porting of the LINUX development tools. IDE means Integrated Development Environment. 9.6.2015 Page 8 Computer Languages for Engineering - SS 15 • Debugger To check an application the executable can be started within a debugger program. The common features of a debugger are supported. The disadvantage of the Open Watcom suite is, that the three regarded components are not integrated within one application. Figure 1.4 shows the one and only source file of the famous hello world appliaction. Double clicking the source file in the manager will start the editor application (see figure 1.4). Figure 1.4: Watcom-IDE-Manager If standard FORTRAN 77 is used, formating could be a problem, because Open Watcom6 does not support any formating highlightings. So we use a formating comment, which marks the first 7 columns of source line. You have to be careful too with the comment columns starting with column 73. This comment columns also are not marked within the coding. If the program sources are compiled, the linker creates the executable linking the object files7 with the libraries. To check the executable the debugger can be started from the Watcom-IDE-Manager. Figure 1.6 shows a debugger session to check our famous startup example hello world. 6 7 Here we talk about the version 1.8. Object files are compiled source files, created by the compiler. E. Baeck 1.5. THE OPEN WATCOM DEVELOPMENT SUITE Page 9 Figure 1.5: Watcom-Editor Session to write the hello world Figure 1.6: Watcom-Debugger Session to check the hello world 9.6.2015 Page 10 1.6 Computer Languages for Engineering - SS 15 The MinGW Package The MinGW package is a porting of the development tools which are by standard available on the Linux/Unix platform. It’s no longer necessary to use a Linux emulating layer like cygwin8 . MinGW contents above all compilers, a linker, a debugger and a make utility. There is no IDE available within the MinGW package. That means if you want to use the pure MinGW package to develop applications you have to start the tool-chain manually. Because some of the tools are spawning other tools as helper application you have setup the correct path access to the MinGW binary folder. If we want to use the MinGW tools from the console window, we have to setup the correct path access. This is shown in figure 1.7. The MinGW binary folder is chained to the actual search path. There are a lot of compiler options which are supported by the MinGW compilers. Unfortunately there is no help page starting the compiler with the usual help option. So you have to read the manual or download some helping informations. Figure 1.7: Setup MinGW Path Access If we want to compile FORTRAN sources, we can use the G77 or the GFortran compiler. The GFortran compiler is the successor of the G77 and supports also the newer FORTRAN languages, FORTRAN 90, FORTRAN 95 and FORTRAN 2003 and parts of FORTRAN 2008. 1.6.1 Running Compilers of the MinGW Package Figure 1.8 shows a simple compiler call of the GFortran compiler, which compiles the FORTRAN source file and links it with the used libraries. Figure 1.8: Create an Executable with one Compiler Call After having called the compiler, we check the existents of the executable with a dir call. Then the 8 More information about cygwin is available on the project site http://www.cygwin.com/. E. Baeck 1.6. THE MINGW PACKAGE Page 11 executable is called and prints it’s legendary hello world to the screen. The most simple call of the compiler is shown in figure 1.8. The file hello.f is selected by the input filter *.f. The output file is set by the option -o. 1.6.2 Installing the MinGW Package If you have an archive file of the MinGW package you can simple extract the package files into a folder. It’s highly recommended to select a folder name, which is free from blanks. A blank usually is used to separate command line parameters, which are passed to programs to control their run-time behavior. Therefore blanks would break path names, if the module call is not set into double quotes. A second kind of installation is given by the MinGW Installation Manager, which you get, if you download newer installations. The installation manager comes with a dialog, where you should select the kind of installation you want (see figure 1.9). For our requirements we should install the following packages from the section Basic Setup. • minqw-32-base, the basic package is used, if we use the code::blocks IDE. • minqw-32-gcc-fortran, this is used for the Fortran part of the lecture. • minqw-32-gcc-g++, this is used for the C++ part of the lecture. You click them right and mark them for installation. After have marked all items to install, you start the function Apply Changes from the menu Installation. Using the MinGW Intallation Manager provides the advantage that we can easily select only the packages we really need. They only will be installed and updated if some newer versions are available. Figure 1.9: MinGW Installation Manager 9.6.2015 Page 12 1.7 Computer Languages for Engineering - SS 15 The G95 Compiler Alternative to the GFortran compiler you can use the G95 compiler, which provides nearly the same features as the GFortran compiler. G95 determines how an input file should be compiled based on its extension. Allowable file name extensions for Fortran source files are limited to .f, .F, .for, .FOR, .f90, .F90, .f95, .F95, .f03 and .F03. The filename extension determines whether Fortran sources are to be treated as fixed form, or free format. Files ending in .f, .F, .for, and .FOR are assumed to be fixed form source compatible with old f77 files. Files ending in .f90, .F90, .f95, .F95, .f03 and .F03 are assumed to be free source form. Files ending in uppercase letters are pre-processed with the C preprocessor by default, files ending in lowercase letters are not pre-processed by default. The basic options for compiling Fortran sources with g95 are: -c Compile only, do not run the linker. -v Show the actual programs invoked by g95 and their arguments. Particularly useful for tracking path problems. -o Specify the name of the output file, either an object file or the executable. An .exe extension is automatically added on Windows systems. If no output file is specified, the default output file is named a.out on unix, or a.exe on Windows systems. Informations about the compiler package are available on the project site http://g95.org. 1.8 The Code::Blocks IDE Code::Blocks was developed as a free C++ cross platform IDE. A FORTRAN version of Code::Blocks was initiated by Darius Markauskas with his project site http://darmar.vgtu.lt/. In 2014 Markauskas’s FORTRAN plugins for the Code::Blocks IDE were implemented and shipped with the standard version of Code::Blocks, which is version 13.13. So this is the version ore a junger version we should use to work on our FORTRAN and C++ projects. Code::Blocks is not written for a specific development package. The IDE provides a very general interface which is able to support a wide variety of compilers. Figure 1.10 shows the starting page of Code::Blocks. You can create a new project or can open an already existing project. You can also select one project from the recent list. To configure the compiler settings the menu command Settings/Compiler and debugger... should be executed. The first line selects the compiler for the file respectively the project. The tab compiler flags shows all supported compiler flags. Within the tab Toolchain executable the compiler, linker and make utility executable are set. Figure 1.12 shows the hello world FORTRAN source file in an editor window. The coloring supports the FORTRAN syntax. There are FORTRAN key-word lists for a fast completion of key words writing source files. Within the left browser window all files of the project are listed in tree mode. Within the right window you find the editor window and below the editor window the output section with information concerning the build tool chain execution. E. Baeck 1.8. THE CODE::BLOCKS IDE Page 13 Figure 1.10: Starting Code::Blocks To check the execution of an application, a breakpoint is set on the first source line 1.13. The debugging session is started with the menu command Debug/Start or with the F8 key. The execution is stopped at the first breakpoint, that is the second line of code. The yellow arrow shows the execution position of the program. The execution of the compiled and linked application is started with the menu command Build/Run or with the Ctrl-F10 key (see figure 1.14). The program is executed within a console window and stops at the end of execution. The window will be closed with any key. 9.6.2015 Page 14 Computer Languages for Engineering - SS 15 Figure 1.11: Setup the Projects Compiler Settings Figure 1.12: Code::Blocks Editor with hello world Source E. Baeck 1.8. THE CODE::BLOCKS IDE Page 15 Figure 1.13: Debugging the hello world Application Figure 1.14: Run the hello world Application 9.6.2015 Page 16 E. Baeck Computer Languages for Engineering - SS 15 2 FORTRAN Basics The description of FORTRAN 77 can be taken from the WATCOM FORTRAN Language Reference[2]. A detailed description of all statements is given there. The description of G95 and GFortran is available on the info.server1 with the documents G95Manual.pdf and GFortran.pdf. 2.1 Structure of a FORTRAN Program A FORTRAN program consists of a mixture of executable and non executable statements, which must occur in a specific order. So the FORTRAN program is divided into three sections. • The Declaration Section The first section is the declaration section with it’s non executable statements of declaration. Here all variables and constants or parameters are declared. Variables are optional initialized. • The Execution Section The second section is the execution section with it’s executable statements. This is the mostly the largest section, because this section contains all the statements, which describe what’s to do. And mostly there is a lot to do. • The Termination Section The termination section consists of statements, which stop the execution of the program and telling the compiler, that the program is complete. All FORTRAN statements should satisfy this requirements. If statements are found in the wrong section, the compiler stops compiling and the executable will not be build. 1 See http://info.statik.uni-due.de module Lehre/CM-CLFE. Actual versions are available on the project’s site. 17 Page 18 2.2 Computer Languages for Engineering - SS 15 Format and Comments Comment lines are denoted by placing a ”C” or ”*” in column one of the line. Blank lines are treated as comment lines too. Comment lines may be placed anywhere in the program source (i.e., they may appear before a FORTRAN statement, they may be intermingled with continuation lines, or they may appear after a statement). There is no restriction on the number of comment lines. Comment lines may contain any characters from the processor character set. Watcom FORTRAN 77 and the standardized with FORTRAN 90 allow end-of-line comments. If a ”!” character appears in column 1 or anywhere in the statement portion of a source line, the remainder of that line is treated as a comment unless the ”!” appears inside quotation marks or in column 6, if it’s a fixed format file. The FORTRAN 77 has by default the format shown in table 2.1. Columns Remarks 01 - 05 Column for label numbers (1 to 99999). Column 1 is also used to set comment lines. 06 Column 6 marks a continuation line for the previous line. The mark character can be every character of the FORTRAN character set but not zero or a blank. By default there are up to 19 continuation lines available. 07 - 72 Column for statements. 73 - 80 Comment Column. And are used in the days of the punch cards as card number field. Table 2.1: Fixed Fortran Format Some source of errors related to the fixed format are discussed below. • You should check the length of the code line. If the code line runs into the numbering field, i.e. into the columns 73 to 80, the code will be truncated and this can produce a very subtle error situation. • If the code is shifted into the header section of a line, i.e. into the columns 1 to 6, the code also will be truncated, but it’s more probable, that the compiler will detect this error. In free format FORTRAN which is introduced by standard with the FORTRAN 90 the sources are no more restricted by some fixed column positions. So the ”C” - comment character is no more available and the continuation column as well. A comment is always starting with the exclamation mark ”!” and a continuation line is introduced by the continuation character ”&” of the previous line which should be continued. The following example shows an implementation of the helloworld application which only will write it’s hello to the screen in a fixed formatted coding. The first line is a comment line, note the ”c” character in the first column. The third line is a continuation line, note the used continuation column 6. So the ”Hello again!” is printed directly behind the ”Hello World”. The last line closes the program code. Listing 2.1: Print one Hello 1 2 3 c234567 this is a comment write ( * , * ) ’Hello World!’ end E. Baeck 2.3. CHARACTER SET Page 19 Listing 2.2: Print a second Hello 1 2 3 4 c234567 this is a comment write ( * , * ) ’Hello World’, & ’ Hello again!’ end The next implementation of the Hello Word application uses the FORTRAN 90 free format. Listing 2.3: Print two Hellos with 90 1 2 3 4 !234567 this is a comment write ( * , * ) ’Hello World’, & ! this is a 2nd comment ’ Hello Again to 90!’ end You see within the first line of code the ”c” character was substituted by an exclamation mark ”!”. All code is shifted to the 1st column. A line end comment is used in the 2nd line using the exclamation mark ”!”. Right before the line end comment the line continuation character ”&” is set. The code position in the 3rd line is arbitrary. The end statement again closes the code. 2.3 Character Set The character set of FORTRAN 77 is the following.2 • upper case letter A - Z • 10 digits 0 - 9 • the 12 special characters: + - * / = ( ) : , . ’ $ and the space character. The character set of FORTRAN 90/95 is the following.3 • upper case letter A - Z • lower case letter a -z • 10 digits 0 - 9 • Miscellaneous common symbols, such as + - * / = ( ) : , . ’ $ ” { } [ ] ! • and any special letter of symbol required by the language, such as a¨ , o¨ , u¨ . 2 Modern FORTRAN 77 compiler support also lower case letters. The ASCII coding system which is used in computing stores one character in one byte. So ASCII is able to code maximal 256 characters. To support languages with more then 256 characters the Unicode coding, a multibyte coding, was developed, which is also supported by the actual FORTRAN 90/95. 3 9.6.2015 Page 20 2.4 Computer Languages for Engineering - SS 15 Available Data Formats The types of data format available on a computer are depending on the hardware. General there are some integer formats for integer values available. Standard is a short type with 2 bytes and a long type with 4 bytes. With the development of the 64 bit operating systems also 8 byte integers are available. Especially this is necessary to access to large files with a size larger then 2 GB or to be able to handle memory blocks with a size above 2 GB. On Windows32 with some provider depending tricks this was already possible, but according to a standard it was not till 64 bit systems are available. For floats general two formats are available, a 4 byte real and a 8 byte real. In the case of complex calculations as solving a linear equation system it is strongly recommended to use the 8 byte float. Within the following table data types are listed with there specific properties. Type Size Restriction INTEGER 2 Bytes Range −215 · · · 215 − 1 = 32767 INTEGER 4 Bytes Range −231 · · · 231 − 1 = 2147483647 ≈ 2.14 · 109 INTEGER 8 Bytes Range −263 · · · 263 − 1 = 9, 22 · 1018 REAL 4 Bytes Exponent range −38 · · · 38, 7 digits precision REAL 8 Bytes Exponent range −308 · · · 308, 15 digits precision CHARACTER 1 Byte one byte character Table 2.2: General Data Types with their Restrictions 2.4.1 Negative Numbers Negative numbers can be represented by complements of the positive number. Each digit of then positive number and it’s complement gives the maximum value of the digit. This would be 1 within the binary system. Table 2.3shows the scheme for the construction of a negative number for 7 using the b-complement. Comments 0000|01112 0716 , the positive number has the value 7 1111|10002 F 816 , number complement of 7 12 b-complement = complement +1, F 916 + 1111|10012 b-complement or negative number Check of the b-complement to be the searched negative number 0000|01112 positive number + 1111|10012 b-complement or negative number 1|0000|00002 The overflowing bit will be truncated and therefor the sum vanishes. Table 2.3: Construction of a negative Number E. Baeck 2.4. AVAILABLE DATA FORMATS 2.4.2 Page 21 Endianness A data type commonly consists of more then one byte, so the order of storing this bytes in memory becomes important. There is not only one order which is used to store this bytes. The terms endian and endianness, refer to how this bytes are stored in memory. Big-endian systems are systems in which the most significant byte is stored in the smallest address given and the least significant byte is stored in the largest address. This is shown in right part of figure 2.1. In contrast to this in the Little Endian order, the most significant byte is stored at the largest address and the least significant byte is stored in the smallest address. This is shown in left part of figure 2.1. Big Endian is like a German time stamp: hour-minute-second and Little Endian is like a German date: day-month-year. Figure 2.1: Little and Big Endian The Little Endian was introduced by Intel’s x86 processors. Therefor this Endian is very popular on all the Intel related computers. On the other hand the Big Endian was introduced by IBM processor architecture and is also used in the network like the IPv6. 9.6.2015 Page 22 2.5 Computer Languages for Engineering - SS 15 Data Types, Variables and Constants Within the following section we will discuss the available data types for the variables of FORTRAN 77 and FORTRAN 90 and above. We will discuss how to declare and initialize variables and constants. 2.5.1 Data Types of FORTRAN 77 Table 2.4 shows the available data types of FORTRAN 77 under standard. Type Comment Example INTEGER Integer number 15, -100, 2500 REAL Floatingpoint number, single precision 3.1415, -5.5, .7E3, 12.5E-5 DOUBLE PRECISION Floatingpoint number, double preci- 3.1415D0, -5.5D0, .7D3, 12.5D-5 sion COMPLEX Complex numbers, (two REAL num- (3.1415, -5.5), (1.4, 7.1E4) bers) LOGICAL Logical values .TRUE., .FALSE. Table 2.4: Fortran 77 Data Types Because the standard FORTRAN 77 only supports this few data types, the de facto standard of the FORTRAN 77 implementation supports some data types more, which are defined by the number of used bytes too. INTEGER*1, LOGICAL*1 INTEGER*2, LOGICAL*2 INTEGER*4, REAL*4, LOGICAL*4, (default on 32 bit platforms) INTEGER*8 REAL*8 (on 32 bit platforms the same as DOUBLE PRECISION) COMPLEX*16 (complex numbers with two DOUBLE-PRECISION numbers) 1 byte 2 bytes 4 bytes 8 bytes 8 bytes 16 bytes In FORTRAN77 the data type for characters or strings are the following. • A character like ’A’ can be stored in the data type CHARACTER or CHARACTER*1. • A string of n characters like ’The End!’ can be stored in the data type CHARACTER*n. E. Baeck 2.5. DATA TYPES, VARIABLES AND CONSTANTS 2.5.2 Page 23 Data Types, Constants and KIND in FORTRAN 90 Table 2.5 shows the available intrinsic buildin data types of FORTRAN 90 under standard. Type Comment Example INTEGER Integer number 15, -100, 2500 REAL Floatingpoint number, single precision 3.1415, -5.5, .7E3, 12.5E-5 COMPLEX Complex numbers, (two REAL numbers) (3.1415, -5.5), (1.4, 7.1E4) LOGICAL Logical values .TRUE., .FALSE. CHARACTER String variable with given length ’This is a string! Hey!’ Table 2.5: Fortran 90 Data Types Because there are de facto so many special data types in FORTRAN 90 a new concept was introduced to setup the desired byte length of a data type with the so called KIND number. The KIND number is depended of the implementation of the FORTRAN 90, that means that the byte length is not generally known, if the KIND number is given. The KIND number of a given value can be evaluated by the function KIND. The data type of a constant in FORTRAN90 can be specified by the extension [number of bytes]. The number of bytes in general are 1,2, 4 and 8 on a 64 bit platform. So we can check this with the following code snippet. Listing 2.4: Constants and it’s Check with KIND 1 2 3 4 5 6 7 ! check the write(*,’(" write(*,’(" write(*,’(" write(*,’(" write(*,’(" write(*,’(" kind of constants created with a literal constant 2_1 kind:",i2)’) kind(2_1) constant 2_2 kind:",i2)’) kind(2_2) constant 2_4 kind:",i2)’) kind(2_4) constant 2_8 kind:",i2)’) kind(2_8) constant 2._4 kind:",i2)’) kind(2._4) constant 2._8 kind:",i2)’) kind(2._8) The code in listing 2.4 will print the following output. 1 2 3 4 5 6 constant constant constant constant constant constant 2_1 2_2 2_4 2_8 2._4 2._8 kind: kind: kind: kind: kind: kind: 1 2 4 8 4 8 9.6.2015 Page 24 Computer Languages for Engineering - SS 15 The following FORTRAN 90 code evaluates the KIND number of 0.0 and 0.0d0, i.e. of single and double precision for the gfortran compiler on 32 bit Windows XP.4 Listing 2.5: Evaluating the Kind Number 90 1 2 3 4 5 6 7 ! Program to evaluate the kind number of ! - single precision (REAL) ! - double precision (DOUBLE PRECISION) program kinds write( * ,’(" Kind of single precision:",i2)’) kind(0.0) write( * ,’(" Kind of double precision:",i2)’) kind(0.0d0) end program ! it’s single ! it’s double Figure 2.2 shows the compile and execution step of the Kind application. And you can see from the output, that single precision gives a KIND number of 4, double precision a KIND number of 8 like real*4 and real*8 in the case FORTRAN 77 standard extension (see above). Figure 2.2: Evaluating the Kind Numbers Constant values can also be set in binary, octal and hexadecimal number system. This we can do using the prefixes b, o and x and appending the number with the universal kind of number representation. Listing 2.6: Constants in Some Number Systems 1 2 3 4 c234567 integer i_b /b’1101’/ integer i_o /o’12’/ integer i_x /x’ff’/ The code in listing 2.6 assigns the value 13 in binary, the number 10 in octal and the number 255 in hexadecimal and will print the following output. We can also see, how a variable in old FORTRAN can be initialized. The initialization value is put in between two slashes. 1 2 3 4 i_b = 13 i_o = 10 i_x = 255 You see, that in FORTRAN 90 a program starts with the key word program <program name> and ends with end program. E. Baeck 2.5. DATA TYPES, VARIABLES AND CONSTANTS 2.5.3 Page 25 Representation of a Float Number Float number can be described by the following equation. x = s · m · be s b e m where (2.1) is one bit for the sign is the base (floats according IEEE 754 use a base of 2) is the exponent (can overflow) is the mantissa Float numbers are stored according to the IEEE Standard for Floating-Point Arithmetic (IEEE 754). This technical standard for floating-point computation was established in 1985 by the Institute of Electrical and Electronics Engineers (IEEE). The representation of some floats are given in the table below5 . Name Common name binary16 Half precision 2 binary32 Single precision binary64 Double precision binary128 Quadruple precision Base Digits E min E max Decimal digits Decimal E max 10+1 −14 +15 3.31 4.51 2 23+1 −126 +127 7.22 38.23 2 52+1 −1022 +1023 15.95 307.95 2 112+1 −16382 +16383 34.02 4931.77 The number of decimal digits, the approximate precision in decimal is given by digits · log10 (base). The maximal decimal exponent is given by Emax · log10 (base). Figure 2.3 shows the usage of the bits in a Single Precision data type according to the IEEE-754. The sign V is stored in the highest bit. The next 8 lower bits are used for the exponent and the rest is used for the mantissa. Figure 2.3: Single Precision Float 2.5.4 Data Ranges If you select a data type for an implementation, then you should know the available data range of that data type. The data range can be depended of the installed operating system and from the processor himself. So an integer of the Windows3 has half the range as the integer from WindowsXP, if the range is not explicit set using the * version of the data type. If a not properly data type is selected, then the code will be inefficient if the data range is larger than needed. If the data range is smaller then needed, the code will not run properly because the data will go out of range, which will result an overflow with all it’s unpredictable consequences6 . 5 E is used for exponent. If you will implement a complex calculation algorithm with a lot of operations like the triangulation of a matrix, then you should use the double precision data types, to avoid a losses of information.) 6 9.6.2015 Page 26 Computer Languages for Engineering - SS 15 The following table shows the ranges of the FORTRAN 77 data types. Daten Type Bytes Data Range INTEGER*1 1 −128 · · · + 127 INTEGER*2 2 −32768 bis +32767 = 215 − 1 INTEGER*4 4 −231 bis +231 − 1 ≈ 2 · 109 REAL*4 4 ±10±38 , precision of 7 digits. REAL*8 8 ±10±308 , precision of 16 digits. 2.5.5 Declaration of Variables and Constants in FORTRAN 77 With an explicit declaration the data type of variables and parameter (constants) is set. The declarations should be made in the header section of a program or a subprogram, which is the section before the first executable statement is made (like an assignment or a write statement). The declaration is given by: Data Type [List of Variables] The list of variables contents one or more variable names, which are separated by commas. The following example shows the declaration of the variables V1,V2 and V3, which are declared as REAL*4. The variables I1,I2 and I3 are declared as INTEGER*2. Listing 2.7: Declaration in 66/77 1 2 3 c234567 REAL*4 V1,V2,V3 INTEGER*2 I1,I2,I3 Constants are called PARAMETER in FORTRAN 77. A constant is declared by an additional PARAMETER statement. The PARAMETER defines the value of the constant. This value is set by compilation and can not changed at run time (therefor it’s called a constant). The parameter is set by: PARAMETER (name = value) The following example declares three real variables and three integer parameters. Listing 2.8: Declaration and Parameters in 66/77 1 2 3 4 5 6 c234567 REAL*4 INTEGER*2 PARAMETER PARAMETER PARAMETER V1,V2,V3 I1,I2,I3 (I1 = 4) (I2 = 8) (I3 = I1+I2) Within the declaration section of a program variable can be initialized with the statement DATA, which follows the declaration of the variables to be initialized. Listing 2.9: Declaration and Initialization in 66/77 1 2 3 4 c234567 REAL*4 V1,V2,V3 INTEGER*2 I1,I2,I3 DATA I1 /4/, I2 /8/, I3 /12/ E. Baeck 2.5. DATA TYPES, VARIABLES AND CONSTANTS 2.5.6 Page 27 Declaration of Variables and Constants in FORTRAN 90 The declaration of variables in FORTRAN 90 is slightly different from the FORTRAN 77 format and is given by the following statement. Data Type :: List of Variables In the following example the three integer variables month, day and year are declared. In the second line the real variable seconds is declared. The character variables which should store some strings should be declared with a properly length. The first declaration shows the declaration with an explicit length parameter, the second declaration uses only the value of the length. Listing 2.10: Declaration in 90++ 1 2 3 4 INTEGER :: day, month, year REAL :: seconds CHARACTER(len = 10) :: prename CHARACTER(20) :: famname In the following example you can see how to initialize the just declared variables. This is like declaration and initialization in the language C simply within one step. So the obsolete statement DATA is no longer needed and can be canceled from modern FOTRAN 90 sources (see also section 2.5.5). Listing 2.11: Declaration and Initialization in 90++ 1 2 3 4 INTEGER :: day = 16, month = 10 , year = 2010 REAL :: seconds = 10.5 CHARACTER(len = 10) :: prename = ’Ernst’ CHARACTER(20) :: famname = ’Baeck’ The following example shows, how to declare parameters (constants). You see the code is very simular compared with the previous code. Only the attribut PARAMETER was added. Within the first code the content of the variables can be changed by an assignment. Within the second code, we have declared parameters and the content of this parameters can not be changed, because a parameter is fixed. So you see, that the FORTRAN 77 statement PARAMETER, now obsolete, is also no longer necessary and can be deleted from modern FORTRAN 90 code (see also section 2.5.5). Listing 2.12: Declaration and Parameters in 90++ 1 2 3 4 2.5.7 INTEGER, PARAMETER :: day = 16, month = 10 , year = 2010 REAL, PARAMETER :: seconds = 10.5 CHARACTER(len = 10), PARAMETER :: prename = ’Ernst’ CHARACTER(20), PARAMETER :: famname = ’Baeck’ Complex Data As we already have seen Fortran provides in contrast to many languages a complex data type. The complex data type is like a vector with two components, the real part and the imaginary part. The parts itself consists of Single Precision floats in the case of COMPLEX*8 and of Double Presision floats in the case of COMPLEX*16. 9.6.2015 Page 28 Computer Languages for Engineering - SS 15 The following code snipped 2.13 shows how to declare a complex variable, initialize it and do some basic calculations. It is shown too, how to get a square root of an arbitrary complex number, which can also be a negative real. Especially this is used in the example to solve the quadratic equation using complex variables. The programs output is shown in figure 2.4. We see in the first step the square root of −1 is calculated giving the imaginary unit i. Listing 2.13: Playing with Some Complex Numbers 1 2 program complexNumbers implicit none 3 4 5 6 complex(8) :: c1,c2,c3 c1 = (-1., 0.) write (*,*) "c1 = ",c1 ! here we use the double precision (kind 8) c2 = cdsqrt(c1) write (*,*) "c2 = ",c2 ! cd: is Complex Double precision c3 = c1 + c2 write (*,*) "c3 = ",c3 ! add two complex numbers c3 = cmplx(1.2,3.4) write (*,*) "c3 = ",c3 ! create a new complex value with the cmplx function c3 = c3/(-c2) write (*,*) "c3 = ",c3 ! division c3 = c3/2. +1 write (*,*) "c3 = ",c3 ! division and subtraction 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 end program Figure 2.4: Console Oupt of the Complex Data Example E. Baeck 2.6. OPERATORS 2.6 Page 29 Operators Fortran offers a set of operators, which are offered from the most computer languages. Fortran uses the same precedence as we know form the mathematics. The power operation has the strongest binding followed by the point operators (products and divisions) followed by the line operators (plus and minus). Unary operators will always be applied first. To change the standard precedence of the operators we use like in mathematics parenthesis to dictate the way a formula should be evaluated. 2.6.1 Unary Operators Unary operators are working only on one value, therefor unary. In all Fortran dialects there are the following unary operators available. Operator Comment 2.6.2 Example + plus operator - minus operator a = 2 >>> x = -a >>> -2 .not. logical inverse a = .false. >>> x = .not.a >>> .true. a = 2 >>> x = +a >>> +2 Arithmetic Operators Fortran offers the following arithmetic operators. You should be careful with the usage of data types especially within divisions. If you use integers, the result generally will be truncated. The following operators are available in all Fortran dialects. Operator Comment Example + sum operator x = 2+3 >>> 5 - subtraction operator x = 4-2 >>> 2 * product operator x = 2*4 >>> 8 / division operator x = 9/2 >>> 4 x = 9./2. >>> 4.5 ** power operator // concatenate of strings x = "hello"//"world" >>> "hello world" x = a**2 9.6.2015 Page 30 2.6.3 Computer Languages for Engineering - SS 15 Comparison Operators Boolean operators are used to branch and to make decisions. The comparing operators of Fortran 90 are now nearly identical to the C comparing operators. The following operators are available in all Fortran dialects. Operator Comment Example .lt. less than x = 23 .lt. 13 >>> .false. .le. less equal x = 23 .le. 23 >>> .true. .gt. greater x = 23 .gt. 13 >>> .true. .ge. left shift of bits x = 23 .ge. 23 >>> .true. .eq. equal x = 23 .eq. 23 >>> .true. .ne. not equal x = 23 .ne. 13 >>> .false. Within Fortran 90 there also the following C like comparing operators available. Operator Comment Example < less than x = 23 < <= less equal x = 23 <= 23 >>> .true. > greater x = 23 > 13 >>> .true. >= left shift of bits x = 23 >= 23 >>> .true. == equal x = 23 == 23 >>> .true. /= non equal x = 23 /= 23 >>> .false. 13 >>> .false. The result of a boolean expression like above are the boolean values .false. or .true.. E. Baeck 2.6. OPERATORS Page 31 To combine comparing expressions the following logical operators can be used.78 This operators are available in all FORTRAN dialects. There are no C-like boolean operators available of this type. Operator Comment Example .and. logical and x = 1 .lt. 2 .and. 2 .lt. 3 >>> .true. .or. logical or x = 1 .lt. 2 .or. .true. .equ. logical or x = .true. >>> y = .false. >>> x .eqv. y >>> .false. x = .true. >>> y = .false. >>> x .neqv. y >>> .true. .nequ. logical or .not. 2 .gr. 3 >>> logical not x = .not. (1 .lt. 2) >>> .false. The truth table for the AND operator ∧ is given as follows. .true. ∧ .true. = .true. (2.2) .true. ∧ .false. = .false. .false. ∧ .true. = .false. .false. ∧ .false. = .false. The truth table for the OR operator ∨ is given as follows. .true. ∨ .true. = .true. (2.3) .true. ∨ .false. = .true. .false. ∨ .true. = .true. .false. ∨ .false. = .false. The precedence of the operators, i.e the order of evaluation is given in the following list. 1 2 3 4 5 6 7 8 9 10 highest precedence | | | | | | | V highest precedence ** * + + // .EQ. .NOT. .AND. .OR. .EQV. / - unary operators .NE. .LT. .LE. .GT. .GE. .NEQV. 7 To make expressions clear parenthesis should be used. A term within a parenthesis is evaluated first and it’s result then is used in further evaluations outside the parenthesis. With parenthesis the order of the evaluation can be set. 8 The logical operator .NEQV. implements the exclusive OR. 9.6.2015 Page 32 2.7 Computer Languages for Engineering - SS 15 File IO, Screen and Keyboard To communicate with the program the two statements read and write are necessary. In the most simple case we read the input data from the keyboard and write the output data into a console window. A simple output gives us the example helloworld, see section 2.2. 2.7.1 Open a File The open statement assigns a file name to a channel number. Every read or write access to this file uses this channel number. An open statement also has to set up the access type with the status parameter. The opening of a file should be checked by the iostat parameter. Some status values of the open statement are discussed in the table below. In Fortran77 only ’new’ and ’old’ are available. status Comment ’old’ The file must exist ’new’ A new file is created. No file with this name must exist. ’replace’ A new file is created. No matter whether there is an old file or not. ’scratch’ A temporary file is created, a file name is not needed. if the file is closed, the file is automatically removed ’unknown’ The same as ’replace’. It’s a pre Fortran90 statement and should be replaced by ’replace’ The access type can be specified with the ’action’ parameter. The available values are discussed in the table below. This is only available in Fortran90. action Comment ’read’ The file is opened only for read access. ’write’ The file is opened only for write access. ’readwrite’ This ’action’ is set if the parameter is not given. The file is opened for reading and writing. With the parameter ’position’ the write position can be specified. By default the write position of a text file is set to the beginning of the file. Sometimes it is needed to append information to a text file, for example, if data should be written into a log file. position Comment ’append’ The file is opened and the write position is set to the file’s end. In example 2.14 shows a snipped to open an existing file. If the file is not available or an open error occur, the error handler will set the iostat variable to a value no equal zero. In this case the program will stop. E. Baeck 2.7. FILE IO, SCREEN AND KEYBOARD Page 33 Listing 2.14: Open a file in 90++ with error Handler 1 2 3 4 5 2.7.2 open(ioin,file=infile,status=’old’,iostat=ioerr) if (ioerr /= 0) then write(*,*) "error: file ’",infile(1:len_trim(infile)),"’ not found." stop endif Writing Texts, write Statement The write statement has two parameters and a data list. 1 WRITE (<parameter1>,<parameter2>) <data list> Parameter 1 sets up the output channel. If we use the value *, the output is written into the console. If we use an integer number, we have to open a channel before for writing and assign this channel with the desired channel number. The following example shows the usage of the console window and the output into a used channel. 1 2 write ( *,*) ’Hallo my Channel’ write (10,*) ’Hallo Channel 10’ ! output into the console ! output into channel 10 Parameter 2 sets up the output format. If we use the value *, a standard format is used. Strings are written in total length and values, integer or floats, are written with full precision. In example 2.15 we see, how to open a new file and write some formated data into it. The file will be created, if not existing, in the current folder, because we have not specified an path. Listing 2.15: Open a file and write some data 1 2 3 4 open(10,file="MyFile.txt",status=’replace’) write(10,100) "Let’s write the number",10," into a file!" 100 format(a,i3,a) close(10) If we opne the file MyFile.txt, we will see the following content. 1 Let’s write the number 10 into a file! Further details to format types are found int the next section. 9.6.2015 Page 34 2.7.3 Computer Languages for Engineering - SS 15 Formatting, FORMAT Statement If data should be formatted, we can use the format statement or we can use the parameter of the format statement as string as parameter 2 of the write statement. The general form of an output format is [m]Tw[.n]. 1. Format Type T • I, integer decimal format. • Z, integer hexadecimal format. • O, integer octal format. • B, integer binary format. • F, float format with fixed decimal point. • E, float format with exponential representation, single precision. • D, float format with exponential representation, double precision. • A, text format, the with of the format is optional. • X, output of spaces • /, a new line, linebreak 2. repeat factor m 3. width of the output field w 4. number of significant digits n within float formatting. If n is used in integer formattings, the leading blanks are filled with zeros. Formats can also be iterated by iteration factors. If more complex formats should be iterated, the format block must be bracketed by round parentheses. Within the following example the output of 3 variable values a, b and c should be written. We use one format with an iteration of 3. The format starts with 2 blanks (2x). Then a text will be written, (a format) and at the end the value of the variable should be written using a fixed float format (f). The output should have a width of 10 with 2 digits after the decimal point. Listing 2.16: Writing some formatted Values to the Screen in F77 1 2 3 4 c W| 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 write(*,8000) ’a=’,a,’b=’,b,’c=’,c 8000 format(3(2x,a,f10.2)) E. Baeck 2.7. FILE IO, SCREEN AND KEYBOARD 2.7.4 Page 35 Read from Keyboard or File The read statement is similar to the write statement. We simply exchange source and destination. The first parameter of the read statement passes the io channel. The second parameter passes the input format. If a * is given, the format is free. That is, the input information units are separated by blanks (white spaces). Optionally a read io error can be handled by the iostat parameter. 1 READ (<channel>, <format> [,<iostat=errvariable>]) <variable list> Example 2.17 shows, how to read data from a text file. We will read two integer numbers using a free formatting, i.e. the numbers in the file have to be separated by white spaces (i.e. blanks or tabs). If an error occur, the error handler will set an error value not equal to zero. An error for example can occur, if there are no numbers in this textfile, so that a format clashing will happen. Listing 2.17: Read two numbers from a text file with a free formatting read(20,*,iostat=ioerr) n1,n2 if (ioerr /= 0) then write(*,*) ’*** error reading numbers from file!’ endif 1 2 3 4 The following line will be read correctly. 1 1 2 If we would read the following line, the error handler will throw an exception, because we can not convert text into numbers. 1 2.7.5 Hello world without any numbers! Close a File If an open file is no longer needed, the file should be closed. Optionally with the status parameter the file can be deleted (’delete’) or saved. If the status parameter is not given, a standard file is kept (’keep’). 1 CLOSE (<channel>, [<status = ’keep’/’delete’>]) Listing 2.18: Deleting a file by open/close calls in F77 1 2 3 4 c W| 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012 open(ioin,file=infile,status=’old’,iostat=ioerr) if (ioerr .eq. 0) close(ioin,status=’delete’) 9.6.2015 Page 36 2.8 Computer Languages for Engineering - SS 15 Loops Within this section we discuss the explicit loop statement. Until the FORTRAN 90 there are no statements to cancel the loop or to cancel a cycle to start the next one immediately. Within the FORTRAN versions up to 77 this statements have to be implemented by explicit jumps using goto statements. 2.8.1 Explicit Loops with Counter in 66, 77 and 90+ A loop that executes a block of statements a specified number of times is called an iterative DO loop or counting loop. A counting loop construct has the following form in FORTRAN 66, where index is the counting variable, which starts with the value istart and iterates up to the value iend. After each cycle the counting variable is incremented by inc. If inc is not given, inc is set to 1. In FORTRAN 66 the DO statement needs an ending label number (in this example 100). The loop executes all lines of code until the line with the given label number. 1 2 3 4 5 c DO loop in FORTRAN66 c234567 do 100 index = istart, iend [, inc] <Block of statements> 100 continue In FORTRAN 77 the labeled DO is substituted by a DO ... END DO construct. Single loop cycles as well as the total loop can be broken by a goto statement jumping to an appropriate label. 1 2 3 4 5 c DO loop in FORTRAN90 c234567 do index = istart, iend [, inc] <Block of statements> end do A loop cycle can be broken by the statement cycle. The cycle statement breaks the actual cycle and starts the next one, if a new cycle should be executed. The total loop can be broken by an exit statement. If we want apply the FORTRAN 90 we use the free format.9 1 2 3 4 ! DO loop in FORTRAN90 do index = istart, iend [, inc] <Block of statements> end do 9 Please note, that there are a lot of non standardized loop constructions, which are added to some Fortran 77 implementations like Watcom77. The Basic like loop structure with LOOP...UNTIL (BOOLEAN EXPRESSION) is not supported by standard and so it should be substituted by the new Fortran 90 loop version. E. Baeck 2.8. LOOPS 2.8.2 Page 37 Simple Nested Loop Example Within the following example all we have a loop from 1 to 10 and an nested loop form 1 to 5. The numbers of the counter variables should be printed and their product to. The implementation ist given in FORTRAN 66, FORTRAN 77 and FORTRAN 90. The FORTRAN 66 is given as follows. The outer loop is labeled by 100 and the inner loop is labeled by 50. Note that each loop needs his own unique counter variable and nested loops must be nested totally, i.e. they must not overlap. Listing 2.19: A nested 66-Loop Example 1 2 3 4 5 6 7 c234567 do 100 i=1,10 do 50 j=1,5 write( * , * ) ’i=’,i,’ j=’,j,’ i * j=’,i * j 50 continue 100 continue end The FORTRAN 77 version differs from the FORTRAN 66 version only by removing the label and closing the loop with an end do statement. Listing 2.20: A nested 77-Loop Example 1 2 3 4 5 6 7 c234567 do 100 i=1,10 do 50 j=1,5 write( * , * ) ’i=’,i,’ j=’,j,’ i * j=’,i * j end do end do end Using an indent of code blocks in loops can be problematic because in FORTRAN 77 we only have the columns from 7 to 72. A much more nicer code can be received, if we use the FORTRAN 90 version with it’s free formating. Now FORTRAN has the look of a realy modern language. We start with the 1st column and we have enough space for indenting. The problem of truncating the code with the comment region is no longer existing. Listing 2.21: A nested 90-Loop Example 1 2 3 4 5 6 7 program loop90 do i=1,10 do j=1,5 write( * , * ) ’i=’,i,’ j=’,j,’ i * j=’,i * j end do end do end loop90 9.6.2015 Page 38 2.8.3 Computer Languages for Engineering - SS 15 Quit a Cycle or a Loop In FORTRAN 66 everything is done using an goto statement, so in this dialect we use goto too to quit a cycle or a the entire loop. Listing 2.22: Breaking a 66-Loop 1 2 3 c234567................................................................ do 100 i = 1,10 write(*,*)"Begin of a cycle!" 4 5 c for i=3 the cycle should not be executed if (i.eq.3) goto 100 c break the entire loop for all i > 5 if (i.gt.5) goto 110 6 7 8 9 10 11 12 13 14 write(*,*)"i=",i 100 continue 110 write(*,*)"End of the loop!" end With FORTRAN 90 two very helpful statements came into picture to avoid explicit jumps using goto. • cycle cancels the current cycle and starts the next • exit cancels the entire loop Listing 2.23: Breaking an explicit 90-Loop 1 2 3 program loopescape do i = 1,10 write(*,*)"Begin of a cycle!" 4 ! for i=3 the cycle should not be executed if (i.eq.3) cycle 5 6 7 ! break the entire loop for all i > 5 if (i.gt.5) exit 8 9 10 write(*,*)"i=",i end do 11 12 13 14 15 write(*,*)"End of the loop!" end program loopescape E. Baeck 2.8. LOOPS 2.8.4 Page 39 Implicit, General Loop without a Control Structure in 90+ In Fortran 90 the do ... end do statement can be used without any control elements too, so the break condition of the look has to be implemented manually by the usage of an if statement (see section 2.9 too). If we compare the listings 2.23 with 2.24, we see, that incrementing the loop counter has to be done manually as well as it’s initialization. Using an implicit (general) loop we have to be careful with the break condition, so that we avoid hanging in an endless loop. Listing 2.24: Breaking a general 90-Loop 1 2 3 4 5 program GeneralLoop i = 0 do i = i+1 write(*,*)"Begin of a cycle!" 6 7 8 ! for i=3 the cycle should not be executed if (i.eq.3) cycle 9 10 11 ! break the entire loop for all i > 5 if (i.gt.5) exit 12 13 write(*,*)"i=",i 14 15 end do 16 17 18 write(*,*)"End of the loop!" end program GeneralLoop 9.6.2015 Page 40 2.8.5 Computer Languages for Engineering - SS 15 Factorial in FORTRAN 90++ The factorial of n is defined iteratively by n! = n Y i (2.4) i=1 so we can translate the product easily into a loop. A Q P product symbol or a sum will be represented in a programming language by a loop. The counter variable of the loop is given by the index range of the product or sum symbol. The result of a product or a sum will be assigned to a result variable. Start p = 1; i = 1; p = p·i A flowchart for the factorial algorithm is given in figure 2.5. After having initialized the product variable p and the loop variable i which is used as a factor in every loop cycle, we multiply the product variable p by the loop variable i and assign this new product value to our product variable. After the product step we have to check whether we have reached the end, i.e. if i is equal to n which is given as the argument of our factorial function. i = i +1 no i =n yes print results End Figure 2.5: The Factorial’s Flowchart Listing 2.25: Implementation of the Factorial 1 2 3 4 ! calculating the factorial ! ! analysing effects of the data type size ! using the factorial 5 6 7 8 9 program factorial90 ! integer(2) :: f ! integer(4) :: f ! real(4) :: f ! result value integer 2 bytes ! result value integer 4 bytes ! result value real 4 bytes 10 11 12 13 real(8) :: f integer(2) :: i integer(2) :: n ! result value real 8 bytes ! loop index variable ! input value 14 15 16 n = 200 write(*,*) ’n = ’,n 17 18 19 f = 1 do i = 1,n 20 21 22 23 24 E. Baeck !new old f = f * i write(*,*) i,’! = ’,f end do 2.8. LOOPS 25 Page 41 end program Within this example you can check the range of the available data types. If p is set to integer 4 bytes are used to store the factorial. The larges integer within 4 bytes can be calculated as follows.10 max = 231 − 1 = 2147483647 ≈ 2, 147 · 109 (2.5) Figure 2.6 shows the problem of an overflow, if the numbers to store exceeds the limit of the data type. The factorial von 12 looks like plausible. But the next must be wrong, because 4 · 13 6= 19. If we check the next factorials, we will see, that the numbers not really exceed 2 · 109 . Thats the 2GB Problem. Furthermore we see, that some numbers become negative. This obvious is also wrong and a consequence of the overflow situation. Figure 2.6: Factorial using INTEGER with 4 Bytes For larger numbers the choice of a real data type is recommended. Within a real data type only the exponent can overflow. The mantissa can not overflow, because it’s independent of the number’s size. If we use the largest available data type double precision11 , we will get the situation of figure 2.7. We see, that with the 8 byte real the factorial reaches 170!. The next result is indicated as Infinity. This happens, if the memory for the exponent overflows. A 4 byte real has an exponent range of −38 · · · + 38, a 8 byte real has an exponent range of −308 · · · + 308. We see, that the largest occurring exponent is +306. 10 This is also known as 2GB problem. This problem occur for example, if a file exceeds the size of 2GB, because a 4 byte integer is used to address the file’s bytes reading or writing them. A similar problem occur, if you want to have more then 256 columns in MS-EXCEL2003, that’s a one byte limit, or if you need more then 65536 rows in MS-EXCEL2003, that’s a 2 byte limit. 11 The data type double precision with FORTRAN 90 is obsolete, see also section 2.5.2, but can be used, if you want. 9.6.2015 Page 42 Computer Languages for Engineering - SS 15 Figure 2.7: Factorial using Real with 8 Bytes 2.9 Branching Branching and decisions can be implemented in Fortran like in the most languages with an if -statement. The application of if constructs will be discussed using the implementation of the general form of a quadratic equation. Within a first approach the implementation in a Fortran66 like program is shown. 2.9.1 if Statement, Fortran 66 like Within the first standardized Fortran version, which will be compiled from modern Fortran compilers too, the if statement is very rough. It’s like a branch in assembler language, that is, the if statement only is able to process one statement. Note the two statements in the macro assembler code of section 1.2.1 resumed below. Within a first step, a registers data is compared. In a second step a conditional jump is performed. 1 2 3 4 ... cmp edx, 2 ja @f ... < comparing the content of the edx with 2 < jump, if equal Listing 2.26: Syntax of the if Statement 1 IF (<logical expression>) <statement> Because of the similar structure of assembler branches and branches in Fortran 66 the developed Fortran code will be very similar to the assembler code comparing their case structures. An example to implement an implicit loop with an if statement and a backward jump to calculate the relative precision is given in section 3.3. E. Baeck 2.9. BRANCHING 2.9.2 Page 43 Implementation of a Quadratic Equation Solver The solver of a generall quadratic equation is a vell known problem, we know from our school days. The implementation whowever requires the solution of a set of subproblems, which indeed is a very good example to show branching. 2.9.2.1 Some Theory The following example implements the solver for a general quadratic equation. a · x2 + b · x + c = 0 (2.6) We have the following cases. • a =0∧b =0∧c =0 Infinite solutions. x can be set arbitrary. • a = 0 ∧ b = 0 ∧ c 6= 0 No solution possible. • a = 0 ∧ b 6= 0 Linear case, x = −c/b. • a 6= 0 Quadratic case, √ 1 (−b + b 2 − 4ac). x1 = 2a √ 1 x2 = 2a (−b − b 2 − 4ac). 9.6.2015 Page 44 Computer Languages for Engineering - SS 15 2.9.2.2 A Flow-Chart of the QuadSolver The following flow chart shows all the case, which we have to handle. The algorithm is given for a real arithmetic, i.e. no complex data types are used. The relevant source code will be developed within the next section. Start yes a=0 c=0 no no yes x1,2 = yes infinit solutions Stop no x = − bc d = b2 − 4 · a · c d <0 yes b=0 Stop √ −b±i −d 2·a no solution Stop Stop no x1,2 = 2.9.2.3 √ −b± d 2·a Stop Quadratic Equation, Solver Implementation Fortran 66 like The first implementation in strict Forteran 66 shows the subsequent solution of the above discussed cases. Because the if statement can only process one statement, the inverse case should be checked to jump over the succeeding code. After the code block a further jump should be performed to the last statement of the program. Listing 2.27: Implementation of a 66-Quad-Solver 1 2 3 4 5 6 7 8 9 10 11 12 c quadratic equation c a*x**2 + b*x + c = c c a,b,c are arbitray input parameters c c explicit declaration should be done c c234567 implicit none c c input parameters double precision a,b,c 13 14 c working variables double precision d,p c output float values double precision x1,x2 15 16 17 18 19 E. Baeck 2.9. BRANCHING 20 Page 45 c output complex values double precision x1r,x2r,x1i,x2i c Initialization p = 1.d-15 a = 1.d0 b = 0.d0 c = 4.d0 21 22 23 24 25 26 27 28 29 c c 30 31 32 33 34 35 36 37 c c 38 39 40 41 input parameters from the keyboard write(*,*) ’input of a:’ read(*,*) a write(*,*) ’input of b:’ read(*,*) b write(*,*) ’input of c’ read(*,*) c list input write(*,*) write(*,*) write(*,*) write(*,*) parameters for checking ’ Input parameters:’ ’ a = ’,a ’ b = ’,b ’ c = ’,c 42 43 c linear, constant branch if (dabs(a) .gt. p) goto 500 c contant branch if (dabs(b) .gt. p) goto 400 44 45 46 47 48 if (dabs(c) .gt. p) 1write(*,*) ’No solution found, constant case.’ 49 50 51 if (dabs(c) .le. p) 1write(*,*) ’Infinit solutions found, constant case.’ 52 53 54 goto 600 55 56 57 c 58 59 60 61 linear branch 400 continue x1 = -c/b write (*,*) ’Linear case: x = ’, x1 goto 600 62 63 64 65 66 c quadratic branch 500 continue c calculating the discriminante d = b**2 -4.d0*a*c 67 68 69 c reel branch if (d .lt. 0.) goto 550 70 71 72 73 d = dsqrt(d) x1 = (-b +d)/(2.e0*a) x2 = (-b -d)/(2.e0*a) 74 75 write(*,*) ’quadratic case, reel values’ 9.6.2015 Page 46 Computer Languages for Engineering - SS 15 write(*,*) ’ x1 = ’,x1 write(*,*) ’ x2 = ’,x2 76 77 78 goto 600 79 80 81 c 82 complex branch 550 continue 83 d x1r x2r x1i x2i 84 85 86 87 88 = = = = = dsqrt(-d) -b/(2.e0*a) x1r d/(2.e0*a) -d/(2.e0*a) 89 write(*,*) ’quadratic case, complex values’ write(*,*) ’ x1 = ’,x1r,’ +i ’,x1i write(*,*) ’ x2 = ’,x2r,’ +i ’,x2i 90 91 92 93 600 continue end 94 95 2.9.2.4 Quadratic Equation, Solver Implementation Fortran 90 like The following code implements the solution of a quadratic equation (see equation 2.6) in a Fortran 90 version. Note, that we are able to implement the case tree without any goto jumps, which were essential in an 66 approach. In contrast to solution 2.27 in the following code 2.28 we use a complex data type for the solutions of the quadratic case. So we can waive the last branching concerning the complex case. Listing 2.28: Implementation of a 90-Quad-Solver 1 2 3 ! Solver for a quadratic equation ! Implementation in Fortran90 program quadequation 4 5 implicit none ! only explicit declarations real(8)::a, b, c real(8)::d real(8)::p real(8)::x complex(8)::c1,c2 ! ! ! ! ! 6 7 8 9 10 11 parameters of the equation discriminant precision for the linear solutions for the quadratic solutions 12 13 14 15 16 17 ! setup a = b = c = p = the parameters for the quadratic equation 1. 0. 4. 1.d-15 18 19 20 ! print parameters values to the screen write(*,’(3(a,F10.3))’) ’a=’,a,’ b=’,b,’ c=’,c 21 22 23 24 E. Baeck ! case a=0: it’s not a quadratic equation! if (dabs(a) < p) then 2.9. BRANCHING 25 26 Page 47 ! subcase b=0: => infinit solutions or no solution if (dabs(b) < p) then 27 ! subcase c=0: Trival case => infinit solutions ! it’s independent of x if (dabs(c) < p) then write(*,*) ’Trivial solution, infinit solutions for x.’ 28 29 30 31 32 ! subcase c!=0: Trivial solution => no solution ! it’s independent of x else write(*,*) ’No solution found.’ 33 34 35 36 37 endif 38 39 40 41 42 43 ! subcase b!=0 => linear case -> one solution else x = -c/b write(*,’(a,f10.4)’) ’Linear case: x=’,x 44 45 endif 46 47 48 49 50 ! a!=0 => quadratic problem -> two solutions ! if we solve the problem with reals, we have to handle ! the real and the complex subcase. else 51 52 53 ! calculate the discriminant to make the case check d = b**2 -4.*a*c 54 55 56 c1 = 1./(2.*a)*(-b +cdsqrt(dcmplx(d,0.d0))) c2 = 1./(2.*a)*(-b -cdsqrt(dcmplx(d,0.d0))) 57 58 59 60 61 ! positive discriminant -> 2 real roots if (d >= 0.) then write(*,’(a,f10.4," +i",f10.4)’) ’Quadratic real case, x1=’,c1 write(*,’(a,f10.4," +i",f10.4)’) ’ x2=’,c2 62 63 64 65 66 67 ! negative discriminant -> 2 complex roots else write(*,’(a,f10.4," +i",f10.4)’) ’Quadratic complex case, x1=’,c1 write(*,’(a,f10.4," +i",f10.4)’) ’ x2=’,c2 endif 68 69 70 endif end program In line 55 the function cdsqrt is applied. The leading character c specifies that it is a complex function. The second character d specifies the kind of the real respectively the imaginary part (see section 2.5.2). 9.6.2015 Page 48 2.10 Computer Languages for Engineering - SS 15 Subroutines and Functions A very important feature of a programming language is the possibility to encapsulate code into reusable packages. Such a package is called in Fortran Subroutine or Function. The only difference between a function and a subroutine is the return value of the function. So a function can be called within an expression like sin(ϕ) or cos(ϕ). A function as well as a subroutine in general receives a list of parameters, which are called formal parameters. A parameter can be used to pass information from the calling program into the function or the subroutine. Then the parameter is called input parameter. A parameter can be used as well to pass information out of the function or the subroutine into the calling program. Then the parameter is called output parameter. 2.10.1 Functions The implementation of a function is given in Fortran66/77 as follows. Listing 2.29: Syntax of a 66/77-Function 1 2 3 4 5 6 7 c234567 <type> FUNCTION <name> ([<Parameter list>]) [<Declaration Block>] [<Code Block>] <name> = <return value> RETURN END The implementation of a function is given in Fortran90+ as follows. Note that the end of a function is set by the statement end function. Listing 2.30: Syntax of a 90-Function 1 2 3 4 5 6 <type> FUNCTION <name> ([<Parameter list>]) [<Declaration Block>] [<Code Block>] <name> = <return value> RETURN END FUNCTION [<name>] The function Test1 in listing 2.31 calculates the function value of a line. The lines parameter and the x-value are passed by the list of the formal parameters. Listing 2.31: A Function and it’s Testing Environment 1 2 3 ! Main program as testing environment for function calls program functions implicit none 4 5 6 7 8 9 10 11 E. Baeck real(8)::Test1 real(8)::p1,p2,x1 real(8)::x0 real(8)::xD real(8)::t integer::i ! ! ! ! ! ! function’s return data type the function’s parameters initial value increment function’s value loop counter 2.10. SUBROUTINES AND FUNCTIONS Page 49 ! initialize x0 and xDel x0 = -2. xD = 0.25 12 13 14 15 ! set function parameters p1 = 2. p2 = -1. x1 = x0 16 17 18 19 20 do i=1,16 t = Test1(p1,p2,x1) ! calculating the function value write(*,’(2(a,f12.6))’)’ x=’,x1,’ f(x)=’,t x1= x1 + xD ! increment the function parameter end do 21 22 23 24 25 26 27 end program 28 29 30 31 32 33 34 ! function to calculate some values real(8) function Test1(a,b,x) implicit none ! we have to declare everything explicitly real(8)::a,b,x ! declaring the parameters of the function Test1 = a*x +b ! calculating and asigning the return value end function Test1 ! the end of function Test1 2.10.2 Subroutines The implementation of a subroutine is given in Fortran66/77 as follows. You see, there is no return value. The only difference between subroutine and function is the keyword subroutine instead of function, the missing return type, and the missing assignment of the return value. Listing 2.32: Syntax of a 66/77-Subroutine 1 2 3 4 5 6 c234567 SUBROUTINE <name> ([<Parameter list>]): [<Declaration Block>] [<Code Block>] RETURN END The implementation of a subroutine is given in Fortran90+ as follows. Note that the end of a subroutine is set by the statement end subroutine. Listing 2.33: Syntax of a 90-Subroutine 1 2 3 4 5 SUBROUTINE <name> ([<Parameter list>]) [<Declaration Block>] [<Code Block>] RETURN END SUBROUTINE [<name>] 9.6.2015 Page 50 2.10.3 Computer Languages for Engineering - SS 15 Functions as Parameters if a function should be used as a parameter, the functions should be declared with the return data type and the attribute external. A typical and nice example is the implementation of Newton’s algorithm to calculate the roots of an arbitrary equation (see section 3.5). The following example shows the implementation of the numerical calculation of a function’s derivative, which is used within the Newton’s algorithm. The function is passed as parameter to the derivative function fs. Within the function’s code the function f is declared as a real(8) function. To distinguish a function from a variable the function’s symbolic name can be declared with an external attribute. Listing 2.34: Passing a Function as a Parameter 1 2 3 4 5 6 function to calculate the deviation of a function real(8) function fs (f,x,h) real(8), external:: f real(8):: x,h fs = (f(x +h/2) - f(x -h/2))/h end function Subroutines may also be passed to procedures as calling arguments. if a subroutine is to be passed as a calling argument, it must be declared in an external statement. The corresponding dummy argument should appear in a call statement in the procedure. 2.11 Arrays An array is a compound of data of the same type. The items of the array are addressed by an index value. A static array is declared by the definition of the data type and the index range of an array. 2.11.1 Static Array An static array is declared in FORTRAN77 with the following statements. One really big problem in FORTRAN77 is, that there are only static arrays, i.e. the developer has to decide about the size of an array. If the array size is to small, the code must be recompiled. So a FORTRAN77 software in general is not able to fit to the problems size. Listing 2.35: Array Declaration in 66/77 1 2 3 4 5 6 <data type> <arrayname> DIMENSION <arrayname> (<range-1>,<range-2>,..,<range-n>) ... ... or ... ... <data type> <arrayname> (<range-1>,<range-2>,..,<range-n>) E. Baeck 2.11. ARRAYS Page 51 Listing 2.36 shows a little example code to calculate the scalar product of the vectors v1 and v2. The vectors are allocated using static arrays. Listing 2.36: Allocating Static Arrays in FORTRAN66/77 1 2 3 4 5 c1234567 real*8 v1(3), v2(3), s c ... initialize the variables do 100 i=1,3 100 s = s + v1(i)*v2(i) In Fortran90 we declare a statical array with the following format. Listing 2.37: Array Declaration in 90 1 <data type>, dimension(<range-1>,<range-2>,...,<range-n>):: <arrayname> Listing 2.38 shows the example 2.36 in a FORTRAN90 style. Listing 2.38: Allocating Static Arrays in FORTRAN90 1 2 3 4 5 6 7 real(8), dimension (3)::v1 ! vector 1 real(8), dimension (3)::v2 ! vector 2 real(8)::s ! result value ! ... initialize the variables do i=1,3 s = s + v1(i)*v2(i) end do 2.11.2 Dynamical Array A dynamical array can be allocated, i.e. created at run time. So first we can evaluate the necessary array size and then we can allocate the used memory for the array. This feature is available starting from FORTRAN90. Listing 2.39: Dynamical Array Declaration only in 90 1 2 3 4 5 6 7 8 9 <data type>, allocatable, dimension(:,:,...,:):: <arrayname> ... ... next step we allocate the array ... allocate(<arrayname>(<range-1>,<range-2>,...,<range-n>) [,stat=<statname>]) ... ... after the usage we deallocate the memory ... deallocate(<arrayname>,[stat=<statname>]) After having declared the array name, the array can be allocated by the allocate statement. After the allocation the array items can be accessed. If an array item is accessed before the array is allocated, the program in general will crash. If the memory of an dynamical array is no longer needed, the array should be deallocated with the deallocate statement. Example 2.40 shows how to allocate and deallocate a double indexed array dynamically. We also see a memory handler, which prevents crashing in the case of allocation problems. 9.6.2015 Page 52 Computer Languages for Engineering - SS 15 Listing 2.40: Array Declaration in 90 1 2 real(8),allocatable,dimension(:,:)::a integer::nDim ! declar a as a dynmical array ! this variable controls the array dimension 3 4 5 6 7 8 9 10 11 ! allocate the array a, an allocation error is handled nDim = 3 allocate(a(nDim,nDim),stat=memstat) if (memstat /= 0) then write (*,*) ’*** Error: array a is not allocatable.’ else deallocate(a,stat=memstat) end if 2.11.3 Automatic Array An automatic array will be created automatic in a function or in a subroutine. If the function or subroutine is exited the automatic array is deallocated automatically. The dimensions of an automatic array are passed into the function or the subroutine as formal parameters. Example 2.41 shows the usage of an automatic array. This feature is only available in FORTRAN90 and newer. If we compare this to an implementation in FORTRAN66/77 we would see, that in the old FORTRAN situation, the array has to be allocated in the main program as a static array and has to be passed by a parameter to the subroutine. In FORTRAN90++ the array is allocated within the scope of the sub program automatically without any allocation code. Leaving the sub program the automatic array automatically is deallocated without any additional code. Listing 2.41: Usage of an Automatic Array in FORTRAN90 1 2 3 4 5 program AutoArray implicit none integer::nDim = 3 call UseAutoArray(nDim) end program 6 7 8 9 subroutine UseAutoArray(nDim) integer::nDim, i real(8),Dimension(nDim)::a; 10 ! set automatic array do i=1,nDim a(i) = i end do ! print content of automatic array write(*,*) (a(i),i=1,nDim) 11 12 13 14 15 16 17 18 end subroutine E. Baeck 2.11. ARRAYS 2.11.4 Page 53 A little Array Example The following code shows how to work with static, dynamic and automatic array in a more complex situation. Listing 2.42: Static, Dynamic and Automatic Arrays in FORTRAN90 1 2 3 4 ! This example shows the 3 types of array available in ! Fortran 90++ program Arrays implicit none 5 6 7 8 9 10 11 integer:: integer:: integer:: integer:: integer:: integer:: i,j nDim memstat ioerr ionr = 10 nDim1,nDim2 ! ! ! ! ! ! loop counter used as matrix dimension used as memory error flag used as file io error flag channel number dimensions of the matrx in file 12 13 14 ! if we use functions, we have to declare their retruns integer:: iwritemat,ireadmatdim,ireadmat 15 16 17 real(8),dimension(3,3)::a real(8),allocatable,dimension(:,:)::b ! static array ! dynmical array 18 19 character(256)::logname ! name of the output file and input file logname = "arrays.log" ! initialize the filename ! note: the file is written into the ! projects folder 20 21 22 23 24 25 26 27 28 29 30 31 ! open the log file as a new blank file open(ionr,file=logname,status=’replace’,iostat=ioerr) ! .ne. if (ioerr /= 0) then write (*,*) ’*** Error: log file not opened!’ stop endif 32 33 34 35 36 37 38 ! allocate the array b, an allocation error ist handled nDim = 3 allocate(b(nDim,nDim),stat=memstat) if (memstat /= 0) then write (*,*) ’*** Error: array b is not allocatable.’ end if 39 40 41 42 43 44 ! allocation check: if b is not allocated, we stop if (.not. allocated(b)) then write (*,*) ’*** Error: array b not allocated’ stop end if 45 46 47 48 ! initialize array a and b with a special number pattern ! - over the rows (1st index) do i=1,3 49 9.6.2015 Page 54 50 51 52 53 54 Computer Languages for Engineering - SS 15 ! - over the columns do j=1,3 a(i,j) = i*10 +j b(i,j) = i*10 +j +100 ! +100, because we want to end do ! know, thats the b 55 56 end do 57 58 59 60 ! print array data of a and b to the sceen write(*,’(3(f10.3,1x))’) ((a(i,j),j=1,3),i=1,3) write(*,’(3(f10.3,1x))’) ((b(i,j),j=1,3),i=1,3) 61 62 63 64 65 ! and ! for ioerr ioerr write the array data of a and b into the log file later readings = iwritemat(ionr,a,3,3) = iwritemat(ionr,b,3,3) 66 67 68 69 70 ! if not longer needed, free the memory of array b deallocate(b,stat=memstat) ! close log file close(ionr) 71 72 73 74 75 76 77 ! open the log file to read the data of the first matrix open(ionr,file=logname,status=’old’,iostat=ioerr) if (ioerr /= 0) then write(*,*) ’*** Error: file ’,logname,’ not found!’ stop endif 78 79 80 ! read the matrix dimension if (ireadmatdim(ionr,nDim1,nDim2) == 0) then 81 82 83 84 ! Check the dimensions: size < 1 is not valid if (nDim1 < 1 .or. nDim2 < 1) then write (*,*) ’*** Error: wrong dimensions ’,nDim1,’ ,’,nDim2 85 86 87 88 89 90 91 ! dimensions ok else ! now we reallocate the array b allocate(b(nDim1,nDim2),stat=memstat) ! and read the matrix data from the file ioerr = ireadmat(ionr,b,nDim1,nDim2) 92 93 endif 94 95 96 97 98 99 else ! wrong format -> close the file and stop it close(ionr) stop endif 100 101 102 ! close the input file close(ionr) 103 104 105 E. Baeck ! now we print the read data into the screen write(*,*) ’Data of the first matrix in file:’,logname 2.11. ARRAYS 106 107 108 Page 55 do i=1,nDim1 write(*,’(10(f10.3,1x))’) (b(i,j),j=1,3) enddo 109 110 111 ! and deallocate the matrix b deallocate(b,stat=memstat) 112 113 114 115 116 117 ! the usage of an automatic array of the dimension 4x5 ! is shown in the next call. Only the dimension of the array ! ist passed, the array is allocated automatically in the ! suboutine call checkautomat(4,5) 118 119 end program Arrays 120 121 122 123 124 ! print matrix data into a file ! integer function iwritemat(io,m,n1,n2) implicit none 125 126 127 128 129 130 integer::io ! io channel number integer::n1 ! number of rows integer::n2 ! number of columns real(8), dimension(n1,n2):: m ! matrix data integer::ioerr, i, j 131 132 133 134 135 136 137 write(io,*,iostat=ioerr) n1,n2 ! write the dimensions if (ioerr /= 0) then write(*,*) ’*** Error: writing not possible’ iwritemat = -1 ! exit, if io error return endif 138 139 140 141 142 ! over the rows do i=1,n1 write(io,*) (m(i,j),j=1,n2) enddo 143 144 145 iwritemat = 0 end function iwritemat ! 0 return: everything ok 146 147 148 ! Read the dimension of a matrix from a file integer function ireadmatdim(io,n1,n2) 149 150 151 152 integer::io integer::n1,n2 integer::ioerr ! io channel ! dimensions of the matrix ! error flag 153 154 155 156 157 158 159 read(io,*,iostat=ioerr) n1,n2 if (ioerr /= 0) then ! if io-error, perhaps a wrong format write(*,*) ’*** Error: wrong file format’ ireadmatdim = -1 return endif 160 161 ireadmatdim = 0 9.6.2015 Page 56 162 Computer Languages for Engineering - SS 15 end function ireadmatdim 163 164 165 166 167 ! read matrix data from a file ! integer function ireadmat(io,m,n1,n2) implicit none 168 integer::io ! io channel number integer::n1 ! number of rows integer::n2 ! number of columns real(8), dimension(n1,n2):: m ! matrix data integer::ioerr, i, j 169 170 171 172 173 174 ! over the rows do i=1,n1 read(io,*,iostat=ioerr) (m(i,j),j=1,n2) if (ioerr /= 0) then ! important to check the read-io write(*,*)’*** Error: format’ ireadmat = -1 return endif enddo 175 176 177 178 179 180 181 182 183 184 185 186 ireadmat = 0 end function ireadmat 187 188 189 ! example for an automic array subroutine checkautomat(nDim1,nDim2) 190 real(8), dimension(nDim1,nDim2)::m 191 ! automatic array 192 ! initialize the array with a number pattern do i=1,nDim1 do j=1,nDim2 m(i,j) = i*10 +j enddo enddo 193 194 195 196 197 198 199 ! print the pattern to the screen write(*,*) ’M:’,nDim1,’,’,nDim2 do i=1,nDim1 write(*,’(10(f10.3,1x))’) (m(i,j),j=1,nDim2) enddo 200 201 202 203 204 205 206 end subroutine 2.11.5 Pseudo Dynamic Arrays in FORTRAN77 If we need a dynamic array using FORTRAN77 the only chance to implement this is, to use a static memory buffer, i.e. a static array which has to be large enough to hold the largest dimension of our pseudo dynamic array. How to implement this we can see in listing 2.43. The dimension of the vectors are read in from a text input file. E. Baeck 2.11. ARRAYS Page 57 Listing 2.43: Dot Product of Vectors using a Pseudo Dynamic Array 1 2 3 c implementing a pseudo dynamic array in FORTRAN 77 c implicit none 4 5 6 7 c 8 integer nDim,i real*8 GetScalProd memory buffer real*8 dBuffer(20) 9 10 c 11 12 13 c c 14 15 16 c c 17 18 19 c c 20 open input file open(10,file=’SkalProd.in’,status=’old’) read the dimension read(10,*) nDim read the first vector read(10,*) (dBuffer(i),i=1,nDim) read the second vector read(10,*) (dBuffer(i),i=nDim+1,nDim*2) 21 22 23 24 25 26 27 28 29 30 c print out write(*,*) ’>> scal product of 2 vectors’ write(*,9000) ’v1 = ’,(dBuffer(i),i=1,nDim) write(*,9000) ’v2 = ’,(dBuffer(i),i=nDim +1,nDim*2) write(*,9000) ’v1*v2 = ’, &GetScalProd(dBuffer(1),dBuffer(nDim+1),nDim) close(10) 9000 format(a,20f10.3) end 31 32 33 34 c calculation of the dot product of two vectors real*8 function GetScalProd(a,b,n) implicit none 35 36 37 integer i,n real*8 a(1),b(1) 38 39 40 41 42 GetScalProd = 0. do 100 i=1,n 100 GetScalProd = GetScalProd + a(i)*b(i) end As we can see from line 20 of listing 2.43 the size of the buffer has to be set large enough. If not, the input data will be read into a memory outside of our buffer, which can produce a lot of ugly side effects. The second problem will occur, if the pointer calculation is not perfect (line 17, 20 and 27). If we use incorrect pointers, internal data can be overwritten without provoking any error situation. Side effects like this are very hard to find and can be avoided using allocatable arrays with FORTRAN90. 9.6.2015 Page 58 2.12 Computer Languages for Engineering - SS 15 Global Data Global data in Fortran are handled with specific access statements. 2.12.1 Classical Fortran and Common Global data in Fortran classically are handled with so called common blocks. A common block is a block of memory, which can be used from all routines, which are permitted to do this. A routine will be permitted to access a common block, if the common block is included into this routine with the statement common. Global data can be initialized with the block data statement. Listing 2.44: common Block and block data 1 2 3 4 c initialization of a common c c | block data’s name block data global 5 6 7 c c 8 9 10 name of the common ! | start with longest datatype common /old77/ dOld,nOld real*8 dOld integer nOld 11 12 13 data nOld /123/ data dOld /3.14/ 14 15 end In listing 2.44 we see, that global data are introduced with a common statement. The name of the common in this case is old77. If this statement and the following declarations (lines 8 to 10) are found within a subroutine ore a function, this common variables are available in terms of global data. The block data statement, which can be only once in a code, will initialize the variables of a common block. In this case we assign a value to nOld and dOld. 2.12.2 Some Aspects of the Module Concept of FORTRAN90 Using FORTRAN90, the classical concept of common blocks should be considered as obsolete. Common blocks can be considered as a source of many possible errors and side effects. With FORTRAN90 common blocks can be substituted by modules. A module in FORTRAN90 is a compound of data and methods according to the object orientated concept of modern languages. So using modules we also can develop software in FORTRAN using modern OOP concepts.12 12 OOP is discussed later in the C++ section, see section 6.11. E. Baeck 2.12. GLOBAL DATA Page 59 In listing 2.45 some global constants are introduced and initialized. Further a method is implemented inside the contains block, to print this constants. Listing 2.45: A module to Handle Some Constants 1 2 3 ! global data in FORTRAN 90 ! making common and block data obsolete module constants 4 implicit none 5 6 ! data section real, parameter::e = 2.7 character(*), parameter::room = "V15-S03-D03" integer ::nrtel= 2613 7 8 9 10 11 ! methodes section contains 12 13 14 subroutine printConstants 15 16 write(*,*) write(*,*) 17 18 "my room..",room "my telnr.",nrtel 19 end subroutine printConstants 20 21 22 end module constants 2.12.3 Using global Data A really big benefit in FORTRAN’s history is, that old FORTRAN code can be used in modern FORTRAN environments with nearly no required changes. This we can see within the next example, which uses the classical common block of listing 2.44 and the modern module of listing 2.45. Listing 2.46: Using commons and modules 1 2 3 ! example to show the usage of common and module ! within one 90 code program GlobalData 4 5 use constants 6 7 implicit none 8 9 10 11 12 13 14 ! insert the common-code here ! name of the common ! | start with longest datatype common /old77/ dOld,nOld real*8 dOld integer nOld 15 16 17 ! print the global common data call printGlobals 18 19 ! change the global common data 9.6.2015 Page 60 Computer Languages for Engineering - SS 15 nOld = 321 dOld = 4.13 20 21 22 ! print the global common data call printGlobals ! print the global module data call printConstants 23 24 25 26 27 28 end program 29 30 31 32 ! subroutine to print the common data ! to this we have to insert the common code subroutine printGlobals 33 common /old77/ dOld,nOld real*8 dOld integer nOld 34 35 36 37 write(*,*) "nOld = ",nOld write(*,*) "dOld = ",dOld 38 39 40 41 end subroutine E. Baeck 2.13. MEMORY MAPPING 2.13 Page 61 Memory Mapping Especially in old FORTRAN66/77 sources the technique of momory mapping using the equivalence statement is applied. The problem which should be solved is a dynamic memory management. In contrast to modern languages like C or C++ there is no possibility of dynamical memory allocation. If a program without dynamical memory management should be efficient the only solution is, to allocate a statical memory buffer and try to map all large arrays onto this buffer. Is we can do with the equivalence statement. In FORTRAN90 sources, where the option of dynamical memory allocation is available, we should not use memory mapping, because it’s highly prone to errors. We easily can overwrite data simple by mapping errors, which are only found in general by an incredibly high effort of testing. The equivalence statement links a list of source memory to their destination memory. Listing 2.47: Equivalence Statement in FORTRAN66/77 1 equivalence (<destination memory>, <source memory>), [(...)] So in example 2.48 we introduce a real*8 buffer, which is used as a memory mapping array. On the other hand we allocate two real*8 arrays v1 and v2 and map them onto the buffer. After having initialized them we calculate the scalar product of this vectors and print them into the console window. After this we initialize the integer*4 vectors i1 and i2, which have the same length as the real*8 vectors and were also allocated inside the start up section of the program. We calculate the scalar product of them and print their values. After this we again print the content of the real*8 vectors and we can see, that their content totally is overwritten by the integer*4 vectors. The values, which are printed don’t make any sense, because it’s a real interpretation of an integer bit pattern. Listing 2.48: Equivalence of Real and Integer Vectors 1 2 3 4 5 c1234567 c allocate the memory integer ndim parameter (ndim = 3) real*8 buff(nDim*2) 6 7 c some real memory real*8 v1(ndim), v2(ndim), s c some integer memory integer*4 i1(ndim*2), i2(ndim*2), is c map v1, v2, s onto the memory buffer equivalence(buff(1),v1), (buff(ndim+1),v2), (buff(ndim*2+1),s) c map i1, i2, is onto the memory buffer equivalence(buff(1),i1), (buff(ndim+1),i2), (buff(ndim*2+1),is) c initialize the buffer do 100 i=1,ndim*2 buff(i) = i 8 9 10 11 12 13 14 15 16 17 18 19 20 100 21 22 23 c print the vector’s content 9.6.2015 Page 62 Computer Languages for Engineering - SS 15 write (*,’(a,3f10.3)’) "v1 =",(v1(i),i=1,ndim) write (*,’(a,3f10.3)’) "v2 =",(v2(i),i=1,ndim) 24 25 26 27 c 28 29 200 30 31 performe the scalar product and print the result s = 0 do 200 i=1,ndim s = s + v1(i)*v2(i) write (*,’(a,f10.3)’) "s =",s 32 33 c 34 35 300 36 set integer vectors do 300 i=1,ndim*2 i1(i) = i*2-1 i2(i) = i*2 37 38 c print the vector’s content write (*,’(a,6i5)’) "i1 =",(i1(i),i=1,ndim*2) write (*,’(a,6i5)’) "i2 =",(i2(i),i=1,ndim*2) c performe the scalar product and print the result s = 0 do 400 i=1,ndim*2 is = is + i1(i)*i2(i) write (*,’(a,i5)’) "is =",is 39 40 41 42 43 44 400 45 46 47 48 49 c c 50 51 print the real vector’s content again we see, that the values are obviously overwritten write (*,’(a,3e10.3)’) "v1 =",(v1(i),i=1,ndim) write (*,’(a,3e10.3)’) "v2 =",(v2(i),i=1,ndim) 52 53 end Figure 2.8 shows the output of example 2.48. Figure 2.8: Output listing E. Baeck 2.14. COMMANDLINE ARGUMENTS 2.14 Page 63 Commandline Arguments In old FORTRAN there is by standard no possibility to have access to the commandlines’s parameter because on the mainframe IBM computers everything concerning the runtime environment was handled by the so called Job Control Language. With the advent of the command shells like BASH (Born Again SHell) on LINUX systems, programs are generally started with a wide set of commandline parameters which are passing input data to the startup of the application (see section A too). Because the FORTRAN90++ is highly influenced by the C language the concept of passing commandline arguments is very close to C’s strategy. In C the commandline parameters are parameters of the main function, i.e. the main program. Unlike C FORTRAN90 provides a function which returns the number of given parameters called IARGC(). With the given number of parameters we can get the value of each parameter by calling the subroutine GETARG(I,ARG), where I is the parameter’s index (starting from 1) and ARG is a CHARACTER variable, which will come with the desired value on return. Like in C we get the program’s name with an index value of 0. The following example 2.49 shows how to access this commandline parameters. Listing 2.49: A Commandline Example 1 2 program main implicit none 3 4 5 6 7 character integer ( integer ( integer ( ( len = 255 ) arg kind = 4 ) i kind = 4 ) iargc kind = 4 ) numarg 8 9 10 11 numarg = iargc ( ) write ( *, ’(a,i8,a)’ ) & ’ Program executed with ’, numarg, ’ commandline options’ 12 13 14 15 16 17 write write write write write ( ( ( ( ( *, *, *, *, *, ’(a)’ ’(a)’ ’(a)’ ’(a)’ ’(a)’ ) ) ) ) ) ’ ’ ’ Found commandline options’ ’ ’ ’ I ARG ’ ’ ’ 18 19 20 21 22 23 do i = 0, numarg call getarg ( i, arg ) write ( *, ’(2x,i3,2x,a20)’ ) i, arg end do end 9.6.2015 Page 64 Computer Languages for Engineering - SS 15 Figure 2.9: Running the Commandline Example We see in figure 2.9, that we run the application commandline with three commandline options called option1, option2 and option3. This options we access within our program and print the found option values into the screen. E. Baeck 3 Some Examples 3.1 Hello World One famous application which does n’t make any sense is the program helloworld. There are only two statements: the first writes the famous text to the screen, the second closes the application. Listing 3.1: A Startup Hello 1 2 3 4 3.2 c234567 c comment write(*,*) "Hello World " end 70 12345 Simple Sum The second example shows the implementation of a simple loop in FORTRAN66 style. The result of the loop (do-loop with labeled end) is the sum of all integers from 1 to 10. Each step is printed to the screen. S= 10 X i (3.1) i=1 Listing 3.2: Sum up all Numbers from 1 to 10 1 2 3 4 5 6 7 c234567 n = 0 ! sum variable, set to zero do 100 i=1,10 ! performing the sum in fortran IV style n = n + i write (*,’(a,i2,a,i4)’) ’ i = ’,i,’ sum = ’,n ! screen dump 100 continue ! end of loop end ! end of application 65 Page 66 Computer Languages for Engineering - SS 15 The screen output running example 3.2 will be the following. i i i i i i i i i i 1 2 3 4 5 6 7 8 9 10 3.3 = 1 = 2 = 3 = 4 = 5 = 6 = 7 = 8 = 9 = 10 sum sum sum sum sum sum sum sum sum sum = = = = = = = = = = 1 3 6 10 15 21 28 36 45 55 Calculation of real*4/8 Precision This example calculates the relative precision of a 4 and 8 byte float arithmetic. In listing 3.3 a strict 66 coding is used, if we forget the line end comment. The idea of this algorithm is, to divide a variable’s value of 1 by 2 as long as the sum of 1 and this reduced value is greater than 1. If we would have an infinite precision, this loop would be an endless loop. Because we only have a few digits, this reduced value will vanish with some cycles. The last visible value than will be our relative precision. In figure 3.1 the algorithm to calculate the relative precision is shown. The first part will calculate the relative precision for a 4 byte arithmetic, the second part will calculate the relative precision for the 8 byte arithmetic. Start x1 = 1. x2 = 1. d = 2. x2 = x2 /d s = x1 + x2 no s = x1 yes result = x2 ∗ d Stop Figure 3.1: Algorithm’s Flowchart Listing 3.3: Calculation of the Arithmetic’s Relative Precision 1 2 3 C234567890 real*4 x14, x24, x34, d4 real*8 x18, x28, x38, d8 ! variables for real*4 analysis ! variables for real*8 analysis 4 5 c 6 7 8 calculation of real*4 relative precision x14 = 1. x24 = 1. d4 = 2. 9 10 11 12 100 x24 = x24 /d4 x34 = x14 + x24 c write (*,1001) x34, x24 E. Baeck ! back jump label and increment ! reduction ! dump is disabled 3.4. RELATIVE PRECISION WITH FUNCTIONS 13 14 15 16 17 18 19 c if (x34 .gt. x14) goto 100 x24 = x24 * d4 output write (*,1000) x24 Page 67 ! if increment still seen next run ! prints result to screen using ! a format statment (1000) 1000 format(’ real*4 relative precision: ’,e10.3) 1001 format(’ x14+x24 = ’,e20.14,’ x24 = ’,e20.14) 20 21 22 23 24 25 26 27 28 29 30 31 32 c calculation of real*8 relative precision x18 = 1. x28 = 1. d8 = 2. ! now the same for real*8 200 x28 = x28 /d8 ! arithmetic x38 = x18 + x28 c write (*,2001) x38, x28 ! dump is disabled if (x38 .gt. x18) goto 200 x28 = x28 * d8 c output write (*,2000) x28 33 34 35 36 2000 format(’ real*8 relative precision: ’,e10.3) 2001 format(’ x18+x28 = ’,e20.14,’ x28 = ’,e20.14) end ! If we run this code, we will get the following screen output. 1 2 real*4 relative precision: real*8 relative precision: 0.119E-06 0.222E-15 We see, that with 4 byte real we nearly get 7 digits, for a 8 byte real we nearly get 16 digits. 3.4 Function to Calculate the Relative Precision The following code consists of two routines, the first is the main program, which calls the evaluation function getRelPrec to get the relative precision. The function is working with one integer parameter. If the parameter is set to 0, the function evaluates the 4 byte relative precision, if the parameter is set to any value but not 0, the function evaluates the 8 byte relative precision. 9.6.2015 Page 68 Computer Languages for Engineering - SS 15 The main program (line 1 to 7) is a testing environments and performs the calls. The code of the function is given starting from line 9. Listing 3.4: Function to Evaluate the Relative Precision for 4 and 8 byte floats 1 2 3 4 5 6 7 ! evaluate the relative precision for ! 4 and 8 byte float arithmetic program getRelPrecMain real(8) :: getRelPrec, eps write(*,*) ’4 byte relative precision: ’, getRelPrec(0) write(*,*) ’8 byte relative precision: ’, getRelPrec(1) end program getRelPrecMain 8 9 10 11 12 13 14 ! function to calculate the relative precision real(8) function getRelPrec(nBytes) ! function interface ! return type is real*8, name is getRelPrec integer :: nBytes ! nByte: 0:4 bytes / 1:8 bytes real(4) :: x14,x24,s4,d4 ! 4 byte data real(8) :: x18,x28,s8,d8 ! 8 byte data 15 ! calculation for 4 byte arithmetic if (nBytes == 0) then x14 = 1. x24 = 1. d4 = 2. do ! implicit loop without a counter x24 = x24/d4 s4 = x14+x24 if (s4 <= x14) exit end do eps = x24*d4 ! last division should be canceld ! by multiplication 16 17 18 19 20 21 22 23 24 25 26 27 28 else x18 = 1. x28 = 1. d8 = 2. do x28 = x28/d8 s8 = x18+x28 if (s8 <= x18) exit end do eps = x28*d8 29 30 31 32 33 34 35 36 37 38 39 ! implicit loop without a counter ! last division should be canceld ! by multiplication endif 40 41 42 43 getRelPrec = eps end function getRelPrec E. Baeck ! assigning the return value ! end of function 3.5. NEWTON’S ALGORITHM TO CALCULATE A ROOT 3.5 Page 69 Newton’s Algorithm to calculate a Root The following example shows, how to pass a function as a functions parameter. Within the Newton’s algorithm a root of an equation should be calculated. So we have to specify the function of interest. The function can be considered as an input parameter. The function’s name is passed to the derivative calculator and to the newton main routine. So, if we want to calculate the roots of an equation f (x ) = 0, we can apply the iteration scheme 3.3. The derivative in the denominator is calculated numerically in equation 3.2. We see that in both equations we need the values of the function f . This problem can be solved by passing the function as a parameter. Figure 3.2: Scheme of the Newton Algorithm The derivative - it’s called fs in the code - is calculated numerical as follows. h h df 0 ≈ f (x + ) − f (x − ) /h f (x ) = dx 2 2 (3.2) The Newton scheme can be described as follows. f (x ) xi+1 = xi − 0 f (x ) (3.3) The same formula we get from the triangle of the slope (see figure 3.2) resolving for xn1 . f 0 (xn ) = f (xn ) xn − xn+1 (3.4) There are three possible situations to handle within the iteration loop. • The function value is vanishing with respect to our selected precision. The iteration loop will be broken and the found result is passed back to the caller. • The slope of the function is vanishing. With a vanishing slope we would divide by zero. With an infinitesimal slope we would get a nearly infinite jump length, which in every case would not be helpful. Therefore in this case we simple jump to the side and try it once more. • During the iteration each cycle is counted. So the iteration loop will be broken, if the maximum available iterations are reached. The actual values and an error message is passed bake to the caller. The code consists of a main program which calls the function newton. Within newton the functions f and fs are called. So we have to implement the following functions. • Myf, the function of our interest. • fs, the function which calculates the slope of a given function numerically. • newton, implements the newton scheme. 9.6.2015 Page 70 Computer Languages for Engineering - SS 15 Figure 3.3 shows a flowchart of Newton’s algorithm. First we set the system’s parameter. Then in the iteration loop we calculate the functions value and it’s derivative. In the simple Newton case we handle the three break conditions in an open loop. Start Set Parameters: x0 ; itmax ; ; ... Initializing: x = x0 ; it = 0; The first checks the function value. If the function value is close enough to zero, we have found a root and the cycles will stop. Function Data: F = f (x ); FS = f 0 (x ); The second checks the function’s slope. We should not divide by zero, therefor we break, if the slope comes close to zero. The third checks the iteration break condition, i.e. the whether the maximal allowed iteration number is reached. In this case the loop breaks, because there is no root found. If we survive all break conditions, the next iteration step is introduce, i.e the new x value is calculated and the iteration counter is incremented. it ≥ itmax yes No Solution no |F | < it = it + 1 yes Solution no The code can be separated in two parts or modules. The first module, which is called NewtonMain.f90, contents the specific code, i.e. the main program and a testing function. The second module, which is called Newton.f90, contents the newton scheme and the numerical calculation of the function’s derivative. x = x −s yes |FS | < no F x = x − FS it = it + 1 Figure 3.3: The Newton’s Flowchart Listing 3.5: Testing Environment to Check the Newton Function 1 2 ! Main program to test the implementation ! of newton’s algorithm 3 4 5 program NewtonMain implicit none 6 7 8 9 10 11 12 ! setup the testing parameters real(8):: x0 ! starting value real(8):: eps ! precision integer:: nmax ! maximum number of iterations real(8), external:: Myf ! declaration of the function integer, external:: newton ! declaration of the newton function 13 14 integer:: nret 15 16 E. Baeck real(8):: x ! return of the newton function ! solution ! root 3.5. NEWTON’S ALGORITHM TO CALCULATE A ROOT 17 18 19 real(8):: f0 real(8):: fs0 integer:: nit ! ! ! x0 = 4. eps = 1.e-6 nmax= 100 ! starting position ! the root’s minimal precison ! available iterations Page 71 function’s value slope at root’s position number of used cycles 20 21 22 23 24 25 26 27 28 29 ! print input values write(*,*) ’ Test program for the newton function’ write(*,’(A,F12.4)’) ’ Starting value..........: ’,x0 write(*,’(A,E12.5)’) ’ Precision...............: ’,eps write(*,’(A,I6)’) ’ Maximum number of cycles: ’,nmax 30 31 32 33 34 ! the newton is implemented as a function, which returns a status ! value. The result values are return by the output parameters ! --- input ----- -- output -nret = newton(Myf,x0,eps,nmax,x,f0,fs0,nit) 35 36 37 38 39 ! solution found ! .eq. (F66/77) if (nret == 0) then ! here we use C-like F90 operators write (*,*) ’ Solution found!’ 40 41 42 43 ! error: vanishing slope, avoid to divide by zero else if (nret == 1) then write (*,*) ’ Vanishing slope, no result found!’ 44 45 46 47 48 ! maximum cycles reached. Break to avoid an infinit loop else write (*,*) ’ No solution found, maximum iterations reached!’ end if 49 50 51 52 53 54 ! output section write (*,’(A,F15.8)’) write (*,’(A,F15.8)’) write (*,’(A,F15.8)’) write (*,’(A,I8)’) ’ " " " Solution value....:’,x Function’s value..:",f0 Function’s slope..:",fs0 Used cycles.......:",nit 55 56 end program 57 58 59 60 61 62 63 64 ! user function is an example to tehst newton’s algorithm ! This function is passed to the newton function. real(8) function Myf(x) real(8)::x Myf = x**2 -1 return end function Myf The second module contents the more general code, i.e. the code of Newton’s scheme and the derivatives calculator. General it’s recommended to encapsulate the general code1 into separate modules. This modules can also be packed into library files2 . 1 2 This is code, which is general applicable and therefore has no dependence with your application A library file contents compiled module code and can be linked without any compilation to an application. 9.6.2015 Page 72 Computer Languages for Engineering - SS 15 Listing 3.6: Simple Newton Function 1 2 3 ! implementation of Newton’s algorithm integer function newton (f,x0,e,ix,x,fx,fsx,nx) implicit none 4 5 6 real(8), external:: f real(8), external:: fs ! user function ! deviation calculator real(8)::x0 real(8)::e real(8)::h real(8)::x real(8)::fx real(8)::fsx real(8)::s = 1.0 integer::ix,nx ! ! ! ! ! ! ! 7 8 9 10 11 12 13 14 15 start value root precision deviation step width result: root value function’ value at x function’s deviation at x jump length 16 ! initialization section x = x0 ! initialize the iteration variable nx = 1 ! iteration counter h = e ! set step width for numerical deviation 17 18 19 20 21 ! iteration loop (note it’s a named loop) mainloop: do 22 23 24 fx = f(x) fsx= fs(f,x,h) 25 26 ! calculating the function’s value ! and the function’s solpe 27 ! check the number of cycles, if exceeded return with error if (nx == ix) then newton = 2 return 28 29 30 31 32 ! check the function value, if success return else if (dabs(fx) < e) then newton = 0 return 33 34 35 36 37 ! check the slope, if vanishing, jump to the side else if (dabs(fsx) < e) then x = x + s 38 39 40 41 ! calculating the next x value else x = x - fx/fsx end if 42 43 44 45 46 ! count the cycle nx = nx +1 47 48 49 end do mainloop 50 51 52 end function newton 53 54 ! function to calculate the deviation of a function E. Baeck 3.5. NEWTON’S ALGORITHM TO CALCULATE A ROOT 55 56 57 58 59 60 Page 73 ! note, that a vanishing step width is not handled real(8) function fs (f,x,h) real(8), external:: f real(8):: x,h fs = (f(x +h/2) - f(x -h/2))/h end function The second version of the Newton program will be extended by a function iwritefunction, which should print the function’s values and the derivative of the function in a given range. We extend the main module by a log file newtonlog.txt and a static array for the iteration path. Before we call the newton function the function iwritefunction will be called to print the function values. The allocated array is passed to the newton function to get the iteration path data. Listing 3.7: Testing Environment to Check the Newton Function 1 2 ! Main program to test the implementation ! of newton’s algorithm 3 4 5 program NewtonMain implicit none 6 7 8 9 10 11 12 ! setup the testing parameters real(8):: x0 ! starting value real(8):: eps ! precision integer:: nmax ! maximum number of iterations real(8), external:: Myf ! declaration of the function integer, external:: newton ! declaration of the newton function 13 14 integer:: nret 15 16 17 18 19 20 21 22 real(8):: x real(8):: f0 real(8):: fs0 integer:: nit character(256)::filename integer::ioerr integer::iwritefunction ! ! ! ! ! ! ! ! ! return of the newton function solution root function’s value slope at root’s position number of used cycles name of the log file return of the write function return value of the function 23 24 25 ! F66/77 version static array real(8), dimension(100)::xp ! iteration path, x values 26 27 28 29 x0 = 0. eps = 1.e-6 nmax= 100 ! starting position ! the root’s minimal precison ! available iterations 30 31 filename = ’newtonlog.txt’ 32 33 34 35 36 37 ! print input values write(*,*) ’ Test program for the newton function’ write(*,’(A,F12.4)’) ’ Starting value..........: ’,x0 write(*,’(A,E12.5)’) ’ Precision...............: ’,eps write(*,’(A,I6)’) ’ Maximum number of cycles: ’,nmax 38 39 40 ! print the function’s values into the log file ioerr = iwritefunction(filename,Myf,-10.D0,10.D0,0.5D0) 9.6.2015 Page 74 Computer Languages for Engineering - SS 15 if (ioerr == 0) then write(*,’(a)’) ’ No problems writing functions values.’ else write(*,’(a,i10)’) " Error writing function values, code=",ioerr endif 41 42 43 44 45 46 ! the newton is implemented as a function, which returns a status ! value. The result values are return by the output parameters ! --- input ----- -- output -nret = newton(Myf,x0,eps,nmax,x,f0,fs0,nit,xp) 47 48 49 50 51 ! solution found ! .eq. (F66/77) if (nret == 0) then ! here we use C-like F90 operators write (*,*) ’ Solution found!’ 52 53 54 55 56 ! error: vanishing slope, avoid to divide by zero else if (nret == 1) then write (*,*) ’ Vanishing slope, no result found!’ 57 58 59 60 ! maximum cycles reached. Break to avoid an infinit loop else write (*,*) ’ No solution found, maximum iterations reached!’ end if 61 62 63 64 65 ! output section write (*,’(A,F15.8)’) write (*,’(A,F15.8)’) write (*,’(A,F15.8)’) write (*,’(A,I8)’) 66 67 68 69 70 ’ " " " Solution value....:’,x Function’s value..:",f0 Function’s slope..:",fs0 Used cycles.......:",nit 71 72 end program 73 74 75 76 77 78 79 ! user function is an example to tehst newton’s algorithm ! This function is passed to the newton function. real(8) function Myf(x) real(8)::x ! Myf = x**2 -1 ! Test 1 80 81 82 83 84 ! Myf = x**2 +1 Myf = x**3 +1 return end function Myf ! Test 2 ! Test 3 The array for the storage of the iteration path data is passed to the newton function. Within the mainloop the positions on the iteration path are saved into the array xpos. At the end of the module the function iwritefunction is added to print the function’s values. E. Baeck 3.5. NEWTON’S ALGORITHM TO CALCULATE A ROOT Page 75 Listing 3.8: Simple Newton Function 1 2 3 ! implementation of Newton’s algorithm integer function newton (f,x0,e,ix,x,fx,fsx,nx,xpos) implicit none 4 5 6 real(8), external:: f real(8), external:: fs ! user function ! deviation calculator real(8)::x0 real(8)::e real(8)::h real(8)::x real(8)::fx real(8)::fsx real(8)::s = 1.0 integer::ix,nx real(8),dimension(ix)::xpos ! ! ! ! ! ! ! 7 8 9 10 11 12 13 14 15 16 start value root precision deviation step width result: root value function’ value at x function’s deviation at x jump length 17 18 19 20 21 ! initialization section x = x0 ! initialize the iteration variable nx = 1 ! iteration counter h = e ! set step width for numerical deviation 22 23 24 ! iteration loop (note it’s a named loop) mainloop: do 25 26 27 fx = f(x) fsx= fs(f,x,h) ! calculating the function’s value ! and the function’s solpe 28 29 30 ! be sure, that the array is dimensioned properly xpos(nx) = x 31 32 33 34 35 36 ! check the number of cycles, if exceeded return with error if (nx == ix) then newton = 2 return 37 38 39 40 41 ! check the function value, if success return else if (dabs(fx) < e) then newton = 0 return 42 43 44 45 ! check the slope, if vanishing, jump to the side else if (dabs(fsx) < e) then x = x + s 46 47 48 49 50 ! calculating the next x value else x = x - fx/fsx end if 51 52 53 ! count the cycle nx = nx +1 54 9.6.2015 Page 76 Computer Languages for Engineering - SS 15 end do mainloop 55 56 57 end function newton 58 59 60 61 62 63 64 65 ! function to calculate the deviation of a function ! note, that a vanishing step width is not handled real(8) function fs (f,x,h) real(8), external:: f real(8):: x,h fs = (f(x +h/2) - f(x -h/2))/h end function 66 67 68 69 ! write function values to a file ! integer function iwritefunction(name,f,xfrom,xto,xstep) 70 71 72 73 74 75 character(256):: name real(8),external::f real(8)::xfrom real(8)::xto real(8)::xstep ! ! ! ! ! files name function start value end value step width 76 77 78 79 integer::ioerror real(8)::x real(8)::h ! return status ! actual position ! step width calculating the derivative 80 81 82 83 84 85 86 ! check the input parameters if (xstep < 1.e-6) then write (*,*) ’ *** Error: xstep not ok!’ iwritefunction = -1 return endif 87 88 89 90 91 92 93 ! open the file open (10,file=name,status=’replace’,iostat=ioerror) if (ioerror .ne. 0) then iwritefunction = ioerror return ! return if there is an error endif 94 95 96 97 98 ! start with the tables header ! 123456789012345678901234567890 write(10,’(a)’)" x f(x) f’(x)" write(10,’(a)’)"------------------------------" 99 100 101 102 103 104 ! write the function’s values h = 1.e-6 x = xfrom do write(10,’(3(F10.4))’,iostat=ioerror) x,f(x),fs(f,x,h) 105 106 107 108 109 110 E. Baeck ! break the loop, if an error occure if (ioerror.ne.0) exit x = x + xstep if (x > xto) exit end do 3.5. NEWTON’S ALGORITHM TO CALCULATE A ROOT Page 77 111 112 113 ! close the output file close(10) 114 115 116 ! return the error code iwritefunction = ioerror 117 118 end function iwritefunction 9.6.2015 Page 78 3.6 Computer Languages for Engineering - SS 15 Matrix Product with 77-Main and 90-Library Within this section the mix of FORTRAN code of version 77 and 90 are discussed. Because in FORTRAN77 only static arrays are available, a buffer array is introduced. This array is used for the program memory management. For the three matrices A, B and C , which we use, index pointers into the buffer array are used for mapping. The dimensions of matrix A and B are read from an input file. The matrix product is calculated according to the following formula. C =A·B (3.5) with an matrix element Ci,j Ci,j = n X Ai,k · Bk ,j = Ai,k · Bk ,j (3.6) k =1 The last expression follows the Einstein notation, which means that we sum over the indices which ocrur in every term. Listing 3.9: 77 Environment to call subsequent FORTRAN90 Routines 1 2 c Fortran77 example to handle a pseudo dynamical memory c manager based on a buffer array 3 4 5 6 7 c234567 integer buffersize parameter (buffersize = 100) real*8 buffer(buffersize) ! we use a parameter to ! allocate the work buffer ! statical allocation of the buffer 8 9 10 11 integer integer integer ipA ipB ipC ! pointer of array A ! pointer of array B ! pointer of array C integer iret ! return code integer integer integer nDimA(2)! Dimension of arra A nDimB(2)! Dimension of arra B nDimC(2)! Dimension of arra C integer integer ionr ioerr integer ireadmatdim,ireadmat,imatmult 12 13 14 15 16 17 18 19 20 ! io channel number ! error parameter 21 22 23 character*32 24 InpFile 25 InpFile = ’MatMult77.inp’ ionr = 10 26 27 ! fixed input file name 28 29 c 30 31 32 33 E. Baeck open the input file write(*,*) ’> open the file:’,InpFile open(ionr,file=InpFile,status=’old’,iostat=ioerr) if (ioerr .ne. 0) then write(*,*) ’*** Error: open file ’,InpFile 3.6. MATRIX PRODUCT WITH 77-MAIN AND 90-LIBRARY Page 79 stop endif 34 35 36 37 c read dimension of the 1st matrix iret = ireadmatdim(ionr,nDimA) if (iret .ne. 0) goto 900 ! jump to the error exit write(*,*) ’> dimension of array 1 read:’,nDimA(1),’,’,nDimA(2) c read data of array 1 ipA = 1 iret = ireadmat(ionr,buffer(ipA),nDimA) if (iret .ne. 0) goto 900 ! jump to the error exit c and list it’s data call listmat(’Data of matrix 1:’,buffer(ipA),nDimA) c read dimension of the 2nd matrix iret = ireadmatdim(ionr,nDimB) if (iret .ne. 0) goto 900 ! jump to the error exit write(*,*) ’> dimension of array 2 read:’,nDimB(1),’,’,nDimB(2) c read data of array 2 ipB = ipA + nDimA(1)*nDimA(2) iret = ireadmat(ionr,buffer(ipB),nDimB) if (iret .ne. 0) goto 900 ! jump to the error exit c and list it’s data call listmat(’Data of matrix 2:’,buffer(ipB),nDimB) 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 c 64 65 66 67 68 69 70 multiply matrix 1 with matrix 2 ipC = ipB +nDimB(1)*nDimB(2) iret = imatmult(buffer(ipA),buffer(ipB),buffer(ipC), & nDimA,nDimB) if (iret .ne. 0) then write(*,*) ’*** Error: wrong dimensions for product’ goto 900 ! jump to the error exit endif 71 72 c print the result nDimC(1) = nDimA(1) nDimC(2) = nDimB(2) call listmat(’Data of product matrix 1x2:’,buffer(ipC),nDimC) c no problems therefore jump to the regular end goto 999 73 74 75 76 77 78 79 80 c error exit 900 write(*,*) ’> Programm canceled due to an error!’ goto 999 ! at last we have to close the input file c close the input file 999 close(ionr) 81 82 83 84 85 86 87 88 stop end 9.6.2015 Page 80 Computer Languages for Engineering - SS 15 Within the FORTAN77 code some array functions are called. This functions are coded in FORTRAN90+. You see that using the new GNU FORTRAN compiler, it is possible to mix Fortran77 with FORTRAN90+ without any problems. Listing 3.10: 90 Library to Perform a Matrix Product 1 2 ! Note: every read statements reads exactly one line ! empty lines are NOT ignored 3 4 5 ! read the matrix dimensions integer function ireadmatdim(io,nDim) 6 integer:: io integer, dimension(2):: nDim 7 8 ! io channel number ! dimension array 9 read(io,*,iostat=ioerr) nDim(1),nDim(2) if (ioerr /= 0) then ! handle io errors write(*,*) ’*** Error: reading dimension data!’ ireadmatdim = -1 return ! if not ok, then return endif 10 11 12 13 14 15 16 ! simple check of the dimensions if (nDim(1) < 1 .or. nDim(2) < 1) then write(*,*) ’*** Error: invalid dimension data:’,nDim(1),’,’,nDim(2) ireadmatdim = -2 return endif 17 18 19 20 21 22 23 ireadmatdim = 0 24 ! return code for ok 25 26 end function ireadmatdim 27 28 29 ! function for reading a matrix integer function ireadmat(io,a,nDim) 30 integer::io integer,dimension(2)::nDim real(8),dimension(nDim(1),nDim(2))::a 31 32 33 ! io channel number ! declare the dimensions ! array 34 ! over the rows do i=1,nDim(1) read(io,*,iostat=ioerr) (a(i,j),j=1,nDim(2)) 35 36 37 38 if (ioerr /= 0) then ! if an error occure, return write(*,*)’*** Error: reading matrix data!’ ireadmat = -1 return endif enddo ireadmat = 0 39 40 41 42 43 44 45 46 47 end function ireadmat 48 49 50 ! print array data to the screen with a comment subroutine listmat (comment,a,nDim) E. Baeck 3.6. MATRIX PRODUCT WITH 77-MAIN AND 90-LIBRARY Page 81 51 52 53 54 character*(*) comment integer, dimension(2)::nDim real(8),dimension(nDim(1),nDim(2))::a ! comment to print ! dimension of the matrix ! array data to print write(*,*) comment ! write the comment line 55 56 57 58 59 60 61 ! over the rows do i=1,nDim(1) write(*,*) (a(i,j),j=1,nDim(2)) enddo 62 63 end subroutine listmat 64 65 66 ! product of 2 matrices integer function imatmult(a,b,c,nDimA,nDimB) 67 68 69 70 71 integer, dimension(2)::nDimA,nDimB real(8),dimension(nDimA(1),nDimA(2))::a real(8),dimension(nDimB(1),nDimB(2))::b real(8),dimension(nDimA(1),nDimB(2))::c ! ! ! ! declare declare declare declare the dimension arrays array a array b array c 72 73 74 75 76 77 ! check the dimension of the matrices if (nDimA(2) /= nDimB(1)) then imatmult = -1 return endif 78 79 80 81 82 83 84 85 86 87 88 89 90 ! calculate the product of the matrices ! - over the rows of C do i=1,nDimA(1) ! - over the columns of c do j=1,nDimB(2) c(i,j) = 0. ! initialize it do k=1,nDimA(2) c(i,j) = c(i,j) + a(i,k)*b(k,j) enddo enddo enddo imatmult = 0 91 92 end function imatmult If the main program is coded in FORTRAN90+, the application can be much more flexible as if it would be coded in FORTRAN77. One very useful addon coming from FORTRAN90+ is a buildin access to the command line parameters. With this extension we can pass a for example the name of an input file through the programs interface. The second extension, which is very important here, is an dynamical memory management. In our FORTRAN77 - version, we have to introduce a static memory buffer and have to map all the variables which should be dynamical onto it, i.e. we have to do the memory management, the OS would provide us by free. 9.6.2015 Page 82 Computer Languages for Engineering - SS 15 Listing 3.11: A Dynamical 90+ Version of our first 77 Environment 1 2 3 ! dynamical arrays in FORTRAN 90+ program MatMult90 implicit none 4 5 6 integer::ioin integer::ioout = 10 = 11 ! Input channel ! Output channel 7 8 9 10 ! helper variables integer::ioerr integer::memerr ! error code for io activities ! error code for allocation ! static integer, integer, integer, ! dimension of matrix A ! dimension of matrix B ! dimension of matrix C 11 12 13 14 15 arrays dimension(2)::nDimA dimension(2)::nDimB dimension(2)::nDimC 16 17 18 19 20 ! allocatable arrays real(8), allocatable, dimension(:,:)::A real(8), allocatable, dimension(:,:)::B real(8), allocatable, dimension(:,:)::C 21 22 23 ! function’s return values integer::ireadmatdim,ireadmat,ilistmat,imatmult 24 25 26 27 ! some strings character(256)::infile character(256)::outfile ! input file name ! output file name 28 29 30 31 ! setup standard names infile = "matmult90.in" outfile = "matmult90.out" 32 33 34 ! start comment write(*,*) ’ open input file...’ 35 36 37 38 39 40 41 ! open the input file open(ioin,file=infile,status=’old’,iostat=ioerr) if (ioerr /= 0) then write(*,*) "*** error: file ’",infile(1:len_trim(infile)),"’ not found!" stop endif 42 43 44 45 46 47 48 ! open the output file open(ioout,file=outfile,status=’replace’,iostat=ioerr) if (ioerr /= 0) then write(*,*) "*** error: file",outfile," not found!" stop endif 49 50 !>>>>> get the matrix A <<<<<<< 51 52 53 54 E. Baeck ! read the dimesion of matrix A if (ireadmatdim(ioin,nDimA) < 0) then write(*,*) "*** error: no dimension for matrix A available!" 3.6. MATRIX PRODUCT WITH 77-MAIN AND 90-LIBRARY 55 56 Page 83 stop endif 57 58 write(*,*) "dimension of A: ",nDimA(1),"/",nDimA(2) 59 60 61 62 63 64 65 ! allocate the matrix A allocate(A(nDimA(1),nDimA(2)),stat=memerr) if (memerr /= 0) then write(*,*) ’*** error: allocation of matrix A’ stop endif 66 67 write(*,*) "read matrix data of A..." 68 69 70 71 72 73 ! read the matrix data for A if (ireadmat(ioin,nDimA,A) /= 0) then write(*,*) ’*** error: reading matrix data!’ stop endif 74 75 write(*,*) "print matrix data of A..." 76 77 78 ! print the data of matrix a ioerr = ilistmat(ioout,’Matrix a’,nDimA,a) 79 80 !>>>>> get the matrix B <<<<<<< 81 82 83 84 85 86 ! read the dimesion of matrix B if (ireadmatdim(ioin,nDimB) < 0) then write(*,*) "*** error: no dimension for matrix B available!" stop endif 87 88 write(*,*) "dimension of B: ",nDimB(1),"/",nDimB(2) 89 90 91 92 93 94 95 ! allocate the matrix B allocate(B(nDimB(1),nDimB(2)),stat=memerr) if (memerr /= 0) then write(*,*) ’*** error: allocation of matrix B’ stop endif 96 97 write(*,*) "read matrix data of B..." 98 99 100 101 102 103 ! read the matrix data for B if (ireadmat(ioin,nDimB,B) /= 0) then write(*,*) ’*** error: reading matrix data!’ stop endif 104 105 write(*,*) "print matrix data of B..." 106 107 108 ! print the data of matrix B ioerr = ilistmat(ioout,’Matrix b’,nDimB,b) 109 110 !>>>>> perform the the product A*B -> C <<<<<<< 9.6.2015 Page 84 Computer Languages for Engineering - SS 15 111 ! check the dimension! if (nDimA(2) /= nDimB(1)) then write (*,*) "*** dimension error! A*B not possible!" stop endif 112 113 114 115 116 117 nDimC(1) = nDimA(1) ! rows of A nDimC(2) = nDimB(2) ! columns of B 118 119 120 ! allocate the matrix C allocate(C(nDimC(1),nDimC(2)),stat=memerr) if (memerr /= 0) then write(*,*) ’*** error: allocation of matrix C’ stop endif 121 122 123 124 125 126 127 ! perform the product ioerr = imatmult(a,b,c,nDimA,nDimB) 128 129 130 ! print the data of matrix C = A*B ioerr = ilistmat(ioout,’Matrix c = a*b’,nDimC,c) 131 132 133 ! deallocate the memory deallocate(a,stat=memerr) deallocate(b,stat=memerr) deallocate(c,stat=memerr) 134 135 136 137 138 write(*,*) "close files..." 139 140 ! close the files close(ioin) close(ioout) 141 142 143 144 145 146 147 end program E. Baeck 4 Linear Algebra, Vectors and Matrices This chapter was written as support for the first lectures only dealing with FORTRAN77 development. Later FORTRAN90 and C++ were added to the curiculum, so that this chapter can be considered as obsolete with respect to our current curiculum. 4.1 4.1.1 Helper Functions Outlines Within the following sections we will discuss some helper functions which we will use to implement the gauss decomposition algorithm and its testing environment. 4.1.2 Reset and List a Matrix To check matrices which are decomposed in a lower and upper triangle for example by GaussLU decomposition its helpful to have some helper functions for checking. The helper function ExtractLU extracts the upper and lower triangle matrix of an arbitrary matrix. If we multiply the upper by the lower matrix we should get the original matrix which was decomposed in triangles. The following new statements are used. • include1 , includes a source code file into a main file. • dimension2 , allocates arrays of items of the same data type. • subroutine3 , declares a subroutine which could be seen as function without return value. • call4 , a subroutine can be called by the use of the call statement. • read5 , the read statement is used to read data from keyboard or file. 1 include statement, Page 108, [2] dimension statement, Page 51, [2] 3 subroutine statement, Page 166, [2] 4 call statement, Page 26, [2] 5 read statement, Page 145, [2] 2 85 Page 86 Name Computer Languages for Engineering - SS 15 Comments ResetMat ResetMat resets the content of a given matrix. The content can be reseted optionally to the values of a zero matrix and a unity matrix ListMat ListMat prints the values of given matrix into the screen window. The values of the matrix can be titled with an arbitrary comment string. ExtractLU ExtractLU extracts the values of a LU-decomposed matrix into a normalized lower triangle matrix and a upper triangle matrix. This function is necessary to check the decomposed matrix automatically. MatMult MatMult performs the multiplication of two arbitrary matrices whose dimensions fit to the multiplication algorithm. We use this function to check the decomposed matrix automatically. DiffMat DiffMat calculates the difference matrix of two matrices and returns the norm of the greatest deference item. This function will be used to check the the decomposed matrix. ReadDim ReadDim is used to read the dimension of the matrices from an input file. ReadMat ReadMat reads the matrix values from a text file. This function is used to import test values into the testing environment of the decomposition application. WriteMat WriteMat writes the matrix values to a text file. This function is used to export test values from the testing environment of the decomposition application. Table 4.1: Helper Functions The global trace flag is part of the common block defined in the header file trace.h. Listing 4.1: Global Data 1 2 3 c234567 common /trace/ ntrace integer*4 ntrace E. Baeck ! global flag ! defined in a header file 4.1. HELPER FUNCTIONS Page 87 The subroutine ResetMat sets the data of an array optional to zero or to unit matrix. Listing 4.2: Reset a Matrix’s Data 1 2 3 c subroutine to initialize a matrix (n1xn2) c if mode = 1 a unit-matrix is set c in all other cases a zero-matrix is set 4 subroutine ResetMat(rm,n1,n2,mode) 5 6 7 8 c c mode = 0 >> zero-matrix mode = 1 >> unit-matrix 9 10 11 12 integer*4 n1, n2 integer mode real*8 rm(n1,n2) ! matrix dimensions ! reset flag ! declaring the passed matrix column-index do i=1,n2 ! in fortran it’s faster to run first ! over all rows then over all columns 13 14 c 15 16 17 c 18 row-index do j=1,n1 19 20 c version 1 with nested ifs if (i.eq.j) then if (mode.eq.1) then rm(j,i) = 1. else rm(j,i) = 0. endif else rm(j,i) = 0. endif c c c c c c version 2 with only one if without nesting if (i.eq.j .and. mode.eq.1) then a(j,i) = 1. else a(j,i) = 0. endif 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 enddo 39 40 enddo 41 42 43 return end 9.6.2015 Page 88 Computer Languages for Engineering - SS 15 The subroutine ListMat prints the data of a matrix to the console window. Besides the data of a matrix the subroutine should print a little title. Listing 4.3: List a Matrix’s Data subroutine ListMat(com,rm,n1,n2) 1 2 character *(*) com integer*4 n1, n2 real*8 rm(n1,n2) 3 4 5 ! *(*) means with variable length ! matrix dimensions ! matrix to print 6 7 c loop over all rows write (*,’(a)’) com do i=1,n1 ! print a little title ! loop over all lines ! implicit loop in write statement write(*,’(10f10.2)’) (rm(i,j),j=1,n2) 8 9 10 11 12 enddo 13 14 return end 15 16 To test the helper subroutines ResetMat and ListMat a main program should be developed. Listing 4.4: Check of previous Routines 1 2 3 c234567 c Main program for step 1 program matrices1 4 5 real*8 a(dim,dim) ! declaring a matrix call ResetMat(a,dim,dim,1) call ListMat(’a-matrix’,a,dim,dim) read (*,*) i ! initialize matrix ! list matrix values to the screen ! read a value 6 7 8 9 10 11 E. Baeck end 4.1. HELPER FUNCTIONS 4.1.3 Page 89 LU-Extract, Product and Matrix Compare In this section we will add some further subroutines and functions to the helper functions library of section 4.1.2. If we want to check the LU decomposition A=L·U (4.1) we have to extract the L and U part of a decomposed matrix Ax , which holds the upper triangle and the diagonal of the U in its upper triangle and it’s diagonal values. The values of the L can extracted form the lower triangle of Ax . Because the L has only 1 values on it’s diagonal, we don’t need to store this values. So we need an extractor subroutine, which creates the L and U matrix. Further we need a subroutine for the multiplication of matrices which is called MatMult. At the end we will need a subroutine which searches for the maximum difference of the elements of two matrices which is called DiffMat. The next coding shows the implementation of the extractor of lower and upper triangle. Listing 4.5: Extract the Triangle Data of a Matrix subroutine ExtractLU (a,l,u,n) 1 2 3 4 5 6 c c c c a: l: u: n: result matrix form LU decomposition lower triagle extracted to n x n upper triagle extracted to n x n dimension of a,l,n >> n x n 7 integer real*8 8 9 n a(n,n),l(n,n),u(n,n) 10 11 c 12 rows do i = 1,n 13 14 c 15 columns do j = 1,n 16 17 18 c c if upper triangle, the lower is set to 0, the upper triangle value is taken if (i.lt.j) then l(i,j) = 0. u(i,j) = a(i,j) c c if lower triangle, the upper is set to 0, the lower triangle value is taken else if (i.gt.j) then u(i,j) = 0. l(i,j) = a(i,j) c c if diagonale the lower is set to 1 the diagonale value is assigned to the upper diagonal else l(i,j) = 1. u(i,j) = a(i,j) 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 9.6.2015 Page 90 Computer Languages for Engineering - SS 15 endif 35 36 enddo 37 38 enddo 39 40 return end 41 42 To calculate the original matrix which was decomposed, we have to calculate the product of equation 4.1. In the first version we only use quadratic matrices. C =A·B (4.2) with an matrix element Ci,j Ci,j = n X Ai,k · Bk ,j (4.3) k =1 Listing 4.6: Product of quadratic Matrices subroutine MatMult(a,b,c,n) 1 2 include ’tracegl.h’ 3 ! give access to common block 4 5 6 7 8 c c c c a: b: c: n: input matrix input matrix a x b matrix dimension of n x n n x n n x n quadratic array 9 real*8 a(n,n),b(n,n),c(n,n) 10 11 ! Trace Code if (ntrace.gt.0) write(*,*) ’> MatMult started...’ 12 13 14 15 c 16 row index do i=1,n 17 18 c 19 column index do j=1,n 20 21 22 23 c c c 24 25 26 c 27 28 performing the scalar product of row vector with column vector but at first we have to initialize the matrix element c(i,j) c(i,j) = 0. do k=1,n ! remember: ’*’-Operator has greater priority ! then ’+’-Operator like in mathematics c(i,j) = c(i,j) + a(i,k)*b(k,j) 29 30 31 32 end do enddo enddo 33 34 E. Baeck ! Trace Code ! end of k-loop ! end of j-loop ! end of i-loop 4.1. HELPER FUNCTIONS Page 91 if (ntrace.gt.0) write(*,*) ’> MatMult ended...’ 35 36 return end 37 38 The function DiffMat calculates the norm of the greatest element of a difference matrix. d = max (|Ai,j − Bi,j |) (4.4) Listing 4.7: Difference Matrix of two Matrices real function DiffMat(a,b,n) 1 2 3 4 integer real*8 n a(n,n),b(n,n) real*8 d ! dimension of the quadratic matrices ! matrices to analyse 5 6 7 8 c starting value, any of them we use the element (1,1) d = dabs(a(1,1) -b(1,1)) c rows do i=1,n 9 10 11 12 13 14 c 15 columns do j=1,n 16 17 c 18 19 20 c if the next is greater then actual take the next if (dabs(a(i,j) -b(i,j)) .gt. d) then dabs(x) calculates the norm of x d = dabs(a(i,j) -b(i,j)) 21 22 endif 23 24 25 enddo enddo 26 27 28 29 DiffMat = d return end ! the greatest difference value is passed back ! to the calling program 9.6.2015 Page 92 4.1.4 Computer Languages for Engineering - SS 15 Matrix Import from Input File In this section we want to read data from a input text file and save it to array variables. Therefor we implement a integer function called ReadMat. We pass the io channel number, the reference to the array variable and its rows and column size. We use the read statement6 . If there happen any error, the function returns an integer of 1. If no error occur then then function returns a zero value. This return value can be used in a calling code to handle the error situation. Listing 4.8: Read Matrix Data from a Text File integer function ReadMat(io,rm,n1,n2) 1 2 integer integer integer real*8 3 4 5 6 io n1 n2 rm(n1,n2) ! ! ! ! io-chanal no. dimension of rm (rows) dimension of rm (columns) matrix 7 8 c row loop do i=1,n1 9 10 read (io,*,err=900) (rm(i,j),j=1,n2) 11 12 ! read row values in an implicit ! loop enddo 13 14 ReadMat = 0 return 15 16 ! no error ! return to calling program 17 18 19 20 900 ReadMat = 1 return end ! reading error The testing environment has the following code. We declare the array variable and the used input function ReadMat. Then we open the input file called matrices2.inp with the open statement7 . This file is a text file and contains the matrix element data separated by spaces. The rows of the matrix data are separated by linefeeds. The second file is an output file. The file will only be created. After having read the data, both files will be closed using the close statement8 . If the io-statements are not executed with success the error handler will be activated and will perform a jump to the specified label. 6 read statement, Page 145, [2] open statement, Page 131, [2] 8 close statement, Page 34, [2] 7 E. Baeck 4.1. HELPER FUNCTIONS Page 93 Listing 4.9: Checking the Matrix IO-Functions 1 2 c234567 program matrices2 ! defining program name for linker 3 4 5 6 integer real*8 real*8 ReadMat a(3,3) r(3) io = 5 io = 6 >>> keyboard >>> screen ! declaration of functions ! test matrix ! test vector 7 8 9 c c 10 io1 = 10 io2 = 11 11 12 ! io-number for input ! io-number for output 13 14 c open files open(io1,file=’matrices2.inp’,status=’old’,err=900)! open an existing file open(io2,file=’matrices2.out’,status=’unknown’) ! create a new file c Input section if (ReadMat(io1,a,3,3) .gt. 0) goto 901 call ListMat(’matrix values of a’,a,3,3) 15 16 17 18 19 20 ! read from file ! list read data 21 22 23 24 c close files 800 close(io1,status=’keep’) close(io2,status=’keep’) ! close input file ! close output file 25 26 27 pause ’press return key...’ stop ! wait for a look 28 29 30 31 900 write(*,*) ’ file matrices2.inp not found!’ pause ’press return key...’ ! wait for a look stop 32 33 34 901 write(*,*) ’ format error, reading matrix a.’ goto 800 35 36 end 9.6.2015 Page 94 4.1.5 Computer Languages for Engineering - SS 15 Memory Manager and Pseudo Dynamical Allocation In this section a quasi dynamical memory management approach will be shown. This approach will overcome the statical declaration of arrays in FORTRAN overlaying them with a statical declared memory block. After having opened the input file matrices3.inp (see figure 4.1) we should read the dimension of the first matrix from line 1. This is done by a new helper function which is called ReadDim. With well-known matrix dimensions we can request the used memory form the memory manager. With the calculated memory block index the matrix can be read form the file. This is done with the helper function ReadMat of Section 4.1.4. With ReadMat the matrix values are read line by line. Figure 4.1: Input data for matrices3 After having read the matrix data ReadDim is called which read the dimension of the vector, i.e. the dimension of a matrix with only one column. Then the values of the vector are read line by line with the function ReadMat. Listing 4.10: Read Matrix’s Dimension from File integer function ReadDim(io,rows,cols) 1 2 integer 3 io, rows, cols 4 read (io,*,err=900) rows,cols ReadDim = 0 return 5 6 7 ! read the dimensionvalues rows and cols ! from the input file 8 9 10 11 900 ReadDim = 1 return end ! error branch, error code set if an ! format error is dedected The main program Marices3 shows an approach to solve the problem of dynamical memory allocation in FORTRAN. Because at last FORTRAN offers only the possibility of statical memory allocation we have to develop our own memory manager. E. Baeck 4.1. HELPER FUNCTIONS Page 95 Therefor a memory block is allocated statically and the used memory of the matrices is overlaid on it. If we want to use the memory we have to calculate the index of each overlay. We start at the beginning of the memory block. So the first array will get the index 1. The index of the second array is calculated as the total length of all memory which is already assigned (that means in our example the length of the first matrix) plus 1. This is the first item of the second matrix. Listing 4.11: Checking Matrix Allocation 1 2 c234567 program Matrices3 ! set the applications name 3 integer maxmem parameter (maxmem = 1000) real*8 mem(maxmem) 4 5 6 ! define a parameter to allocate ! the memory block statically ! memory block 7 integer integer integer integer 8 9 10 11 n1,n2,n3,n4 np1 np2 ReadDim, ReadMat ! ! ! ! matrix/vector dimension position of 1st array position of 1st vector declaring the return data type of functions 12 io1 = 10 np1 = 1 13 14 ! input channel for input file ! position of 1st matrix on memoryblock 15 16 c open the input file open (io1,file=’matrices3.inp’,status=’old’,err=900) c reading the dimension of the matrix if (ReadDim(io1,n1,n2) .gt. 0) 17 18 19 20 goto 901 21 22 c reading the dimension of the matrix if (ReadMat(io1,mem(np1),n1,n2) .gt. 0)goto 902 c write matrix content to output window call ListMat(’content of 1st matrix’,mem(np1),n1,n2) c reading the dimension of the vector if (ReadDim(io1,n3,n4) .gt. 0) 23 24 25 26 27 28 29 goto 901 30 31 c reading the dimension of the matrix np2 = n1*n2 +1 if (ReadMat(io1,mem(np2),n3,n4) .gt. 0)goto 902 c write vector content to output window call ListMat(’content of 1st vector’,mem(np2),n3,n4) c halt a little bit pause ’ press return...’ goto 999 32 33 34 35 36 37 38 39 40 ! jump to the end 41 42 c 43 44 45 error branch for missing input file 900 write(*,*) ’ file matrices3.inp not found!’ pause ’ press return...’ stop 46 47 c error branch for incorrect format of dimension line 9.6.2015 Page 96 Computer Languages for Engineering - SS 15 901 write(*,*) ’ format error reading the dimension’ pause ’ press return...’ stop 48 49 50 51 52 c error branch for incorrect format of matrix value line 902 write(*,*) ’ format error reading the matrix’ pause ’ press return...’ stop c "this is the end" 999 stop end 53 54 55 56 57 58 59 The output of the program Matrices3 is shown in figure 4.2. Figure 4.2: Output screen In figure 4.3 the memory overlay for our little example are shown. Left-aligned we see the memory of the 3x3 matrix and at the right of the matrix we see the memory block of the vector. In our case the parameter memmax must be greater equal 12 to have enough memory to allocate the examples data. Figure 4.3: Memory overlays E. Baeck 4.1. HELPER FUNCTIONS 4.1.6 Page 97 Automatic Allocation of a Set of Matrices In this section the main module of the last example is extended. The error handling is changed to a more variable one. The matrices are read within a loop. The memory pointer are calculated step by step based on the read matrix dimensions. We have added two further vectors to our input data file. Figure 4.4: Input data for matrices3, version 2 The code of the main program is given below. Listing 4.12: Checking Automatic Memory Manager 1 2 c234567 program Matrices3 ! set the applications name 3 4 include ’tracegl.h’ ! access to global data, common block integer maxmem parameter (maxmem = 1000) real*8 mem(maxmem) ! define a parameter to allocate ! the memory block statically ! memory block integer maxmat parameter (maxmat = 7) ! number of matrices ! we want to handle 7 matrices 5 6 7 8 9 10 11 12 13 14 15 integer integer integer nrow(maxmat) ncol(maxmat) np(maxmat) ! matrix/vector row dimension ! matrix/vector column dimension ! position of matrix "i" 16 17 18 19 integer nCode character*32 cCode integer ReadDim, ReadMat, ! error-Code of io functions ! error text ! declaring the return data type of functions 9.6.2015 Page 98 Computer Languages for Engineering - SS 15 & 20 21 22 23 WriteMat io1 = 10 io2 = 11 np1 = 1 ! input channel number for input file ! output cannel number for output file ! position of 1st matrix on memoryblock initalization of tracing ntrace = 0 ! tracing disabled 24 25 c 26 27 28 c open the input file cCode = ’*** error: input file not found!’ open (io1,file=’matrices3.inp’,status=’old’,err=900) c c c c read trace information from input file ntrace = 0: tracing disabled ntrace = 1: tracing level one (start/stop of routines) ntrace = 2: tracing level two (al trace data) cCode = ’*** error: traceinfo not found!’ read(io1,*,err=900) ntrace ! c open the output file cCode = ’*** error: output file could not be created!’ open (io2,file=’matrices3.out’,status=’unknown’,err=900) c reading the matrix and vector information of 4 matrices do i =1,4 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 c setup 1st memory pointer if (i.eq.1) then np(i) = 1 c setup momory pointers starting from the 2nd ... else previous + length of previous matrix np(i) = np(i-1) + nrow(i-1)*ncol(i-1) endif 47 48 49 50 51 52 c 53 54 55 56 c 57 58 59 & 60 61 62 reading the dimension of the matrix nCode = ReadDim(io1,nrow(i),ncol(i)) if (nCode .gt. 0) then write (cCode,’(a,i1,a)’) ! setup error information ’*** error: dimension format of ’,i,’. matrix’ goto 900 endif 63 64 c 65 66 67 & 68 69 70 reading the dimension of the matrix nCode = ReadMat(io1,mem(np(i)),nrow(i),ncol(i)) if (nCode .gt. 0) then write (cCode,’(a,i1,a)’) ! setup error information ’*** error: matrix data format of ’,i,’. matrix’ goto 900 endif 71 72 c 73 74 75 E. Baeck write matrix content to output window write(cCode,’(a,i1,a)’) ’content of ’,i,’. matrix’ call ListMat(cCode,mem(np(i)),nrow(i),ncol(i)) 4.1. HELPER FUNCTIONS 76 c 77 Page 99 write matrix content to output file nCode = WriteMat(io2,cCode,mem(np(i)),nrow(i),ncol(i)) 78 end do 79 80 81 c 82 83 84 85 86 1. product np(5) = np(4) +nrow(4)*ncol(4) ! setup memory for product matrix nrow(5)= nrow(1) ! number of rows and columns are given ncol(5)= ncol(2) ! by the numbers of the other matrices call MatMult(mem(np(1)),mem(np(2)),mem(np(5)), & nrow(1),ncol(1),ncol(2)) 87 88 c write matrix content to output window write(cCode,’(a)’) ’product matrix of a * r1’ call ListMat(cCode,mem(np(5)),nrow(5),ncol(5)) nCode = WriteMat(io2,cCode,mem(np(5)),nrow(5),ncol(5)) c halt a little bit pause ’ press return...’ goto 999 89 90 91 92 93 94 95 ! jump to the end 96 97 c only one error branch, because we have an error text 900 write(*,*) cCode pause ’ press return...’ stop c "this is the end" 999 close(io1) close(io2) stop end 98 99 100 101 102 103 104 105 106 9.6.2015 Page 100 Computer Languages for Engineering - SS 15 Figure 4.5 shows the output screen of the discussed program. In the do loop the matrices are read from the input file and after that the matrix values are written to the screen. The last matrix output shows the value of the product of the first and the second matrix. Because the second matrix is the unity vector in the first component direction, the result matrix is equal to the first column vector of the first matrix. Figure 4.5: Output Screen Whats new? • The write statement can also be used to perform formated output to a string. If we want a formated output to a string the first parameter - which is usually the io channel number - is used to pass the destination string. Listing 4.13: Writing into a String 1 write(cCode,’(a,i1,a)’) ’content of ’,i,’. matrix’ • As an index value for the memory pointer we use a indexed array value, mem(np(i)). The pointer value np(i) is stored for every used array and is used as an index value to access the memory block mem. Listing 4.14: Allocating with Memory Pointers 1 E. Baeck nCode = ReadMat(io1,mem(np(i)),nrow(i),ncol(i)) 4.1. HELPER FUNCTIONS 4.1.7 Page 101 Implementing Tracing In this section we will implement tracing in our example of section 4.1.7. Therefore we introduce a new data line in our input file. Its only one value with the following meaning. 0: Tracing is disabled. 1: Tracing is enabled. The start and the end of a subroutine or function call is logged on screen. 2: Tracing is enabled. An extended Tracing is activated. Values of the called subroutines or functions are logged on screen too. We have added one further line at the beginning of our input data file. Figure 4.6: Input data for matrices3, version 3 So we have to introduce a common block as a possibility to access to a global trace flag which is called ntrace. This common block is stored in the file tracegl.h and is used by an include in all subroutines and functions which should use tracing. Listing 4.15: Global Data tracegl.h 1 2 3 c234567 common /trace/ ntrace integer ntrace So we insert the following line of code in all subroutines and functions which should have access to the trace common. Listing 4.16: Include the global Data Header 1 include ’tracegl.h’ 9.6.2015 Page 102 Computer Languages for Engineering - SS 15 To initialize tracing we set a default trace value. We have used disabled tracing, that means a value of 0. Then the first line should be read from the input file, which contents the trace value. The following code is added to our main program. Listing 4.17: Implementing Trace Functionality, 77 like 1 2 3 4 5 6 7 8 9 10 11 12 13 ... c initalization of tracing ntrace = 0 ! tracing disabled c open the input file cCode = ’ *** error: input file not found!’ open (io1,file=’matrices3.inp’,status=’old’,err=900) c read trace information from input file c ntrace = 0: tracing disabled c ntrace = 1: tracing level one (start/stop of routines) c ntrace = 2: tracing level two (all trace data) cCode = ’ *** error: traceinfo not found!’ read(io1, * ,err=900) ntrace ! ... To implement tracing in our function ReadDim we have to give access to the common block including the common code. Then we can optionally write trace information to the screen. If ntrace is greater 0, we log the beginning and the end of the routine. If ntrace is greater 1, we log the read dimension values as well. Listing 4.18: Read Matrix’s Dimension from File integer function ReadDim(io,rows,cols) 1 2 include ’tracegl.h’ 3 ! give access to common block 4 integer 5 io, rows, cols 6 ! Trace Code if (ntrace.gt.0) write(*,*) ’> ReadDim started...’ 7 8 9 read (io,*,err=900) rows,cols 10 ! read the dimensionvalues rows and cols 11 ! Trace Code if (ntrace.gt.1) &write(*,’(a,i3,a,i3)’) ’> 12 13 14 rows: ’,rows,’ columns: ’, cols 15 ReadDim = 0 16 ! from the input file 17 ! Trace Code if (ntrace.gt.0) write(*,*) ’> ReadDim ended, no errors...’ return 18 19 20 21 22 900 ReadDim = 1 ! error branch, error code set if an 23 24 25 26 27 E. Baeck ! Trace Code if (ntrace.gt.0) write(*,*) ’> ReadDim ended with read error...’ return ! format error is dedected end 4.1. HELPER FUNCTIONS Page 103 Like in ReadDim we implement tracing also in the routine ReadMat which is given in the following code. Listing 4.19: Read Matrix’s Data from File integer function ReadMat(io,rm,n1,n2) 1 2 include ’tracegl.h’ 3 ! give access to common block 4 integer integer integer real*8 5 6 7 8 io n1 n2 rm(n1,n2) ! ! ! ! io-chanal no. dimension of rm (rows) dimension of rm (columns) matrix 9 ! Trace Code if (ntrace.gt.0) write(*,*) ’> ReadMat started...’ 10 11 12 13 14 c row loop do i=1,n1 15 16 read (io,*,err=900) (rm(i,j),j=1,n2) 17 18 ! read row values in an implicit ! loop enddo 19 20 ReadMat = 0 ! no error 21 22 23 24 ! Trace Code if (ntrace.gt.0) write(*,*) ’> ReadMat ended without errors...’ return ! return to calling program 25 26 27 28 29 30 900 ReadMat = 1 ! reading error ! Trace Code if (ntrace.gt.0) write(*,*) ’> ReadMat ended with read error...’ return end 9.6.2015 Page 104 Computer Languages for Engineering - SS 15 Figure 4.7 shows the trace output in the output window. Figure 4.7: Screen Output with Tracing Level 2 E. Baeck 4.2. GAUSS-LU-ALGORITHM 4.2 4.2.1 Page 105 Gauss-LU-Algorithm Gauss Decomposition with Pivot Search The following linear equation system is given in matrix notation. A·x =b (4.5) If we write it in components we will get: a1,1 a1,2 . . . a1,n a2,1 a2,2 . . . a2,n .. .. · .. .. . . . . an,1 an,2 . . . an,n x1 b1 x2 b2 = ... ... bn xn (4.6) or n X aik · xk − bi = 0 (i = 1, 2, ..., n) (4.7) k =1 To decompose the matrix A we can use the following theorem.9 For a regular matrix A there is a permutation matrix P so that the product P · A can be decomposed into L a lower and U a upper triangle matrix. P ·A=L·U (4.8) If we insert 4.8 into 4.5 after having multiplied with P from the left hand side. P ·A·x =P ·b (4.9) To calculate the solution vector x we introduce the helper vector c. L·U ·x =P ·b (4.10) L·c =P ·b (4.11) With the forward substitution we get the helper vector c.10 The permutation matrix has to be considered on the right hand side. L·c =P ·b (4.12) With the backward substitution we get the solution vector x .11 R·x =c (4.13) After having decomposed the matrix A, forward- backward substitution can be done for as many as desired right sides. So we can calculate the inverse of A column by column with a respective unity vector as right side. A·X =1 (4.14) 9 The proof is given in H.R. Schwarz [4] Satz 2.4. Forward substitution is possible, because the the vector c is multiplied form the left by a lower triangle matrix. 11 backward substitution is possible, because the the vector x is multiplied form the left by an upper triangle matrix. 10 9.6.2015 Page 106 4.2.2 Computer Languages for Engineering - SS 15 Step 1: Memory Manager and Main Program To handle different data types using only one memory block we apply the equivalence statement12 , which overlays a set of variables. In our example we overlay the real*8 array dmem and the integer*4 array imem. Listing 4.20: Setup Memory Manager for GaussLU, 77 like 1 2 3 4 c234567 program GaussLU c include ’tracegl.h’ 5 6 c 7 memory allocation integer*4 mem8X 8 parameter (mem8X = 1000) 9 10 11 12 13 c 14 15 real*8 dmem(mem8X) integer*4 imem(1) integer*4 imem(mem8X*2) equivalence(dmem(1),imem(1)) c 16 17 18 character*32 integer*4 integer*4 cCode nDim nLC ! dimension of matrix a ! number of load cases (vector b) integer integer integer nPA nPS nPI ! pointer for matrix a ! pointer for scaling vector ! pointer for permutation vector (integer*4) integer ReadMat ! function return values 19 20 21 22 23 24 25 26 c 27 initialization ntrace = 0 ! no tracing io - channels io1 = 10 io2 = 11 ! input ! output 28 29 c 30 31 32 33 c Read input file cCode = ’*** error: input file not found!’ open (io1,file = ’gausslu.inp’,status = ’old’,err=900) c read 1st line cCode = ’*** error: format error in 1st line!’ read (io1,*,err=900) nDim, nLC, ntrace c memory pointer values nPA = 1 nPS = nPA + nDim*nDim nPI = (nDim*nDim + nDim)*2 +1 c read the matrix from input file 34 35 36 37 38 39 40 41 42 43 44 45 46 12 equivalence statement, Page 84, [2] E. Baeck 4.2. GAUSS-LU-ALGORITHM Page 107 nCode = ReadMat(io1,dmem(nPA),nDim,nDim) if (nCode .gt. 0) then cCode = ’*** error: format error reading matrix’ goto 900 endif 47 48 49 50 51 52 53 c close input stream close(io1) c list input values call ListMat(’Matrix values of A:’,dmem(nPA),nDim,nDim) 54 55 56 57 58 write(*,’(a,i5)’)’ Pointer of matrix..............:’, nPA write(*,’(a,i5)’)’ Pointer of scaling vector......:’, nPS write(*,’(a,i5)’)’ Pointer of permutation vector..:’, nPI 59 60 61 62 goto 999 63 64 65 c error branch 900 write(*,*) cCode c end of program 999 continue 66 67 68 69 70 71 72 73 pause ’press enter....’ stop end The first input line sets the dimension of the matrix, the number of the right hand sides and the trace flag. Figure 4.8: Input data 9.6.2015 Page 108 Computer Languages for Engineering - SS 15 Figure 4.9 shows the screen output of the first steps main program. After having read the matrix values the matrix values are listed. Then the real*8 pointer index values of the matrix (1), the scaling vector (3 · 3 + 1 = 10) and the integer*4 pointer index value ((3 · 3 + 3) ∗ 2 + 1 = 25) are listed. Figure 4.9: Screen output of step 1 4.2.3 Step 2: Decomposition with Pivot Search We encapsulate the Gauss LU decomposition in a subroutine of our mathlib module. Listing 4.21: GaussLU Decomposition, 77 like 1 2 3 4 5 6 7 8 c c c c c c c c 9 10 c 11 12 13 Interface description a : matrix to decompose n : dimension of a ip : permutation vector d : scaling vector flag: 0: singular matrix 1: positiv sign -1: negativ sign subroutine gaussLU(a,n,ip,d,flag) c c 14 15 16 17 include ’tracegl.h’ ! give access to common block declaration integer*4 n real*8 a(n,n) integer*4 ip(n) real*8 d(n) ! matrix to decompose ! permutation vector ! scaling vector 18 real*8 real*8 real*8 real*8 19 20 21 22 23 28 ! precision deps = 1.e-7 flag = 1 24 26 ! helper variables c 25 27 s dh dF deps c c c 29 (1) build up the scaling vector for pivot search - initialization for row i do i=1,n 30 31 c 32 33 E. Baeck initialize with original row index ip(i) = i 4.2. GAUSS-LU-ALGORITHM 34 c 35 36 37 38 39 40 c c 41 Page 109 calcation of sum of all values in column j of a s = 0. do j=1,n s = s + dabs(a(i,j)) end do s = 0 ? if yes, then no regular matrix if (s .lt. deps) 42 43 c 44 45 matrix singular return with error flag flag = 0 return 46 47 c else 48 49 50 51 c c 52 scaling with respect to the column sum value of column i d(i) = 1./s 53 endif 54 55 end do 56 57 58 59 c c 60 (2) decomposition - loop over all rows do i =1,n-1 61 62 63 c c 64 65 66 67 68 c c c 69 70 71 72 73 74 75 76 77 c c 78 pivot search - initalization: dpvt = dabs(a(i,i)*d(i)) ipvt = i ! scaled diagonal value ! and its position - no we look at all elements in the matrix below the actual diagonal element do j=i+1,n dh = dabs(a(j,i)*d(j) if (dh .gt. dpvt) then ! scaled element is greater then dpvt = dh ! actual pivot, we take this value ipvt = j ! and its position endif end do check the pivat if (dpvt .lt. deps) then ! if pivot is less then precision ! we have non regular matrix 79 80 c 81 82 83 84 85 86 87 c c c matrix is singular => return with error flag = 0 return end if if neccesary we swap the matrix lines i to pivot-line if (i .ne. ipvt) then 88 89 flag = -flag ! sign information for determinant 9.6.2015 Page 110 90 91 92 c c Computer Languages for Engineering - SS 15 ! calculation to swap data, we need a third variable to avoid over writing! 93 94 c swap permutation values j = ip(i) ip(i) = ip(ipvt) ip(ipvt)= j c swap scaling factors dh = d(i) d(i) = d(ipvt) d(ipvt) = dh c swap matrix do j = 1,n dh a(i,j) a(ipvt,j) end do 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 elements of row i and ipvt = a(i,j) = a(ipvt,j) = dh 110 endif 111 112 113 114 c c 115 elemination step for all rows below actual diagonal element do j = i+1,n 116 117 c 118 119 120 121 c c 122 123 124 elimination factor a(j,i) = a(j,i)/a(i,i) dF = a(j,i) for all values at the right hand side of the actual diagonal elements column do k = i+1,n a(j,k) = a(j,k) -dF*a(i,k) end do ! end of k-loop column-loop of elemeination step 125 126 127 end do end do 128 129 130 E. Baeck return end ! end of j-loop ! end of i-loop row-loop of elemination step global row-loop 4.2. GAUSS-LU-ALGORITHM 4.2.4 Page 111 Step 3: Tracing and Memory Pointers With this step we implement tracing to be able to check the implemented algorithm. 4.2.4.1 Tracing the Decomposer To be able to check the decomposer algorithm we introduce some tracing code into the decomposer routine gaussLU. We implement the following add ons. • We will access to the global variable ntrace of the common-block trace. • We trace the start of gaussLU. • We trace the break of gaussLU, if a singular matrix is found. • We trace the sum of the norm of all items of a matrix row and it’s inverse, the row scaling factor. • Before starting the elimination step we trace the pivot value and the corresponding matrix row. • If row swapping is needed, we trace the original row index and the pivot row index. • We trace the modified matrix elements during elimination. Listing 4.22: GaussLU Decomposition with Tracing, 77 like 1 2 3 4 5 6 7 8 c c c c c c c c 9 10 c 11 12 13 Interface description a : matrix to decompose n : dimension of a ip : permutation vector d : scaling vector flag: 0: singular matrix 1: positiv sign -1: negativ sign subroutine gaussLU(a,n,ip,d,flag) c c 14 15 16 17 18 include ’tracegl.h’ ! give access to common block declaration integer*4 n real*8 a(n,n) integer*4 ip(n) real*8 d(n) integer flag ! ! ! ! matrix to decompose permutation vector scaling vector permutation flag 19 real*8 real*8 real*8 real*8 20 21 22 23 24 25 26 s dh dF deps ! helper variables ! precision c deps = 1.e-7 flag = 1 27 28 ! Trace Code 9.6.2015 Page 112 Computer Languages for Engineering - SS 15 if (ntrace.gt.0) write(*,*) ’> GaussLU started...’ 29 30 31 32 c c c (1) build up the scaling vector for pivot search - initialization for row i do i=1,n 33 34 35 c initialize with original row index ip(i) = i c calcation of sum of all values in column j of a s = 0. do j=1,n s = s + dabs(a(i,j)) end do 36 37 38 39 40 41 42 43 44 c c s = 0 ? if yes, then no regular matrix if (s .lt. deps) then 45 46 47 c matrix singular return with error flag flag = 0 48 49 50 & 51 52 53 if (ntrace.gt.0) write(*,*) ’> GaussLU ended with vanishing column...’ return c else 54 55 56 57 c c scaling with respect to the column sum value of column i d(i) = 1./s 58 59 60 & & 61 62 if (ntrace.gt.1) write (*,’(a,i5,a,e12.5,a,e12.5)’) ’>> row:’,i,’ sum: ’,s,’ scaling value: ’, d(i) 63 endif 64 65 end do 66 67 68 69 c c 70 (2) decomposition - loop over all rows do i =1,n-1 71 72 73 c c 74 75 76 77 78 c c c 79 80 81 82 83 84 E. Baeck pivot search - initalization: dpvt = dabs(a(i,i)*d(i)) ipvt = i ! scaled diagonal value ! and its position - no we look at all elements in the matrix below the actual diagonal element do j=i+1,n dh = dabs(a(j,i)*d(j)) if (dh .gt. dpvt) then ! scaled element is greater then dpvt = dh ! actual pivot, we take this value ipvt = j ! and its position endif 4.2. GAUSS-LU-ALGORITHM end do 85 86 87 Page 113 c c check the pivat if (dpvt .lt. deps) then 88 ! if pivot is less then precision ! we have non regular matrix 89 90 c 91 92 & 93 94 95 96 matrix is singular => return with error flag = 0 if (ntrace.gt.0) write(*,*) ’> GaussLU ended with vanishing pivot...’ return end if c 97 & & 98 99 if (ntrace.gt.1) write(*,’(a,i5,a,i5,a,e12.5)’) ’>> row:’,i,’ pivot row:’,ipvt,’ pivot value:’,dpvt 100 101 102 103 c c c if neccesary we swap the matrix lines i to pivot-line if (i .ne. ipvt) then 104 105 flag = -flag ! sign information for determinant ! calculation to swap data, we need a third variable to avoid over writing! 106 107 108 109 c c 110 111 c swap permutation values j = ip(i) ip(i) = ip(ipvt) ip(ipvt)= j c swap scaling factors dh = d(i) d(i) = d(ipvt) d(ipvt) = dh c swap matrix do j = 1,n dh a(i,j) a(ipvt,j) end do 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 elements of row i and ipvt = a(i,j) = a(ipvt,j) = dh 127 128 & 129 130 if (ntrace.gt.1) write(*,’(a,i5,a,i5)’) ’>> swap of ’,i,’ with ’,ipvt endif 131 132 133 c c 134 elemination step for all rows below actual diagonal element do j = i+1,n 135 136 c elimination factor a(j,i) = a(j,i)/a(i,i) dF = a(j,i) c for all values at the right hand side of the actual 137 138 139 140 9.6.2015 Page 114 141 Computer Languages for Engineering - SS 15 c diagonal elements column do k = i+1,n a(j,k) = a(j,k) -dF*a(i,k) end do ! end of k-loop column-loop of elemeination step 142 143 144 145 146 147 148 & & if (ntrace.gt.1) write(*,’(a,i3,a,i3,2x,10e12.3)’) ’ el.step:’,i,’ row:’,j,(a(j,k),k=j,n) 149 150 151 152 end do end do ! end of j-loop ! end of i-loop row-loop of elemination step global row-loop 153 154 155 if (ntrace.gt.0) & write(*,*) ’> GaussLU ended with decomposition...’ 156 157 158 4.2.4.2 return end Preparing the Test Environment To decompose a matrix using gaussLU and check the result automatically we need memory for the following arrays. • real*8 array(n,n) AC, the copy of the matrix to decompose. • real*8 array(n,n) A, the matrix to decompose. • integer*4 vector(n) Ip, to store the permutation information. • real*8 vector(n) S, to store the scaling information for pivot search. • real*8 array(n,n) L, the lower triangle matrix with zero values in the upper triangle. • real*8 array(n,n) U, the upper triangle matrix with zero values in the lower triangle. • real*8 array(n,n) AR, the product matrix of lower and upper triangle. After having read the input data from the file gausslu.inp we copy the matrix values to a backup matrix using the subroutine MemCpyR8. The helper function MemCpyR8 will be implemented later. After having decomposed the matrix A into a lower and upper triangle, the matrices L and U will be extracted using the helper function extractLU. We now can recalculate the original matrix multiplying L and U and swapping the matrix rows to the original positions. In a final step we calculate the greatest difference value of the last and the first matrix. If this value vanishes with respect to the used arithmetic (precision) the decomposed matrix is checked. E. Baeck 4.2. GAUSS-LU-ALGORITHM Page 115 Listing 4.23: Main Program for GaussLU Decomposition with Tracing 1 2 3 4 c234567 program GaussLU c include ’tracegl.h’ ! access to trace common 5 6 c 7 8 memory allocation integer*4 mem8X parameter (mem8X = 1000) ! length of memory block in *8 units 9 10 11 12 c 13 14 c 15 16 17 real*8 dmem(mem8X) ! memory block integer*4 imem(1) ! integer*4 memory integer*4 imem(mem8X*2) ! the same as a above equivalence(dmem(1),imem(1))! dmem and imem, the same memory ! character*32 cCode ! error code string integer*4 nDim ! dimension of matrix a integer*4 nLC ! number of load cases (vector b) 18 19 20 21 integer integer integer nPA nPS nPI ! pointer for matrix a ! pointer for scaling vector ! pointer for permutation vector (integer*4) integer ReadMat ! function return values 22 23 24 25 c 26 initialization ntrace = 0 ! no tracing io - channels io1 = 10 io2 = 11 ! input ! output 27 28 c 29 30 31 32 c Read input file cCode = ’*** error: input file not found!’ open (io1,file = ’gausslu.inp’,status = ’old’,err=900) c read 1st line cCode = ’*** error: format error in 1st line!’ read (io1,*,err=900) nDim, nLC, ntrace c memory pointer values nPA = 1 nPS = nPA + nDim*nDim nPI = (nDim*nDim + nDim)*2 +1 33 34 35 36 37 38 39 40 41 42 43 ! memory pointer of matrix a ! memory pointer of scaling vector ! memory pointer of permutation vector 44 45 c read the matrix from input file nCode = ReadMat(io1,dmem(nPA),nDim,nDim) if (nCode .gt. 0) then cCode = ’*** error: format error reading matrix’ goto 900 endif c close input stream close(io1) 46 47 48 49 50 51 52 53 54 9.6.2015 Page 116 55 c Computer Languages for Engineering - SS 15 list input values for checking call ListMat(’Matrix values of A:’,dmem(nPA),nDim,nDim) 56 57 write(*,’(a,i5)’)’ Pointer of matrix..............:’, nPA write(*,’(a,i5)’)’ Pointer of scaling vector......:’, nPS write(*,’(a,i5)’)’ Pointer of permutation vector..:’, nPI 58 59 60 61 62 c c save original matrix: AC = A call MemCpyR8(dmem(1),dmem(nPAC),nDim*nDim) 63 64 65 c c LU-decomposition....: A = (L,U) call GaussLU(dmem(1),nDim,imem(nPI),dmem(nPS),flag) 66 67 68 69 c c c Extract.L,U-martices: L = extract(L,U), U= extract(L,U) A L U call ExtractLU(dmem(1),dmem(nPL),dmem(nPU),nDim) 70 71 72 c c Product.............: AR = L * U call MatMult(dmem(nPL),dmem(nPU),dmem(nPAR),nDim,nDim,nDim) 73 74 75 c c difference of AC and AR dif = DiffMat(dmem(nPAC),dmem(nPAR),nDim) 76 77 if (dif .lt. dEpsR) then write (*,*) ’ LU-Decomposition ok!’ else write (*,*) ’ Error in LU-Decomposition!’ endif 78 79 80 81 82 83 goto 999 84 85 86 c error branch 900 write(*,*) cCode c end of program 999 continue 87 88 89 90 91 92 93 94 E. Baeck pause ’press enter....’ stop end 4.2. GAUSS-LU-ALGORITHM Page 117 Figure 4.10: Output of gaussLU Tracing 4.2.5 Forward-Backward Substitution After the decomposition (equation 4.8) the solution of a linear equation system can be calculated using the so called forward-backward substitution (equations 4.10, 4.12 and 4.13). Equation 4.12 describes the forward substitution, which calculates the intermediate vector b. Equation 4.13 describes the backward substitution, which calculates the result vector x . The code of the forward-backward substitution is given below. Listing 4.24: Forward Backward Substitution, 77 like 1 2 3 4 5 6 7 8 c c c c c c c c 9 10 forward / backward substituation a : decomposed n : dimension m : number of load cases ip: permutation vector b : right hand side matrix x : soluton matrix subroutine SubstFB(a,n,m,ip,b,x) c integer n,m real*8 a(n,n),b(n,m),x(n,m) integer*4 ip(n) 11 12 13 14 15 c 16 17 18 19 20 c c c 21 helpers real*8 s integer ipvt ! original row of load item over all right hand sides do k=1,m 22 23 24 25 26 c spectial case dim = 1 if (n.eq.1) then x(1,k) = b(1,k)/a(1,1) goto 100 ! next load case 9.6.2015 Page 118 endif 27 28 29 30 Computer Languages for Engineering - SS 15 c c c forward substituation: L*c = P*b ipvt = ip(1) x(1,k) = b(ipvt,k) 31 32 33 do i=2,n s = 0. do j=1,i-1 s = s + a(i,j)*x(j,k) enddo 34 35 36 37 38 39 ipvt = ip(i) x(i,k) = b(ipvt,k) -s enddo 40 41 42 43 44 45 c c c backward substituation: U*x = c x(n,k) = x(n,k)/a(n,n) 46 47 do i=n-1,1,-1 s = 0. do j=i+1,n s = s + a(i,j)*x(j,k) enddo 48 49 50 51 52 53 x(i,k) = (x(i,k) -s)/a(i,i) enddo 54 55 56 57 58 100 continue enddo ! end of load case loop 59 return end 60 61 4.2.6 Linear Solver The solution of a linear equation system is performed in two steps. In a first step the equation matrix is decomposed by a Gauss-LU decomposer. The second step will be performed for every right hand side calculating the result vector for the selected right hand side. The input data consist of the system dimension, the matrix data and the data of the right hand sides. An example is given below. Listing 4.25: Inputdata to Read from a Text File 1 2 3 4 5 6 7 3 4 2.0 0.7 1.3 1.0 0.0 0.0 E. Baeck 2 1.1 1.3 3.0 -0.7 -3.0 1.5 0.0 1.0 3.0 1.0 2.0 2.0 0.0 3.0 1.0 >> dimension , number of "load cases" >> start of matrix data >> 4 vectors of right hand side 4.2. GAUSS-LU-ALGORITHM Page 119 The main program implementing the solution of a linear equation system is given below. Listing 4.26: Implementation of a Linear Solver Environment 1 2 C1234567 program LinSolve 3 include ’tracegl.h’ 4 5 6 c memory block integer*4 mem8X parameter (mem8X = 1000) c Memory real*8 integer*4 7 8 9 10 11 12 13 c equivalence(dmem,imem) 14 15 16 c c 17 18 19 20 21 c c 22 23 24 25 26 27 28 35 45 51 52 io1 = 10 io2 = 11 ntrace = 0 ! input channel ! output channel nLA nLB nLS 44 50 ! status code c 43 49 character*32 cCode cCode = ’*** error: format error!’ read(io1,*,err=900) nDim, nLC,ntrace1 41 48 ! declare the functions c 40 47 ReadMat cCode = ’*** error: input file not found!’ open (io1,file = ’linsolve.inp’,status = ’old’,err=900) 38 46 Pointer of A right hand side matrix pointer solution matrix pointer scaling vector pointer pointer for the product A*X pointer for permutation information c 37 42 ! ! ! ! ! ! c 34 39 adress pointer integer nPA integer nPB integer nPX integer nPS integer nPAX integer nPI integer 33 36 ! Systemdimension ! Number of load cases ! Permutation sign c 31 32 system parameters integer*4 nDim integer*4 nLC integer*4 nFlag c 29 30 dmem(mem8X) imem(1) c c = nDim*nDim = nDim*nLC = nDim ! number of items in A ! number of items in B,X and A*X ! number of items in scaling vector calculation of memory pointers nPA = 1 nPB = nLA +1 nPX = nLA +nLB +1 nPS = nLA +2*nLB +1 nPAX = nLA +2*nLB +nLS +1 ! ! ! ! ! pointer pointer pointer pointer pointer of of of of of A B X scaling vector the testing matrix A*X 9.6.2015 Page 120 nPI 53 54 55 Computer Languages for Engineering - SS 15 c c = (nLA +3*nLB +nLS)*2 +1 ! pointer of permutation vector reading input data of A nCode = ReadMat(io1,dmem(nPA),nDim,nDim) if (nCode .gt. 0) then cCode = ’*** error: format of matrix A’ goto 900 endif 56 57 58 59 60 61 62 c reading input data of B nCode = ReadMat(io1,dmem(nPB),nDim,nLC) if (nCode .gt. 0) then cCode = ’*** error: format of matrix B’ goto 900 endif 63 64 65 66 67 68 69 c c close input file close(io1) 70 71 72 c list of decomposed call ListMat(’Input matrix:’,dmem(nPA),nDim,nDim) c list of right hand side matrix call ListMat(’Right hand sides:’,dmem(nPB),nDim,nLC) 73 74 75 76 77 78 c c decompose A using Gauss-LU-Decomposition call GaussLU(dmem(nPA),nDim,imem(nPI),dmem(nPS),nflag) if (nflag.eq.0) then write(*,*) ’*** error: matrix A is singular!’ goto 999 endif 79 80 81 82 83 84 85 c c list of decomposed call ListMat(’Decomposed matrix:’,dmem(nPA),nDim,nDim) 86 87 88 c c forward / backward substituation call SubstFB (dmem(nPA),nDim,nLC,imem(nPI),dmem(nPB),dmem(nPX)) 89 90 91 c c list of solution vectors call ListMat(’Solution vectors:’,dmem(nPX),nDim,nLC) 92 93 94 c c jump to the end goto 999 95 96 97 c 900 98 error output write(*,*) cCode 99 100 101 c 999 102 103 104 E. Baeck the end continue pause ’press enter...’ stop end 4.2. GAUSS-LU-ALGORITHM Page 121 Figure 4.11 shows the output window of the LinSolve program. To check the execution of the program some test strings were written to the output screen. Figure 4.11: Outputwindow of LinSolve 9.6.2015 Page 122 E. Baeck Computer Languages for Engineering - SS 15 Part II C/C++ 123 Page 125 C++ is a statically typed, free-form, multi-paradigm, compiled, general-purpose programming language. It is regarded as a ”middle-level” language, as it comprises a combination of both high-level and low-level language features. It was developed by Bjarne Stroustrup starting in 1979 at Bell Labs as an enhancement to the C language and originally named C with Classes. It was renamed C++ in 1983. As one of the most popular programming languages ever created, C++ is widely used in the software industry. Some of its application domains include systems software, application software, device drivers, embedded software, high-performance server and client applications, and entertainment software such as video games. Several groups provide both free and proprietary C++ compiler software, including the GNU Project, Microsoft, Intel and Embarcadero Technologies. C++ has greatly influenced many other popular programming languages, most notably C# and Java. A nice C++ documentation in the web you get with the link [6]. 9.6.2015 Page 126 E. Baeck Computer Languages for Engineering - SS 15 5 Development Tools 5.1 The Code::Blocks IDE We have talked a lot about the Code::Blocks IDE in section 5.1. For our C/C++ codes we use the Code::Blocks IDE as well. We only need to install the Fortran version of Code::Blocks, because it’s a simple add on to the original one, which is written for C/C++ developments. The only thing we should change developing C/C++ projects is the selection of the project type and the selection of the compiler. This is discussed below. If we create a new project, to implement C/C++ coding, we have to select the Consol Application from the project creation wizzards template list (see figure 5.1). Figure 5.1: Selecting a Console Application, a C/C++-Project If the template is selected, we have to select the standard type of Consol Application, that means we have to select a C project or a C++ project. Because the C features are also available in a C++ development, we select the C++ console application type (see figure 5.2). If we close the console application type, we get the well known start up form for our project (see figure 5.3). We have to specify the project’s root folder and the project’s name. Within the next step we have to setup the targets. It’s very important to select now the proper compiler. 127 Page 128 Computer Languages for Engineering - SS 15 Figure 5.2: Selecting of a C++ Console Application Figure 5.3: Selecting of a C++ Console Application In this case it is the GNU C++ compiler, which is called GCC (see figure 5.4). If we click on next, we have completed the creation procedure of a console application. If we open the project browser an load the code file main.cpp, we can study the startup code (see figure 5.5). The HelloCpp node is inserted into the project tree, like we have seen in the case of Fortran projects. In the case of C++ project, the code generator creates a file with the name main.cpp. cpp is used as extension of C++ files. This extension is important to select the GCC compiler building the project’s executable. main.cpp is used, because a C/C++ application needs a main function. This main function will be called, if the executable is started. E. Baeck 5.1. THE CODE::BLOCKS IDE Page 129 Figure 5.4: Specifying the project targets Figure 5.5: Startup Code of a C++ Project 9.6.2015 Page 130 E. Baeck Computer Languages for Engineering - SS 15 6 Basics of C/C++ A detailed description of C and C++ is available in literature or can be downloaded from several web pages. 6.1 The Preprocessor C and C++ are native languages which are compiled into the native code. Before a code file (cpp-file) is compiled by the compiler a preprocessing step is performed. The preprocessor includes further files into a file or substitutes textually so-called macros, which are defined in terms of preprocessor commands. A further feature of the preprocessor is, the possibility of macro driven optional compiling. Preprocessor commands in general starts with a # character. In the following box in the first line a file of the development package is included. If we use the parenthesis < and > the compiler searches for the files first in the folders of the compiler packages and then in the actual folder. In the second line the file name is bracketed with the quotes ". Therefore the compiler searches first in the current folder and then in the folders of the development package. With #define a macro is introduced. With #undefine the macro can be removed from the macro list. With #ifdef, #elif, #else and #endif optional compiling can be performed. This preprocessor commands are working like the if, else, endif of the FORTRAN compiler (section 2.9). Listing 6.1: Working with Macros 1 2 #include <a_library_file.h> #include "a_user_file.h" 3 4 5 6 7 8 9 #define MYMACRO #ifdef MYMACRO ... code, which will be compiled if MYMCRO was defined #else ... code, which will be compiled if MYMCRO was not defined #endif 131 Page 132 6.2 Computer Languages for Engineering - SS 15 Comments The C++ language offers to different comment types. With the classical block comment, which comes form the C language, we can bracket a part of a line or some lines and make them comments. We start the comment with /* and close the comment with */. We can put an arbitrary number of characters and lines in between this brackets. The second type of comment is a line end comment. This line end comment is set with the character //. In the following box we see at the beginning a block comment followed by several line end comments. The one and only statement of this code includes the <stdio.h>, which is the declaration header of the printf function, which comes from C and is used by standard to print to the screen. Listing 6.2: Block- and Line End Comment 1 2 3 4 /* A little comment example (this is a block comment) */ 5 6 // ... and this is a line end comment 7 8 9 10 11 // declare the printf function #include <stdio.h> // <..> compiler searches first in // the lib folder // ".." compiler searches in actual folder 6.3 Data Types C/C++ offers the following data types, which are also related the the platform. This data types can be grouped into character, integral, floating point and void types. type char bytes comment 2 character or small integer range signed: -128 to 127 unsigned: 0 to 255 short 2 short integer signed: -32768 to 32767 unsigned: 0 to 65535 int 4 standard integer signed: -2147483648 to 2147483647 unsigned: 0 to 4294967295 long 4 long integer signed: -2147483648 to 2147483647 unsigned: 0 to 4294967295 long long 8 64 bit integer signed: −263 to +263 − 1 ≈ 9, 2218 unsigned: 0 to 264 − 1 float 4 floating point number. +/- 3.4e +/- 38 ( 7 digits) double 8 double precision floating point number. +/- 1.7e +/- 308 ( 15 digits) E. Baeck 6.4. OPERATORS Page 133 The integer data type can be modified with the unsigned key. 6.4 Operators C/C++ offers the following operators. 6.4.1 Assignment Operator operator comment Assignment = 6.4.2 example i = 5; Arthmetic Operators operator comment example + addition i = 9 +5; ⇒ 14 - substraction i = 9 -5; ⇒ 4 * multiplication i = 9 *5; ⇒ 45 / division i = 9 /5; ⇒ 1 % modulo i = 9 %5; ⇒ 4 6.4.3 Compound Arithmetic Assignment operator comment example += addition i = 2; i += 9; ⇒ 11 -= substraction i = 3; i -= 1; ⇒ 2 *= multiplication i = 2; i *= 4; ⇒ 8 /= division i = 6; i /= 2; ⇒ 3 %= modulo i = 9; i %= 2; ⇒ 1 6.4.4 Increment - Decrement Operators operator comment example ++ incrementation i = 2; i++; ⇒ 3 -- decrementation i = 3; i--; ⇒ 2 9.6.2015 Page 134 6.4.5 Computer Languages for Engineering - SS 15 Relational and Equality Operators operator comment example == equal to 3 == 2; ⇒ 0 != not equal to 3 != 2; ⇒ 1 > greater 3 > 2; ⇒ 1 < less than 3 < 2; ⇒ 0 >= greater equal 3 >= 3; ⇒ 1 <= less equal 3 <= 4; ⇒ 0 6.4.6 Logical Operators operator comment example ! not !(3 == 2); ⇒ 1 && and (3 == 2) && (1==1); ⇒ 0 || or (3 == 2) || (1==1); ⇒ 1 The following table shows the truth values of the && and the || operator. Truth tabel of the && operator 6.4.7 Truth tabel of the || operator a b a && b a b a || b true true true true true true true false false true false true false true false false true true false false false false false false Bitwise Operators operator comment example & and 0x20 & 0xff; ⇒ 0x20 | or 0x20 | 0xff; ⇒ 0xff ˆ exclusive or 0x20 ˆ 0xff; ⇒ 0xdf ˜ complement ˜0x20; ⇒ 0xffffffdf << shift left 0x20<<1; ⇒ 0x40 >> shift right 0x20>>1; ⇒ 0x10 E. Baeck 6.4. OPERATORS 6.4.8 Page 135 Compound Bitwise Assignment operator comment example &= and i = 0x20; i &= 0xff; ⇒ 0x20 |= or i = 0x20; i |= 0xff; ⇒ 0xff ˆ= exclusive or i = 0x20; i ˆ= 0xff; ⇒ 0xdf <<= shift left i = 0x20; i <<=1; ⇒ 0x40 >>= shift right i = 0x20; i >>=1; ⇒ 0x10 6.4.9 Explicit Type Casting Operator The type casting operator converts one type into another type. The destination type is set in between two round brackets. In the following example an integer -1 is casted into an unsigned character, so the sign bit is no more interpreted as a sign and becomes part of the one byte number information. Every bit now is set in the one byte variable and therefore we get 255 as value. Listing 6.3: Casting an integer to a unsigned char 1 2 3 4 unsigned char i; int j = -1; k = (unsigned char)j; printf("k= %d\n",k); 5 6 7 ... output: k= 255 6.4.10 sizeof Operator The sizeof operator determines the length of a variable in bytes. In the following example the length is determined of an unsigned char, an int, a float and a double. So we get the numbers 1, 4, 4 and 8. Listing 6.4: Evaluating the Size with the sizeof Operator 1 2 3 4 5 unsigned char i int j float f double d printf("%d, %d, = 1; = -1; = 1.2; = 3.4; %d, %d\n",sizeof(i),sizeof(j),sizeof(f),sizeof(d)); 6 7 8 ... output: 1, 4, 4, 8 9.6.2015 Page 136 6.4.11 Computer Languages for Engineering - SS 15 Address and Value Operator The address operator & determines the address of a variable in memory. The value operator * determines the value of data with a given address. In our example we determine the address of a variable i and then the value of the second item of an integer array is determined. The address of an array is given by the arrays name. The address of the second item of an array therefore is given by the name of the array plus one. operator comment example & address printf("%X\n",&i); ⇒ 22FF40 * value int s[2] = {1,2}; *(s+1); ⇒ 2 6.4.12 C++ operator synonyms For some operators there are macros which provides some synonyms for C/C++ operators. Some of them are given within the following table 6.5 operator synonym example && and (3 == 2) and (1==1); ⇒ 0 || or (3 == 2) or (1==1); ⇒ 1 ! not not (3 == 2); ⇒ 1 != not_equ 3 not_equ 2; ⇒ 1 & bitand 0x20 bitand 0xff; ⇒ 0x20 | bitor 0x20 bitor 0xff; ⇒ 0xff ˆ xor 0x20 xor 0xff; ⇒ 0xdf ˜ compl compl 0x20; ⇒ 0xffffffdf &= and_equ i = 0x20; i and_eq 0xff; ⇒ 0x20 |= or_equ i = 0x20; i or_eq 0xff; ⇒ 0xff ˆ= xor_equ i = 0x20; i xor_eq 0xff; ⇒ 0xdf Taking about the Hello We have seen, if we create a new project within Code::Blocks IDE a Hello code is created automatically. From this code we can see how a C/C++ code is buildup. We see from the following example code, that a C/C++ application has to have a function, which is called main. A function in C/C++ starts the the type of return, in this case an int, which is an integer. The type of return is followed by the functions name. In this case the name should be main, because it’s the main routine. A function, like in FORTRAN, will have some formal parameters, which are listed in between a pair of brackets. In this case we don’t use parameters, therefore the brackets are empty. The definition of the function’s interface is followed by a code block. In C/C++ a code block is enclosed within curled brackets {...code block...}. Every statement, which is an executable or declaring code, has to be closed by a semicolon. Therefore in general more than one statement can be written into E. Baeck 6.6. LINE OUTPUT WITH PRINTF Page 137 one line. Note the line end comments too, which are listed with a green highlighting. We also can see, that a function can be exited with a return statement. If a function should return a value, then the return statement is followed by an expression. In this case a value of zero is returned to the operating system, which is calling the main function. Listing 6.5: Our First Hello 1 2 // declare the iostream #include <iostream> 3 4 5 // setup the namespace using namespace std; 6 7 8 9 10 11 12 6.6 int main() { // stream the string to the screen cout << "Hello world of C++!" << endl; return 0; } Line Output with printf The C function printf to print something into the console window will get at least one parameter. This parameter is the text to print. In general there are a lot of data to be filled into this text with so called escapes. For every escape there is one data item needed. The data items to be filled into the text have to be passed in the order of their escapes starting from the second argument of the printf function. So the following printf will print four values into the console window, therefore we have to pass five parameters, the text with the escapes and the four values which have to be filled into the escapes. Listing 6.6: printf Example 1 2 3 4 5 int i = 1; char c = ’A’; float f = 1.2; double e = 2.4e-4; printf("Let’s print 4 values: i=%d, c=%c, f=%10.3f and e=%12.3e\n",i,c,f,e); The consol output is the following. Listing 6.7: Consol Output of above’s Example 1 2 3 4 5 int i = 1; char c = ’A’; float f = 1.2; double e = 2.4e-4; Let’s print 4 values: i=1, c=A, f= 1.2 and e= 2.400e-004 9.6.2015 Page 138 Computer Languages for Engineering - SS 15 There is a wide set of escapes. The most used escapes are described below. An escape starts with the % character and ends with the type specifier. The type specifier selects the data type to print. 1 %[flags][width][.precision][length]specifier • flags specify the kind of output, for example the justification • width specify field width, which should be used for the output • precision specify the number of digits in the case of a floating point item • length specify data type length, for example a long float The following table shows the most important specifiers we use in an escape.1 1 specifier output example d or i signed decimal integer -123 u unsigned decimal integer 123 o unsigned octal 771 x unsigned hexadecimal integer ff X unsigned hexadecimal integer (uppercase) FF f decimal floating point 392.65 e scientific notation (mantissa/exponent), lowercase 3.9265e+2 E scientific notation (mantissa/exponent), uppercase 3.9265E+2 g use the shortest representation: %e or %f 392.65 G use the shortest representation: %E or %F 392.65 c character a s string of characters sample Hello % % followed by another % will write a single % % A complete description of all escape specifiers si given in http://www.cplusplus.com/reference/cstdio/printf/. E. Baeck 6.7. A FOR LOOP 6.7 Page 139 A For Loop The most used loop statement in C/C++ is the for statement. The for statement executes a statement or a code block. The execution of the loop is controlled by tree expressions. Listing 6.8: Syntax of a for Loop 1 2 3 4 for ([<initialization>];[<break condition>];[<iterator condition>]) { [<code block>] } • The initialization is executed before the loop is starting. • The break condition is evaluated after each loop cycle. If the value of the expression is vanishing the loop will be broken, if not a next cycle will be started. • The iteration condition is a statement, which is executed after the execution of a cycle. The following example shows how to iterate an integer starting from 5 up to 50 with a step width of 5. The first step is to declare and initialize the used variables. We declare them all as int to get 4 byte integers. After the data type - in this case int - the variable optionally can be initialized using an = operator for the assignment. After we have declared and initialized the variables, a little header line is printed onto the screen with the printf function. Here we see, that a constant string in C/C++ is bracketed with double quotes ". A line break character \n uses the so-called escape character \. Then the for loop is started with the for key word and it’s control fields. Within the first field the i loop counter variable is set to the starting value, which we have stored in the variable iFrom. The second field contents the boolean expression, which should control the cycles. If the loop counter i is less equal the end value, which is stored in the variable iTo a next cycle is started. Here we use the operator <=, which is also introduced in FORTRAN90+. The third field adds the value of iStep to the loop counter i. This is done with the += operator. We also can write instead of the third expression i=i+iStep. Within the cycle the value of i is printed to the screen using the printf function. The string, which should be printed to the screen contents the escape %d this is an escape to fill in an integer value. In this case it’s the value of i, which follows the string as first data parameter. The loop is completed by the closing parenthesis of the code block. Listing 6.9: Simple loop Example 1 2 3 4 /* A little loop example (this is a block comment) */ 5 6 // ... and this is a line end comment 7 8 9 10 11 // declare the printf function #include <stdio.h> // <..> compiler searches first in // the lib folder // ".." compiler searches in actual folder 9.6.2015 Page 140 Computer Languages for Engineering - SS 15 12 13 // note: every statement has to be closed by an ’;’ character 14 15 16 17 18 19 20 int main() { int i; int iFrom = 5; int iTo = 50; int iStep = 5; // note: there is no implicit declaration // here we use an integer datatype 21 // "\n" : line break // printf is the standard C print routine, which is declared // within the stdio.h file of the standard C library printf("This is a Loop Application\n"); 22 23 24 25 26 // run a loop for (i=iFrom;i<=iTo;i+=iStep) { // %d: integer escape printf("i = %3d\n",i); } 27 28 29 30 31 32 33 } Figure 6.1 shows the screen output of the example ALittleLoop. Figure 6.1: Screen Output of ALittleLoop E. Baeck 6.8. STATIC ARRAYS 6.8 Page 141 Static Arrays In this section we will see how to allocate static arrays in C/C++. A static array is like in FORTRAN a sequence of data of the same data type. The data within the array are accessed by an index. A very important difference between C/C++ and FORTRAN is, that the first index of an array in C/C++ is zero and a standard array in FORTRAN starts with the index one. So if we mix FORTRAN and C/C++ code we should be very careful with the array indexing. To show the usage of static arrays we will discuss a little example, which calculates the scalar product of two vectors, which are represented in the code by two arrays. An array in C/C++ is declared by the data type followed by the arrays name and the arrays dimension. Listing 6.10: Declaration of a Static Array 1 <data type> <name> [ <index 1>[,<index 2>]...[,<index n>] ]; In our example we use two array with one index and it’s dimension of 3. The data type should be a large float which in C/C++ is called double. Within the first part our example declares and initializes the vectors and the result variable sp. If an array should be initialized, a list of values separated by commas and bracketed into curled brackets is assigned to the arrays name. Then the result is calculated within a for loop, which runs over all vector items. See also equation 6.1, which describes how to calculate the scalar product of two vectors ~a and ~b. The loop counter is initialized with zero, the first array index. The loop runs up to the index 2, which is the last valid array index. After each cycle the counter i is incremented with the incrementation operator ++. s = ~a · ~b = n X ai · bi (6.1) i=1 After having calculated the scalar product the result is printed to the screen with the printf function. In this case we use the escape %10.3f. The first number sets the width of the output field, in this case 10 character. The second number sets the decimal places. Listing 6.11: Multiplying Vectors by a Dot Product 1 2 // declare the printf function #include <stdio.h> 3 4 5 6 7 8 9 int main() { int double double double i; s1[3] = {1.1,1.2,1.3}; s2[3] = {2.1,2.2,2.3}; sp = 0.; // // // // loop counter vector 1 vector 2 scalar product 10 11 12 13 14 15 16 // calculate the scalar product sp = 0.; // you should initialize everything for (i=0;i<3;i++) // over all components of the vector { sp += s1[i]*s2[i]; // array access by index } 17 9.6.2015 Page 142 Computer Languages for Engineering - SS 15 // print the result printf("s1*s2 = %10.3f\n",sp); 18 19 // print the result 20 // print the addresses of the arrays. If we do this, the address // should be casted to an unsigned integer value // - this is an cast operator: new = (datatype)old // it converts the old into the new printf("address of s1: %X\n",(unsigned int)s1); printf("address of s2: %X\n",(unsigned int)s2); 21 22 23 24 25 26 27 // calculate the scalar product // using item pointers, i.e. the address of the data in memory // // s1 is the address of the first item in the array s1 // so we add up the index value [0..2] and get all the values // of the array. The * - operator gets the value which is stored // at the given address sp = 0.; for (i=0;i<3;i++) { // s1[i] s2[i] // ----------------sp += (*(s1+i)) * (*(s2+i)); } 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 // print the result (we hope it’s the same as above ;) printf("s1*s2 = %10.3f\n",sp); 43 44 45 } Figure 6.2 shows the screen output of the example ScalarProduct. Figure 6.2: Screen Output of ScalarProduct E. Baeck 6.9. BRANCHING 6.9 Page 143 Branching In C/C++ we have the if statement, which is very close to FORTRAN’s if. On the other hand we have the switch statement, which is not so general because it is only working on comparing variable values of the type int and char to constant values. 6.9.1 if Branching Like in FORTRAN too C/C++ provides some statements for branching. The most used and most common branching statement is the if..else if..else statement. The if is followed by an expression. If the expression’s value is zero the assigned code block is not executed. If the expression’s value is not vanishing the code block is executed. Optionally we can extend the first if by further else if conditions. At the end an optional else can be introduced with it’s own code block. Listing 6.12: Syntax of the if Statement 1 2 3 4 5 if ([expression 1]) code block 1 else if ([expression 2]) code block 2 else if ([expression 3]) code block 3 ... else code block n The following example shows the usage of the if..else if..else statement and the usage of some operators, which are discussed in section 6.4. Listing 6.13: Simple Operator Example Using if 1 2 3 /* Example to discuss some C/C++ operators */ 4 5 6 // this we need for printf #include <stdio.h> 7 8 9 10 11 int main() { int int i=2, j=-2; flag; // two number // and a flag 12 13 14 15 16 17 18 19 20 21 22 23 // compare two numbers // 1st if if (i == j) printf("%d is equal %d\n",i,j); // 2nd if else if (i == -j) printf("%d is equal minus %d\n",i,j); // ... and the else branch else { printf("no if executed!\n"); printf("i = %d, j = %d\n",i,j); } 24 25 26 // here we use the not equal operator if (i != j) printf("%d is not equal %d\n",i,j); 9.6.2015 Page 144 Computer Languages for Engineering - SS 15 27 // now there are some very useful bit operators flag = 0; 28 29 30 // set the 20th bit, coded in a hexdecimal number flag |= 0x00080000; // 4 byte number 31 32 33 // we print the number first in hex and then in decimal printf("flag = 0x%8.8X - %d\n",flag,flag); 34 35 36 // the 20th bit is set and all bits of the lower 2 bytes. // (if we set all bits we can simple use the F digit) flag |= 0x0008FFFF; // 4 byte number 37 38 39 40 // to clear the 2nd bit we first code it in hex 0x02 // and then invert it with the ˜ operator. Then the inverse // of the second is overlayed bit by bit with the and operator & flag &= ˜0x2; 41 42 43 44 45 // the result is printed to the screen printf("flag = 0x%8.8X - %d , inverted 2nd bit: 0x%8.8X\n",flag,flag,˜0x2); 46 47 48 6.9.2 } switch Branching An alternative branching can be implemented using the switch statement. This statement preferable is used to check a variables value’s. We will find this statement for example in message filters of windowed environments. The statement starts with the switch key followed by a series of case blocks. A case block is executed either if the case value is found in the switch variable or if the previous case block is not closed by a break statement. The break statement at the end of a case block is optional and exits the switch statement. Listing 6.14: Syntax of the switch Statement 1 2 3 4 5 6 7 8 9 10 11 12 switch (variable) { case value1: code block 1 [break] case value2: code block 2 [break] ... default: dfault code block } E. Baeck 6.9. BRANCHING Page 145 In example 6.15 a variable i is declared. The value 2 is checked inside the switch statement. Because there is a unclosed 2 block, the code of the 2 and the following 3 block will be executed. The result is shown below. Listing 6.15: Simple switch example 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 // usage of a SWITCH-CASE statement int var = 2; .... switch (var) { case 1: printf("this is case 1: var = 1\n"); break; case 2: printf("this is case 2 and 3: we will enter next case block\n"); case 3: printf("this is case 3: var = %d\n",var); break; default: printf("sorry no case available for var = %d\n",var); } Listing 6.16: Console Output of above’s Example 1 2 this is case 2 and 3: we will enter next case block this is case 3: var = 2 9.6.2015 Page 146 6.10 Computer Languages for Engineering - SS 15 Exceptions The classical error handling is working with return codes, i.e. a function is called and will return a code (mostly a number), to show whether the the function had success or not. If we only have a very flat calling hierarchy, this is not a big work to implement. But if we have a lot of functions to call, the way back to the master call can be very difficult. So in modern languages so called error handlers have been implemented, which are working similar to event handlers, i.e. an error event is created, the function is canceled and the error handler is searching for a goal to continue with the work. An error event like this is called an exception. The exception will be thrown in the case on an error. The error event handler is enabled, if we put the code, which can produce exceptions into an try block. So, if an error occur in this try block, the error event handler is searching for a snipped of code, which should treat this event. This codes are put into so called catch blocks, which should catch a thrown exception. So the syntax of this error handling is as follows. Listing 6.17: Syntax of the if Statement 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 // enabling error handling try { ... // this code should be executed with an error handler } // catch block with exception data catch(<type> <variable>) { ... // error code for a specified exception } // at the end all unspecified exceptions are handled catch(...) { ... // error code for all unspecified exceptions } The following example shows how to throw an exception. In line 8 an integer exception is thrown with the code 123. So the try code is aborted in this line and the catch block in line 23 is executed. If we would set the first throw an comment, the second would throw an const char* exception, which passes a simple C string. In this case the catch block starting at line 17 would be executed. Listing 6.18: Simple Exception Example 1 2 3 4 5 6 7 8 #include <stdio.h> int main() { // try to do something try { // create an integer exception throw 123; 9 10 11 E. Baeck // create a string exception throw "*** kill the startup!"; 6.10. EXCEPTIONS Page 147 12 printf("some unreachable statements....\n"); 13 } 14 15 // catch the exception, if something is wrong#include "node.h" catch(const char* str) { printf("*** Exception: ’%s’\n",str); } 16 17 18 19 20 21 // catch an int exception catch(int nCode) { printf("*** Exception: %d\n",nCode); } // catch everything catch(...) { printf("typeless exception\n"); 22 23 24 25 26 27 28 29 30 31 } return 0; 32 33 34 } The program’s output is shown in figure 6.3. We can see the format of the integer catch block with the data of the throwing statement 123. Figure 6.3: Output of the Exception Example 9.6.2015 Page 148 6.11 Computer Languages for Engineering - SS 15 OOP with Classes C++ is an object orientated programing language2 . So what is the concept of an object or a class? A class or an object combines data, called attributes, with functions, called methods. This can be described by so called UML3 An instance of a class, that is the realization of the class in memory, is created simply by assigning the class’s name followed by the constructors parameter list to a symbolic name, which than is the pointer or the reference to this instance. To access member attributes or methods of a class the dot notation is used, i.e. <instance>.<member>. The instance will be deleted if it runs out of scope. 6.11.1 Some UML Diagrams UML structure diagrams of the emphasize the things that must be present in the system being modeled. Since structure diagrams represent the structure they are used extensively in documenting the architecture of software systems. In our description of the examples we want to implement we use the Class Diagram which describes the structure of a system by showing the system’s classes, their attributes, and the relationships among the classes. A UML class diagram (see figure 6.4) consists of a rectangular box, which is divided into three sections. The fist section contents the class’s name. This name is written centered in bold letters. The second section contents the attribute’s names of the class and the third section contents the method’s names. Class Name attribute 1 attribute 2 method 1 method 2 Figure 6.4: A UML Class Diagram A UML note diagram (see figure 6.5) consists of a stylized note sheet which is filled with some information. Class Name attribute 1 attribute 2 method 1 method 2 Figure 6.5: A UML Note Diagram This is Class 1 A UML note diagram (see figure 6.6) will be assigned to an other component of the diagram scene with a simple line. Figure 6.6: A UML Note Diagram Assignment Figure 6.7 shows how to draw diagrams for inheriting classes. An arrow with a white filled arrowhead points from the inheriting class, the special class, to the inherited class, the base class. The attributes and the methods of the Base class are now available in the name space of the inheriting class, i.e. the special class now has the attributes attributB1, attributB2, attributS1 and attributS2. 2 This is only a simple Note Base Class attributeB1 attributeB2 methodB1 methodB2 Special Class attributeS1 attributeS2 methodS1 methodS2 Object orientated Programming is often used with the abbreviation OOP. Figure 6.7: A UML InThe Unified Modeling Language includes a set of graphic notation techniques to create visual models of software-intensive heritance Diagram systems. The Unified Modeling Language is an international standard see [7], UML 2.3 was formally released in May 2010. 3 E. Baeck 6.11. OOP WITH CLASSES Class 1 List A List B Page 149 * method 1 1 2..* Class A attribute A1 attribute A2 1..* Class B attribute B1 attribute B2 Figure 6.8 shows a aggregation and a composition. An aggregation is drawn by a white filled rhombus. An composition is drawn by a black filled rhombus. Aggregation and compositions describe a container or a list of several instances of an object, which are members of a main class. If for example a profile consists of several parts, the parts can be described as an composition, if a part only exists within a profile. If a part exists also without a profile, the parts are described within the profile with an aggregation. At the ends of the connecting lines the multiplicities are noted. The multiplicity gives the range of referFigure 6.8: A UML Diagram for a Composition and enced instances in the form from..to. For the Class A an Aggregation we have 2 up to infinite instances in an composition, therefor at the end of the line we can not have a multiplicity of zero. In our example we have exactly one instance of the class 1. On the other hand Class B is referred to Class 1 within an aggregation. In our example on instance of Class B can be reverenced by an undefined number of instances of Class 1. This is shown by the * icon. On the other hand the class 1 references at least on instance of the Class B. Otherwise the number of references is arbitrary. This is also shown by the * icon. method A 6.11.2 method B C++ Class A C++ class consists of a declaration and an implementation part. The declaration part often lives in a header file with the class’s name as prefix and h as suffix. The implementation of a class lives often in a code file which uses the classes name as a prefix and cpp as suffix. So for example if we want to implement a class with the name MyClass we put the declaration code into the file MyClass.h and the implementation code into the file MyClass.cpp.4 6.11.2.1 Declaration A class declaration is introduced by the key word class followed by the class’s name. If a class inherits base classes, the list of the base classes to inherit follows a colon. The declaration code is bracketed into curled parenthesizes. With the keys public:, protected: and private: the access permissions are set. Every item (attribute or method) following an access permission will get this permission. Within the declaration blocks members are declared like variables and methods are declared like functions. 4 This concept is rigorously applied in the implementation of the Microsoft Foundation Classes MFC. 9.6.2015 Page 150 Computer Languages for Engineering - SS 15 The access permissions to an item (attribute or method) are discussed as follows. If it is set to • public, all permissions are assigned, i.e. the item is accessible from outside of the class. • protected, only the class itself and classes, which are derived from this class, are allowed to access. • private, only the class itself is allowed to access the item. The following example shows a little class declaration with 2 public attributes and methods and 2 protected attributes and methods as well. The declaration starts with the classes’ constructor which is declared without return type. The name of the constructor is given by the classes’ name. In general a constructor can be used with some parameters for special initializations. The declaration of the constructor is followed by the declaration of the destructor. A destructor is declared without a return type. The name of the destructor is given by the classes’ name with a leading tilde ˜ character. The class does not inherit a base class. Listing 6.19: Declaring a Class 1 2 3 class MyClass: public MyBaseClass { public: 4 MyClass(parameter1, parameter2); ˜MyClass(); 5 6 // constructor // destructor 7 8 9 int int publicAttribut1; publicAttribut2; int int publicMethode1(); publicMethode2(); 10 11 12 13 protected: int protectedAttribut1; int protectedAttribut2; 14 15 16 17 int int 18 19 20 protectedMethode1(); protectedMethode2(); } if an attribute is declared as static, the attribute is an object attribute and not an attribute of the instance of an object, that means, that this attribute only exists once for an object. A non static attribute will be created however for each single instance of a class. E. Baeck 6.11. OOP WITH CLASSES 6.11.2.2 Page 151 Implementation If a class should be implemented in an cpp file, the header file with the class declarations has to be included in the header section of the cpp file. To avoid multiple inclusion of a header file, the header file should be protected with a macro #define discussed in section ??. The implementation file contents in general all implementations of the class’s methods and the implementation of the object attributes, i.e. the attributes with a static property. Every name of a part of a class, which should be implemented starts with the class name followed by the part’s name. The class name and the part name are separated by two colons. So an implementation of the above discussed declaration could be as follows. Listing 6.20: Implementing a Class 1 #include "MyClass.h" // include the declaration 2 3 4 5 6 7 // implement the constructor code MyClass::MyClass(parameter1, parameter2) { ... lines of code ... } 8 9 10 11 12 13 // implement the destructor code MyClass::˜MyClass() { ... lines of code ... } 14 15 16 17 18 19 // implement one of MyClass’s methods int MyClass::publicMethod1() { ... lines of code ... } 20 21 ... 9.6.2015 Page 152 6.11.2.3 Computer Languages for Engineering - SS 15 Structures and Classes In example 6.21 we discuss the relation between a struct and a class. We can say, that a struct is the simplest kind of a class. A struct only has public items and only attributes no methods. Therefor we can use a struct as a base class of a our class. In this example the struct is introduced using a typedef statement. With typedef a new name is introduced as an alias (lines 11-16). In our case we use MYSTRUCT instead of struct tagMYSTRUCT. Then we inherit from MYSTRUCT the class MyClass. This class comes with a public method f1 and a private method f2. In addition we introduce a public attribute lValue. The constructor performs all initializations. We initialize all inherited attributes and the additional class attribute. Inside the methods we throw an exception and stop the execution of the program. The attributes of a struct as well as the attributes of a class can be initialized using the fast and very low leveled function memset. A given byte, in this case the zero byte, is copied into the struct’s attributes. The main routine shows how to initialize the struct and how to copy an instance’s data into a struct’s data using the socalled cast operator. In this case this operator by default is available, because the class is inherited from the struct. The programs execution is stopped by an exception throw, which is done in one of the called methods. Listing 6.21: Stuctures and Classes 1 2 3 4 /* Example to show the relation between a structure and a class. The program is stopped by an exception */ 5 6 7 #include <stdio.h> #include <memory.h> 8 9 10 11 12 13 14 15 16 // declaring a structure with some // attributes typedef struct tagMYSTRUCT { int nValue; double dValue; float fArray[2]; } MYSTRUCT; 17 18 19 20 // inhereting this structure by a class class MyClass : public MYSTRUCT { 21 22 // methodes 23 24 25 26 27 // visible, accessable from outside public: // constructor is called if the instance is created MyClass(); // no return type 28 29 E. Baeck // destructor is called if the instance is deleted 6.11. OOP WITH CLASSES ˜MyClass(); 30 Page 153 // no return type 31 // method 1 int f1(); 32 33 34 35 36 37 38 // not accessable from outside private: // method 2 int f2(); 39 40 41 42 43 // accessable attributs public: long lValue; }; 44 45 46 47 48 49 // classname :: methode name // this is called the constructor MyClass::MyClass() { printf("create MyClass instance...\n"); 50 // initialization of the attributes nValue = 1; dValue = 2.; fArray[0] = 3.1; fArray[1] = 3.2; lValue = 4; 51 52 53 54 55 56 57 } 58 59 60 61 62 63 // this is called the destructor MyClass::˜MyClass() { printf("delete MyClass instance...\n"); } 64 65 66 67 68 69 // method 1 calls private method f2 // and throws an int exception int MyClass::f1() { f2(); 70 // integer exception (error!) throw (123); 71 72 73 } 74 75 76 77 78 79 // method 2 ony throws a char* exception int MyClass::f2() { throw("error form f2"); } 80 81 82 83 84 85 // this is the main routine, which is called // from the OS int main() { // type name 9.6.2015 Page 154 Computer Languages for Engineering - SS 15 MYSTRUCT m; 86 87 // print the struct’s content printf("1> m: %d %lf %f %f\n",m.nValue,m.dValue,m.fArray[0],m.fArray[1]); 88 89 90 // initialize a struct // address initvalue memset(&m , 0, sizeof(MYSTRUCT)); 91 92 93 94 // print the struct’s content printf("2> m: %d %lf %f %f\n",m.nValue,m.dValue,m.fArray[0],m.fArray[1]); 95 96 97 try { 98 99 // create an instance of the class MyClass MyClass c; 100 101 102 // print the instance’s content printf("3> c: %d %lf %f %f\n",c.nValue,c.dValue,c.fArray[0],c.fArray[1]); 103 104 105 // assign the instances data to the struct m = (MYSTRUCT)c; 106 107 108 // print the struct’s content printf("4> m: %d %lf %f %f\n",m.nValue,m.dValue,m.fArray[0],m.fArray[1]); 109 110 111 // call f1 to throw an exception c.f1(); 112 113 114 // print the struct’s content printf("3> %d %lf %f %f %ld\n",c.nValue,c.dValue,c.fArray[0],c.fArray[1],c.lValue); 115 116 } 117 118 // error handler for integer exceptions catch(int nError) { printf("*** int Error: %d\n",nError); } 119 120 121 122 123 124 // handle unspecified exceptions catch(...) { printf("*** unhandled error\n"); } 125 126 127 128 129 130 return 0; 131 132 } E. Baeck 6.11. OOP WITH CLASSES Page 155 Figure 6.9 shows the console output of program 6.21. We see, that printing the not initialized struct, we get some arbitrary data pattern from the stack, i.e. data patterns which by accident are found int the uninitialized variables. After that the struct is initialized, so that now the items values are zero. Then the instance of the class is allocated and we assign their values to the struct. At the end the program stops after a const char* exception is thrown, which obviously not is handled, so that the catch for unspecified exceptions will be used. Figure 6.9: Output of the Stuct and Classes Example 9.6.2015 Page 156 E. Baeck Computer Languages for Engineering - SS 15 7 Profile Example In this chapter we discuss the implementation of little software project in C++ which implements the thin-walled approach for some profile types. 7.1 Class Concept for the Tin-Walled Approach If we want to build up a class library for modeling and describing a problem, it’s recommended to start with a general base class, which should content all the common features of our classes. One of this features could be a general logging code, which should be implemented in every class of the library. The logging code should write informations into a log file. This feature should be available in every class of the class library. A second feature is an instance counter. Every instance, which is created should be counted by the base class. 7.2 Implementation Figure 7.1 shows the class tree of our profile class library. The common base class Base is inherited by a general profile class Profile. The Profile class contents all common features of a profile like in our case the model of the thin-walled approach with it’s nodes and elements. A special profile, in our example an U, H or L profile, then is implemented in the frame of it’s special class, i.e. UProfile, HProfile or LProfile class. To show the inheritance we also inherit from the Profile a specialized profile like the UProfile class. The differences between the UProfile class and the LProfile class for example are obviously the input parameters. A second difference between the specialized profiles are the methods to create the geometry of the profile, i.e. the method to create nodes and elements. If now the geometry of the specialized profile is created, the general method, to calculate the section values in the frame of the TWA is called from each element and from the base class Profile.1 The thin-walled profile approach is given by a set of nodes and elements which describes the profile part as lines with a constant thickness. The nodes, described in terms of a Node class, and the elements, described in terms of a Element class, are created dynamically. The pointer of the nodes and elements 1 Base classes are also called superclass or parent class 157 Page 158 Computer Languages for Engineering - SS 15 Node 2..* 1 Base 2..* 1 Element 1..* 1 Profile UProfile HProfile LProfile Figure 7.1: Class Hierarchy of the Profile Implementation are stored in the profile class Profile in a pointer array m_pNC and m_pNE. Because the element should be able to calculate it’s section values, we make the element a collection of it’s nodes. So a direct access to the node’s data is possible. If we would only store the node numbers in the element, an element would not be able to calculate it’s section values, because a direct access to the node’s coordinate values would not be possilble. E. Baeck 7.2. IMPLEMENTATION 7.2.1 Page 159 Base, the Base Class of all Classes A UML class diagram in figure 7.2 shows the concept of the common class Base. The class contents three object attributes, the name of the common log file, a print buffer and the instance counter. This attributes are declared with an static property. Besides the constructor and the destructor, the class has a method which will write the content of the message buffer to the screen and/or into a log file. The object method ResetLog deletes an allready existing log file. The following code gives the declaration header code of the class Base. Base − logFile[256]: char − msg[256]: char − counter:int + Base(): − − ˜Base(): − + appendLog(..): int + resetLog(..): int Figure 7.2: TWD’s Base Class Listing 7.1: Base Class of All Classes 1 2 #ifndef BASE_H_INCLUDED #define BASE_H_INCLUDED // protect against multiple // inclusion 3 4 5 6 // class declaration class Base { 7 8 public: // can be accessed from everywhere 9 10 11 12 13 // attributes // ========== // instance counter static int counter; // class attribute // logfile’s name static char logFile[256]; // old c string 14 15 16 17 18 19 // string buffer to log static char msg[256]; 20 21 22 23 24 25 26 27 28 // methodes // ======== Base(); ˜Base(); int appendLog(char* str); static void resetLog(); // // // // constructor to initalize destructor to free memory or to close files ... logging methode reset log }; #endif // BASE_H_INCLUDED 9.6.2015 Page 160 Computer Languages for Engineering - SS 15 The implementation code of the class Base is given below. Note, that object attributes must be implemented like methods. Therefore in line 6 to 8 we implement the object attributes like global variables. Listing 7.2: Base Class Implementation 1 2 3 #include "base.h" // specify the interface #include <stdio.h> // io functions of the c library #include <string.h> // string functions of the c library 4 5 6 7 8 // implement class attributes int Base::counter = 0; char Base::logFile[256] = {0}; char Base::msg[256] = {0}; // counter // log file // message buffer 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 // constructor // |class name // | | method, the contructor Base::Base() { #ifdef _COMMENT // initialization: it is sufficient to set the first item to 0 logFile[0] = 0; // ... as number logFile[0] = ’\0’; // ... or as character code #endif // set the default filename // |destination // | |source if (!logFile[0]) strcpy(logFile,"Base.log"); 24 // count the new instance counter++; 25 26 27 // create the output string // - print into an char array sprintf(msg,"> %d instance(s) created.\n",counter); 28 29 30 31 // - print the message into the file appendLog(msg); 32 33 34 } 35 36 37 38 39 40 // destructor Base::˜Base() { // create message sprintf(msg,"> Instance %d deleted.\n",counter); 41 // print message appendLog(msg); 42 43 44 // decrement the counter, because an instance is deleted counter--; 45 46 47 } 48 49 50 51 // AppendLog prints a message to the log and to the screen int Base::appendLog(char* pMsg) { E. Baeck 7.2. IMPLEMENTATION Page 161 // print to the screen printf("%s",pMsg); 52 53 54 // open the log // |file structure (like a channel number // |file name // | | open mode FILE* pHnd = fopen(logFile,"a"); 55 56 57 58 59 60 // check the file return // (same as if (pHnd == 0) if (!pHnd) return 0; 61 62 63 64 // print the message // | pointer to the file structure FILE fprintf(pHnd,"%s",msg); 65 66 67 68 // close the log fclose(pHnd); 69 70 71 return 1; 72 73 } 74 75 76 77 78 79 7.2.2 // ResetLog deletes an allready existing log-file void Base::resetLog() { remove(logFile); } Node Class for Model Nodes A UML class diagram in figure 7.3 shows the concept of the class Node. The class contents the node’s number and a double array to hold the instance’s coordinate data. Besides the constructor and the destructor, the class has a method List which will write the instances data to the log using the inherited method of the class Base::appendLog. To keep it simple all attributes and methods of the class are declared as public. In line 6 we can see, that the class inherits it’s base class Base (see section 7.2.1), therefore we have to include the Base class header file in line 4. Node + nNo: int + dx[2]: double + Node(..): − − ˜Node(): − + listData(..): void Figure 7.3: The Class Node Listing 7.3: Node Class’s Header 1 2 #ifndef NODE_H_INCLUDED #define NODE_H_INCLUDED 3 4 5 6 7 #include "Base.h" // include Base class’s header // | inherit from Base class Node: public Base { 8 9.6.2015 Page 162 9 Computer Languages for Engineering - SS 15 public: 10 // attribute int nNo; double dx[2]; 11 12 13 // node number // x,y values (vector) 14 15 16 17 public: // constructor: x,y optional Node(int nNo, double x = 0., double y = 0.); 18 // destructor ˜Node(); 19 20 21 // list the node’s data void listData(); 22 23 24 25 }; #endif // NODE_H_INCLUDED The implementation code of the class Node is given below. Note, that constructor is used with it’s parameters to initialize the node coordinates. If coordinate values are passed by the constructor, the values are assigned to the attributes. The constructor at the end will print the nodes data into the log file. Listing 7.4: Node Class’s Implementation 1 2 #include "node.h" #include <stdio.h> 3 4 5 6 7 8 9 10 11 // constructor // parameters call Base class Node::Node(int nNo, double x, double y) : Base() { // assign the coordinates nNo = nNo; dx[0] = x; dx[1] = y; 12 // print the data listData(); 13 14 15 } 16 17 18 // destructor (nothing to do) Node::˜Node() { } 19 20 21 22 23 24 25 26 27 // List method, prints Node’s data void Node::listData() { sprintf(msg, "> node: no = %2d, x = %10.3f y = %10.3f\n", nNo,dx[0],dx[1]); AppendLog(msg); } E. Baeck 7.2. IMPLEMENTATION 7.2.3 Page 163 Checking the Node Class To check the node class, we implement a little testing frame, which simple creates some Node instances, assign some data, print the Node data using it’s List method and deleting the Node instance at the end. Because we inherit the Base class features, the creation and the deletion of an instance is logged to the screen, so we can simple check this event. The testing code is given within the following main steps. • If we want to use the Node class, we have to declare it with the inclusion of it’s declaration header, this is done in the first line. • A simple Node instance is created by declaring a variable N1 of the type Node. • To initialize with special parameters we assign the return of the constructor to a variable of type Node, called N2. • Then we print the data of the created nodes N1 and N2 by calling their method List. This is done using the dot access (<variable name>.<method>). • After this we declare an address pointer pN32 of the to a Node instance using the Node* type. Because there is no valid instance address in N3 we initialize it with a zero value.3 • To assign a valid address of a Node instance to pN3, we create a Node instance dynamically by the usage of the new operator. The return of the new Node is an address and this address is assigned to our third Node* variable N3. • To print the data of the third node we have to use the arrow access because it’s a pointer variable (<variable name>-><method>). • After the printing of the content of the third Node instance we remove the instance from memory using the delete operator. To avoid the access to an invalid address, we check the address pointer against zero. The code of the testing frame is given below. Listing 7.5: Checkenvironment for the Node Class 1 2 #include "Node.h" #include <stdio.h> // load Node header // used for printf int main() { Node Node // the main routine 3 4 5 6 7 n1 = Node(1); n2 = Node(2,1.1,1.2); // standard coordinates used // use parameters to initialize 8 9 10 n1.listData(); n2.listData(); // list the data of node 1 // list the data of node 2 11 12 // creating 3rd Node instances dynamically 2 A prefix in the variable name of p shows, that it’s a pointer variable, which means, that the variable contents the address of the assigned data and not the data itself. 3 A zero address value means, that no memory is referenced to this address pointer. 9.6.2015 Page 164 Computer Languages for Engineering - SS 15 Node* 13 14 pN3 = 0; // address of a Node instance // 0: no memory assigned 15 // create an instance dynamically with the new operator printf("now create the instance for pN3\n"); pN3 = new Node(3,2.1,2.2); 16 17 18 19 // list data of node 3. note, we use the ’->’ operator // if the variable contents an adress an not the data itself pN3->listData(); 20 21 22 23 // delete the Node instance using the delete operator if (pN3) delete pN3; 24 25 26 } The Output of the above discussed code is given in figure 7.4. If an instance is created, we see in the log the message of the Base class constructor, which counts the instances. This output is followed by the output of the Node class constructor, which prints the coordinates of the Node instance. Here we can see, that the first node will get the standard zero values and the second instance will get the coordinates, which are passed to the constructor. Figure 7.4: Output of the Node Check Program Then the created instances’ data are printed using their List methods. Then we create the third instance and their data are printed like in the previous cases. After that the third Node instance is deleted explicitly and we see that the instance counter in the Base class destructor is decremented. Closing the program the destructor of the Base class of the statically created instances N1 and N2 is automatically executed. So we can see the decrementation of the instance counter down to zero. E. Baeck 7.2. IMPLEMENTATION 7.2.4 Page 165 Element Class for Model Elements A UML class diagram in figure 7.5 shows the concept of the class Element. The Element class get’s the pointer of it’s Node instances. If this pointers (the addresses of the Node instances) are known, the Element will be able to calculate all it’s sections values. So we can encapsulate this features into the Element. At the end to get the profiles total section values we simply have to add up all it’s Element’s values. Element + nNoX : int + nNo : int + dt : double + pN : Node** + dL : double To get a general implementation, we introduce a simple container for the Node instance pointers, which is a simple C array. In the general case we can have elements with an arbitrary number of nodes, so we introduce the number of Node pointers m_nNoX, which in our case will be 2. The Node instance pointers will be stored in the array m_pN, which we have to allocate dynamically. + dA : double Besides the Node instance pointers, our element need to know it’s thickness. The thickness is stored in the attribute m_dt. + initResults() : void The following attributes will hold the element’s section values. + dS[2] : double + dI[3] : double + Element(..) : − − ˜Element(): − + listDagatList() : void + setData() : int Figure 7.5: The Class Element • dL, the length of the element. This value is used to calculate the section values. • dA, the area of the element (see section C.1.1). • dS, the static moments of the element, [0]: Sy , [1]: Sz (see section C.1.2). • dI, the moments of inertia of the element, [0]: Iyy , [1]: Izz , [2]: Iyz , (see section C.1.3). Besides the constructor and the destructor, the class provides a method which will write the instances data to the log using the inherited method of the class Base::AppendLog. The method InitResults will initialize all internal and public items of the element, which are used to get and hold the result values. The method SetData will calculate the elements section values. The Element class provides only one constructor, i.e only one interface. The constructor comes with four parameters which will describe a linear element.4 • int nNo, the element’s number. • Node* pN1, the instance pointer to the element’s starting node. • Node* pN2, the instance pointer to the element’s terminating node. • double dt, the element’s thickness. We have to note, that all element’s methods suppose, that the input data are correct. Error checking will be performed later, if we create an general TWA profile. 4 Introducing further interfaces, the library can be extended by other element types. 9.6.2015 Page 166 Computer Languages for Engineering - SS 15 As we have seen in figure 7.5 the Element class may be considered as composition of Node instances. One element will have 2 nodes or more. Element Node 1 2..* Figure 7.6: UML-Diagram of the Element Composition Listing 7.6: Element Class’s Header 1 2 #ifndef ELEMENT_H_INCLUDED #define ELEMENT_H_INCLUDED 3 4 5 #include "Base.h" class Node; // include Base-Class’s header // we need a Node pointer 6 7 8 9 class Element: public Base { public: 10 // attributes int nNoX; // number of Nodes int nNo; // element number double dt; // thickness // adress of the Node instances pointer array (therefore 2 *) Node** pN; 11 12 13 14 15 16 17 // section attributes double dL; // double dA; // double dS[2]; // double dI[3]; // 18 19 20 21 22 element length element area element’s static moments element’s moments of inertia 23 // methods // constructor: only with all element data Element(int no,Node* n1, Node* n2, double t); 24 25 26 27 // destructor ˜Element(); 28 29 30 // list the data void listData(); 31 32 33 // initialize the results void initResults(); 34 35 36 // calculate results int setData(); 37 38 39 40 }; #endif // ELEMENT_H_INCLUDED The implementation code of the class Element is given below. Note, that constructor is used with it’s parameters to initialize the element data. E. Baeck 7.2. IMPLEMENTATION Page 167 Listing 7.7: Element Class’s Implementation 1 2 3 4 #include #include #include #include "element.h" <stdio.h> <memory.h> <math.h> // // // // declare the element class for printing used for memory access (memset) to calculate the root 5 6 #include "node.h" // we should know the Node too 7 8 9 10 11 // constructor Element::Element(int nNo, Node* pN1, Node* pN2, double dt) { nNoX = 2; // we only have elements with two Nodes 12 // assigning attributes this.nNo = nNo; this.dt = dt; 13 14 15 16 // Node address array // !!! we have allocate it !!! pN = new Node*[nNoX]; 17 18 19 20 // assign the Node pointers pN[0] = pN1; pN[1] = pN2; 21 22 23 24 // initialize the result attributs initResults(); 25 26 27 printf("> element %d created.\n",nNo); 28 29 } 30 31 32 33 34 35 // destructor: note: the Node instances have to be freed outside Element::˜Element() { printf("> element %d deleted.\n",nNo); } 36 37 38 39 40 41 // initialize the result attributs void Element::initResults() { dL = 0.; dA = 0.; 42 // initialize an array // | destination (array’s address) // | | byte to copy // | | | number of byte to copy memset((void*)&dS[0],0,2*sizeof(double)); // dS == &dS[0] memset((void*)&dI[0],0,3*sizeof(double)); 43 44 45 46 47 48 49 } 50 51 52 53 54 // list the elements data void Element::listData() { // print header 9.6.2015 Page 168 Computer Languages for Engineering - SS 15 sprintf(msg,"element %2d, t = %10.2f\n",nNo,dt); appendLog(msg); 55 56 57 // print element’s Node data for (int i=0;i<nNoX;i++) { pN[i]->listData(); } 58 59 60 61 62 63 // list result values sprintf(msg," L..........: appendLog(msg); sprintf(msg," A..........: appendLog(msg); sprintf(msg," Sx.........: appendLog(msg); sprintf(msg," Sy.........: appendLog(msg); sprintf(msg," Ixx........: appendLog(msg); sprintf(msg," Iyy........: appendLog(msg); sprintf(msg," Ixy........: appendLog(msg); 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 %12.3f mm\n",dL); %12.3f cmˆ2\n",dA/1.e2); %12.3e cmˆ3\n",dS[0]/1.e3); %12.3e cmˆ3\n",dS[1]/1.e3); %12.3e cmˆ4\n",dI[0]/1.e4); %12.3e cmˆ4\n",dI[1]/1.e4); %12.3e cmˆ4\n",dI[2]/1.e4); } 80 81 82 83 84 85 86 // calculate all supported the section values of an element int Element::setData() { // introduce some helper variables double dxc[2]; // element center coordinates double dLp[2]; // projected length of the element 87 88 89 90 91 92 // ... and calculate them for (int i=0;i<2;i++) { // center coordinates dxc[i] = (pN[1]->dx[i] + pN[0]->dx[i])/2.; 93 // projected length dLp[i] = pN[1]->dx[i] - pN[0]->dx[i]; 94 95 96 } 97 98 99 // calculate the length of the element dL = sqrt(dLp[0]*dLp[0] + dLp[1]*dLp[1]); 100 101 102 // calculate the area dA = dL * dt; 103 104 105 106 107 108 // calculate the static moment for (int i=0;i<2;i++) { dS[i] = dxc[(i+1)%2]*dA; } 109 110 E. Baeck // calculate the moment of inertia 7.2. IMPLEMENTATION Page 169 // Ixx, Iyy int j; for (int i=0;i<2;i++) { j = (i+1)%2; // index on right side dI[i] = ( (pow(dLp[j],2)/12. + pow(dxc[j],2)) * dA); } // Ixy dI[2] = (dLp[0]*dLp[1]/12. + dxc[0]*dxc[1])*dA; 111 112 113 114 115 116 117 118 119 120 121 return 1; 122 123 7.2.5 // 1==true -> ok! } Checking the Element Class To check the element class, we implement a little testing frame, which simple creates an element to calculate the section values of a flat steel Fl200x4. The origin of the used coordinate system is in the center of the element. The element is orientated into the vertical direction, i.e. in the direction of the second coordinate. The program executes the following steps. • Allocate the Node instances. We create 2 node 1 at the position (0, −100) and node 2 at (0, 100). • Allocate the Element instances. Element 1 gets the Node instance pointer 1 and 2 and the element’s thickness. • Calculate the section values. • List the element’s data. • Delete all created instances. The code of the testing frame is given below. Listing 7.8: Check Environment for the Node Class 1 2 3 #include <stdio.h> #include "node.h" #include "element.h" // what is a Node // what is an Element 4 5 6 7 8 9 10 11 12 int main() { // create a flat steel Fl 200x4 // we need 2 Nodes for one Element // the are allocated at the heap // No x y Node* pN1 = new Node(1, 0., -100.); Node* pN2 = new Node(2, 0., 100.); 13 14 15 // create one Element on the heap Element* pE1 = new Element(1,pN1,pN2,4.); 9.6.2015 Page 170 Computer Languages for Engineering - SS 15 16 // calculate the values pE1->setData(); 17 18 19 // list element’s data pE1->listData(); 20 21 22 // clear the memory delete pN1; delete pN2; delete pE1; return 0; 23 24 25 26 27 28 } The Output of the above discussed code is given in figure 7.4. If an instance is created, we see in the log the message of the Base class constructor, which counts the instances. This output is followed by the output of the Node class constructor, which prints the coordinates of the Node instance. After this the Element instance is created. The constructor of the element shows us the element’s number. After the calculation of the section values Figure 7.7: Output of the Element Check Program is done, the element’s List method is called, i.e. we see the element’s number and it’s thickness. Then all node data of this element are listed and at the end we see the calculated section values. Because the element is vertical we only get an area of A = L · t = 200 · 4 mm 2 = 8cm 2 and a moment of inertia of Ixx = E. Baeck 1 1 · L3 · t = · 2003 · 4 mm 2 = 266, 7cm 4 12 12 7.2. IMPLEMENTATION 7.2.6 Page 171 Profile Class for Model Profiles A UML class diagram in figure 7.8 shows the concept of the class Profile. The class contents the profile’s name and a container for the profile’s nodes and profile’s elements. The containers are build by a simple dynamical array and an integer, which holds the length of the array. The methods of the class creates the containers with a specified length. A Node instance and an Element instance can be added by an specific Add function. To check the node and element information Check methods are implemented, which checks the existents of the referenced node and element number. Besides the constructor and the destructor, the class has a method which will write the instances data to the log using the inherited method of the class Base::appendLog. To keep it simple all attributes and methods are set to public. We have the following class attributes. The meaning of the section values is discussed in section C.1. • pName[256], the profile’s name • dA, total area • dS[2], total static moment’s Profile + pName[256] : char + dA : double + dS[2] : double + de[2] : double + dIu[3] : double + dIc[3] : double + dIm[2] : double + dAlpha : double + pNC : Node** + nNC : int + pEC : Element** + nEC : int + Profile(..) : − − ˜Profile(): − + addNodeContainer(..) : int + addElementContainer(..) : int + addNode(Node* pN) : int + addElement(..) : int + listData() : void • de[2], center of mass coordinates • dIu[3], moment of inertia in user coordinates + deleteNodes() : int + deleteElements() : int + checkNode(..) : int • dIc[3], moment of inertia in center of mass coordinates + checkElement(..) : int • dIm[2], moment of inertia in main coordinates + resetSectionValues() : void + getSectionValues() : int • dAlpha, main axis angle + listSectionValues() : void • pNC, pointer to Node instance array + pTrans() : double • nNC, length of node container Figure 7.8: The Class Profile • pEC, pointer to Element instance array • nEC, length of element container 9.6.2015 Page 172 Computer Languages for Engineering - SS 15 The profile class has the following methods. methode parameter type comment constructor Profile cName const char profile’s name ˜Profile destructor addNodeContainer create the node container nLength int length of the array create the element container addElementContainer nLength int length of the array adds a Node pointer into it’s array slot, error checking is done addNode pN Node* instance pointer to store adds an Element pointer into it’s array slot, error checking is done addElement pE Element* instance pointer to store listData prints all data of the profile deleteNodes deletes all nodes and their container deleteElements deletes all elements and their container checkNode check the data of a node pN Node* Node instance to check check the data of an element checkElement pE Element* Element instance to check resetSectionValues initialize all section values getSectionValues print the section values pTrans calculate the principal values of the moment of inertia and the rotation angle In this version of the implementation the input data is checked only by the Profile class, because in this context it does not sence to access nodes and elements directly. So all input data runs through the profile’s input functions. If errors are detected we give up running backward using error codes. If an item is checked and an error is detected an exception is thrown. So the applying routine has to catch the exception using the code of the profile within a try block.. E. Baeck 7.2. IMPLEMENTATION Page 173 The declaration of the class Profile is given by the following listing. Listing 7.9: Profile Class’s Header 1 2 #ifndef PROFILE_H_INCLUDED #define PROFILE_H_INCLUDED 3 4 5 6 #include "Base.h" class Node; class Element; // we should know something about the Base // we use Node pointers // we use Element pointers 7 8 9 10 11 // a profile’s class class Profile: public Base { public: 12 13 14 // interface - constructor Profile(const char* cName); 15 16 17 // destructor ˜Profile(); 18 19 20 // create a Node container int addNodeContainer(int nLength); 21 22 23 // create a Element container int addElementContainer(int nLength); 24 25 26 // add a Node int addNode(Node* pN); 27 28 29 // add a Element int addElement(Element* pE); 30 31 32 // list all values void listData(); 33 34 35 // delete/clear all Nodes int deleteNodes(); 36 37 38 // delete/clear all Nodes int deleteElements(); 39 40 41 // check the Node instance int checkNode(Node* pN); 42 43 44 // check the Element instance int checkElement(Element* pE); 45 46 47 48 // methods to calculate the section values // - reset void resetSectionValues(); 49 50 51 // - calculate the section values int getSectionValues(); 52 9.6.2015 Page 174 Computer Languages for Engineering - SS 15 // - list them void listSectionValues(); 53 54 55 // - main axis transformation double mTrans(); 56 57 58 // ----------------------// attributes of the class // - profile’s name char pName[256]; 59 60 61 62 63 64 65 66 67 68 69 70 71 // - section values double dA; double dS[2]; double de[2]; double dIu[3]; double dIc[3]; double dIm[2]; double dAlpha; // // // // // // // // Node container Node** pNC; int nNC; // Node instance array // array’s dimension // Element container Element** pEC; int nEC; // Element instance array // array’s dimension area static moment center of mass coordinates M o I in user coordinates M o I in main CS M o I main values rotation angle 72 73 74 75 76 77 78 79 80 81 }; #endif // PROFILE_H_INCLUDED The implementation code of the class Profile is given below. Note, that constructor is used with it’s parameters to initialize the Profile data. Listing 7.10: Profile Class’s Implementation 1 2 3 4 5 6 /* Implementation of the Profile class */ #include "profile.h" #include "node.h" #include "element.h" 7 8 9 10 #include <stdio.h> #include <string.h> #include <math.h> // standard io (printing) // to use string functions // to use math functions 11 12 13 14 15 16 17 // constructor Profile::Profile(const char* pName): Base() { // copy name // dest. source strcpy(pName,pName); 18 19 20 21 E. Baeck // initialize the container pNC = 0; // for Nodes nNC = 0; 7.2. IMPLEMENTATION pEC = 0; nEC = 0; 22 23 Page 175 // for Elements 24 // reset the results resetSectionValues(); 25 26 27 } 28 29 30 31 32 33 34 // destructor Profile::˜Profile() { // first delete the content deleteNodes(); deleteElements(); 35 // delete the containers if (pNC) delete [] pNC; if (pEC) delete [] pEC; 36 37 38 39 // for Nodes // for Elements } 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 // reset the section values void Profile::resetSectionValues() { dA = 0.; // destination pointer // | | byte to copy // | | | number of bytes to copy memset((void*)dS, 0, sizeof(double)*2); memset((void*)de, 0, sizeof(double)*2); memset((void*)dIu, 0, sizeof(double)*3); memset((void*)dIc, 0, sizeof(double)*3); memset((void*)dIm, 0, sizeof(double)*2); dAlpha = 0.; } 55 56 57 58 59 60 // delete all Nodes int Profile::deleteNodes() { // check the container if (!pNC) return 0; 61 // delete the Node instances for (int i=0;i<nNC;i++) { Node* pN = pNC[i]; // get the instance "i" if (pN) delete pN; // delete instance "i", if available } 62 63 64 65 66 67 68 delete [] pNC; pNC = 0; nNC = 0; 69 70 71 // delete the container // NO conaitner // NO length 72 return 1; 73 74 } 75 76 77 // delete all Elements int Profile::deleteElements() 9.6.2015 Page 176 78 Computer Languages for Engineering - SS 15 { // check the container if (!pEC) return 0; 79 80 81 // delete the Node instances for (int i=0;i<nEC;i++) { Element* pE = pEC[i]; // get the instance "i" if (pE) delete pE; // delete instance "i", if available } 82 83 84 85 86 87 88 delete [] pEC; pEC = 0; nEC = 0; 89 90 91 // delete the container // NO conaitner // NO length 92 return 1; 93 94 } 95 96 97 98 99 100 101 // list all profile attributs void Profile::listData() { sprintf(msg,"Profile ’%s’ (node space: %d, element space: %d)\n", pName,nNC,nEC); AppendLog(msg); 102 // list the nodes data for (int i=0;i<nNC;i++) { Node* pN = pNC[i]; if (pN) pN->listData(); } 103 104 105 106 107 108 109 // list the elements data for (int i=0;i<nEC;i++) { Element* pE = pEC[i]; if (pE) pE->listData(); } 110 111 112 113 114 115 116 // list the result values listSectionValues(); 117 118 119 } 120 121 122 123 124 125 126 127 128 129 130 131 132 133 // list the section values void Profile::listSectionValues() { appendLog((char*)"Section values:\n"); sprintf(msg," area.....................: appendLog(msg); sprintf(msg," static moment S_y........: appendLog(msg); sprintf(msg," static moment S_z........: appendLog(msg); sprintf(msg," center of mass y.........: appendLog(msg); sprintf(msg," center of mass z.........: E. Baeck %10.2f cmˆ2\n",dA/1.e2); %10.2f cmˆ3\n",dS[0]/1.e3); %10.2f cmˆ3\n",dS[1]/1.e3); %10.2f mm\n",de[0]); %10.2f mm\n",de[1]); 7.2. IMPLEMENTATION appendLog(msg); appendLog((char*)"Moment of Inertia in user cooordinates:\n"); sprintf(msg," I_yy.....................: %10.2f cmˆ4\n",dIu[0]/1.e4); appendLog(msg); sprintf(msg," I_zz.....................: %10.2f cmˆ4\n",dIu[1]/1.e4); appendLog(msg); sprintf(msg," I_yz.....................: %10.2f cmˆ4\n",dIu[2]/1.e4); appendLog(msg); appendLog((char*)"Moment of Inertia in centroid cooordinates:\n"); sprintf(msg," I_yy.....................: %10.2f cmˆ4\n",dIc[0]/1.e4); appendLog(msg); sprintf(msg," I_zz.....................: %10.2f cmˆ4\n",dIc[1]/1.e4); appendLog(msg); sprintf(msg," I_yz.....................: %10.2f cmˆ4\n",dIc[2]/1.e4); appendLog(msg); appendLog((char*)"Moment of Inertia in main cooordinates:\n"); sprintf(msg," I_eta....................: %10.2f cmˆ4\n",dIm[0]/1.e4); appendLog(msg); sprintf(msg," I_zeta...................: %10.2f cmˆ4\n",dIm[1]/1.e4); appendLog(msg); sprintf(msg," alpha....................: %10.2f ˆ A ◦ \n",dAlpha*45./atan(1.)); appendLog(msg); 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 Page 177 } 157 158 159 160 161 162 // allocate the Node space int Profile::addNodeContainer(int nLength) { // delete the old container deleteNodes(); 163 // create the Node array pNC = new Node* [nLength]; if (!pNC) return 0; // no memory available 164 165 166 167 // initialize the memory with Null (0) // destination address // | | byte to copy memset((void*)pNC,0,sizeof(Node*)*nLength); 168 169 170 171 172 // save the length nNC = nLength; 173 174 175 return nLength; 176 177 } 178 179 180 181 182 183 // allocate the Element space int Profile::addElementContainer(int nLength) { // delete the old container deleteElements(); 184 185 186 187 // create the Element array pEC = new Element* [nLength]; if (!pEC) return 0; // no memory available 188 189 // initialize the memory with Null (0) 9.6.2015 Page 178 Computer Languages for Engineering - SS 15 // destination address // | | byte to copy memset((void*)pEC,0,sizeof(Element*)*nLength); 190 191 192 193 // save the length nEC = nLength; 194 195 196 return nLength; 197 198 } 199 200 201 202 203 204 // Add an element and it’s Nodes int Profile::addElement(Element* pE) { // container available if (!pEC) throw "*** Error: no element container!\n"; 205 // check element address if (!pE) throw "*** Error: no element pointer\n"; 206 207 208 // check the instance and throw an exception, // if there is an error checkElement(pE); 209 210 211 212 // add the element pEC[pE->nNo -1] 213 214 = pE; 215 // add the element’s Nodes pNC[pE->pN[0]->nNo -1] = pE->pN[0]; pNC[pE->pN[1]->nNo -1] = pE->pN[1]; 216 217 218 219 return 1; 220 221 } 222 223 224 225 226 227 // check an element instance int Profile::checkElement(Element* pE) { // check instance pointer if (!pE) throw "*** Error: invalid element pointer!"; 228 229 230 if (pE->nNo < 1 || pE->nNo > nEC) throw "*** Error: invalid element number!"; 231 232 233 234 235 236 237 238 239 240 // check the Node instances for (int i=0;i<2;i++) { Node* pN = pE->pN[i]; if (!pN) throw "*** Error: Node instance not found!"; if (pN->nNo < 1 || pN->nNo > nNC) throw "*** Error: Node number invalid!"; } 241 242 243 244 245 E. Baeck // check Node numbers if (pE->pN[0]->nNo == pE->pN[1]->nNo) throw "*** Error: Invalid Node numbers"; 7.2. IMPLEMENTATION return 0; 246 247 Page 179 } 248 249 250 251 252 253 // calculate the section values int Profile::getSectionValues() { // initialization resetSectionValues(); 254 // sum over all elements for (int i=0;i<nEC;i++) { // get element Element* pE = pEC[i]; 255 256 257 258 259 260 // element exists? if (!pE) continue; 261 262 263 // calculate Element results pE->setData(); 264 265 266 // sum up the values dA += pE->dA; for (int j=0;j<2;j++) for (int j=0;j<3;j++) 267 268 269 270 dS[j] += pE->dS[j]; dIu[j] += pE->dI[j]; } 271 272 // calculate the center of mass de[0] = dS[1]/dA; de[1] = dS[0]/dA; 273 274 275 276 // calculate the principal values pTrans(); 277 278 279 return 1; 280 281 } 282 283 284 285 286 287 288 289 290 // principal axis transformation // return: angle double Profile::pTrans() { // M o I in CCS (center of mass) dIc[0] = dIu[0] - de[1]*de[1]*dA; dIc[1] = dIu[1] - de[0]*de[0]*dA; dIc[2] = dIu[2] - de[0]*de[1]*dA; 291 292 293 294 295 // helper values double dIdel = dIc[0] - dIc[1]; double dIsum = dIc[0] + dIc[1]; double dIsqr = sqrt(dIdel*dIdel +4.*dIc[2]*dIc[2]); 296 297 298 299 // M o I in principal coordinate system dIm[0] = 0.5*(dIsum + dIsqr); dIm[1] = 0.5*(dIsum - dIsqr); 300 301 // calcualate the rotation angle 9.6.2015 Page 180 Computer Languages for Engineering - SS 15 dAlpha = 0.5*atan2(2.*dIc[2],dIdel); 302 303 return dAlpha; 304 305 7.2.7 } Checking the Profile Class To check the profile class, we implement a little testing frame, which simple creates an profile to calculate the section values of a flat steel Fl200x4. The origin of the used coordinate system is in one endpoint and the profile is rotated by 45◦ . We have used this example already to check the element in section 7.2.5. The program executes the following steps. • Allocate the Profile instances. We pass the name of the profile. • Allocate the containers. We need a node space of two and an element space of one. • Create the Node instances. To check it flexible we use macros to enable or disable this part of code. Because in our case there is no macro definition, we get the code within the #else branch. • Create the Element instance. We pass the element number, the Node instance pointers and the element’s thickness. • Add the element to the profile. Adding the element to the profile, the element’s data are checked, so if there will be detected an error, the routine will be canceled by an exception, we have to handle. • Calculate the section values and • list the profiles data. • At the end we should not forget to clear the memory, deleting the Profile instance. All the above discussed steps should be done within a try block, so that we can handle detected errors in the following catch blocks. • The first catch block will handle our exceptions, because we pass a const char* to the exception throwing them. • The second catch block will handle all other exceptions. So, if we divide by zero or if we have forgotten to check one case, the program is not crashing, an unspecified exception is thrown. E. Baeck 7.2. IMPLEMENTATION Page 181 The code of the testing frame is given below. Listing 7.11: Check Environment for the Profile Class 1 2 3 4 5 #include #include #include #include #include <stdio.h> <math.h> "profile.h" "node.h" "element.h" // we start with the profile // and will need some nodes // and some elements 6 7 8 9 10 11 12 13 14 15 // #define _CENTERED // disabled, to get the #else branch int main() { // run the code using an exception handler, to // handle errors detected by the error checker try { // create the profile Profile* pProf = new Profile("Fl200x4"); 16 // add containers pProf->addNodeContainer(2); pProf->addElementContainer(1); 17 18 19 20 21 22 23 24 // create the Nodes (double symmetric) #ifdef _CENTERED Node* pN1 = new Node(1,0., 100.); Node* pN2 = new Node(2,0.,-100.); 25 26 #elif 27 28 29 _SHIFTED // create the Nodes, shifted Node* pN1 = new Node(1,0., 0.); Node* pN2 = new Node(2,0.,-200.); 30 31 #else // create the Nodes, shifted and rotated Node* pN1 = new Node(1,0., 0.); Node* pN2 = new Node(2, 200./sqrt(2.),-200./sqrt(2.)); 32 33 34 35 #endif // create the Elements Element* pE1 = new Element(1,pN1,pN2,4.); 36 37 38 // add element pProf->addElement(pE1); 39 40 41 // calculate the section values pProf->getSectionValues(); 42 43 44 // list profile data pProf->listData(); 45 46 47 // delete the profile delete pProf; 48 49 50 } 51 52 // handle the errors throwing string exceptsions 9.6.2015 Page 182 Computer Languages for Engineering - SS 15 catch(const char* str) { printf("Exception: %s\n",str); } 53 54 55 56 57 // handle all unspecified exceptions catch(...) { printf("Unknow exception!"); } return 0; 58 59 60 61 62 63 64 } The Output of the above discussed code is given in figure 7.9. We see, that we create four instances (two nodes, one element and one profile). We see the test printing of there constructors. Then after having assembled the profile the section values are calculated and printed. We see, that we get the same area as in the case of the element check (section 7.2.5). The center of mas we get at the center of the element, which is 200 1 = 70.71 ex = · √ 2 2 At the end we see, that we will get the same values for the moment of inertia in the main coordinate system as we have calculated in the case of the element check. The calculated rotation angle is -45◦ as expected. E. Baeck Figure 7.9: Output of the Element Check Program 7.2. IMPLEMENTATION 7.2.8 Page 183 H-Profile Class for Model Profiles A UML class diagram in figure 7.10 shows the concept of the class HProfile. The class is derived from the class Profile, which itself is derived from the class Base. The only attributes the class HProfile gets, are the parameters to describe the profile’s geometry. Base Profile • dh, the height, • dw, the width HProfile • dt, the flange thickness + dh : double + dw : double • ds, the web thickness + dt : double The HProfile instance is created with the call of the constructor, so we put the data checking and the creation of the profile’s geometry data into the constructor. Therefore with one statement all the things to do are done. We pass the name of the profile and it’s geometry parameter to the constructor. + ds : double + HProfile(..) : − − ˜HProfile(): − + check(): int + create() : int + listData() : void The HProfile class therefore will have the following methods. methode parameter type Figure 7.10: The Class HProfile comment constructor HProfile cName const char profile’s name, send to the base class dh double profile’s height dw double profile’s width dt double profile’s flange thickness dw double profile’s weg thickness ˜Profile destructor check checks the parameter passed by the constructor call create create the geometry of the profile in terms of nodes and elements listData list the profile’s data calling the List method of the base class too If errors are detected we give up running backward using error codes. If an item is checked and an error is detected an exception is thrown. So the applying routine has to catch the exception using the code of the profile within a try block.. 9.6.2015 Page 184 Computer Languages for Engineering - SS 15 The declaration of the class HProfile is given by the following listing. Listing 7.12: HProfile Class’s Header 1 2 #ifndef HPROFILE_H_INCLUDED #define HPROFILE_H_INCLUDED 3 4 #include "profile.h" // we have to know something about Profile 5 6 7 8 class HProfile : public Profile { public: 9 HProfile(const char* pName, double dh, double dw, double dt, double ds); 10 11 12 13 14 // // // // // profile’s name height width thickness of the flanges thickness of the web 15 int check(); int create(); void listData(); 16 17 18 // check the parameters // create the profile // list the data 19 // attributes: profile parameters double dh; // height double dw; // width double dt; // thickness of the flanges double ds; // thickness of the web 20 21 22 23 24 25 26 }; #endif // HPROFILE_H_INCLUDED The implementation code of the class HProfile is given below. Note, that constructor is used with it’s parameters to initialize the HProfile data. Listing 7.13: HProfile Class’s Implementation 1 2 3 4 #include #include #include #include "hprofile.h" "node.h" "element.h" <stdio.h> // we need the HProfile header // we will create Nodes // and Elements 5 6 7 8 9 10 11 12 13 14 15 // constructor don’t forget to call the base classes constructor HProfile::HProfile(const char* pName, double dh, double dw, double dt, double ds) : Profile(pName) { // assign the input data this.dh = dh; this.dw = dw; this.dt = dt; this.ds = ds; 16 17 18 // check the data check(); 19 20 21 E. Baeck // create the profile create(); 7.2. IMPLEMENTATION 22 Page 185 } 23 24 25 26 27 28 29 30 31 32 33 // Check the H-profile’s data, throw an exception if something is not ok int HProfile::check() { double dEps = 0.5; if (dt < dEps) throw "error: dt invalid!"; if (ds < dEps) throw "error: ds invalid!"; if (dw < 2.*ds) throw "error: dw invalid!"; if (dh < 3.*dt) throw "error: dh invalid!"; return 1; } 34 35 36 37 38 39 40 // create the geometry int HProfile::create() { // add node and element space addNodeContainer(6); // for 6 nodes addElementContainer(5); // for 5 elements 41 // create nodes Node* pN[6]; pN[0] = new Node(1,-dw/2., (dh-dt)/2.); pN[1] = new Node(2, 0., (dh-dt)/2.); pN[2] = new Node(3, dw/2., (dh-dt)/2.); pN[3] = new Node(4,-dw/2.,-(dh-dt)/2.); pN[4] = new Node(5, 0.,-(dh-dt)/2.); pN[5] = new Node(6, dw/2.,-(dh-dt)/2.); 42 43 44 45 46 47 48 49 50 // create elements Element* pE[5]; pE[0] = new Element(1,pN[0],pN[1],dt); pE[1] = new Element(2,pN[1],pN[2],dt); pE[2] = new Element(3,pN[3],pN[4],dt); pE[3] = new Element(4,pN[4],pN[5],dt); pE[4] = new Element(5,pN[1],pN[4],ds); 51 52 53 54 55 56 57 // bottom flange // top flange // web 58 // add elements to the profile for (int i=0;i<5;i++) addElement(pE[i]); 59 60 61 return 1; 62 63 } 64 65 66 67 68 69 70 71 72 73 74 75 76 77 // list all profile data void HProfile::listData() { // list input data sprintf(msg,"profile name........: appendLog(msg); sprintf(msg," height.............: appendLog(msg); sprintf(msg," width..............: appendLog(msg); sprintf(msg," flange thickness...: appendLog(msg); sprintf(msg," web thickness......: %s\n",pName); %10.2f mm\n",dh); %10.2f mm\n",dw); %10.2f mm\n",dt); %10.2f mm\n",ds); 9.6.2015 Page 186 Computer Languages for Engineering - SS 15 appendLog(msg); 78 79 // call the base class’s methode! Profile::listData(); 80 81 82 7.2.9 } Checking the HProfile Class To check the profile class, we implement a little testing frame, which simple creates an profile to calculate the section values of a HEA100 standard profile. The program executes the following steps. • Allocate the HProfile instances, specifying all profile parameters. • Calculate the section values and • list the profiles data. • At the end we should not forget to clear the memory, deleting the HProfile instance. All the above discussed steps should be done within a try block, so that we can handle detected errors in the following catch blocks. • The first catch block will handle our exceptions, because we pass a const char* to the exception throwing them. • The second catch block will handle all other exceptions. So, if we divide by zero or if we have forgotten to check one case, the program is not crashing, an unspecified exception is thrown. The code of the testing frame is given below. Listing 7.14: Check Environment for the HProfile Class 1 2 3 #include <stdio.h> #include <math.h> #include "hprofile.h" // we start with the hprofile 4 5 6 7 8 9 10 11 12 int main() { // run the code using an exception handler, to // handle errors detected by the error checker try { // check the H-Profile class HProfile* pProf = new HProfile("HEA-100",96.,100.,8.,5.); 13 14 15 // calculate the section values pProf->getSectionValues(); 16 17 18 19 E. Baeck // list profile data pProf->listData(); 7.2. IMPLEMENTATION Page 187 delete pProf; 20 } 21 22 // handle the errors throwing string exceptsions catch(const char* str) { printf("Exception: %s\n",str); } 23 24 25 26 27 28 // handle all unspecified exceptions catch(...) { printf("Unknow exception!"); } return 0; 29 30 31 32 33 34 35 } The Output of the above discussed code is given in figure 7.11. We see the input data used and the coordinates of the created nodes. Because the output given is rather long, we have split this output into the input section (upper picture) and the result section (lower picture). The lower figure shows the calculated section values. Because the origin is in the center of the H-profile, we will not get any eccentricity, i.e. static moments. If we use this symmetric origin, we see that the moments of inertia with respect to our three coordinaten systems are the same. Figure 7.11: Output of the Element Check Program 9.6.2015 Page 188 Computer Languages for Engineering - SS 15 In the following table we compare the calculated values with the values we will find in the book of standard profiles.5 We see, that the values are on the secure side, i.e. the calculated values are smaller than the exact value and the error in this case is less than 4%. value exact TWA error [cm 2 ] [cm 2 ] area 21.2 20.4 −3.4 [cm 4 ] [cm 4 ] 5 [%] [%] big moment of inertia 349 338 −3.3 small moment of inertia 134 133 -0.8 Areas are given in cm 2 and moments of inertia are given in cm 4 , according to German standard table books. E. Baeck Part III Appendix 189 Appendix A The Console’s Friends If you want to work with a console window you should know the console’s best friends, the commands to navigate through the folder tree, the commands to create, delete and chance directories, the commands to setup paths and environment variables and commands to copy and delete files. And if you want to be happy using the console window, it’s recommended to know something about assembling so-called batch files, which are in the most simple kind only a list of commands which should be executed without typing them a dozen times. The console window can be created with the command cmd from the execution input field in the start menu of windows. A list of a lot of console commands is given by the command help1 . If you need a specific information related to a special command, you will get this information calling the help with the command’s name as parameter. If the command is not part of the command line help you can run the command with the option /?. An A-Z Index of the Windows CMD command line ADDUSERS Add or list users to/from a CSV file ADmodcmd Active Directory Bulk Modify ARP Address Resolution Protocol ASSOC Change file extension associations ASSOCIAT One step file association AT Schedule a command to run at a specific time ATTRIB Change file attributes b BCDBOOT Create or repair a system partition BCDEDIT Manage Boot Configuration Data BITSADMIN Background Intelligent Transfer Service BOOTCFG Edit Windows boot settings BROWSTAT Get domain, browser and PDC info c CACLS CALL CERTREQ CERTUTIL CD CHANGE CHKDSK 1 Change file permissions Call one batch program from another Request certificate from a certification authority Utility for certification authority (CA) files and services Change Directory - move to a specific Folder Change Terminal Server Session properties Check Disk - check and repair disk problems The list of commands is given in the language of the computer. 191 Page 192 CHKNTFS CHOICE CIPHER CleanMgr CLEARMEM CLIP CLS CLUSTER CMD CMDKEY COLOR COMP COMPACT COMPRESS CON2PRT CONVERT COPY CSCcmd CSVDE Computer Languages for Engineering - SS 15 Check the NTFS file system Accept keyboard input to a batch file Encrypt or Decrypt files/folders Automated cleanup of Temp files, recycle bin Clear memory leaks Copy STDIN to the Windows clipboard Clear the screen Create or configure a cluster Start a new CMD shell Manage stored usernames/passwords Change colors of the CMD window Compare the contents of two files or sets of files Compress files or folders on an NTFS partition Compress one or more files Connect or disconnect a Printer Convert a FAT drive to NTFS Copy one or more files to another location Client-side caching (Offline Files) Import or Export Active Directory data d DATE Display or set the date DEFRAG Defragment hard drive DEL Delete one or more files DELPROF Delete user profiles DELTREE Delete a folder and all subfolders DevCon Device Manager Command Line Utility DIR Display a list of files and folders DIRUSE Display disk usage DISKPART Disk Administration DISKSHADOW Volume Shadow Copy Service DISKUSE Show the space used in folders DOSKEY Edit command line, recall commands, and create macros DriverQuery Display installed device drivers DSACLs Active Directory ACLs DSAdd Add items to active directory (user group computer) DSGet View items in active directory (user group computer) DSQuery Search for items in active directory (user group computer) DSMod Modify items in active directory (user group computer) DSMove Move an Active directory Object DSRM Remove items from Active Directory e ECHO Display message on screen ENDLOCAL End localisation of environment changes in a batch file ERASE Delete one or more files EVENTCREATE Add a message to the Windows event log EXIT Quit the current script/routine and set an errorlevel EXPAND Uncompress CAB files EXTRACT Uncompress CAB files f FC FIND FINDSTR FOR /F FOR /F FOR FORFILES E. Baeck Compare two files Search for a text string in a file Search for strings in files Loop command: against a set of files Loop command: against the results of another command Loop command: all options Files, Directory, List Batch process multiple files Page 193 FORMAT FREEDISK FSUTIL FTP FTYPE Format a disk Check free disk space File and Volume utilities File Transfer Protocol File extension file type associations GETMAC GOTO GPRESULT GPUPDATE Display the Media Access Control (MAC) address Direct a batch program to jump to a labeled line Display Resultant Set of Policy information Update Group Policy settings g h HELP Online Help HOSTNAME Display the host name of the computer i iCACLS IF IFMEMBER IPCONFIG INUSE Change file and folder permissions Conditionally perform a command Is the current user a member of a group Configure IP Replace files that are in use by the OS LABEL LOGEVENT LOGMAN LOGOFF LOGTIME Edit a disk label Write text to the event viewer Manage Performance Monitor Log a user off Log the date and time in a file MAKECAB MAPISEND MBSAcli MEM MD MKLINK MODE MORE MOUNTVOL MOVE MOVEUSER MSG MSIEXEC MSINFO32 MSTSC Create .CAB files Send email from the command line Baseline Security Analyzer Display memory usage Create new folders Create a symbolic link (linked) Configure a system device Display output, one screen at a time Manage a volume mount point Move files from one folder to another Move a user from one domain to another Send a message Microsoft Windows Installer System Information Terminal Server Connection (Remote Desktop Protocol) NET NETDOM NETSH NETSVC NBTSTAT NETSTAT NOW NSLOOKUP NTBACKUP NTDSUtil NTRIGHTS Manage network resources Domain Manager Configure Network Interfaces, Windows Firewall & Remote access Command-line Service Controller Display networking statistics (NetBIOS over TCP/IP) Display networking statistics (TCP/IP) Display the current Date and Time Name server lookup Backup folders to tape Active Directory Domain Services management Edit user account rights l m n o OPENFILES Query or display open files p 9.6.2015 Page 194 Computer Languages for Engineering - SS 15 PATH Display or set a search path for executable files PATHPING Trace route plus network latency and packet loss PAUSE Suspend processing of a batch file and display a message PERMS Show permissions for a user PERFMON Performance Monitor PING Test a network connection POPD Return to a previous directory saved by PUSHD PORTQRY Display the status of ports and services POWERCFG Configure power settings PRINT Print a text file PRINTBRM Print queue Backup/Recovery PRNCNFG Display, configure or rename a printer PRNMNGR Add, delete, list printers set the default printer PROMPT Change the command prompt PsExec Execute process remotely PsFile Show files opened remotely PsGetSid Display the SID of a computer or a user PsInfo List information about a system PsKill Kill processes by name or process ID PsList List detailed information about processes PsLoggedOn Who’s logged on (locally or via resource sharing) PsLogList Event log records PsPasswd Change account password PsPing Measure network performance PsService View and control services PsShutdown Shutdown or reboot a computer PsSuspend Suspend processes PUSHD Save and then change the current directory q QGREP Query Query Query Query Search file(s) for lines that match a given pattern Process / QPROCESS Display processes Session / QWinsta Display all sessions (TS/Remote Desktop) TermServer /QAppSrv List all servers (TS/Remote Desktop) User / QUSER Display user sessions (TS/Remote Desktop) r RASDIAL Manage RAS connections RASPHONE Manage RAS connections RECOVER Recover a damaged file from a defective disk REG Registry: Read, Set, Export, Delete keys and values REGEDIT Import or export registry settings REGSVR32 Register or unregister a DLL REGINI Change Registry Permissions REM Record comments (remarks) in a batch file REN Rename a file or files REPLACE Replace or update one file with another Reset Session Delete a Remote Desktop Session RD Delete folder(s) RMTSHARE Share a folder or a printer ROBOCOPY Robust File and Folder Copy ROUTE Manipulate network routing tables RUN Start | RUN commands RUNAS Execute a program under a different user account RUNDLL32 Run a DLL command (add/remove print connections) s SC E. Baeck Service Control Page 195 SCHTASKS Schedule a command to run at a specific time SCLIST Display Services SET Display, set, or remove session environment variables SETLOCAL Control the visibility of environment variables SETX Set environment variables SFC System File Checker SHARE List or edit a file share or print share ShellRunAs Run a command under a different user account SHIFT Shift the position of batch file parameters SHORTCUT Create a windows shortcut (.LNK file) SHOWGRPS List the groups a user has joined SHOWMBRS List the Users who are members of a group SHUTDOWN Shutdown the computer SLEEP Wait for x seconds SLMGR Software Licensing Management (Vista/2008) SOON Schedule a command to run in the near future SORT Sort input START Start a program, command or batch file SUBINACL Edit file and folder Permissions, Ownership and Domain SUBST Associate a path with a drive letter SYSTEMINFO List system configuration t TAKEOWN TASKLIST TASKKILL TELNET TIME TIMEOUT TITLE TLIST TOUCH TRACERT TREE TSDISCON TSKILL TSSHUTDN TYPE TypePerf Take ownership of a file List running applications and services End a running process Communicate with another host using the TELNET protocol Display or set the system time Delay processing of a batch file Set the window title for a CMD.EXE session Task list with full path Change file timestamps Trace route to a remote host Graphical display of folder structure Disconnect a Remote Desktop Session End a running process Remotely shut down or reboot a terminal server Display the contents of a text file Write performance data to a log file VER VERIFY VOL Display version information Verify that files have been saved Display a disk label WAITFOR WEVTUTIL WHERE WHOAMI WINDIFF WINMSDP WINRM WINRS WMIC WUAUCLT Wait for or send a signal Clear event logs, enable/disable/query logs Locate and display files in a directory tree Output the current UserName and domain Compare the contents of two files or sets of files Windows system report Windows Remote Management Windows Remote Shell WMI Commands Windows Update XCACLS XCOPY Change file and folder permissions Copy files and folders v w x 9.6.2015 Page 196 :: Computer Languages for Engineering - SS 15 Comment / Remark Within the following sections some of the most important commands are discussed which are very helpful, if you are working within the console window. The command’s details should be studied from the original literature too. Within this short scratch only the most important command options are discussed. A.1 Directory Commands A.1.1 Selecting a Drive Selecting the active drive, the name of the drive should be given as command, like c: to select the standard drive or d: to select the cd-rom or z: to select the user-drive on a computer of the computerpool. Informations of the drive are listet with the command vol. A.1.2 Listing the Content of a Directory To list the content of an directory we can use the command dir. dir lists all files and subdirectories of the actual directory. You can also call the dir command with some wildcards filtering. Figure A.1 shows a directory list using a wildcard2 filtering of *.pdf. Note, that the volume information of the actual drive is also listed. We will get the same listing, if we use the absolute path of the desired directory and start the command from anywhere. dir c:\CM\Cm-CLFE\BookOfExamples\*.pdf Figure A.1: Create a Directory List of all pdf Files A.1.3 Creating and Removing a Directory A directory can simply be created with the command mkdir. Figure A.2 shows in the first step a directory list of the directory c:\cm\commands. Then a new directory is created with the name commands. After that the creation is checked with a further dir call. You can remove the directory with the inverse command rmdir. 2 A wildcard is a joker or a filter definition for a command. * means everthing within one section of a file name in between dots. A ? character is a joker for only one character within a string. So the wildcards t?st.* would filter a file with the name test.pdf or tost.nothing . E. Baeck A.2. FILE COMMANDS Page 197 Figure A.2: Create a Directory A.1.4 Browsing through Directories To browse through directories you can use the command cd, which is also called change directory. With the command cd .. you step one level up to the root directory. You can specify the directories name relative then you will jump out from the actual directory to the specified. You also can specify the directories name absolute. Then you will jump from the roots directory into the specified. Figure A.3 shows, that we start in the root directory. Then we clime up with relative jumps into the directory c:\cm, c:\cm\commands and at last c:\cm\commands\test. After that we jump back with one jump into the root with cd \ and then back in our test directory with an absolute jump cd \cm\commands\test. Figure A.3: Browsing through the Directories A.2 File Commands How can we check the content of a file with a simple command. You can use the type command. Figure A.4 shows how to pipe3 a screen stream into a file. This can be done by using a pipe character >. Within a first step the directory list of the actual directory is created with the dir command. This list is piped into a file with the name dir.lst using the command dir > dir.lst. After having created the text file with the directory list this list is viewed by the command type. If you have to list larger files with a lot of lines, you can use the command more to give a list page by page. So you have to pipe the output of the type command into the postprocessing more command by using the command type longfile.txt | more. 3 Pipe means, that a output stream of one process is used as an input stream for a following process, | character, or is used as input stream for a file, > character. 9.6.2015 Page 198 Computer Languages for Engineering - SS 15 Figure A.4: Create a Directory To delete a file the command del can be used. Figure A.5 shows in the first step a directory list. Two files are found in the directory, dir.lst and helloworld.f90. In the second step the file dir.lst is deleted by the command del dir.lst. This is shown in the third step giving an actual directory list with the command dir. You can also use wildcards within the del command, so you can delete all files from the directory with the wildcard *.*, so in this case we use the command del *.*. Figure A.5: Create a Directory A.3 Environment Commands One of the most important environment commands is the command path which is used to specify the search path for the executables. If you want to execute a program from the command line, it should be accessible by the command executer. Therefore the command path should be used to extend the standard path by the access path of the desired program. If we want to use the compiler gfortran.exe, which lives in the folder c:\programs\mingw\bin we should use the following command to extend the search path. path = %path%;c:\programs\mingw\bin %path% sets the still active path, which should not be overwitten by the extension. If you want to check the actual path, then the command path can be given without parameters (see figure A.6). The figure shows that the MinGW\bin is set and that the installed compiler g95.exe was found. E. Baeck A.3. ENVIRONMENT COMMANDS Page 199 Figure A.6: Checking the System Path Please note, that no quotes are used to specify the search path of the MinGW package, even if there are space characters inside the path name. It is astonishing, that the compiler executable is found, if quotes are used, but the secondary processes can obviously not be executed, so that you will get the following error message (see figure A.7), if the compiler should compile a source file. Figure A.7: Compiler Error due to wrong Path Settings 9.6.2015 Page 200 E. Baeck Computer Languages for Engineering - SS 15 Appendix B Code::Blocks’s first Project Within this chapter the creation of a project with the Code::Blocks IDE for Fortran is discussed as well as the settings which are essential (see also section 5.1). A new project can be created following the steps discussed below. 1. Start the IDE To create a project we start the Code::Blocks application (see figure 1.10). 2. Check the Toolchain Because the Code::Blocks IDE is a general IDE for a set of compilers, we have to setup the parameters for the toolchain executable. Especial the compiler’s root directory has to be fitted (see figure 1.11). 3. Start new Project Click on the link Create a new project. 4. Select the Fortran Template Figure B.1 shows the selection of the template to initialize the new project. 5. Setup the Project’s Name and it’s Folder To setup the project’s name and folder, we have to fill in the project’s name into the first edit control (see figure B.2). The second edit contents the root folder of our project. We can use the browser control - the button with the three dots - to browse the folder tree. The third edit contents the project file. The file name is created automatically. The whole project file with the total path is created within the fourth edit. 6. Setup the Project’s Configuration Clicking on next we will see the form to specify the project’s configuration (see figure B.3). Note, that is very important to select the proper compiler in the first combo. Please select the GNU Fortran Compiler. If the proper compiler is not selected, the build chain will crash with a strange error message. It’s recommended to use the standard configurations release and debug. If the software will work properly, you can build the release version of your software for shipping. If the software is working faultily or bad the debug version can be used to find the software bugs with the debugger tool. 201 Page 202 Computer Languages for Engineering - SS 15 Figure B.1: Select the Fortran Template Figure B.2: Setup the Project’s Name and Folder 7. Project’s Wizard finished If we click on finish, the wizard is closed and we will see the project within the IDE’s project browser (see figure B.4). If we open the node Fortran Sources, we find a source module called hello.f90. 8. Rename the hello by a meaningful Name To use a meaningful name for our main module we have to rename the hello.f90 by simply clicking right on it. You select the rename item from the context menu and will get a small form to overwrite the hello.f90 (see figure B.5). So we overwrite the default name by LittleProjectMain.f90. E. Baeck Page 203 Figure B.3: Setup the Project’s Configuration Figure B.4: The Project now is created 9. Open the Main Module A module’s source can be loaded into the editor by double clicking it in the source folder tree. Figure B.6 shows the content of the default source of the renamed module. 9.6.2015 Page 204 Computer Languages for Engineering - SS 15 Figure B.5: Overwriting the default source filename Figure B.6: Open the Main Module E. Baeck Page 205 10. Overwrite the Module Source After having opend the main module, we can overwrite it’s source. In line 6 we implement a call to the subroutine MySubModule. This source lives in a second module, which should be created as follows. Figure B.7 shows the overwritten source of the main module. Figure B.7: Overwriting the Main Module’s Source 11. Create a new Module for the Subroutine To create a new source file, a new module we use the command file/New/empty file from the main menu. The new source file is initialized with a standard name, in our case with untitled1. So we are asked, whether we want to add this new source to our project. The answer should be yes (see figure B.8). 12. Save the new Module using a meaningful Name After having added the new module to our project, we should specify the name of the new module within the standard save as dialog (see figure B.9). 13. Setup the Configuration for the new Module Before we can start to write the new module’s source we have to set up it’s configuration. We select both configurations, the release and the debug configuration (see figure B.10). 14. Writing the new Module’s Source After having installed the new module for our project we can double click it’s node within the source module tree and load it into the IDE’s editor. Because we want to show how to work with more than one module our main program should call the subroutine MySubModule which only should print the content of it’s input string parameter. Figure B.11 shows the source of the sub module. 15. Build the Executable The executable now can be build by the command Build/Build from the main menu or by the acceleration key Ctrl-F9. If you have executed the build you should see a build log like in figure 9.6.2015 Page 206 Computer Languages for Engineering - SS 15 Figure B.8: Starting with a new Module Figure B.9: Save the new Module specifying it’s name B.12 in the build log window. Before a total build is executed it is recommended to check and compile each module for it’s own by the acceleration key Ctrl-Shift-F9. 16. Executing the Executable The executable can be started from the IDE with the command Build/Run or by usage of the acceleration key Ctrl-F10. The programm starts within a command window and will print it’s output (see figure B.13). E. Baeck Page 207 Figure B.10: Select the Configurations to support Figure B.11: Writing the Sub Module’s Source 9.6.2015 Page 208 Computer Languages for Engineering - SS 15 Figure B.12: Check the Output int the Build Log Window Figure B.13: Result of our little Project E. Baeck Appendix C Some Theory C.1 Section Properties Within this chapter the formulas for the section properties of a thin walled model are given. A thin walled model for a profile section consists of a set of lines which describes the profile section geometry at the centerline. C.1.1 The Area of a Profile Section The Area is approximately the sum of the areas of the lines of the thin walled model. Z A= eµ · dA ≈ A with: Li ti eµ,i C.1.2 n X eµ,i · Li · ti (C.1) i=1 the length of line i the thickness of line i the relative elasticity of line i (1 for only one material) First Moments of an Area The first moments of an area are the area integrals given below. The (y,z) values are related to an given coordinate system. Z Sy = eµ · z · dA ≈ A Z Sz = eµ · y · dA ≈ A with: Ai yi zi n X i=1 n X eµ,i · z i · Ai eµ,i · y i · Ai i=1 the area of a line i the y coordinate of the center of line i the z coordinate of the center of line i 209 (C.2) Page 210 C.1.3 Computer Languages for Engineering - SS 15 Second Moments of an Area or Moments of Inertia The moments of inertia can be calculated with the formulas below. If we have a given arbitrary coordinate system in general we have three values of inertia the Iy , the Iz and the mixed Iyz . If we use the main coordinate system, the mixed moment of inertia is vanishing, so we use the symbols Iξ and Iη . Z 2 eµ · z · dA ≈ Iy = A Z A (zb,i − za,i )2 /12) + z 2i · Ai n X eµ,i · (yb,i − ya,i )2 /12) + y 2i · Ai i=1 n X Z eµ · y · z · dA ≈ Iyz = A with: Ai yi zi ya,i za,i yb,i zb,i eµ,i · i=1 eµ · y 2 · dA ≈ Iz = n X eµ,i · (((yb,i − ya,i )(zb,i − za,i )/12) + y i · z i ) · Ai ) (C.3) i=1 the area of a line i the y coordinate of the center of line i the z coordinate of the center of line i the y coordinate of the first point of line i the z coordinate of the first point of line i the y coordinate of the second point of line i the z coordinate of the second point of line i R To solve an integral like Iy = A z 2 · dA for a polyline we can split up the integral into the sum of integrals over the polyline segments. Z Iy = A z 2 · dA = n Z X i=1 z 2 · dA (C.4) Ai To solve an integral for a polyline segment we simple calculate it for the center of mass, because a simple shift only will give us an additional term, the Steiner term. If we now want to calculate the polyline integral at the center of mass we rotate the coordinate system by an angle ϕ into the line’s longitudinal direction, because the transversal dimension, the thickness, is constant and so the respective integral will be trivial. E. Baeck C.1. SECTION PROPERTIES Page 211 Thus we make the following substitution. (y, z ) ⇒ (η, ξ) (C.5) z = ξ/cos(ϕ) (C.6) With this substitution we will get the following integral. Z η=+t Z ξ=+L/2 Iy,i = ξ2 · dη · dξ cos(ϕ)2 η=−t ξ=−L/2 Z ξ=+L/2 ξ2 =t· · dξ cos(ϕ)2 ξ=+L/2 1 ξ3 =t· · 3 cos(ϕ)2 ξ=−L/2 ξ=−L/2 =t· L3 1 · 12 cos(ϕ)2 (zb,i − za,i )2 = · Ai 12 C.1.4 with t · L = Ai (C.7) Center of Mass The coordinates of the center of mass are calculated with the arithmetic mean. Because the numerator of the arithmetic mean is identical with the first moment of the area (see section C.1.2) and the denominator is identical with the area of the profile, which is calculated in section C.1.1 we can use this values. R y · dA Sz = yc = AR A dA R A z · dA Sy = zc = AR (C.8) A A dA C.1.5 Moments of Inertia with Respect to the Center of Mass If we know the center of mass coordinates given in section C.1.4 we can calculate the moments of inertia with respect to the center of mass using Steiner’s Theorem as follows. Iy,c = Iy − zc2 · A Iz ,c = Iz − yc2 · A Iyz ,c = Iyz − yc · zc · A (C.9) 9.6.2015 Page 212 C.1.6 Computer Languages for Engineering - SS 15 Main Axis Transformation To get the moments of inertia Iξ and Iη we have to transform the moments of inertia into the main coordinate system. Using this coordinate system the mixed moment of inertia is vanishing. The main axis transformation is given with equation C.10.1 Idel = Iy,c − Iz ,c Isum = Iy,c + Iz ,c q 2 + 4 · I2 Isqr = IDel yz ,c 2 · Iyz ,c 1 · arctan( ) 2 Idel 1 Iξ = · (Isum + Isqr ) 2 1 Iη = · (Isum − Isqr ) 2 ϕ= 1 (C.10) The rotation angle ϕ should be shifted into the intervall [−π/2... + π/2]. To avoid a zero division calculating the rotation angle ϕ a special version of the atan function should be used, which is able to handle the pole problem. In Python like in C this function is called atan2(y, x ), which calculates the atan( xy ). E. Baeck Appendix D Conventions D.1 The Java Code Conventions The following code convention [8] is published by Oracle (successor of Sun Microsystems, Inc). We apply this convention to choose names for our software items. 1. Classes Class names should be nouns, in mixed case with the first letter of each internal word capitalized. Try to keep your class names simple and descriptive. Use whole words-avoid acronyms and abbreviations (unless the abbreviation is much more widely used than the long form, such as URL or HTML). 2. Methods Methods should be verbs, in mixed case with the first letter lowercase, with the first letter of each internal word capitalized. 3. Variables Except for variables, all instance, class, and class constants are in mixed case with a lowercase first letter. Internal words start with capital letters. Variable names should not start with underscore _ or dollar sign $ characters, even though both are allowed. Variable names should be short yet meaningful. The choice of a variable name should be mnemonic- that is, designed to indicate to the casual observer the intent of its use. One-character variable names should be avoided except for temporary ”throwaway” variables. Common names for temporary variables are i, j, k, m, and n for integers; c, d, and e for characters. 4. Constants The names of variables declared class constants and of ANSI constants should be all uppercase with words separated by underscores ("_"). (ANSI constants should be avoided, for ease of debugging.) 213 Page 214 E. Baeck Computer Languages for Engineering - SS 15 Bibliography [1] Photo: Lawrence Livermore National Laboratory (Link: http://www.columbia.edu/acis/history/704.html) [2] Watcom FORTRAN 77 Language Reference Edition 11.0c [3] Stefen J. Chapman Fortran 90/95 for Scientists and Engineers, Second Edition McGraw-Hill, 2004 [4] H.R. Schwarz, N. K¨ockler Numerische Mathematik BI Wissenschaftsverlag Mannheim/Wien/Z¨urich, 1988 [5] Wikipedia, The Free Encyclopedia [6] cplusplus.com - The C++ Resources Network http://www.cplusplus.com/ [7] ISO/IEC 19501:2005 Information technology – Open Distributed Processing – Unified Modeling Language (UML) Version 1.4.2 [8] Java Code Conventions Oracle Inc., Sun Microsystems, Inc., September 12, 1997 215 Index .AND., 31 .EG., 30 .EQU., 31 .FALSE., 31 .GE., 30 .GT., 30 .LE., 30 .LT., 30 .NE., 30 .NEQU., 31 .NOT., 31 .OR., 31 .TRUE., 31 2GB, 20 A format, 34 ACTION, 32 area, 209 arithmetic mean, 211 array, 50 automatic, 52 dynamical, 51 static, 50 assembler, 4 backward substitution, 105 base class, 157 bash, 63 binary numbers, 20 block data, 58 break, 144 build, 13, 205 bytecode, 4 C, 3, 5, 131 C++, 3, 131 C/C++ projects, 127 C#, 3 CALL, 85 card punch, 6 cast, 152 catch, 146, 172, 183 cd, 197 cdsqrt, 47 center of mass, 211 channel, 33 char, 132 CLOSE, 92 cmd, 191 Code::Blocks, 7, 12, 201 columns, 18 COMMON, 101 common, 58 compiler, 3 complement, 20 complex, 20 console window, 191 constants, 23 contains, 59 CONTINUE, 36 CYCLE, 36 CYGWIN, 10 debugger, 8 del, 198 derivative, 69, 71 dir, 196, 197 DO, 36 double, 132, 141 DOUBLE PRECISION, 22 E format, 34 editor, 7, 12 encapsulating, 108 END DO, 36 end function, 48 end subroutine, 49 endian, 21 endianness, 21 equivalence, 61 216 INDEX exception, 146, 172, 183 EXIT, 36 exponent, 41 extension, 12 F format, 34 factorial, 40 FB substitution, 117 Fibunacci numbers, 4 first moment, 209 fixed format, 18 float, 20, 132 formal parameters, 48 FORMAT, 34 Forth, 4 FORTRAN 2008, 5 FORTRAN I, 5 FORTRAN II, 5 FORTRAN IV, 5 Fortran66, 36, 37, 42, 65 Fortran77, 18, 19, 22, 36, 37 Fortran90/95, 18, 19, 22, 37 forward substitution, 105 free format, 18 free FORTRAN compiler, 12 free FORTRAN tools, 7 function, 48 G95, 7, 17, 198 GCC C++ compiler, 128 GFortran, 17, 198 Grace Hopper, 5 I format, 34 IBM Type704, 5 IDE, 7 IEEE 754, 25 IF, 42 if, 143 INCLUDE, 85 info.server, 17 input parameters, 48 integer, 20 INTEGER*2, 22 INTEGER*4, 22 interpreter, 3 IPv6, 21 Page 217 ISML, 5 iteration, 69 Java Code Conventions, 213 Job Control Language, 63 kind, 23, 47 label, 36 linear equation system, 105, 118 LINUX, 3, 10, 63 LOGICAL*1, 22 LOGICAL*2, 22 long, 132 long long, 132 machine code, 4 main axis, 212 main function, 136 memory manager, 106 memset, 152 MFC, 149 MinGW, 7, 10, 198 mkdir, 196 module, 58 moment of inertia, 210 MS Excel, 41 NAG, 5 negative numbers, 20 new, 163 Newton, 50 Newton’s Algorithm, 71 OOP, 148 OPEN, 92 output parameters, 48 parent class, 157 path, 198 permissions, 150 piping, 197 POSITION, 32 preprocessor, 131 private, 150 program structure, 17 protected, 150 public, 150 9.6.2015 Page 218 punch card and 80, 6 Python, 4 Computer Languages for Engineering - SS 15 Watcom, 7, 17 Windows32, 20 WRITE, 32 quadratic equation, 43 X format, 34 READ, 32, 85, 92 REAL*8, 22 return, 137 return code, 146 rmdir, 196 run, 206 second moment, 210 short, 132 Smalltalk, 4 stack, 155 STATUS, 32 struct, 152 SUBROUTINE, 85 subroutine, 48 suffix cpp, 149 h, 149 superclass, 157 switch, 144 the most famous code, 65 throw, 146, 172, 183 try, 146, 172, 183 type, 197 typedef, 152 UML, 148 aggregation, 149 class diagram, 148 composition, 149 inheritance diagram, 148 note diagram, 148 note diagram assignment, 148 unicode, 19 UNIX, 10 unsigned, 132 upper lower extractor, 85 VBA, 3 virtual machine, 3 void, 132 vol, 196 E. Baeck
© Copyright 2024