```
$>cat hello.f
program my
write (*, *) , ' hello'
stop
end program my
$> f95 hello.f
$> ./a.out
hello
$> gfortran -o myhello hello.f
$> ./myhello
hello
```

```
Col. 1 : blank | "c" for comments | "*" for comments
Col. 1-5 : statement label (optional)
Col. 6 : continuation of previous line (optional)
Col. 7-72 : statements
```

the exclamation mark (!) for comments is not a standard part of Fortran77, but is explicitly allowed in Fortran90. the exclamation mark may appear anywhere on a line

```
integer list of variables
real list of variables
double precision list of variables
complex list of variables
logical list of variables
character list of variables
```

the list of variables should consist of variable names separated by commas. you cannot have a variable with the same name as the program
Some constants appear many times in a program. It is then often desirable to define them only once, in the beginning of the program. This is what the *parameter*
statement is for. It makes programs more readable. For example, the circle area program should rather have been written like this:

```
real r, area, pi
parameter (pi = 3.14159)
```

The rules for the parameter statement are:
For constants that are larger than the largest real allowed, or that requires high precision, double precision should be used. The notation is the same as for real constants except the "E" is replaced by a "D":

```
2.0D-1
1D99
```

logical constants:

```
.TRUE.
.FALSE.
```

Note that the dots enclosing the letters are required
An array of characters, called a string, consist of an arbitrary sequence of characters enclosed in apostrophes (single quotes):

```
'ABC'
'Anything goes!'
'It is a nice day'
```

Strings and character constants are case sensitive. A problem arises if you want to have an apostrophe in the string itself. In this case, you should double the apostrophe:
```
'It''s a nice day'
```

The precedence of arithmetic operators in Fortran77 are (from highest to lowest):
```
** exponentiation
*,/ multiplication, division
+,- addition, subtraction
```

If you want to change the default evaluation order, you can use parentheses
When different data types occur in the same expression, type conversion has to take place, either explicitly or implicitly. Fortran will do some type conversion implicitly. For example,

```
real x
x = x + 1
```

will convert the integer '1' to the real number, and has the desired effect of incrementing x by one. However, in more complicated expressions, it is good programming practice to force the necessary type conversions explicitly. For numbers, the following functions are available:
```
int
real
dble
ichar ! takes a character and converts it to an integer
char ! takes a integer and converts it to a character
```

Logical expressions can only have the value .TRUE. or .FALSE.. A logical expression can be formed by comparing arithmetic expressions using the following relational operators:
```
.LT. <
.LE. <=
.GT. >
.GE. >=
.EQ. =
.NE. /=
```

logical operators
```
.AND.
.OR.
.NOT.
```

Example:
```
logical a, b
a = .TRUE.
b = a .AND. 3 .LT. 5/2
```

The order of precedence is important. The rule is that arithmetic expressions are evaluated first, then relational operators, and finally logical operators
The logical if statement: "if (logical expression) executable statement". This has to be written on one line. This example finds the absolute value of x:

```
if (x .LT. 0) x = -x
```

If more than one statement should be executed inside the if, then the following syntax should be used:
```
if (logical expression) then
statements
endif
```

The most general form of the if-statement has the following form:
```
if (logical expression) then
statements
elseif (logical expression) then
statements
:
:
else
statements
endif
```

if statements can be nested in several levels:
```
if (x .GT. 0) then
if (x .GE. y) then
write (*, *) 'x is positive and x >= y'
else
write (*, *) 'x is positive but x
endif
```

The do-loop is used for simple counting. Here is a simple example that prints the cumulative sums of the integers from 1 through n (assume n has been assigned a value elsewhere):
```
integer i, n, sum
sum = 0
do 10 i = 1, n
sum = sum + i
write (*, *) 'i =', i
write (*, *) 'sum =', sum
10 continue
```

The number 10 is a statement label. The programmer is responsible for assigning a unique number to each label in each program (or subprogram). Recall that column positions 1-5 are reserved for statement labels. The numerical value of statement labels have no significance, so any integers can be used, in any order. Typically, most programmers use consecutive multiples of 10
The variable defined in the do-statement is incremented by 1 by default. However, you can define the step to be any number but zero. This program segment prints the even numbers between 1 and 10 in decreasing order:

```
integer i
do 20 i = 10, 1, -2
write (*, *) 'i =', i
20 continue
```

The general form of the do-loop is as follows:
```
do label var = expr1, expr2, expr3
statements
label continue
```

The general form of while-loop is as follows:
```
label if (logical expr) then
statements
goto label
endif
```

The general form of until-loop is as follows:
```
label continue
statements
if (logical expr) goto label
```

The simplest array is the one-dimensional array, which is just a sequence of elements stored consecutively in memory. For example, the declaration

` real a(20) `

declares a as a real array of length 20. That is, a consists of 20 real numbers stored contiguously in memory. By convention, Fortran arrays are indexed from 1 and up. Thus the first number in the array is denoted by a(1) and the last by a(20)
The type of an array element can be any of the basic data types:

