Matrizes Em Fortran 90

See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/265283632

Matrizes em Fortran 90  Article

CITATIONS

READS

0

4,684

2 authors:

Flavio Aristone

Juvenal Muniz

Federal University of South Mato Grosso, Brazil

Universidade Federal de Mato Grosso do Sul

37 PUBLICATIONS 126 CITATIONS

1 PUBLICATION 0 CITATIONS

SEE PROFILE

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Electronic transport properties measurements for low-dimensional semiconductor heterostructure devices View project

PLANO DE INOVAÇÃO DA CADEIA DA BOCAIUVA DO ESTADO DO MATO GROSSO DO SUL View project

Matrizes em Fortran 90 Flavio Aristone Juvenal Muniz 1 de maio de 2010 Resumo

Discutiremos neste documento a implementa¸c˜ ao de opera¸co ˜es b´ asicas sobre matrizes na linguagem de programa¸c˜ ao Fortran 90: adi¸ c˜ ao, multiplica¸c˜ ao, transposta e inversa. Calcularemos o determinante de uma matriz pela sua decomposi¸ca ˜o LU.

Sum´ ario 1 Introdu¸ c˜ ao

1

2 Teoria

1

2.1 Arrays, vetores e matrizes . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Aloca¸ca˜o dinˆamica de arrays . . . . . . . . . . . . . . . . . . . . . . 2.2 Fun¸c˜oes e subrotinas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Implementa¸ ca ˜o

4

Referˆ encias

1

2 3 3

11

Introdu¸ c˜ ao

Uma matriz ´e um arranjo retangular de n´ umeros. Por exemplo,

8 3 A =  4

  

7 2 5 10 11 3 1 15 7

As opera¸co˜es b´asicas que podem ser aplicadas a matrizes s˜ ao:   adi¸c˜  ao,  multiplica¸cao, ˜  transposi¸c˜  ao e invers˜  ao. Nossa preocupa¸ca˜o maior neste documento ´e a implementa¸ca˜o destas opera¸c˜oes na linguagem de programa¸ca˜o Fortran 90. Para uma discuss˜ ao te´ orica sobre matrizes, sugerimos o livro [3].

2

Teoria

Antes de implementarmos as opera¸co˜es matriciais, vamos fornecer um breve resumo de alguns conceitos da linguagem Fortran 90 que iremos utilizar mais adiante neste documento.

1

2.1

Arrays, vetores e matrizes

Para implementar matrizes na linguagem Fortran 90 iremos utilizar uma array1 . Uma array ´e uma estrutura de dados que consiste de uma cole¸c˜ aveis ou ao   de valores (vari´ constantes) de um mesmo tipo que s˜ao referidos por um mesmo nome. Um valor individual dentro da array ´e chamado um  elemento da array  e ´e identificado pelo nome da array  juntamente com um ´ındice   indicando a localiza¸c˜ao do elemento dentro da array. Por exemplo, a seguinte cole¸ca˜ o de 5 n´ umeros inteiros

{5 7 9 1 8} ,

,

,

,

poder´a ser representada em Fortran 90 por i n t e g e r ,   dimension ( 5)

: :   numeros   = ( / 5 , 7 , 9 , 1 , 8 / )

O trecho de c´ odigo acima declara uma array chamada “numeros” contendo 5 n´ umeros inteiros. No c´ odigo, a array “numeros” ´e inicializada j´ a na sua pr´ opria declara¸c˜ao pelo construtor  de arrays (/ . . . /). Outra forma de declarar e inicializar esta array seria i n t e g e r ,   dimension ( 5) numeros ( 1 ) numeros ( 2 ) numeros ( 3 ) numeros ( 4 ) numeros ( 5 )

= = = = =

: :   numeros

5 7 9 1 8

O n´umero entre os parˆenteses ´e o ´ındice que localiza a posi¸ c˜a o dentro da array; o ´ındice 1 localiza a primeira posi¸ca˜o na array, o 2 a segunda e assim em diante. Portanto, o c´odigo acima atribui o valor 5 `a posi¸c˜ao 1 da array, o valor 7 a` posi¸ca˜o 2, etc. O sinal de igual “=” ´e a forma de atribuir valores a vari´ aveis em Fortran 90. As arrays em Fortran 90 podem ter mais de um ´ındice, podendo ser organizadas em m´ultiplas dimens˜oes. Essas arrays s˜ ao convenientes para representar dados organizados, por exemplo, em forma matricial (linhas e colunas). Utilizaremos arrays bidimensionais (dois ´ındices) para representar matrizes. Chamaremos de  matrizes  as arrays bidimensionais e de  vetores  as arrays unidimensionais. Por exemplo, se quis´essemos representar a matriz identidade de ordem 3 Id3

