Tutorial de Fortran

[1. Introducción ] [Tutorial de Fortran] [3. LAPACK ]

BLAS

BLAS son las siglas de "Basic Linear Algebra Subroutines". Como su nombre lo indica, éste contiene subprogramas para operaciones básicas con vectores y matrices. BLAS fue diseñado para ser usado como un tabique dentro de otros códigos, por ejemplo en LAPACK. El código fuente se encuentra disponible a través de Netlib. Sin embargo, muchos vendedores de computadoras tienen versiones particulares de BLAS sintonizadas para tener un desempeño y eficiencias altos en sus equipos. Esta es una de las ventajas principales de BLAS: las secuencias de llamadas estan estandarizadas por los programas que llamen a BLAS trabajarán en cualquier computadora donde BLAS este instalado. Si se tiene una versión rápida de BLAS, también se tendrá un alto rendimiento en todos los programas que llamen a BLAS. Por lo que BLAS proporciona una forma simple y portable para llevar a cabo alto rendimiento para cálculos donde se involucre álgebra lineal. LAPACK es paquete de alto nivel construído con las misma ideas.

Niveles y convenciones para nombres

Las subrutinas de BLAS pueden ser divididas en tres niveles:

Nivel 1: Operaciones vector a vector. Datos en O(n) y trabajo en O(n).
Nivel 2: Operaciones matriz-vector. Datos en O(n^2) y trabajo en O(n^2).
Nivel 3: Operaciones matriz-matriz. Datos en O(n^2) y trabajo en O(n^3).

Cada rutina en BLAS y LAPACK se encuentra disponible en varias versiones, una para cada tipo de dato. La primera letra del nombre del subprograma indica la precisión usada:

      S      Precisión simple Real.
      D      Precisión doble Real.
      C      Precisión simple Complex.
      Z      Precisión doble Complex.

La precisión doble Comples no esta estrictamente definida en Fortran 77, pero muchos compiladores aceptarán alguna de las siguientes declaraciones:

      double complex lista_de_variables
      complex*16     lista_de_variables

BLAS 1

Algunos subprogramas de BLAS 1 son:

xCOPY - copia un vector en otro
xSWAP - intercambia dos vectores
xSCAL - escala un vector por una constante
xAXPY - agregar un múltiplo de un vector en otro
xDOT - producto interno
xASUM - norma-1 de un vector
xNRM2 - norma-2 de un vector
IxAMAX - encontrar la entrada máxima en un vector

La primera letra (x) puede ser una de las letras S,D,C,Z dependiendo de la precisión. Una referencia rápida a BLAS 1 puede ser revisada en http://www.netlib.org/blas/blasqr.ps

BLAS 2

Algunos de los subprogramas de BLAS 2 son:

xGEMV - multiplicación general de una matriz por un vector
xGER - actualización general rango-1
xSYR2 - actualización simétrica rango-2
xTRSV - resuelve un sistema triangular de ecuaciones

Una descripción detallada de BLAS 2 puede ser encontrada en http://www.netlib.org/blas/blas2-paper.ps.

BLAS 3

Algunos de los subprogramas de BLAS 3 son:

xGEMM - general matrix-matrix multiplication
xSYMM - symmetric matrix-matrix multiplication
xSYRK - symmetric rank-k update
xSYR2K - symmetric rank-2k update

The more advanced matrix operations, like solving a linear system of equations, are contained in LAPACK. A detailed description of BLAS 3 can be found at http://www.netlib.org/blas/blas3-paper.ps.

Ejemplos

Veamos primero un ejemplo sencillo de una rutina de BLAS, SSCAL. La secuencia de llamada es

      call SSCAL ( n, a, x, incx )

Donde x es el vector, n es la longitud (número de elementos en x que se quieren modificar), y a es un escalar por el que se quiere multiplicar a x. El último argumento incx es el incremento. Usualmente, incx=1 y el vector x corresponde directamente a un arreglo unidimensional de Fortran. Para incx>1 este indica cuantos elementos en el arreglo serán "brincados" entre cada elemento del vector x. Por ejemplo, si incx=2 significa que se deberán escalar los elementos nones (nota: la dimensión física del arreglo x deberán ser al menos 2n-1). Considerar los siguientes ejemplos donde x ha sido declarado como real x(100).

      call SSCAL(100, a, x, 1)
      call SSCAL( 50, a, x(50), 1)
      call SSCAL( 50, a, x(2), 2)

La primera línea escalará los 100 elementos de x por a. La siguiente línea escalará solamente los últimos 50 elementos de x por a. La última línea escalará todos los de índice par de x por a.