```
integer i(10)
logical b(2)
double precision x(100)
```

Each element of an array can be thought of as a separate variable. You reference the i'th element of array a by a(i). Here is a code segment that stores the 10 first square numbers in the array sq:
```
integer i, sq(20)
do 100 i = 1, 10
sq(i) = i**2
100 continue
```

A COMMON BUG IN FORTRAN IS THAT THE PROGRAM TRIES TO ACCESS ARRAY ELEMENTS THAT ARE OUT OF BOUNDS OR UNDEFINEDTHIS IS THE RESPONSIBILITY OF THE PROGRAMMER, AND THE FORTRAN COMPILER WILL NOT DETECT ANY SUCH BUGS!

Matrices are very important in linear algebra. Matrices are usually represented by two-dimensional arrays. For example, the declaration

```
real A (3,5)
```

defines a two-dimensional array of 3*5=15 real numbers. It is useful to think of the first index as the row index, and the second as the column index
It is quite common in Fortran to declare arrays that are larger than the matrix we want to store. (This is because Fortran does not have dynamic storage allocation.) This is perfectly legal. Example:

```
real A (3,5)
integer i, j
c
c We will only use the upper 3 by 3 part of this array
c
do 20 j = 1, 3
do 10 i = 1, 3
a(i,j) = real(i)/real(j)
10 continue
20 continue
```

The elements in the submatrix A (1:3,4:5) are undefined. Do not assume these elements are initialized to zero by the compiler (some compilers will do this, but not all)
Fortran stores higher dimensional arrays as a contiguous sequence of elements.

It is important to know that 2-dimensional arrays are stored by columnEach array will have some memory address assigned to the beginning of the array, that is element (1,1). The address of element (i,j) is then given by

```
addr[A (i,j)] = addr[A (1,1)] + (j-1)*lda + (i-1)
```

where

Note thatldais in general different from the actual matrix dimension. Many Fortran errors are caused by this, so it is very important you understand the distinction!

There is an alternate way to declare arrays in Fortran77. The statements

```
real A, x
dimension x(50)
dimension A (10,20)
```

are equivalent to
```
real A (10,20), x(50)
```

The dimension statement is considered old-fashioned style

Fortran subprogram calls are based on "call by reference"

the calling parameters are not copied to the called subprogram, but the addresses of the parameters (variables) are passed and subroutine changes value of this variable

program mytestref real x x = 10.2 write (*, 100) x call byref (x) write (*, 100) x 100 format ("fm main :", F8.3) end program mytestref subroutine byref (x) write (*, 200) x x = x + 11.7 write (*, 200) x return 200 format ("fm sub :", F8.3) end subroutine byref

*Fortran77 has no global variables*. The only way to pass information between subroutines is to use the subroutine parameter list

EXP | returns the exponential of a real number |

LOG | returns the natural logarithm of a real number |

LOG10 | returns the base 10 logarithm of a real number |

SQRT | returns the square root of a real number |

MAX | returns the maximum value in a list |

MIN | returns the minimum value in a list |

MOD | returns the remainder after division |

DPROD | returns the product of real values X and Y, using double precision |

ABS | returns the absolute value of a numeric quantity |

AINT | returns the value of a real number rounded towards zero |

NIN T | returns the nearest integer to a given real number |

ANINT | returns (as a real value) the nearest integer to a given real number |

AIMAG | returns the imaginary part of a complex number |

CMPLX | returns the complex number formed by real and imaginary parts |

CONJG | returns the conjugate of a complex number |

REAL | returns the real part of a complex number |

DBLE | converts a complex, integer or real value to double precision real |

INT | converts a complex, double precision or real value to integer |

SIN | returns the sine of a real or complex number |

COS | returns the cosine of a real or complex number |

TAN | returns the tangent of a real number |

ASIN | returns the inverse sine of a value between -1 and 1 |

ACOS | returns the inverse cosine of a value between -1 and 1 |

ATAN | returns the inverse tangent of a value |

ATAN2 | returns the inverse tangent of a value |

COSH | returns the hyperbolic cosine of a real number |

SINH | returns the hyperbolic sine of a real number |

TANH | returns the hyperbolic tangent of a real number |

DIM | returns the maximum of X-Y or 0 |

NOT | returns the bitwise NOT of an integer |

BTEST(I,P) | is TRUE if bit P of I is 1 |

IAND | returns the bitwise AND of two integers |

IOR | returns the bitwise inclusive OR of two integers |

ISHFT | shifts the bits in an integer word |

ISHFTC | circular-shifts the bits in an integer word |

IBSET | sets a given bit to one in an integer word |

IEOR | returns the bitwise exclusive OR of two integers |

IBCLR | sets a given bit to zero in an integer word |

MVBITS | extracts bits from one place and copies them elsewhere |

IBITS | extracts LEN bits from word I start at position POS |