1 = 0

 

0 0 1 0 0 0 1

poder´ıamos utilizar o seguinte c´ odigo em Fortran 90 i n t e g e r : : i , j i n t e g e r ,   dimension  ( 3 , 3 ) : :   identidade do i   = 1 , 3 do j   = 1 , 3 i f  ( i /= j )   then identidade ( i ,  j ) = 0 else identidade ( i ,  j ) = 1 e nd i f   e nd d o e nd d o 1

Por falta de uma melhor tradu¸c˜ ao da palavra, iremos utilizar o termo em inglˆes neste documento.

2

Neste trecho de c´ odigo utilizamos as vari´aveis i e j  para representar os valores das linhas e colunas, respectivamente. O la¸co externo percorre as linhas da matriz, equanto que o la¸co interno percorre as colunas. Para cada elemento da matriz ´e atribu´ıdo o valor 0, caso a posi¸ca˜o deste elemento esteja fora da diagonal (i  = j), ou o valor 1, caso o elemento esteja na diagonal (i = j). 2.1.1

Aloca¸ c˜ ao dinˆ amica de arrays

Em todos os trechos de c´ odigo acima, as declara¸co˜es das arrays s˜ ao do tipo denominado aloca¸ c˜ ao est´ atica da mem´ oria, isso porque o tamanho de cada array ´e escolhido no momento da declara¸c˜a o e n˜ ao poder´ a mudar. Desta maneira, o tamanho da array dever´ a ser grande o suficiente para que todos os valores possam ser armazenados sem problemas. A aloca¸c˜ao est´ atica de mem´oria pode apresentar s´erias limita¸ co˜es. Por exemplo, suponha que vocˆe declare uma array de 1000 posi¸ c˜oes. Essa array ir´ a ocupar parte da mem´ oria do computador para armazenar 1000 n´ umeros, mesmo que vocˆe s´ o utilize 2 posi¸co˜es dela, ou seja, vocˆe estaria desperdi¸cando 998 posi¸c˜oes de mem´ oria. Uma solu¸ca˜o para este problema ´e a chamada  aloca¸c˜ ao dinˆ amica da mem´ oria, exemplificada no seguinte c´ odigo i n t e g e r ,   a l l o c a t a b l e ,   dimension ( : , : )

: :   identidade

O c´odigo acima declara a matriz identidade com aloca¸c˜ao dinˆamica da mem´ oria. Observe que os sinais de dois pontos, :, entre os parˆenteses, s˜ ao obrigat´ orios, assim como a palavra “allocatable”. Observe tamb´em que a v´ırgula ap´ os o primeiro sinal de dois pontos foi necess´ aria, pois estamos declarando uma array bidimensional (uma matriz). Com este tipo de declara¸ca˜o, nenhuma mem´ oria ´e reservada. Para utilizarmos a matriz identidade, temos que alocar (reservar) mem´ oria antes. Veja o seguinte c´ odigo i n t e g e r ,   a l l o c a t a b l e ,   dimension ( : , : )

: :   identidade

a l l o c a t e ( identidade  (3 ,3) )

Para alocarmos mem´ oria, utilizamos a fun¸c˜ao   allocate   com o nome da matriz e o 2 tamanho que quisermos .   Feito isto, podemos atribuir valores aos elementos da matriz normalmente. Ap´ os utilizarmos a matriz, poderemos liberar a mem´ oria para ela reservada, utilizando a fun¸ca˜o  deallocate, como a seguir d e a l l o c a t e ( identidade )

2.2

Fun¸ c˜ oes e subrotinas

Utilizamos certos procedimentos em Fortran 90, tais como  fun¸ c˜ oes e  subrotinas, para organizar melhor nossos programas. Essa organiza¸ c˜ao facilita a corre¸ca˜o de erros, estrutura¸ca˜o do programa em tarefas e subtarefas, etc. Iremos apenas mostrar como utilizar fun¸co˜es e subrotinas em Fortran 90. Para mais informa¸co˜es sobre este assunto, sugerimos o livro [1]. Basicamente, subrotinas s˜ ao procedimentos que s˜ ao executados por uma instru¸c˜ao chamada  CALL. Subrotinas podem retornar m´ ultiplos resultados atrav´es da sua lista de argumentos de chamada. Fun¸ c˜oes s˜ao procedimentos que s˜ ao executados pelos pr´ oprios 2

