Fortran: LAPACK

<< zur Fortran-Startseite
< FGSL	OpenMP >

Wikipedia hat einen Artikel zum Thema:

LAPACK

Allgemeines[Bearbeiten]

LAPACK steht für "Linear Algebra Package". LAPACK ist eine Bibliothek zwecks Lösung von

linearen Gleichungssystemen
LLS-Aufgaben
Eigenwertproblemen
Singulärwertproblemen

LAPACK ist in FORTRAN 77 geschrieben und nutzt weitgehend Funktionen von BLAS Level 3. Falls keine für einen bestimmten Prozessor optimierte Version von BLAS, z.B. ATLAS, bereits installiert ist, kann LAPACK mit der eigenen BLAS-Implementierung kompiliert werden. Aus FORTRAN 77 resultierende Namensbeschränkung auf eine maximale Länge von 6 Zeichen führt zu sehr kryptischen Unterprogrammbezeichnungen, z.B.

	D	G	E	T	R	F
Zeichenposition:	D1	M1	M2	O1	O2	O3

Erläuterung der Zeichenpositionen:

D1

Datentyp:

S	...	real
D	...	double precision
C	...	complex
Z	...	double complex

M1, M2

Matrixtyp, z.B.

GE	...	generell
DI	...	diagonal
OR	...	orthogonal (reelle Zahlen)
SY	...	symmetrisch

O1, O2, (O3)

Operation, z.B.

TRF

...

faktorisiere

Eine detailliertere und umfassendere Beschreibung des Funktionsumfanges und der Anwendungsmöglichkeiten der LAPACK-Bibliothek bietet der LAPACK Users' Guide. Die einzelnen Subroutinen inklusive Unterprogrammparameter sind zudem auch in den LAPACK-Sourcecode-Dateien ausführlich dokumentiert.

Beispiel: Lösen eines einfachen Gleichungssystems[Bearbeiten]

Gegeben ist folgendes Gleichungssystem:

$2x+y\ =\ 5$

$3x+y\ =\ 6$

bzw. in Matrixschreibweise:

${\begin{bmatrix}2&1\\3&1\\\end{bmatrix}}{\begin{bmatrix}x\\y\\\end{bmatrix}}={\begin{bmatrix}5\\6\\\end{bmatrix}}\Leftrightarrow A\mathbf {x} =B$

Gesucht sind die Unbekannten x und y:

Fortran 90/95-Code (free source form)

program bsp
  implicit none
 
  integer :: info
  real, dimension(2,2) :: a = reshape ( (/2.,3.,1.,1./), (/2,2/) )
  real, dimension(2)   :: b = (/5., 6./), ipiv
 
!  SUBROUTINE SGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
!  s ... real, ge ... general matrix type, sv ... solver
   call sgesv(2, 1, a, 2, ipiv, b, 2, info)
 
   write(*,*) "Lösung (x, y): ", b
 
   if(info == 0) then
     write(*,*) "Ergebnis OK"
   else
     write(*,*) "Ergebnis NOK"
   end if   
 
!  Ausgabe: Lösung (x, y):     1.000000       3.000000 
!  Ausgabe: Ergebnis OK  
end program bsp

Kompilieren, Linken:

gfortran bsp.f95 -llapack -lblas

Beispiel: Inverse Matrix[Bearbeiten]

Gegeben ist eine 3x3-Matrix

$A={\begin{pmatrix}3&-2&1\\-3&5&0\\2&-1&2\end{pmatrix}}$

die invertiert werden soll. Die Zahlenwerte dieser Matrix A entsprechen einem Beispiel aus Bartsch: Mathematische Formeln, 21. Auflage, VEB Fachbuchverlag Leipzig, 1986, Seite 109, ebenfalls zum Thema "Inverse Matrix".

Verwendet werden hierzu die beiden LAPACK-Subroutinen:

SGETRF
- S ... Datentyp: real
- GE ... Matrixtyp: general
- TRF ... Operation: LU-Faktorisierung (Dreiecksform)
SGETRI
- S ... Datentyp: real
- GE ... Matrixtyp: general
- TRI ... Operation: Invertierung einer LU-faktorisierten Matrix

Fortran 90/95-Code (free source form)

program bsp
  implicit none 

  real, dimension( 3, 3 ) :: A
  integer, dimension( 3 ) :: ipiv
  real, dimension( 3 )    :: work 
  integer                 :: m = 3, n = 3, lda = 3, lwork = 3, info 
  
				
  A =  reshape( (/  3.0, -3.0,  2.0, -2.0,  5.0, -1.0, 1.0,  0.0,  2.0  /),  &
                shape( A ) )
  		
! LU-Faktorisierung (Dreieckszerlegung) der Matrix A		
  call sgetrf( m, n, A, lda, ipiv, info )		

! Inverse der LU-faktorisierten Matrix A	
  call sgetri( n, A, lda, ipiv, work, lwork, info )
  
  write( *, * ) "Inverse Matrix Ai =", A
  write( *, * ) "Testweise wie im Bartsch-Beispiel, Ai = 1/11 * (", A * 11, ")"
  
  if( info == 0 ) then
    write( *, * ) "OK"
  else
    write( *, * ) "Nicht OK"  
  end if  

! Ausgabe:
!   Inverse Matrix Ai = 0.909091 0.54545456 -0.63636374 0.2727273    
!     0.36363637 -0.090909116 -0.45454553 -0.2727273 0.81818193
!   Testweise wie im Bartsch-Beispiel, Ai = 1/11 * ( 10.000001 6. -7.000001  
!     3.0000005 4. -1.0000002 -5.000001 -3.0000005 9.000001 )
!   OK

end program bsp

Weblinks[Bearbeiten]

<< zur Fortran-Startseite
< FGSL	OpenMP >