IBSET | sets a given bit to one in an integer word |

LGE | returns "S1 >= S2" where S1 and S2 are strings |

LGT | returns "S1 > S2" where S1 and S2 are strings |

LLE | returns "S1 <= S2" where S1 and S2 are strings |

LLT | returns "S1 < S2" where S1 and S2 are strings |

LEN_TRIM | returns the "used" length of a string variable |

LEN | returns declared length of a string variable (or string constant) |

INDEX(S1,S2) | returns the location of the first occurrence of substring S2 in string S1 |

ICHAR | returns the character index of a given character |

CHAR | returns the character of given character index |

"read" is used for input, while "write" is used for output. A simple form is

```
read (unit no, format no) list-of-variables
write (unit no, format no) list-of-variables
```

The unit number can refer to either standard input, standard output, or a file. The format number refers to a label for a format statement
Examples

```
integer m, n
real x, y, z(10)
read (*, *) m, n
read (*, *) x, y
read (*, *) z
```

A data file consists of records according to traditional Fortran terminology. In our
example, each record contains a number (either integer or real). Records are separated by
either blanks or commas. Hence a legal input to the program above would be:
```
-1 100
-1.0 1e+2
1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
```

Or, we could add commas as separators:
```
-1, 100
-1.0, 1e+2
1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0
```

Note that Fortran77 input is line sensitive, so it is important not to have extra input elements (fields) on a line (record).

For simple list-directed I/O it is possible to use the alternate syntax

```
read *, list-of-variables
print *, list-of-variables
```

which has the same meaning as the list-directed read and write statements described earlier. This version always reads/writes to standard input/output so the * corresponds to the format
In Fortran each file is associated with a unit number, an integer between 1 and 99. Some unit numbers are reserved: 5 is standard input, 6 is standard output

Before you can use a file you have to open it. The command is

```
open (list-of-specifiers)
```

where the most common specifiers are:
```
[UNIT=] u
IOSTAT= ios
ERR= err
FILE= fname
STATUS= sta
ACCESS= acc
FORM= frm
RECL= rl
```

- u is a number in the range 1-99 that denotes this file (the programmer may chose any number but he/she has to make sure it is unique).
- ios is the I/O status identifier and should be an integer variable. Upon return, ios is zero if the statement was successful and returns a non-zero value otherwise.
- err is a label which the program will jump to if there is an error.
- fname is a character string denoting the file name.
- sta is a character string that has to be either NEW, OLD or SCRATCH. It shows the prior status of the file. A scratch file is a file that is created when opened and deleted when closed (or the program ends).
- acc must be either SEQUENTIAL or DIRECT. The default is SEQUENTIAL.
- frm must be either FORMATTED or UNFORMATTED. The default is UNFORMATTED.
- rl specifies the length of each record in a direct-access file.

```
read ([UNIT=]u, [FMT=]fmt, IOSTAT=ios, ERR=err, END=s)
write ([UNIT=]u, [FMT=]fmt, IOSTAT=ios, ERR=err, END=s)
```

where most of the specifiers have been described above. The END=s specifier defines which statement label the program jumps to if it reaches end-of-file

```
write (*, label) list-of-variables
label format format-code
```

A simple example demonstrates how this works. Say you have an integer variable you want to print in a field 4 characters wide and a real number you want to print in fixed point notation with 3 decimal places
```
write (*, 900) i, x
900 format (I4, F8.3)
```

The code I4 stands for an integer with width four, while F8.3 means that the number should be printed using fixed point notation with field width 8 and 3 decimal places
The most common format code letters are:

- A - text string
- D - double precision numbers, exponent notation
- E - real numbers, exponent notation
- F - real numbers, fixed point format
- I - integer
- X - horizontal skip (space)
- / - vertical skip (newline)

For horizontal spacing, the nX code is often used. This means n horizontal spaces. If n is omitted, n=1 is assumed

Instead of specifying the format code in a separate format statement, one can give the format code in the read/write statement directly:

```
write (*, '(A, F8.3)') 'The answer is x = ', x
```

Sometimes text strings are given in the format statements, e.g. the following version is also equivalent:

```
write (*, 999) x
999 format ('The answer is x = ', F8.3)
```

Often a format statement involves repetition, for example

```
950 format (2X, I3, 2X, I3, 2X, I3, 2X, I3)
```

There is a shorthand notation for this:
```
950 format (4(2X, I3))
```

program test_io_01 integer x, y open (unit = 10, file = "aaa.dat") open (unit = 11, file = "bbb.dat") read (10, *) x, y write (11, 100) (x + y) 100 format (I4) end program test_io_01