Isto ´ e, caso haja mem´oria suficiente para alocar a matriz.

3

nomes em express˜ oes e s´o podem retornar um u´nico resultado a ser usado na express˜ ao. Os exemplos a seguir d˜ao a forma geral de uma subrotina e de uma fun¸ c˜ao, respectivamente s u b r o u t i n e   minha_subrotina (   l i st a d e a r g u m en t o s ) ... (declara¸c˜ oes de vari´ aveis) ... (c´ odigo a ser executado) ... return e nd s u b r o u t i n e

f u n c t i o n   minha_funcao (   l i s ta d e a r g u me n t o s ) ... (declara¸c˜ oes de vari´ aveis) ... (c´ odigo a ser executado) ... ao  mi nh a_ fu nc ao   = express˜ return e nd f u n c t i o n

Tanto para subrotinas quanto para fun¸co˜es, ´e poss´ıvel estabelecer uma lista de argumentos que ser˜ ao utilizados pela subrotina ou fun¸ca˜o. Vale destacar que na fun¸ca˜o, ´e necess´ario que haja uma atribui¸ca˜o de valor ao nome da fun¸c˜ao, indicando o valor que ser´a retornado pela mesma.

3

Implementa¸ c˜ ao

O c´odigo a seguir ´e a implementa¸ca˜o em Fortran 90 das opera¸co˜es b´asicas sobre matrizes e o c´alculo do determinante de uma matriz pela sua decomposi¸c˜ao em matrizes triangulares inferior e superior, conhecida como decomposi¸ca˜o LU. O algoritmo utilizado ´e conhecido como   algoritmo de Doolittle. Utilizamos um algoritmo simples para o c´ alculo da inversa de uma matriz. Esse algoritmo ´e baseado no processo de elimina¸ c˜ao gaussiana e tem um car´ ater mais did´ atico, pois seu desempenho ´e baixo. Tentamos, exageradamente, comentar o c´ odigo para facilitar a compreens˜ a o do material. Utilizamos v´ arias t´ecnicas e constru¸co˜es sint´aticas para mostrar v´ arias formas de implementa¸c˜ao. O c´odigo ´e extenso devido ao seu car´ ater did´ atico, mas vocˆe poder´ a modific´ a-lo `a vontade, buscando melhores formas de expor os algoritmos ou, at´e mesmo, utilizar algoritmos mais eficientes e precisos! ´ aconselh´ E avel que vocˆe digite (ao inv´ es de copiar-e-colar) todo o c´ odigo e tente execut´ a-lo, pois assim, estar´a treinando suas habilidades para detectar e corrigir erros, caso ocorram. Vale ressaltar que este ´e um c´ odigo com car´ ater did´ atico e n˜ao deve ser utilizado em outros ambientes que n˜ ao os de sala de aula. Algumas instru¸c˜oes em Fortran 90 possuem op¸co˜es at´e agora n˜ ao utilizadas, como por exemplo, a op¸c˜ao   advance=‘no’   na instru¸ca˜o  write. N˜ao detalharemos o significado dessas op¸co˜es, mas todas elas s˜ao facilmente encontradas em [1].

4

! A vari´ avel seed ser´a compartilhada com a subrotina que cria os elementos das matrizes A e B ! aleatoriamente. Isto ´e necess´ ario para que os valores das matrizes n˜ ao sejam iguais. ! ! M´ odulos devem sempre vir definidos antes do seu uso no programa. module   global i m p l i c i t no ne i n t e g e r : :   seed end module

! In´ıcio do programa program   matrizes use   global ! ! Prop´ osito: ! Criar procedimentos para somar, multiplicar, transpor, inverter e decompor matrizes e ! calcular o determinante de uma matriz. ! ! Data: 1.5.2010 ! Programador: Juvenal Muniz ! ! Revis˜ oes ! 1.5.2010 : C´ odigo original i m p l i c i t no ne ! Declara¸ c˜ oes das vari´ aveis e das subrotinas utilizadas ! Descri¸ c˜ ao das vari´ aveis e/ou matrizes ! n  Dimens˜ao da matriz ! i  Contador !   condicao  Indica se houve erro ao abrir o arquivo !   opcao  Armazena a opcao selecionada pelo usu´ario !   determinante  Valor do determinante da matriz ! a,  b ,  c  Matrizes a, b e c !   arquivo  Constante utilizada para acessar o arquivo !   resposta  Armazenar´ a a resposta a ` pergunta (S ou N) !   nome arquivo Nome do arquivo contendo os elementos das matrizes ! Descri¸ c˜ ao das subrotinas !   matriz aleatoria  Gera uma matriz quadrada aleat´ oria de ordem n !   matriz adicao  Efetua a adi¸ca ˜o de duas matrizes !   matriz mul  Efetua a multiplica¸ca ˜o de duas matrizes !  matriz transposta  Obt´em a transposta de uma matriz !   matriz inversa  Obt´em a inversa de uma matriz !   matriz lu  Decomp˜ oe a matriz em matrizes triangulares inferior e superior !   matriz imprime  Exibe os elementos da matriz i n t e g e r : : n , i , j ,   condicao  ,   opcao r e a l : :   determinante r e a l ,   a l l o c a t a b l e ,   dimension  ( : , : ) : : a , b , c i n t e g e r ,   parameter : :   arquivo = 17 c h a r a c t e r ( 1 ) : :   resposta c h a r a c t e r (20 ) : :   nome_arquivo ! Valor inicial para a fun¸ c˜ ao RAN() seed   = 3 2 7 6 4 5 6 1

