09/11/11
http://jean-pierre.moreau.pagesperso-orange.fr/f_matrices.html
(1)!
! LU decomposition routines used by test_lu.f90
!
! F90 version by J-P Moreau, Paris
! ---------------------------------------------------
! Reference:
!
! "Numerical Recipes By W.H. Press, B. P. Flannery,
! S.A. Teukolsky and W.T. Vetterling, Cambridge
! University Press, 1986" BIBLI 08,.
!
!
MJDULE LU
CJNTAINS
!
! Given an N x N matrix A, this routine replaces it by the LU
! decomposition of a rowwise permutation of itself. A and N
! are input. INDX is an output vector which records the row
! permutation effected by the partial pivoting; D is output
! as -1 or 1, depending on whether the number of row inter-
! changes was even or odd, respectively. This routine is used
! in combination with LUBKSB to solve linear equations or to
! invert a matrix. Return code is 1, if matrix is singular.
!
Subroutine LUDCMP(A,N,INDX,D,CJDE)
PARAMETER(NMAX=100,TINY=1.5D-16)
REAL8 AMAX,DUM, SUM, A(N,N),VV(NMAX)
INTEGER CJDE, D, INDX(N)
D=1; CJDE=0
DJ I=1,N
AMAX=0.d0
DJ J=1,N
IF (DABS(A(I,J)).GT.AMAX) AMAX=DABS(A(I,J))
END DJ ! j loop
IF(AMAX.LT.TINY) THEN
CJDE = 1
RETURN
END IF
VV(I) = 1.d0 / AMAX
END DJ ! i loop
DJ J=1,N
DJ I=1,J-1
SUM = A(I,J)
DJ K=1,I-1
SUM = SUM - A(I,K)A(K,J)
END DJ ! k loop
A(I,J) = SUM
END DJ ! i loop
AMAX = 0.d0
DJ I=J,N
SUM = A(I,J)
DJ K=1,J-1
SUM = SUM - A(I,K)A(K,J)
END DJ ! k loop
A(I,J) = SUM
DUM = VV(I)DABS(SUM)
IF(DUM.GE.AMAX) THEN
IMAX = I
AMAX = DUM
END IF
END DJ ! i loop
IF(J.NE.IMAX) THEN
DJ K=1,N
DUM = A(IMAX,K)
A(IMAX,K) = A(J,K)
A(J,K) = DUM
END DJ ! k loop
D = -D
VV(IMAX) = VV(J)
END IF
INDX(J) = IMAX
IF(DABS(A(J,J)) < TINY) A(J,J) = TINY
IF(J.NE.N) THEN
DUM = 1.d0 / A(J,J)
DJ I=J+1,N
A(I,J) = A(I,J)DUM
END DJ ! i loop
END IF
END DJ ! j loop
RETURN
END subroutine LUDCMP
!
! Solves the set of N linear equations A . X = B. Here A is
! input, not as the matrix A but rather as its LU decomposition,
! determined by the routine LUDCMP. INDX is input as the permuta-
! tion vector returned by LUDCMP. B is input as the right-hand
! side vector B, and returns with the solution vector X. A, N and
! INDX are not modified by this routine and can be used for suc-
! cessive calls with different right-hand sides. This routine is
! also efficient for plain matrix inversion.
!
Subroutine LUBKSB(A,N,INDX,B)
REAL8 SUM, A(N,N),B(N)
INTEGER INDX(N)
II = 0
DJ I=1,N
LL = INDX(I)
SUM = B(LL)
B(LL) = B(I)
IF(II.NE.0) THEN
DJ J=II,I-1
SUM = SUM - A(I,J)B(J)
END DJ ! j loop
ELSE IF(SUM.NE.0.d0) THEN
II = I
END IF
B(I) = SUM
END DJ ! i loop
DJ I=N,1,-1
SUM = B(I)
IF(I < N) THEN
DJ J=I+1,N
SUM = SUM - A(I,J)B(J)
END DJ ! j loop
END IF
B(I) = SUM / A(I,I)
END DJ ! i loop
RETURN
END subroutine LUBKSB
END MJDULE LU
! end of file lu.f90
!
! Inversion of a symmetric matrix by Cholesky decomposition.
! The matrix must be positive definite.
! --------------------------------------------------------------
! REFERENCE:
! From a Java Library Created by Vadim Kutsyy,
! "http://www.kutsyy.com".
! --------------------------------------------------------------
! SAMPLE RUN:
!
! Inversion of a square real symetric matrix by Cholevsky method
! (The matrix must positive definite).
!
! Size = 4
!
! Matrix A:
! 5.000000 -1.000000 -1.000000 -1.000000
! -1.000000 5.000000 -1.000000 -1.000000
! -1.000000 -1.000000 5.000000 -1.000000
! -1.000000 -1.000000 -1.000000 5.000000
!
! Determinant = 432.000000
!
! Matrix Inv(A):
! 0.250000 0.083333 0.083333 0.083333
! 0.083333 0.250000 0.083333 0.083333
! 0.083333 0.083333 0.250000 0.083333
! 0.083333 0.083333 0.083333 0.250000
!
! F90 Release By Jean-Pierre Moreau, Paris.
! --------------------------------------------------------------
! Release 1.1 : added verification Inv(A) A = I.
!
Program Cholesky
real8, pointer :: A(:,:), A1(:,:), B(:,:), C(:,:)
integer i,j, n, istat, Check_Matrix
character1 answer
n=4
allocate(A(0:n-1,0:n-1),B(0:n-1,0:n-1),stat=istat)
allocate(A1(0:n-1,0:n-1),C(0:n-1,0:n-1),stat=istat)
print ,' Inversion of a square real symetric matrix by Cholesky method'
print ,' (The matrix must positive def.).'
print ,' '
print ,' Size = ', n
! define lower half of matrix
A(0,0)= 5;
A(1,0)=-1; A(1,1)=5
A(2,0)=-1; A(2,1)=-1; A(2,2)= 5
A(3,0)=-1; A(3,1)=-1; A(3,2)=-1; A(3,3)= 5
! define upper half by symmetry
do i=0, n-1
do j=i+1, n-1
A(i,j)=A(j,i)
end do
end do
call MatPrint('Matrix A:',n,A)
print ,' '
if (Check_Matrix(n,A).eq.1) then
A1 = A
print ,' Determinant = ', choldet(n,A)
call cholsl(n,A,B)
call MatPrint('Matrix Inv(A):',n,B)
else
print ,' This matrix is not positive definite !'
end if
print ,' '
write(,10,advance='no'); read , answer
if (answer.eq.'y') then
call MatMult(n,A1,B,C)
call MatPrint('Verification A Inv(A) = I:',n,C)
end if
print ,' '
deallocate(A,A1,B,C)
stop
10 format(' Do you want a verification (y/n) . ')
END
! ------------------------------------------------
! Cholesky decomposition.
! input n size of matrix
! input A Symmetric positive def. matrix
! output aa lower deomposed matrix
! uses choldc1(int,MAT,VEC)
! ------------------------------------------------
Subroutine choldc(n,A,aa)
integer n
real8 A(0:n-1,0:n-1), aa(0:n-1,0:n-1)
integer i,j, ialloc
real8, pointer :: p(:)
allocate(p(0:n-1),stat=ialloc)
aa = A
call choldc1(n, aa, p)
do i = 0, n-1
aa(i,i) = p(i)
do j = i + 1, n-1
aa(i,j) = 0.d0
end do
end do
deallocate(p)
return
End
! -----------------------------------------------------
! Inverse of Cholesky decomposition.
! input n size of matrix
! input A Symmetric positive def. matrix
! output aa inverse of lower decomposed matrix
! uses choldc1(int,MAT,VEC)
! -----------------------------------------------------
Subroutine choldcsl(n,A,aa)
integer n
real8 A(0:n-1,0:n-1), aa(0:n-1,0:n-1)
integer i,j,k, ialloc
real8 sum
real8, pointer :: p(:)
allocate(p(0:n-1),stat=ialloc)
aa = A
call choldc1(n, aa, p)
do i = 0, n-1
aa(i,i) = 1.d0 / p(i)
do j = i + 1, n-1
sum = 0.d0
do k = i, j-1
sum = sum - aa(j,k) aa(k,i)
end do
aa(j,i) = sum / p(j)
end do
end do
deallocate(p)
return
End
! ----------------------------------------------------------------------
! Computation of Determinant of the matrix using Cholesky decomposition
! input n size of matrix
! input a Symmetric positive def. matrix
! return det(a)
! uses choldc(int,MAT,MAT)
! ----------------------------------------------------------------------
real8 Function choldet(n,a)
integer n
real8 a(0:n-1,0:n-1)
real8, pointer :: c(:,:)
real8 d
integer i, ialloc
allocate(c(0:n-1,0:n-1),stat=ialloc)
d=1.d0
call choldc(n,a,c)
do i = 0, n-1
d = d c(i,i)
end do
choldet = d d
deallocate(c)
return
End
! ---------------------------------------------------
! Matrix inverse using Cholesky decomposition
! input n size of matrix
! input A Symmetric positive def. matrix
! output aa inverse of A
! uses choldc1(int,MAT,VEC)
! ---------------------------------------------------
Subroutine cholsl(n,A,aa)
integer n
real8 A(0:n-1,0:n-1), aa(0:n-1,0:n-1)
integer i,j,k
call choldcsl(n,A,aa)
do i = 0, n-1
do j = i + 1, n-1
aa(i,j) = 0.d0
end do
end do
do i = 0, n-1
aa(i,i) = aa(i,i) aa(i,i)
do k = i + 1, n-1
aa(i,i) = aa(i,i) + aa(k,i) aa(k,i)
end do
do j = i + 1, n-1
do k = j, n-1
aa(i,j) = aa(i,j) + aa(k,i) aa(k,j)
end do
end do
end do
do i = 0, n-1
do j = 0, i-1
aa(i,j) = aa(j,i)
end do
end do
return
End
! -------------------------------------------------
! main method for Cholesky decomposition.
!
! input n size of matrix
! input/output a matrix
! output p vector of resulting diag of a
! author: <Vadum Kutsyy, kutsyy@hotmail.com
! -------------------------------------------------
Subroutine choldc1(n,a,p)
integer n
real8 a(0:n-1,0:n-1), p(0:n-1)
integer i,j,k
real8 sum
do i = 0, n-1
do j = i, n-1
sum = a(i,j)
do k = i - 1, 0, -1
sum = sum - a(i,k) a(j,k)
end do
if (i.eq.j) then
if (sum <= 0.d0) &
print ,' the matrix is not positive definite!'
p(i) = dsqrt(sum)
else
a(j,i) = sum / p(i)
end if
end do
end do
return
End
! print a square real matrix A of size n with caption s
! (n items per line).
Subroutine MatPrint(s,n,A)
character() s
integer n
real8 A(0:n-1,0:n-1)
integer i,j
print ,' '
print ,' ', s
do i=0, n-1
write(,10) (A(i,j),j=0,n-1)
end do
return
10 format(10F10.6)
End
Integer Function Check_Matrix(n,A)
integer n
real8 A(0:n-1,0:n-1)
integer i,j,k
real8 sum
Check_Matrix=1
do i = 0, n-1
do j = i, n-1
sum = A(i,j)
do k = i - 1, 0, -1
sum = sum - a(i,k) a(j,k)
end do
if (i.eq.j) then
if (sum <= 0.d0) Check_Matrix=0
end if
end do
end do
return
End
!
! MULTIPLICATIJN JF TWJ SQUARE REAL
! MATRICES
! ---------------------------------------
! INPUTS: A MATRIX NN
! B MATRIX NN
! N INTEGER
! ---------------------------------------
! JUTPUTS: C MATRIX NN PRJDUCT AB
!
!
Subroutine MATMULT(n,A,B,C)
integer n
real8 A(0:n-1,0:n-1), B(0:n-1,0:n-1), C(0:n-1,0:n-1)
real8 SUM
integer I,J,K
do I=0, n-1
do J=0, n-1
SUM= 0.d0
do K=0, n-1
SUM=SUM+A(I,K)B(K,J)
end do
C(I,J)=SUM
end do
end do
return
End
!end of file choles.f90
%%http]]wwwsa||hnetfreeserverscom]eng|neer|ng]fortran_codes]power_methodhtm|
POWER METHOD FOR EIGENVALUES
PJWER METHJD
PRJGRAM FJR DETERMINATIJN JF HIGHEST AND LJWEST EIGENVALUES
(SPECTRAL RADIUS)JF A SQUARE MATRIX
program PJWER
implicit doubleprecision(a-h,o-z)
parameter(nd=40)
doubleprecision AM(nd,nd),SM(nd),CM(nd),EM(nd),HM(nd),DM(nd)
integer count
parameter (eps = 1d-6)
print, 'Enter the order of square matrix:'
read(,) n
print, 'Enter the matrix row-wise:'
read(,) ((AM(i,j), j = 1,n), i = 1,n)
do i = 1, n
SM(i) = 1
enddo
DM = SM ! copy SM to DM
count = 1
maxim = 200
error = 1d9
epsil = 1d-5
do i = 1, n
HM(i) = 0
enddo
10 if ((error epsil).and.(count < maxim)) then
CALL MATMUL(AM,SM,CM,n,n,1)
b = BIG(CM,n)
do i = 1,n
if (b.ne.0) then
SM(i) = CM(i) /b
endif
EM(i) = SM(i) - HM(i)
enddo
error = dabs(BIG(EM,n))
HM = SM ! copy SM to HM
count = count + 1
go to 10
endif
if (count = maxim) then
write(6,) 'Tolerance is not met with 200 iterations:'
endif
write(6,)'-------------------------------------------------------
$----------------------'
write(6,101) b
101 format('',' The dominant Eigenvalue (spectral radius) is:',
$3x,f12.4)
write(6,)'-------------------------------------------------------
$----------------------'
write(6,) ' The Eigenvector corresponding to dominant Eige
$nvalue is:'
write(6,)
write(6,102) (SM(i), i = 1,n)
102 format('',29x,f11.4)
write(6,)'-------------------------------------------------------
$----------------------'
sum = 0
do i = 1, n
sum = sum + AM(i,i)
enddo
do i = 1, n
AM(i,i) = AM(i,i) - b
enddo
if (dabs(sum-b) = eps) then
count = 1
maxim = 250
error = 1d9
epsil = 1d-5
do i = 1, n
HM(i) = 0
enddo
20 if ((error epsil).and.(count < maxim)) then
CALL MATMUL(AM,DM,CM,n,n,1)
bb = BIG(CM,n)
do i = 1, n
if (bb.ne.0d0) then
DM(i) = CM(i) /bb
endif
EM(i) = DM(i) - HM(i)
enddo
error = dabs(BIG(EM,n))
HM = DM ! copy DM to HM
count = count + 1
go to 20
endif
if (count = maxim) then
write(6,) 'Tolerance is not met with 250 iterations:'
endif
bb = bb + b
write(6,)'--------------------------------------------------
$---------------------'
write(6,103) bb
103 format('',' The other extreme Eigenvalue is:',4x,f12.4)
write(6,)'--------------------------------------------------
$---------------------'
write(6,) ' The Eigenvector corresponding to this Eigen
$value is:'
write(6,)
write(6,102) (DM(i), i = 1,n)
elseif (dabs(sum-b) < eps) then
write(6,)'-----------------------------------------------------
$---------------------'
write(6,103) 0
write(6,)'-----------------------------------------------------
$------------------------'
write(6,) ' The Eigenvector corresponding to this Eigen
$value is:'
write(6,)
write(6,102) (DM(i), i = 1,n)
endif
write(6,)'-------------------------------------------------------
$-------------------'
stop
end
FUNCTIJN FJR SELECTING THE LARGEST ABSJLUTE VALUE - (1D ARRAY)
function BIG(A,n)
implicit doubleprecision(a-h,o-z)
parameter(nd=40)
doubleprecision A(nd),BIG
BIG = A(1)
do i = 2, n
if (dabs(BIG) < dabs(A(i))) then
BIG = A(i)
endif
enddo
return
end
MATRIX MULTIPLICATIJN
subroutine MATMUL(A,B,PRJD,m,l,n)
implicit doubleprecision(a-h,o-z)
parameter(nd=40)
doubleprecision A(nd,nd),B(nd,nd),PRJD(nd,nd)
do i = 1,m
do j = 1,n
PRJD(i,j) = 0
do index = 1,l
PRJD(i,j) = PRJD(i,j) + A(i,index) B(index,j)
enddo
enddo
enddo
return
end
