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