! Leitura da dimens˜ ao das matrizes quadradas w r i t e ( ∗ , ' (A) '  , advance= ' no  '  ) ' Q ua l e a d im en sa o d as m a t r iz e s ? read ( ∗ , ' ( I ) ' ) n

'

! Verifica se a dimens˜ ao ´e maior ou igual a 2. Caso contr´ ario, interrompe o programa. i f  ( n >= 2)   then a l l o c a t e ( a ( n ,  n ) , b ( n ,  n ) , c ( n ,  n ) ) else w r i t e ( ∗  , ∗ ) ' As m a t r iz e s devem p o s s u i r d im en sa o m ai or ou i g u a l a 2 ! stop e nd i f  

'

! Pede confirma¸ca ˜o para gerar as matrizes aleatoriamente w r i t e ( ∗ , ' (A) '  , advance= ' no  '  ) ' D e se j a q ue a s m a t r iz e s a e b s ej am g e r ad a s ? ( S o u N) read ( ∗ , ' (A1) ' )   resposta ! Verifica a resposta dada pelo usu´ ario i f  ( ( resposta == ' S  '  ) . or  . ( resposta == ' s ' ) )   then ! Usu´ ario prefere que as matrizes sejam geradas aleatoriamente c a l l   matriz_aleatoria  ( a , n )

5

'

c a l l   matriz_aleatoria  ( b , n ) else ! Usu´ ario prefere fornecer as matrizes em um arquivo w r i t e ( ∗ , ' (/ , x ,A) ' ) 'ATENCAO! Termine o arq uiv o com uma li nh a em branco . ' w r i t e ( ∗ , ' (x ,A) '  , advance= ' no  '  ) ' Nome d o a r q u i v o c o nt e nd o a s m a t r i z e s ( max . 2 0 cara cter es ) : ' read ( ∗ , ' (A20) ' )   nome_arquivo

←

! Abrindo o   nome arquivo  para leitura. !   status=‘old’  for¸ca que nome arquivo j´ a exista !   iostat=condicao  armazena informa¸c˜ oes sobre erros ao abrir o arquivo open ( u n i t=arquivo  , f i l e=nome_arquivo  ,   s t a t u s= ' o l d '  ,   a c t i o n= ' read '  ,   i o s t a t=condicao ) ! O nome arquivo abriu sem problemas? i f  ( condicao   == 0)   then ! Arquivo aberto ! Agora vamos ler os valores das matrizes a e b read ( arquivo  , ∗ ) ( ( a ( i ,  j ) , j = 1 , n ) , i read ( arquivo  , ∗ ) ( ( b ( i ,  j ) , j = 1 , n ) , i w r i t e ( ∗ , ' (/ , x ,A) ' ) ' M a t ri z es l i d a s com else ! Erro ao abrir o arquivo. Terminar o programa. w r i t e ( ∗ , ' (/ , x ,A) ' ) ' ∗ ∗   E rr o ao a b r i r o c l o s e ( arquivo ) stop e nd i f  

= 1 , n) = 1 , n) s u c e s so !

