Portando rotinas Fortran para a GêBR

GêBR Project

22/06/2009

A GêBR atua como interface para programas (preferencialmente) de linha de comando capazes de serem controlados por parâmetros de linha de comando.  Em princípio, não é possível incluir na GêBR diretamente apenas uma subrotina. Para isto seria necessário um programa principal que a evocasse. Este tutorial explica como portar uma subrotina, escrita em Fortran, para dentro da GêBR. A estratégia para isto será criar um programa principal, escrito em C, capaz de obter todos os parâmetros necessários à subrotina Fortran, diretamente da linha de comando, e então executar a subrotina. Discutiremos também como realizar a leitura e escrita dos dados sísmicos no formato nativo do Seismic Un*x, de maneira a poder integrar o programa gerado com todos os outros menus que já se valem deste formato padrão.

A rotina hipotética que estamos portando, depende dos parâmetros ns (número de amostras de um traço), np (um parâmetro inteiro), y (um parâmetro real de precisão simples), e xin (um vetor de números reais de precisão simples que armazena o traço de entrada). O resultado é escrito em xout (também um vetor de números reais de precisão simples). O protótipo desta rotina é:

subroutine proc (ns, np, y, xin, xout)

Apresentaremos a seguir o programa principal, por blocos, comentando-o linha por linha. Em amarelo estão destacadas as porções do programa que devem ser escritas de acordo com sua aplicação.

Cabeçalho

O programa se inicia com:

#include <stdlib.h>

#include <stdio.h>

#include <string.h>

#include <glib.h>

#include "header.h"

As três primeiras linhas carregam bibliotecas padrão do C e não devem ser removidas a menos que você saiba o que está fazendo. A seguir, o header da biblioteca Glib é incluído, pois está será a responsável pelo parse dos parâmetros de linha de comando. A última linha inclui a descrição do header de um dado no formato do pacote Seismic Un*x. Salve o arquivo header.h na mesma pasta do programa C e da subrotina Fortran.

Como a subrotina Fortran está implementada num arquivo .f90 e não dispões de header, é necessário declará-la no programa principal, para poder utiliza-la adequadamente. Note que em Fortran todos os parâmetros são passados por referência. Além disso, se o nome da rotina Fortran for proc o nome do objeto gerado para ser lincado com o código C, tem um “_” acrescido ao final. Assim, dentro do programa C, a rotina será chamada proc_.

extern proc_(int*, float*, float*, int*, float*, float*);

Leitura de parâmetros da linha de comando

O programa principal iniciará, criando as estruturas necessárias para o parse dos parâmetros passados em linha de comando. Isso é feito utilizando a biblioteca Glib, que provê suporte para esta tarefa.