```
c write integer matrix A to file
c m - rows (lines in the file) n - cols
subroutine out_Matrix (filename, m, n, X)
implicit none
character(len=*) filename
integer i, j, m, n, fd
real X (100, 100)
parameter (fd = 21)
open (unit = fd, file = filename)
do j = 1, m ! rows
write (fd, *) ( X (j, i) , i = 1, n ) ! cols
enddo
close (fd)
end subroutine out_Matrix
c read integer matrix A from file
c m - rows (lines in the file) n - cols
subroutine in_Matrix (filename, m, n, X)
implicit none
character(len=*) filename
integer i, j, m, n, fd
real X (100, 100)
parameter (fd = 20)
open (unit = fd, file = filename, status = 'old')
do j = 1, m ! rows
read (fd, *) ( X (i, j) , i = 1, n ) ! cols
enddo
close (fd)
end subroutine in_Matrix
```

Libraries have file names starting with lib and ending in .a. Some libraries have already been installed by your system administrator, usually in the directories /usr/lib and /usr/local/lib

For example, the BLAS library MAY be stored in the file /usr/local/lib/libblas.a. You use the -l option to link it together with your main program, e.g

```
f77 main.f -lblas
```

You can link several files with several libraries at the same time if you wish:
```
f77 main.f mysub.f -llapack -lblas
```

The order you list the libraries is significant. In the example above -llapack should be listed before -lblas since LAPACK calls BLAS routines
If you want to create your own library, you can do so by compiling the source code to object code and then collecting all the object files into one library file. This example generates a library called my_lib:

```
f77 -c *.f
ar rcv libmy_lib.a *.o
ranlib libmy_lib.a
rm *.o
```

Level 1: y ← α * x + y

Level 2: y ← α * A * x + β * y

Level 3: B ← α * A1 * A2 + β * B

## LEVEL 1

## copy

call scopy (n, x, incx, y, incy) call dcopy (n, x, incx, y, incy) call ccopy (n, x, incx, y, incy) call zcopy (n, x, incx, y, incy) copies vector to another vector routines perform a vector-vector operation defined as y = x, where x and y are vectors n INTEGER. Specifies the number of elements in vectors x and y x y REAL for scopy DOUBLE PRECISION for dcopy COMPLEX for ccopy DOUBLE COMPLEX for zcopy Array, size at least (1 + (n-1)*abs(inc.)) inc. INTEGER. Specifies the increment for the elements of x|y Output : y contains a copy of the vector x if n is positive## scal

call sscal (n, a, x, incx) call dscal (n, a, x, incx) call cscal (n, a, x, incx) call zscal (n, a, x, incx) call csscal (n, a, x, incx) call zdscal (n, a, x, incx) Computes the product of a vector by a scalar The routines perform a vector operation defined as x = a*x where: a is a scalar, x is an n-element vector n INTEGER. Specifies the number of elements in vector x a REAL for sscal and csscal DOUBLE PRECISION for dscal and zdscal COMPLEX for cscal DOUBLE COMPLEX for zscal Specifies the scalar a x REAL for sscal DOUBLE PRECISION for dscal COMPLEX for cscal and csscal DOUBLE COMPLEX for zscal and zdscal Array, size at least (1 + (n -1)*abs(incx)) incx INTEGER. Specifies the increment for the elements of x Output: x - updated vector x## nrm2

res = snrm2 (n, x, incx) res = dnrm2 (n, x, incx) res = scnrm2 (n, x, incx) res = dznrm2 (n, x, incx) computes the Euclidean norm of a vector. The routines perform a vector reduction operation defined as res = ||x||, where: x is a vector, res is a value containing the Euclidean norm of the elements of x n INTEGER. Specifies the number of elements in vector x x REAL for snrm2 DOUBLE PRECISION for dnrm2 COMPLEX for scnrm2 DOUBLE COMPLEX for dznrm2 Array, size at least (1 + (n -1)*abs (incx)) incx INTEGER. Specifies the increment for the elements of x res REAL for snrm2 DOUBLE PRECISION for dnrm2 REAL for scnrm2 DOUBLE PRECISION for dznrm2 Contains the Euclidean norm of the vector x## dot

res = sdot (n, x, incx, y, incy) res = ddot (n, x, incx, y, incy) res = dot (x, y) computes a vector-vector dot product n INTEGER. Specifies the number of elements in vectors x and y x,y REAL for sdot DOUBLE PRECISION for ddot Array, size at least (1+(n-1)*abs(incx)) inc. INTEGER. Specifies the increment for the elements of x|y## LEVEL 2

## gemv

