fortran


  • data types
  • if and loops
  • arrays
  • subroutines
  • IO
  • libraries
  • BLAS
  • LAPACK

  • 
    
    $>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


    Data types

    
          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:
  • The name defined here is a constant - you cannot change its value at a later point in the program
  • A name can appear in at most one parameter statement
  • The parameter statements must come before the first executable statement
  • 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

    If and Loops

    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
    

    Arrays

    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 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)
    
    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 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)
    
    are equivalent to
    
          real A (10,20), x(50)
    
    The dimension statement is considered old-fashioned style


    Subroutines

    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

    intrinsic subroutines

    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


    IO

    "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
    

    
          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:

    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))
    

    examples

          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

    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

    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 )

    examples of BLAS usage

    
          
          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
    
          
      

    LAPACK

    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)
    "the minimum norm solution" (underdetermined 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.

    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