int main(int argc, char **argv){

   /* Variables to hold command-line parameters */

   /* and their default values                  */

   static gint    np = -1;

   static gdouble  y =  0;

Primeiramente definimos duas variáveis para armazenar os valores dos parâmetros que serão lidos da linha de comando, sendo uma inteira (np) e outra real (y). Perceba que a rotina Fortran, tem outros parâmetros, como ns (número de amostras em um traço), xin e xou (traço de entrada e traço de saída), porém estas outras informações não serão obtidas da linha de comando, mas sim do header do traço (no caso de ns) e via entrada e saída padrão (no caso dos traços a serem processados).

   /* Variables to read header and data trace */

   su_header_t  hdr;

   float       *xin;

   float       *xout;

   /* Other variables */

   int           ns;

   double        dt;

   float         y_single;

A seguir definimos mais três variáveis, hdr que armazenará o header de um traço lido, xin e xout que armazenarão as amostras do traço de entrada e do traço computado de saída. Repare que tanto xin como xout não tem dimensão definida. Estes vetores terão sua dimensão descoberta em tempo de execução, após a leitura do header do traço, o que nos permitirá conhecer a quantidade de amostras. Após isto, teremos como alocar exatamente a memória necessária para estas duas variáveis.

   /* Command-line parameters definition */

   static GOptionEntry entries[] =    {

         { "npar", 0, 0, G_OPTION_ARG_INT,   &np, "some integer parameter", "n" },

         { "pary", 0, 0, G_OPTION_ARG_DOUBLE, &y, "some real number parameter", "x" },

         { NULL }

     };

O bloco acima define a estrutura necessária para a captura dos parâmetros em linha de comando. Cada linha define um parâmetro, e os campos são, nesta ordem, um string que será utilizado como palavra chave na linha de comando, dois campos mantidos no default (0), a especificação do tipo de parâmetro, o endereço da variável onde o parâmetro lido será armazenado, um string contendo uma descrição sucinta do parâmetro, e um último string que será utilizado de exemplo para um possível valor do parâmetro.

Abaixo está o bloco de programa que realiza o parse dos parâmetros em linha de comando. Como a biblioteca utilizada, automaticamente, cria um help para o programa, é possível entrar com textos descritivos do programa.

   GError *error = NULL;

   GOptionContext *context;

   /* Set a short description for the program */

   context = g_option_context_new ("- Put here an one-line description of the program");

   /* Summary */

   g_option_context_set_summary (context,

       "Brief summary of program functionalities comes here.\n"

       "Indeed, it can take as many lines as you want."

                                );

   /* Description */

   g_option_context_set_description (context,

        "Not so brief description, which comes after\n"

        "the parameter list. Usualy, authors and email\n"

        "addresses are included here, as well as references."

                                    );

   g_option_context_add_main_entries (context, entries, NULL);

   /* Complain about unknown options */

   g_option_context_set_ignore_unknown_options (context, FALSE);

   /* Parse command line */

   if (g_option_context_parse (context, &argc, &argv, &error) == FALSE){

    fprintf(stderr, "%s: syntax error\n", argv[0]);

    fprintf(stderr, "Try %s --help\n", argv[0]);

    return EXIT_FAILURE;

   }

   g_option_context_free (context);

Neste ponto, os parâmetros de linha de comando, np e y foram lidos através das palavras chave “npar” e “pary”.

Bloco principal

   /*-------------------------------------------------------------*/

   /* Main code */

   y_single = (float) y;

   /* Read just the first header */

   fread(&hdr, SIZEOF_SEGYHDR, 1, stdin);

   ns = hdr.ns;

   dt = hdr.dt / 1.0e6;

O comando fread, lê o header do primeiro traço, a partir da entrada padrão (stdin), armazenando-o em hdr. Na variável ns guardamos o número de amostras do traço, obtido do campo ns do header. Na variável dt guardamos a razão de amostragem em segundos (o padrão do SU é armazenar este número em micro-segundos).

De posse do número de amostras de um traço, alocamos memória para armazenar os traços de entra e de saída.

   /* Memory allocation */

   xin  = (float *) malloc(sizeof(float) * ns);

   xout = (float *) malloc(sizeof(float) * ns);

Agora já possível completar a leitura do primeiro traço do dado, lendo suas amostras para o vetor xin. Finalmente, iniciamos o loop de ler as amostras do traço, processa-lo, salvar o resultado (header e amostras), até que não haja mais traços a serem lidos.

   do{

     fread(xin, sizeof(float), ns, stdin);

     proc_(&np, &ysingle, &dt, &ns, xin, xout);

     fwrite (&hdr, SIZEOF_SEGYHDR, 1, stdout);

     fwrite (xout, sizeof(float), ns, stdout);

   }while (fread (&hdr, SIZEOF_SEGYHDR, 1, stdin));

Quando sairmos deste loop, todos os traços terão sido processados. Resta apenas liberar toda a memória alocada dinamicamente (em tempo de execução) e encerrar o programa.

  /* Memory release */

   free (xin);

   free (xout);

   return EXIT_SUCCESS;

}

Compilando

Abaixo apresentaremos o arquivo de Makefile que compila e linca o código Fortran com o programa principal em C.

# Defina na variavel PROGF o nome do arquivo com

# a rotina Fortran (sem a extensão .f90)

PROG = progf

# Compiladores

CC = gcc

F90 = gfortran

# Flags para a biblioteca GLib

GLIB_CFLAGS = `pkg-config --cflags glib-2.0`

GLIB_LIBS = `pkg-config --libs glib-2.0`

# Flag para biblioteca Fortran necessaria ao link com C

FORTRAN_LIBS = `gfortran -print-file-name=libgfortran.so`

CFLAGS = -g -O2 $(GLIB_CFLAGS)

LIBS = $(GLIB_LIBS) $(FORTRAN_LIBS)

$(PROGF): $(PROGF).o main.o

        $(CC) $(LIBS) main.o $(PROGF).o -o $@

# Troque a extensão de .f90 para .f, dependendo do seu

# arquivo fortran

$(PROGF).o: $(PROGF).f90

        $(F90) -O2 -fPIC -c $<

main.o: main.c header.h

        $(CC) $(CFLAGS) -c $<

clean:

        rm -f $(PROGF).o main.o $(PROGF) *~

.PHONY: clean

