$>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.
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, subtractionwhen 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
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/2the 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 = -xif more than one statement should be executed inside the if, then the following syntax should be used:
if (logical expression) then
statements
endifthe 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 here 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) 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 UNDEFINED
this 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)
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 column
each 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 lda is the leading (row) dimension of array A
note that lda is 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)
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
the difference between subroutine and function is that function must return value and should't use side-effect code (IO or changing of varibables' values in main program)
so you should declare its return value before keyword function in function
declaration part. after the keyword function must be its name. the last statement in function
body must be assignment of value to variable with function name
example :
program example_fun
real k , t
print * , 'put your num : '
read (*, 200) t
k = addTwo (t)
print * , '2.0 has been added'
write (*, 200) k
200 format (F8.2)
end program example_fun
real function addTwo (x)
addTwo = x + 2.0
end function
| 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
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)
[UNIT=] u
IOSTAT= ios
ERR= err
FILE= fname
STATUS= sta
ACCESS= acc
FORM= frm
RECL= rl
read ([UNIT=]u, [FMT=]fmt, IOSTAT=ios, ERR=err, END=s)
write ([UNIT=]u, [FMT=]fmt, IOSTAT=ios, ERR=err, END=s)
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:
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
BLAS (Basic Linear Algebra Subprograms) are a set of low-level subroutines that perform common linear algebra operations such as copying, vector scaling, vector dot products, linear combinations, and matrix multiplication. BLAS are used as a building block in higher-level math programming languages and libraries, including LAPACK and R
BLAS functionality is divided into three levels: 1, 2 and 3
Level 1: y ← α * x + y
vector-vector operations that are linear (O(n)) in data and linear (O(n)) in work
Level 2: y ← α * A * x + β * y
matrix-vector operations that are quadratic (O(n^2)) in data and quadratic (O(n^2)) in work
Level 3: B ← α * A1 * A2 + β * B
operations that are quadratic (O(n^2)) in data and cubic (O(n^3)) in work
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
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
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
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
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
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) by matrix (2x2) multiplication with scalar factor plus
c scalar multiplied 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 * 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.
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.
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.
====================================================================
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