'

a r qu iv o ! V e r i f iq u e o a r qu iv o de e nt ra da . '

! Fechar o arquivo, pois j´ a temos os valores armazenados nas vari´ aveis c l o s e ( arquivo ) e nd i f   ! Menu para sele¸ ca ˜o das opera¸c˜ oes matriciais implementadas ! Este la¸co se repetir´ a at´e que o usu´ ario escolha a op¸c˜ a o 0. do w r i t e ( ∗ , ' (3 / ,x ,A) ' ) ' ( 1 ) Somar a s m a t r i ze s A e B ' w r i t e ( ∗  , ∗ ) ' ( 2) M u lt i pl i ca r a s m at ri ze s A e B ' w r i t e ( ∗  , ∗ ) ' ( 3 ) O bt er a t r a n s po s t a d e A ' w r i t e ( ∗  , ∗ ) ' ( 4 ) O bt er a i n v e rs a de A ' w r i t e ( ∗  , ∗ ) ' ( 5 ) C a l c ul a r o d e te r mi n an t e de A ' w r i t e ( ∗  , ∗ ) ' ( 6 ) E x ib i r a s m a tr i ze s A e B ' w r i t e ( ∗  , ∗ ) ' ( 0) S a ir ' w r i t e ( ∗ , ' (/ ,A) '  , advance= ' no  '  ) ' D i g i t e a o pc ao d e s ej a d a ( 0− 6) : read ( ∗ , ∗ )   opcao i f  ( opcao   == 1)   then ! C = A +B c a l l   matriz_soma ( a , b , c , n ) w r i t e ( ∗ , ' (/ , x ,A) ' ) ' A soma da m a tr i z A com B e i g u a l a c a l l   matriz_imprime ( c , n ) e l s e i f   ( opcao   == 2)   then ! C = A *B c a l l   matriz_mul ( a , b , c , n )

'

w r i t e ( ∗ , ' (/ , x ,A) ' ) 'O p r od ut o da m a tr i z A com B e i g u a l a c a l l   matriz_imprime ( c , n ) e l s e i f   ( opcao   == 3)   then ! Matriz transposta de A w r i t e ( ∗ , ' (/ , x ,A) ' ) ' A m at ri z t r a ns p os t a de A e i g u a l a ' w r i t e ( ∗ , ' (/ , x ,A) ' ) 'MATRIX A ' c a l l   matriz_imprime ( a , n ) c a l l   matriz_transposta  ( a , n ) w r i t e ( ∗ , ' (/ , x ,A) ' ) c a l l   matriz_imprime e l s e i f   ( opcao   == 4)   ! Matriz inversa de A w r i t e ( ∗ , ' (/ , x ,A) ' ) w r i t e ( ∗ , ' (/ , x ,A) ' ) c a l l   matriz_imprime

'MATRIX A ( tr an sp os ta ) ' ( a , n ) then 'A

m at ri z i n ve r sa de A e i g u al a A' ( a , n ) 'MATRIX

c a l l   matriz_inversa ( a , b , n ) w r i t e ( ∗ , ' (/ , x ,A) ' ) 'MATRIX B ( i nv e r sa ) ' c a l l   matriz_imprime ( b , n ) c a l l   matriz_inversa_verifica  ( a , b , n )

6

'

'

'

e l s e i f   ( opcao   == 5)   then ! Decomposi¸ c˜ a o de A em LU c a l l   matriz_lu ( a , n ) ! C´ alculo do determinante para uma matriz LU determinante   = 1 . 0 do i = 1 , n determinante =   determinante

∗

a ( i ,  i )

e nd d o w r i t e ( ∗ , ' (/ , x ,A) ' ) 'MATRIX A ' c a l l   matriz_imprime ( a , n ) w r i t e ( ∗ , ' (/ , x ,A,F) ' ) 'O d e t e r m i na n t e d a m a t r i z A v a l e : e l s e i f   ( opcao   == 6)   then ! Exibe os elementos das matrizes A e B w r i t e ( ∗ , ' (/ , x ,A) ' ) 'MATRIX A ' c a l l   matriz_imprime ( a , n ) w r i t e ( ∗ , ' (/ , x ,A) ' ) 'MATRIX B ' c a l l   matriz_imprime ( b , n ) e l s e i f   ( opcao   == 0)   then ! Interrompe o programa w r i t e ( ∗ , ' (/ ,x ,A,/ ) ' ) ' Programa terminado . Ate mais ! exit e nd i f   e nd d o