Compile o programa digitando make. O programa executável criado será progf. Executando o programa com a opção “--help”, será exibida a tela com a descrição dos parâmetros e demais detalhes, como abaixo:

Usage:

   progf [OPTION...] - Put here an one-line description of the program

 

Brief summary of program functionalities comes here.

Indeed, it can take as many lines as you want.

 

Help Options:

   -?, --help       Show help options

 

Application Options:

   --npar=n         some integer parameter

   --pary=x         some real number parameter

 

Not so brief description, which comes after

the parameter list. Usualy, authors and email

addresses are included here, as well as references.

Para utilizar este programa, rode-o como:

./progf --npar=3 --pary=2.5 saida.su

Requisitos

Poucos requisitos devem ser cumpridos para que o esquema apresentado funcione. No exemplo construído, foi utilizado o compilador Fortran livre, parte da suíte de compiladores GNU, gfortran, disponível para diversas arquiteturas diferentes. Não é obrigatório que este seja o compilador empregado, porém uma escolha diferente pode demandar ajustes no arquivo Makefile.

É necessário ainda que as bibliotecas de desenvolvimento da GLib estejam instaladas, uma vez que ela foi utilizada para a leitura dos parâmetros passados em linha de comando. A GLib é uma biblioteca com várias ferramentas auxiliares, implementada em C e pode ser facilmente encontrada. De fato, ela faz parte de praticamente todas as distribuições Linux, pois é base para muitos programas.

Abaixo apresentamos todo o código continuamente, para que fique claro sua inspeção.

#include <stdlib.h>

#include <stdio.h>

#include <string.h>

#include <glib.h>

#include "header.h"

extern proc_(int*, float*, float*, int*, float*, float*);

int main(int argc, char **argv){

   /* Variables to hold command-line parameters */

   /* and their default values                  */

   static gint    np = -1;

   static gdouble  y =  0;

   /* Variables to read header and data trace */

   su_header_t  hdr;

   float       *xin;

   float       *xout;

   /* Other variables */

   int           ns;

   float         dt;

   float         y_single;

   /* Command-line parameters definition */

   static GOptionEntry entries[] =    {

         { "npar", 0, 0, G_OPTION_ARG_INT,   &np, "some integer parameter", "n" },

         { "pary", 0, 0, G_OPTION_ARG_DOUBLE, &y, "some real number parameter", "x" },

         { NULL }

     };

   GError *error = NULL;

   GOptionContext *context;

   /* Set a short description for the program */

   context = g_option_context_new ("- Put here an one-line description of the program");

   /* Summary */

   g_option_context_set_summary (context,

        "Brief summary of program functionalities comes here.\n"

        "Indeed, it can take as many lines as you want."

                                );

   /* Description */

   g_option_context_set_description (context,

         "Not so brief description, which comes after\n"

        "the parameter list. Usualy, authors and email\n"

        "addresses are included here, as well as references."

                                    );

   g_option_context_add_main_entries (context, entries, NULL);

   /* Complain about unknown options */

   g_option_context_set_ignore_unknown_options (context, FALSE);

   /* Parse command line */

   if (g_option_context_parse (context, &argc, &argv, &error) == FALSE){

     fprintf(stderr, "%s: syntax error\n", argv[0]);

     fprintf(stderr, "Try %s --help\n", argv[0]);

     return EXIT_FAILURE;

   }

   g_option_context_free (context);

   /*-------------------------------------------------------------*/

   /* Main code */

   y_single = (float) y;

   /* Read just the first header */

   fread(&hdr, SIZEOF_SEGYHDR, 1, stdin);

   ns = hdr.ns;

   dt = hdr.dt / 1.0e6;

   /* Memory allocation */

   xin  = (float *) malloc(sizeof(float) * ns);

   xout = (float *) malloc(sizeof(float) * ns);

   do{

      fread(xin, sizeof(float), ns, stdin);

      proc_(&np, &ysingle, &dt, &ns, xin, xout);

      fwrite (&hdr, SIZEOF_SEGYHDR, 1, stdout);

      fwrite (xout, sizeof(float), ns, stdout);

   }while (fread (&hdr, SIZEOF_SEGYHDR, 1, stdin));

   /* Memory release */

   free (xin);

   free (xout);

   return EXIT_SUCCESS;

}

Referências

O exemplo apresentado deve ser suficiente para a maioria dos usuários. Caso precise, maiores informações sobre a leitura e o tratamento de parâmetros de linha de comando estão disponíveis na documentação da GLib, a biblioteca utilizada para este fim.