Observar que el arreglo x será sobreescrito por los nuevos valores. Si se necesita preservar una copia de los valores anteriores de x, se tiene que hacer primero una copia, y para eso se puede usar SCOPY.

Ahora se considerara un ejemplo complicado. Suponer que se tienen arreglos bidimensionales A y B, y se quiere encontrar la entrada (i,j) del producto A*B. Esto es hecho fácilmente haciendo el producto interno del renglón i de A y la columna j de B. Se puede usar la subrutina SDOT de BLAS 1. La única dificultad es imaginar los índices correctos y sus incrementos. La llamada a la secuencia SDOT es

      SDOT ( n, x, incx, y, incy )

Suponiendo que las declaraciones de los arreglos fueron

      real A(lda,lda)
      real B(ldb,ldb)

pero en el programa se sabe que el tamaño actua de A es m*p y para B es p*n. El i-ésimo renglón de A inicia en el elemento A(i,1). Pero como Fortran almacena los arreglos bidimensionales en columnas descendentes, el siguiente elemento del renglón A(i,2) será guardado lda elementos después en la memoria (ya que lda es la longitud de una columna). Por lo tanto se pone incx = lda. Para la columna en B no hay tal problema, los elementos son guardados consecutivamente por lo tanto incy = 1. La longitud del calculo del producto interno es p. Por lo que la respuesta es

      SDOT ( p, A(i,1), lda, B(1,j), 1 )

Como obtener BLAS

Lo primero que se debe hacer es revisar si se tiene BLAS en el sistema. Si no, se puede encontrar en Netlib en http://www.netlib.org/blas.

Documentación

Las rutinas de BLAS son casi auto explicatorias. Una vez que se sabe cual es la rutina que se ocupa, se obtiene y se lee la sección de la cabezera que explica los parámetros de entrada y salida en detalle. Por ejemplo, se muestra la cabezera de la subrutina SGEMV de BLAS 2:

      SUBROUTINE SGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
     $                   BETA, Y, INCY )
*     .. Scalar Arguments ..
      REAL               ALPHA, BETA
      INTEGER            INCX, INCY, LDA, M, N
      CHARACTER*1        TRANS
*     .. Array Arguments ..
      REAL               A( LDA, * ), X( * ), Y( * )
*     ..
*
*  Purpose
*  =======
*
*  SGEMV  performs one of the matrix-vector operations
*
*     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
*
*  where alpha and beta are scalars, x and y are vectors and A is an
*  m by n matrix.
*
*  Parameters
*  ==========
*
*  TRANS  - CHARACTER*1.
*           On entry, TRANS specifies the operation to be performed as
*           follows:
*
*              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
*
*              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
*
*              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
*
*           Unchanged on exit.
*
*  M      - INTEGER.
*           On entry, M specifies the number of rows of the matrix A.
*           M must be at least zero.
*           Unchanged on exit.
*
*  N      - INTEGER.
*           On entry, N specifies the number of columns of the matrix A.
*           N must be at least zero.
*           Unchanged on exit.
*
*  ALPHA  - REAL            .
*           On entry, ALPHA specifies the scalar alpha.
*           Unchanged on exit.
*
*  A      - REAL             array of DIMENSION ( LDA, n ).
*           Before entry, the leading m by n part of the array A must
*           contain the matrix of coefficients.
*           Unchanged on exit.
*
*  LDA    - INTEGER.
*           On entry, LDA specifies the first dimension of A as declared
*           in the calling (sub) program. LDA must be at least
*           max( 1, m ).
*           Unchanged on exit.
*
*  X      - REAL             array of DIMENSION 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.
*           Unchanged on exit.
*
*  INCX   - INTEGER.
*           On entry, INCX specifies the increment for the elements of
*           X. INCX must not be zero.
*           Unchanged on exit.
*
*  BETA   - REAL            .
*           On entry, BETA specifies the scalar beta. When BETA is
*           supplied as zero then Y need not be set on input.
*           Unchanged on exit.
*
*  Y      - REAL             array of DIMENSION at least
*           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
*           and at least
*           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
*           Before entry with BETA non-zero, the incremented array Y
*           must contain the vector y. On exit, Y is overwritten by the
*           updated vector y.
*
*  INCY   - INTEGER.
*           On entry, INCY specifies the increment for the elements of
*           Y. INCY must not be zero.
*           Unchanged on exit.
*
*
*  Level 2 Blas routine.
*
*  -- Written on 22-October-1986.
*     Jack Dongarra, Argonne National Lab.
*     Jeremy Du Croz, Nag Central Office.
*     Sven Hammarling, Nag Central Office.
*     Richard Hanson, Sandia National Labs.
*

[1. Introducción ] [Tutorial de Fortran] [3. LAPACK ]