'

! Libera a mem´ oria ocupada pelas matrizes A, B, C d e a l l o c a t e ( a , b , c ) end program

! Esta subrotina exibe os valores dos elementos de uma matriz ! ! Parˆ ametros: ! a  Matriz cujos elementos ser˜ ao exibidos ! n  Dimens˜ ao da matriz (apenas matrizes quadradas) s u b r o u t i n e   matriz_imprime ( a , n ) i m p l i c i t no ne i n t e g e r : : n , i , j r e a l : : a ( n ,  n ) do i = 1 , n do j = 1 , n w r i t e ( ∗ , ' ( f , 2 x ) '  , advance= ' no  '  ) a ( i ,  j ) e nd d o print ∗ e nd d o e nd s u b r o u t i n e

! Esta subrotina gera uma matriz com valores aleat´ orios para seus elementos ! ! Parˆ ametros: ! a  Matriz que armazenar´ a os elementos gerados ! n  Dimens˜ ao da matriz (apenas matrizes quadradas) s u b r o u t i n e   matriz_aleatoria  ( a , n ) use   global i m p l i c i t no ne i n t e g e r : : n , i , j r e a l : : a ( n ,  n ) do i = 1 , n do j = 1 , n ! Calcula valores aleat´ orios de -5 `a +5 ! A fun¸ca ˜o ran() retorna valores entre 0.0 e 1.0 a ( i ,  j ) = ( − 1.0) ∗ ∗ ( i+j ) ∗ 5 . 0 ∗ ra n ( seed ) e nd d o e nd d o e nd s u b r o u t i n e

! Esta subrotina faz a adi¸ ca ˜o de duas matrizes quadradas (a + b) ! ! Parˆ ametros: ! a  e b  Matrizes de entrada

7

'  ,

  determinante

! c  Matriz que armazenar´ a a soma de a  com  b ! n  Dimens˜ ao da matriz (apenas matrizes quadradas) s u b r o u t i n e   matriz_soma ( a , b , c , n ) i m p l i c i t no ne i n t e g e r : : n , i , j r e a l : : a ( n ,  n ) , b ( n ,  n ) , c ( n ,  n ) do i = 1 , n do j = 1 , n c ( i ,  j ) = a ( i ,  j ) + b ( i ,  j ) e nd d o e nd d o e nd s u b r o u t i n e

! Esta subrotina faz a multiplica¸ c˜ ao de duas matrizes quadradas (a * b) ! ! Parˆ ametros: ! a  e b  Matrizes de entrada ! c  Matriz que armazenar´ a o produto de a  por  b ! n  Dimens˜ ao da matriz (apenas matrizes quadradas) s u b r o u t i n e   matriz_mul ( a , b , c , n ) i m p l i c i t no ne i n t e g e r : : n , i , j , k r e a l : : a ( n ,  n ) , b ( n ,  n ) , c ( n ,  n ) c  = 0.0

do i = 1 , n do j = 1 , n do k = 1 , n c ( i ,  j ) = c ( i ,  j ) + a ( i ,  k ) e nd d o e nd d o e nd d o e nd s u b r o u t i n e

∗

b ( k ,  j )

! Esta subrotina cria a transposta de uma matriz ! ! Parˆ ametros: ! a  Matriz de entrada ! n  Dimens˜ ao da matriz (apenas matrizes quadradas) ! !   AVISO: A matriz a  ser´ a modificada no processo s u b r o u t i n e   matriz_transposta  ( a , n ) i m p l i c i t no ne i n t e g e r : : n , i , j r e a l : : a ( n ,  n ) ,   temp do i = 1 , n do j = 1 , n i f  ( i < j )   then temp = a ( i ,  j ) a ( i ,  j ) = a ( j ,  i ) a ( j ,  i ) =   temp e nd i f   e nd d o e nd d o e nd s u b r o u t i n e

! Esta subrotina calcula a inversa de uma matriz por elimina¸ c˜ ao gaussiana ! ! Parˆ ametros: !   matriz  Matriz de entrada !   inversa  Matriz que armazenar´ a a inversa ! n  Dimens˜ ao da matriz (apenas matrizes quadradas) ! !   AVISO: A matriz a  ser´ a modificada no processo ! Referˆencia : ! http://math.uww.edu/˜mcfarlat/inverse.htm ! http://www.tutor.ms.unimelb.edu.au/matrix/matrix inverse.html s u b r o u t i n e   matriz_inversa ( ma tr iz  ,   inversa  , n ) i m p l i c i t no ne integer : : n

8

r e a l ,   dimension ( n ,  n ) : :   matriz r e a l ,   dimension ( n ,  n ) : :   inversa l o g i c a l : :   invertivel = . true . i n t e g e r : : i , j , k , l real :: m r e a l ,   dimension ( n , 2 ∗ n ) : :   matriz_aumentada ! Matriz aumentada com uma matriz identidade do i = 1 , n do j   = 1 , 2 ∗ n i f  ( j <= n )   then  m at ri z_ au me nt ad a  ( i ,  j ) =   matriz ( i ,  j ) e l s e i f   ( ( i+n ) == j )   then  m at ri z_ au me nt ad a  ( i ,  j ) = 1 else  m at ri z_ au me nt ad a  ( i ,  j ) = 0 endif  e nd d o e nd d o ! Reduzir a matriz aumentada a uma matriz triangular superior pela elimina¸ c˜ ao gaussiana do k = 1 , n − 1 ! Verifica se algum elemento da diagonal ´ e zero i f  ( ab s (  m at ri z_ au m en ta da  ( k ,  k ) ) <= 1 . 0 E −6)   then invertivel = .  false  . do i = k + 1 , n ! Verifica se os elementos s˜ ao maiores que zero i f  ( ab s (  m at ri z_ au m en ta da  ( i ,  k ) ) > 1 . 0 E −6)   then do j   = 1 , 2 ∗ n  m at ri z_ au m en ta da  ( k ,  j ) =   matriz_aumentada ( k ,  j ) +   matriz_aumentada  ( i ,  j ) e nd d o invertivel = .  true  . exit endif  ! Se algum elemento da diagonal for zero, n˜ ao podemos calcular a inversa i f  ( invertivel == . false  . )   then w r i t e ( ∗ , ' (/ ,x ,A,/ ) ' ) ' ∗ ∗   A m at ri z nao e i n v e r s t i v e l ! inversa = 0 stop endif  e nd d o endif 

∗∗ '

! Elimina¸ ca ˜o gaussiana do j = k + 1 , n  m =   matriz_aumentada  ( j ,  k ) /   matriz_aumentada  ( k ,  k ) do i = k , 2

∗ n  m at ri z_ au me nt ad a  ( j ,  i ) =   matriz_aumentada  ( j ,  i )

−

m

∗

  matriz_aumentada  ( k ,  i )

e nd d o e nd d o e nd d o ! Teste para invertibilidade do i = 1 , n ! Elementos da diagonal n˜ ao podem ser zero i f  ( ab s (  m at ri z_ au m en ta da  ( i ,  i ) ) <= 1 . 0 E −6)   then w r i t e ( ∗ , ' (/ ,x ,A,/ ) ' ) ' ∗ ∗   A m at ri z nao e i n v e r s t i v e l ! inversa = 0 stop endif  e nd d o ! Elementos da diagonal iguais a 1 do i = 1 , n  m =   matriz_aumentada  ( i ,  i ) do j = i , 2 ∗ n  m at ri z_ au me nt ad a  ( i ,  j ) =   matriz_aumentada  ( i ,  j ) / m e nd d o e nd d o ! Reduzir o lado esquerdo da matriz aumentada a matriz identidade do k = n − 1 , 1 , −1 do i = 1 , k

9

∗∗ '

 m =   matriz_aumentada  ( i ,  k +1) do j = k , 2 ∗ n  m at ri z_ au me nt ad a  ( i ,  j ) =   matriz_aumentada  ( i ,  j )

−

  matriz_aumentada ( k +1 , j )

∗

m

e nd d o e nd d o e nd d o ! Armazene o resultado do i =1 , n do j = 1 , n inversa ( i ,  j ) =   matriz_aumentada ( i , j + n ) e nd d o e nd d o e nd s u b r o u t i n e

! Esta subrotina verifica se a matriz inversa ´e v´ alida (a * b = identidade) ! ! Parˆ ametros: ! a  Matriz original ! b  Matriz inversa de a ! n  Dimens˜ ao da matriz (apenas matrizes quadradas) s u b r o u t i n e   matriz_inversa_verifica  ( a , b , n ) i m p l i c i t no ne i n t e g e r : : n , i , j r e a l ,   dimension ( n ,  n ) : : a , b ,   identidade l o g i c a l : : ok ok = .  true  .

c a l l   matriz_mul ( a , b ,   identidade  , n ) do i = 1 , n do j = 1 , n ! Elementos fora da diagonal principal n˜ ao podem ser diferentes de zero i f  ( ( i /= j ) . a nd  . ( ab s ( identidade ( i ,  j ) ) >= 1 . 0 E −6) )   then ok = .  false  . ! Elementos na diagonal principal n˜ ao podem ser diferentes de um e l s e i f   ( ( i == j ) . a nd  . ( ab s ( identidade ( i ,  i ) −   1 . 0 ) >= 1 . 0 E −6) )   then ok = .  false  . e nd i f   e nd d o e nd d o i f  ( ok == . false  . )   then w r i t e ( ∗ , ' (/ ,x ,A,/ ) ' ) ' A v e r f i c a c a o da m at ri z i n v er s a f a lh o u ' else w r i t e ( ∗ , ' (/ ,x ,A,/ ) ' ) ' A v e r f i c a c a o da m at ri z i n v er s a f o i c o nc l ui d a c om s u ce s so ! e nd i f   e nd s u b r o u t i n e

! Esta subrotina decomp˜ oe uma matriz em matrizes triangulares inferior e superior (LU) ! pelo algoritmo de Doolittle (diagonal principal da matriz L igual a 1). Para uma ! descri¸c˜ ao te´ orica sobre o algoritmo de Doolittle, veja [2]. ! ! Parˆ ametros: ! a  Matriz que armazenar´ a as matrizes LU ! n  Dimens˜ ao da matriz (apenas matrizes quadradas) ! !   AVISO: A matriz a  ser´ a modificada no processo s u b r o u t i n e   matriz_lu ( a , n ) i m p l i c i t no ne i n t e g e r : : n , i , j , k r e a l : : a ( n ,  n ) , s , ss ! Os elementos da diagonal n˜ ao podem ser zero i f  ( ab s ( a (1 ,1) ) <= 1 . 0 E −6)   then w r i t e ( ∗ , ' (/ ,x ,A,/ ) ' ) ' ∗ ∗   I m po s si v e l f a t o r a r a m at ri z stop e nd i f   do i = 2 , n a (i , 1 ) = a(i , 1) / a (1 ,1)

10

∗∗ '

'

e nd d o do i = 2 , n

−

1

s  = 0.0 do k = 1 , i − 1 s = s − a ( i ,  k ) e nd d o

∗

a ( k ,  i )

a ( i ,  i ) = ( a ( i ,  i ) + s )

! Os elementos da diagonal n˜ ao podem ser zero i f  ( ab s ( a ( i ,  i ) ) <= 1 . 0 E −6)   then w r i t e ( ∗ , ' (/ ,x ,A,/ ) ' ) ' i m p o ss i v el f a t o r a r a m at ri z ' stop e nd i f   do j = i + 1 , n ss   = 0 . 0 s  = 0.0 do k = 1 , i − 1 ss = ss − a ( i ,  k ) ∗ a ( k ,  j ) s = s − a ( j ,  k ) ∗ a ( k ,  i ) e nd d o a ( i ,  j ) = a ( i ,  j ) + ss a ( j ,  i ) = ( a ( j ,  i ) + s ) / a ( i ,  i )

e nd d o e nd d o s  = 0.0

do k = 1 , n − 1 s = s − a ( n ,  k ) e nd d o

∗

a ( k ,  n )

a ( n ,  n ) = a ( n ,  n ) + s

e nd s u b r o u t i n e

Referˆ encias [1] Stephen J. Chapman, Fortran 95/2003 for scientists and engineers , McGrall-Hill Companies, New York, NY, 2007. [2] R.L. Burden e J.D. Faires, Numerical analysis , Brooks Cole, California, CA, 2000. ´  [3] J.L. Boldrini et al., Algebra linear , Editora Harbra, Ltda., S˜ ao Paulo, SP, 1980.

11