call sgemv (trans, m, n, alpha, a, lda, x, incx, beta, y, incy) call dgemv (trans, m, n, alpha, a, lda, x, incx, beta, y, incy) call cgemv (trans, m, n, alpha, a, lda, x, incx, beta, y, incy) call zgemv (trans, m, n, alpha, a, lda, x, incx, beta, y, incy) call scgemv (trans, m, n, alpha, a, lda, x, incx, beta, y, incy) call dzgemv (trans, m, n, alpha, a, lda, x, incx, beta, y, incy) Computes a matrix-vector product using a general matrix The gemv routines perform a matrix-vector operation defined as y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, or y := alpha*conjg(A')*x + beta*y, where: alpha and beta are scalars, x and y are vectors, A is an m-by-n matrix trans CHARACTER*1. Specifies the operation: if trans= 'N' or 'n', then y := alpha*A*x + beta*y; if trans= 'T' or 't', then y := alpha*A'*x + beta*y; if trans= 'C' or 'c', then y := alpha *conjg(A')*x + beta*y m INTEGER. Specifies the number of rows of the matrix A. The value of m must be at least zero n INTEGER. Specifies the number of columns of the matrix A. The value of n must be at least zero alpha REAL for sgemv DOUBLE PRECISION for dgemv COMPLEX for cgemv, scgemv DOUBLE COMPLEX for zgemv, dzgemv Specifies the scalar alpha a REAL for sgemv, scgemv DOUBLE PRECISION for dgemv, dzgemv COMPLEX for cgemv DOUBLE COMPLEX for zgemv Array, size (lda, n) Before entry, the leading m-by-n part of the array a must contain the matrix of coefficients lda INTEGER. Specifies the leading dimension of a as declared in the calling (sub)program The value of lda must be at least max(1, m) x REAL for sgemv DOUBLE PRECISION for dgemv COMPLEX for cgemv, scgemv DOUBLE COMPLEX for zgemv, dzgemv Array, size at least (1+(n-1)*abs(incx)) when trans= 'N' or 'n' and at least (1+(m - 1)*abs(incx)) otherwise Before entry, the incremented array x must contain the vector x incx INTEGER. Specifies the increment for the elements of x The value of incx must not be zero beta REAL for sgemv DOUBLE PRECISION for dgemv COMPLEX for cgemv, scgemv DOUBLE COMPLEX for zgemv, dzgemv Specifies the scalar beta. When beta is set to zero, then y need not be set on input y REAL for sgemv DOUBLE PRECISION for dgemv COMPLEX for cgemv, scgemv DOUBLE COMPLEX for zgemv, dzgemv Array, size at least (1 +(m - 1)*abs(incy)) when trans= 'N' or 'n' and at least (1 +(n - 1)*abs (incy)) otherwise Before entry with non-zero beta, the incremented array y must contain the vector y incy INTEGER. Specifies the increment for the elements of y The value of incy must not be zero Output: y - updated vector y## LEVEL 3

## gemm

call sgemm (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc) call dgemm (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc) call cgemm (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc) call zgemm (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc) call scgemm (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc) call dzgemm (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc) Computes a matrix-matrix product with general matrices gemm routines compute a scalar-matrix-matrix product and add the result to a scalar-matrix product, with general matrices. The operation is defined as C := alpha*op(A)*op(B) + beta*C, where: op(X) is one of op(X) = X, or op(X) = X^T, or op(X) = X^H, alpha and beta are scalars, A, B and C are matrices: op(A) is an m-by-k matrix, op(B) is a k-by-n matrix, C is an m-by-n matrix transa CHARACTER*1. Specifies the form of op(A) used in the matrix multiplication: if transa = 'N' or 'n', then op(A) = A; if transa = 'T' or 't', then op(A) = A^T; if transa = 'C' or 'c', then op(A) = A^H transb CHARACTER*1. Specifies the form of op(B) used in the matrix multiplication: if transb = 'N' or 'n', then op(B) = B; if transb = 'T' or 't', then op(B) = B^T; if transb = 'C' or 'c', then op(B) = B^H m INTEGER. Specifies the number of rows of the matrix op(A) and of the matrix C The value of m must be at least zero n INTEGER. Specifies the number of columns of the matrix op(B) and the number of columns of the matrix C. The value of n must be at least zero k INTEGER. Specifies the number of columns of the matrix op(A) and the number of rows of the matrix op(B). The value of k must be at least zero alpha REAL for sgemm DOUBLE PRECISION for dgemm COMPLEX for cgemm, scgemm DOUBLE COMPLEX for zgemm, dzgemm Specifies the scalar alpha a REAL for sgemm, scgemm DOUBLE PRECISION for dgemm, dzgemm COMPLEX for cgemm DOUBLE COMPLEX for zgemm Array, size lda by ka, where ka is k when transa = 'N' or 'n', and is m otherwise Before entry with transa = 'N' or 'n', the leading m-by-k part of the array a must contain the matrix A, otherwise the leading k-by-m part of the array a must contain the matrix A lda INTEGER. Specifies the leading dimension of a as declared in the calling (sub)program When transa = 'N' or 'n', then lda must be at least max(1, m), otherwise lda must be at least max (1, k) b REAL for sgemm DOUBLE PRECISION for dgemm COMPLEX for cgemm, scgemm DOUBLE COMPLEX for zgemm, dzgemm Array, size ldb by kb, where kb is n when transa = 'N' or 'n', and is k otherwise Before entry with transa = 'N' or 'n', the leading k-by-n part of the array b must contain the matrix B, otherwise the leading n-by-k part of the array b must contain the matrix B ldb INTEGER. Specifies the leading dimension of b as declared in the calling (sub)program When transb = 'N' or 'n', then ldb must be at least max(1, k), otherwise ldb must be at least max (1, n) beta REAL for sgemm DOUBLE PRECISION for dgemm COMPLEX for cgemm, scgemm DOUBLE COMPLEX for zgemm, dzgemm Specifies the scalar beta When beta is equal to zero, then c need not be set on input c REAL for sgemm DOUBLE PRECISION for dgemm COMPLEX for cgemm, scgemm DOUBLE COMPLEX for zgemm, dzgemm Array, size ldc by n. Before entry, the leading m-by-n part of the array c must contain the matrix C, except when beta is equal to zero, in which case c need not be set on entry ldc INTEGER. Specifies the leading dimension of c as declared in the calling (sub)program The value of ldc must be at least max(1, m) Output: c - overwritten by the m-by-n matrix ( alpha*op(A)*op(B) + beta*C )

program myblas01
c inner product of two vectors, dimension=3
c copy of vector, dimension=3
c scalar multiplication of two vectors, dimension=3
c norm L2 (euclidean) of vector, dimension=3
real x(3), y(3), z(3)
read *, x ; read *, y
print *, sdot (3, x, 1, y, 1) ! inner product
call scopy (3, x, 1, z, 1) ! copy
print *, z
print *, snrm2 (3, z, 1) ! norm2
call sscal (3, 12.0, z, 1) ! scalar mult
print *, z
end program myblas01
! gfortran bla01.f -L/home/user/.local/lib -lblas
=========================================================================
program myblas02
c matrix (2x2) vector (2) multiplication with scalar factor plus
c scalar mult of the second vector (2)
200 format (2('F8,3'))
real K(4), x(2), y(2), alpha, beta
print *, 'koef a' ; read *, alpha
print *, 'koef b' ; read *, beta
print *, 'rows of your matrix'
read *, K(1), K(3) ! matrix stored by cols
read *, K(2), K(4) ! the first col, the second col
print *, 'vector x' ; read *, x(1), x(2)
print *, 'vector y' ; read *, y(1) , y(2)
call sgemv ('n', 2, 2, alpha, K, 2, x, 1, beta, y, 1) ! update y with result
print *, 'result is : '
write (*,'(F8.3)') y(1)
write (*,'(F8.3)') y(2)
end program myblas02
! gfortran bla02.f -L/home/user/.local/lib -lblas
==================================================================
program myblas03
c matrix (2x2) matrix (2x2) multiplication with scalar factor plus
c scalar mult of the third matrix (2x2)
200 format (2('F8,3'))
real A (4), B(4), C(4), alpha, beta
print *, 'koef a' ; read *, alpha
print *, 'koef b' ; read *, beta
print *, 'rows of your A matrix'
read *, A (1), A (3)
read *, A (2), A (4)
print *, 'rows of your B matrix'
read *, B(1), B(3)
read *, B(2), B(4)
print *, 'rows of your C matrix'
read *, C(1), C(3)
read *, C(2), C(4)
call sgemm ('n', 'n', 2, 2, 2, alpha, A, 2, B, 2, beta, C, 2) ! update C with result
print *, 'result is : '
write (*, *) C(1) , C(3)
write (*, *) C(2) , C(4)
end program myblas03
! gfortran bla03.f -L/home/user/.local/lib -lblas

## a general system of linear equations

A * X = B sgesv real dgesv real double precision cgesv complex zgesv complex double precision SGESV (N, NRHS, A, LDA, IPIV, B, LDB, INFO) DGESV (N, NRHS, A, LDA, IPIV, B, LDB, INFO) CGESV (N, NRHS, A, LDA, IPIV, B, LDB, INFO) ZGESV (N, NRHS, A, LDA, IPIV, B, LDB, INFO) N INTEGER The number of linear equations, i.e., the order of the matrix A. N >= 0. NRHS INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0. A (input/output) array, dimension (LDA,N) On entry, the N-by-N coefficient matrix A. On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored. LDA INTEGER The leading dimension of the array A. LDA >= max(1,N). IPIV (output) INTEGER array, dimension (N) The pivot indices that define the permutation matrix P; row i of the matrix was interchanged with row IPIV(i). B (input/output) array, dimension (LDB,NRHS) On entry, the N-by-NRHS matrix of right hand side matrix B. On exit, if INFO = 0, the N-by-NRHS solution matrix X. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,N). INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed.## "least squares solution" (overdetermined system of linear equations)

A * X = B where A is a general rectangular matrix of full rank sgels real dgels real double precision cgels complex zgels complex double precision SGELS (TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO) DGELS (TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO) CGELS (TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO) ZGELS (TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO) TRANS CHARACTER*1 = 'N': the linear system involves A; = 'T': the linear system involves A**T. M INTEGER The number of rows of the matrix A. M >= 0. N INTEGER The number of columns of the matrix A. N >= 0. NRHS INTEGER The number of right hand sides, i.e., the number of columns of the matrices B and X. NRHS >=0. A (input/output) array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, if M >= N, A is overwritten by details of its QR factorization as returned by SGEQRF; if M < N, A is overwritten by details of its LQ factorization as returned by SGELQF. LDA INTEGER The leading dimension of the array A. LDA >= max(1,M). B (input/output) array, dimension (LDB,NRHS) On entry, the matrix B of right hand side vectors, stored columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS if TRANS = 'T'. On exit, if INFO = 0, B is overwritten by the solution vectors, stored columnwise: if TRANS = 'N' and m >= n, rows 1 to n of B contain the least squares solution vectors; the residual sum of squares for the solution in each column is given by the sum of squares of elements N+1 to M in that column; if TRANS = 'N' and m < n, rows 1 to N of B contain the minimum norm solution vectors; if TRANS = 'T' and m >= n, rows 1 to M of B contain the minimum norm solution vectors; if TRANS = 'T' and m < n, rows 1 to M of B contain the least squares solution vectors; the residual sum of squares for the solution in each column is given by the sum of squares of elements M+1 to N in that column. LDB INTEGER The leading dimension of the array B. LDB >= MAX(1,M,N). WORK (workspace/output) array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK INTEGER The dimension of the array WORK. LWORK >= max( 1, MN + max( MN, NRHS ) ). For optimal performance, LWORK >= max( 1, MN + max( MN, NRHS )*NB ). where MN = min(M,N) and NB is the optimum block size. INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, the i-th diagonal element of the triangular factor of A is zero, so that A does not have full rank; the least squares solution could not be computed.

"the minimum norm solution" (underdetermined system of linear equations)## eigenvalues and left and right eigenvectors of a general matrix

sgeev real dgeev real double precision cgeev complex zgeev complex double precision ZGEEV (JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO) CGEEV (JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO) DGEEV (JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK, INFO) SGEEV (JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK, INFO) JOBVL = 'N': left eigenvectors of A are not computed; = 'V': left eigenvectors of A are computed. JOBVR = 'N': right eigenvectors of A are not computed; = 'V': right eigenvectors of A are computed. N INTEGER The order of the matrix A. N >= 0. A (input/output) array, dimension (LDA,N) On entry, the N-by-N matrix A. On exit, A has been overwritten. LDA INTEGER The leading dimension of the array A. LDA >= max(1,N). WR (output) array, dimension (N) WI (output) array, dimension (N) WR and WI contain the real and imaginary parts, respectively, of the computed eigenvalues. Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first. VL (output) array, dimension (LDVL,N) If JOBVL = 'V', the left eigenvectors u(j) are stored one after another in the columns of VL, in the same order as their eigenvalues. If JOBVL = 'N', VL is not referenced. If the j-th eigenvalue is real, then u(j) = VL(:,j), the j-th column of VL. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and u(j+1) = VL(:,j) - i*VL(:,j+1). LDVL INTEGER The leading dimension of the array VL. LDVL >= 1; if JOBVL = 'V', LDVL >= N. VR (output) array, dimension (LDVR,N) If JOBVR = 'V', the right eigenvectors v(j) are stored one after another in the columns of VR, in the same order as their eigenvalues. If JOBVR = 'N', VR is not referenced. If the j-th eigenvalue is real, then v(j) = VR(:,j), the j-th column of VR. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and v(j+1) = VR(:,j) - i*VR(:,j+1). LDVR INTEGER The leading dimension of the array VR. LDVR >= 1; if JOBVR = 'V', LDVR >= N. WORK (workspace/output) array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK. LWORK INTEGER The dimension of the array WORK. LWORK >= max(1,3*N), and if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good performance, LWORK must generally be larger. RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = i, the QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed; elements i+1:N of WR and WI contain eigenvalues which have converged.## examples of LAPACK usage

program LinearEquations c solving the matrix equation A * x = b using LAPACK Implicit none c declarations Real*4 A (3, 3), b(3) integer i, j, p(3), ok c define matrix A A (1, 1) = 3.1 ; A (1, 2) = 1.3 ; A (1, 3) = -5.7 A (2, 1) = 1.0 ; A (2, 2) = -6.9 ; A (2, 3) = 5.8 A (3, 1) = 3.4 ; A (3, 2) = 7.2 ; A (3, 3) = -8.8 c define vector b b(1) = -1.3 ; b(2) = -0.1 ; b(3) = 1.8 c find the solution using the LAPACK routine SGESV call sgesv (3, 1, A, 3, p, b, 3, ok) c c parameters c order of matrix A, c number of right hand sides (b), c matrix A, c leading dimension of A, c array p that records pivoting, c vector b on entry, result vector on exit, c leading dimension of b c return value of the function c c print the result vector if (ok .eq. 0) then do i = 1, 3 write (*, *) b(i) end do else write (*, *) "error" endif end program LinearEquations ! gfortran lapa_lin.f -L/home/user/.local/lib -llapack -lblas ==================================================================== program LinearOverEstimate c solving the fat (thin) matrix equation A * x = b using LAPACK, LSS Implicit none c declarations Real*4 A (3, 2), b(3), w(6) integer i, ok c define matrix A A (1, 1) = 1 ; A (1, 2) = 2 A (2, 1) = 1 ; A (2, 2) = 1 A (3, 1) = 3 ; A (3, 2) = -1 c define vector b b(1) = 5 ; b(2) = 3 ; b(3) = 1 c find the solution using the LAPACK routine call sgels ('n', 3, 2, 1, A, 3, b, 3, w, 6, ok) c Parameters: c 'n' - the linear system involves A; c #row of the matrix A c #col of the matrix A c #rhs, i.e., #col of the matrices B and X c array, on entry, the #row-by-#col matrix A c on exit, the details of factorization c lda of A. >= max(1,#row) c (input/output) array, dimension (ldb,#rhs) c on entry, the matrix B of #rhs vectors c on exit, if res = 0, B is overwritten by the solution vectors c ldb of the array B. >= max (#row, #col) c (workspace/output) array WORK, dim = LWORK c on exit, if res = 0, WORK(1) returns the optimal LWORK c dim of the array WORK c LWORK >= #col + max (#col, #rhs) c result c = 0: successful exit c < 0: if = -i, the i-th argument had an illegal value c > 0: if = i, the i-th diagonal element of the triangular factor is zero, c so A does not have full rank and could not be computed c c print the result vector if (ok .eq. 0) then do i = 1, 2 write (*, '(F8.3)') b(i) end do else write (*, *) "error" endif end program LinearOverEstimate ! gfortran lapa_over_sys.f -L/home/user/.local/lib -llapack -lblas ==================================================================== program LinearUnderEstimate c solving the (thick) matrix equation A * x = b using LAPACK, LNS Implicit none c declarations Real*4 A (2, 3), b(3), w(6) integer i, ok c define matrix A A (1, 1) = 1 ; A (1, 2) = 2 ; A (1, 3) = 3 A (2, 1) = 1 ; A (2, 2) = 1 ; A (2, 3) = 1 c define vector b b(1) = 6 ; b(2) = 3 c find the solution using the LAPACK routine call sgels ('n', 2, 3, 1, A, 2, b, 3, w, 6, ok) c c Parameters: c 'n': the linear system involves A c #row of the matrix A c #col of the matrix A c #rhs, i.e., #col of the matrices B and X c array, on entry, the #row-by-#col matrix A c on exit, the details of factorization c lda of A >= #row c (input/output) array, dimension (ldb, #rhs) c on entry, the matrix B of rhs vectors, stored columnwise c on exit, if res = 0, B is overwritten by the solution vectors c ldb of the array B >= #col c (workspace/output) array WORK, dim = LWORK c on exit, if res = 0, WORK(1) returns the optimal LWORK c dim of the array WORK c LWORK >= #row + max (#row, #rhs) c result c = 0: successful exit c < 0: if = -i, the i-th argument had an illegal value c > 0: if = i, the i-th diagonal element of the triangular factor is zero, c so A does not have full rank and the could not be computed c c print the result vector if (ok .eq. 0) then do i = 1, 3 write (*, '(F8.3)') b(i) end do else write (*, *) "error" endif end program LinearUnderEstimate ! gfortran lapa_under_sys.f -L/home/user/.local/lib -llapack -lblas ==================================================================== Program Eigenvalues c finding the eigenvalues of a complex matrix using LAPACK Implicit none c declarations complex*16 A (3,3), b(3), DR(3), DL(3), W(6) double precision RW(6) integer i, ok c define matrix A A (1, 1) = ( 3.1, -1.8) A (1, 2) = ( 1.3, 0.2) A (1, 3) = (-5.7, -4.3) A (2, 1) = ( 1.0, 0.0) A (2, 2) = (-6.9, 3.2) A (2, 3) = ( 5.8, 2.2) A (3, 1) = ( 3.4, -4.0) A (3, 2) = ( 7.2, 2.9) A (3, 3) = (-8.8, 3.2) call zgeev ('N', 'N', 3, A, 3, b, DL, 1, DR, 1, W, 6, RW, ok) c c parameters: c no left eigenvectors, c no right eigenvectors, c order of input matrix A, c input matrix A, c leading dimension of A, c array for eigenvalues, c array DL for left eigenvalues, c leading dimension of DL, c array DR for right eigenvalues, c leading dimension of DR, c workspace array W (dim >= 2 * order of A), c dimension of W c workspace array RW (dim = 2 * order of A), c return value c output of eigenvalues if (ok .eq. 0) then do i = 1, 3 write (*, *) b(i) enddo else write (*, *) "error" endif end program Eigenvalues ! gfortran lapa_eig.f -L/home/user/.local/lib -llapack -lblas