Consider the linear system of equations Ax=b with the block structure:
We know that the LU decomposition of A is:
Make a program in fortran to compute the LU decomposition and to find the solution.
_________________________________________________________________
Program ludecomposition
integer INFO1,INFO2, i, j, N, N1,N2,N3
integer, ALLOCATABLE, dimension (:) :: IPIVZ1,IPIVZ2
double precision, ALLOCATABLE, DIMENSION(:,:) :: A11,A13,A22,
c A23,A31,A32,A33,S,Z3,Z4,A11C,A22C
double precision, ALLOCATABLE, DIMENSION(:) ::B1,B2,B3
c we let the user to enter the dimensions N1,N2 and N3.
open (unit=1, file="solution.dat")
write(*,*) 'Input the 3 values n1,n2,n3.After each press enter.'
read (*,*) N1,N2,N3
c The dimension of A is N. Let's compute N.
N=N1+N2+N3
C A is a square matrix. That's why we have to compute M=N.
ALLOCATE(A11(N1,N1),A13(N1,N3),A22(N2,N2),A23(N2,N3),
c A31(N3,N1),A32(N3,N2),A33(N3,N3),S(N3,N3),Z3(N3,N1),
c Z4(N3,N2),A11C(N1,N1),A22C(N2,N2))
ALLOCATE(B1(N1),B2(N2),B3(N3),IPIVZ1(N1),IPIVZ2(N2))
C Let's full the matrixes AII with random numbers.
do i=1,N1
do j=1,N1
A11(i,j)=rand()
enddo
do j=1,N3
A13(i,j)=rand()
enddo
enddo
do i=1,N2
do j=1,N2
A22(i,j)=rand()
enddo
do j=1,N3
A23(i,j)=rand()
enddo
enddo
do i=1,N3
do j=1,N1
A31(i,j)=rand()
enddo
do j=1,N2
A32(i,j)=rand()
enddo
do j=1,N3
A33(i,j)=rand()
enddo
enddo
C Now the matrixs AII are full of random numbers. In order to see if
C it works correcly, we can full the matrix b with the elements of
C each row added.
do i=1,N1
B1(i)=0
do j=1,N1
B1(i)=B1(i)+A11(i,j)
enddo
do j=1,N3
B1(i)=B1(i)+A13(i,j)
enddo
enddo
do i=1,N2
B2(i)=0
do j=1,N2
B2(i)=B2(i)+A22(i,j)
enddo
do j=1,N3
B2(i)=B2(i)+A23(i,j)
enddo
enddo
do i=1,N3
B3(i)=0
do j=1,N1
B3(i)=B3(i)+A31(i,j)
enddo
do j=1,N2
B3(i)=B3(i)+A32(i,j)
enddo
do j=1,N3
B3(i)=B3(i)+A33(i,j)
enddo
enddo
C now, we call our subroutine in order to compute the LU decomposition
C of the matrix A composed by our AII matrixes
call ludec(A11,A13,A22,A23,A31,A32,A33,
C N1,N2,N3,S,Z3,Z4,A11C,A22C,IPIVZ1,IPIVZ2)
C finally, we call our subrouine in order to get the values of x
call solve(A11,A13,A22,A23,A31,A32,A33,
C N1,N2,N3,S,Z3,Z4,B1,B2,B3,A11C,A22C,IPIVZ1,IPIVZ2)
C let's write the result of the computation. If we write B, we are
C writting the solution of the system.
write(1,*)'The solution of the x-array you wanted to know is'
do i=1,N1
write(1,*)B1(i)
enddo
do i=1,N2
write(1,*)B2(i)
enddo
do i=1,N3
write(1,*)B3(i)
enddo
write(*,*)'The solution can be found in the file solution.dat.'
C in our example, x-array should be a vector of ones.
stop
end
C==================================================
subroutine ludec(A11,A13,A22,A23,A31,A32,A33,
C N1,N2,N3,S,Z3,Z4,A11C,A22C,IPIVZ1,IPIVZ2)
integer INFO,i,j,k
CHARACTER(1) TRANSN,TRANST
integer IPIVZ1(N1),IPIVZ2(N2)
double precision A11,A13,A22,A23,A31,A32,A33,
C A31Z1,A32Z2,S,A31T,A32T,Z3,Z4,Z1,Z2,A11C,A22C
dimension A11(N1,N1),A13(N1,N3),A22(N2,N2),A23(N2,N3),
C A31(N3,N1),A32(N3,N2),A33(N3,N3),A(N1+N2+N3,N1+N2+N3),
C A31Z1(N3,N3),A32Z2(N3,N3),S(N3,N3),A31T(N1,N3),Z3(N3,N1),
C A32T(N2,N3),Z4(N3,N2),Z1(N1,N3),Z2(N2,N3),
C A11C(N1,N1),A22C(N2,N2)
TRANSN='N'
TRANST='T'
N=N1+N2+N3
C let's make a copy of A11,A22,A13 and A23
do i=1,N1
do j=1,N1
A11C(i,j)=A11(i,j)
enddo
do j=1,N3
Z1(i,j)=A13(i,j)
enddo
enddo
do i=1,N2
do j=1,N2
A22C(i,j)=A22(i,j)
enddo
do j=1,N3
Z2(i,j)=A23(i,j)
enddo
enddo
c let's find the value of Z1=(A11-1)*A13. We know that A11*Z1=A13
c let's get the LU decomposition of A11
call DGETRF( N1, N1, A11C, N1, IPIVZ1, INFO)
c let's solve the system A11*Z1=A13 to find Z1.
call DGETRS(TRANSN,N1,N3,A11C,N1,IPIVZ1,Z1,N1,INFO)
c=====
c let's find the value of Z2=(A22-1)*A23. We know that A22*Z2=A23
c let's get the LU decomposition of A11
call DGETRF( N2, N2, A22C, N2, IPIVZ2, INFO)
c let's solve the system A22*Z2=A23 to find Z2.
call DGETRS(TRANSN,N2,N3,A22C,N2,IPIVZ2,Z2,N2,INFO)
C=====
c we want to fins S, That's why we have to make some products:
c A31Z1=A31*Z1
call prod(N3,N1,N3,A31,Z1,A31Z1)
c A32Z2=A32*Z2
call prod(N3,N2,N3,A32,Z2,A32Z2)
c and now we can calculate S as:
do i=1,N3
do j=1,N3
S(i,j)=A33(i,j)-A31Z1(i,j)-A32Z2(i,j)
enddo
enddo
C===============================
C we need also to find the products: Z3=A31*(A11-1) and Z4=A32*(A22-1)
c we know that Z3*A11=A31===>(A11)'*(Z3)'=(A31)'
c the first thing we have to do is transposing A31:
call transp(N3,N1,A31,A31T)
c and we solve the system making the LU decomposition before.
call DGETRS(TRANST,N1,N3,A11C,N1,IPIVZ1,A31T,N1,INFO)
c now, A31T is Z3 trasposed. Let's find Z3:
call transp(N1,N3,A31T,Z3)
c====================
c let's do the same with Z4
call transp(N3,N2,A32,A32T)
c and we solve the system making the LU decomposition before.
call DGETRS(TRANST,N2,N3,A22C,N2,IPIVZ2,A32T,N2,INFO)
c now, A32T is Z4 trasposed. Let's find Z4:
call transp(N2,N3,A32T,Z4)
C we have all the elements of the LU decomposition.
end
C===========================================
subroutine solve(A11,A13,A22,A23,A31,A32,A33,
C N1,N2,N3,S,Z3,Z4,B1,B2,B3,A11C,A22C,IPIVZ1,IPIVZ2)
integer INFO,i,j,k
CHARACTER(1) TRANSN
integer IPIVZ3(N3),IPIVZ1(N1),IPIVZ2(N2)
double precision A11,A13,A22,A23,A31,A32,A33,A,A11C,A22C
C A31Z1,A32Z2,S,A31T,A32T,Z3,Z4,Z4B2,Z3B1,B1,B2,B3,A13X3,
C A23X3
dimension A11(N1,N1),A13(N1,N3),A22(N2,N2),A23(N2,N3),
C A31(N3,N1),A32(N3,N2),A33(N3,N3),A(N1+N2+N3,N1+N2+N3),
C A31Z1(N3,N3),A32Z2(N3,N3),S(N3,N3),Z3(N3,N1),Z4(N3,N2),
C Z4B2(N3),Z3B1(N3),B1(N1),B2(N2),B3(N3),A13X3(N1),A23X3(N2),
c A11C(N1,N1),A22C(N2,N2)
TRANSN='N'
N=N1+N2+N3
c let's find the solutions:
c we know that S*X3=B3-Z3*B1-Z4*B2
c we have to do the LU decomposition of S
c then we have to solve the system. Be compute B3=B3-Z3*B1-Z4*B2 and
c after calling DGETRS, B3 will be the solution X3.
call prod(N3,N1,1,Z3,B1,Z3B1)
call prod(N3,N2,1,Z4,B2,Z4B2)
do i=1,N3
B3(i)=B3(i)-Z3B1(i)-Z4B2(i)
enddo
call DGETRF(N3,N3,S,N3,IPIVZ3,INFO)
call DGETRS(TRANSN,N3,1,S,N3,IPIVZ3,B3,N3,INFO)
c now, we have the solutions of X3 in B3
c to calculate x1, we can compute: A11*X1=B1-A13*X3
call prod(N1,N3,1,A13,B3,A13X3)
do i=1,N1
B1(i)=B1(i)-A13X3(i)
enddo
call DGETRS(TRANSN,N1,1,A11C,N1,IPIVZ1,B1,N1,INFO)
c now, the solution X1 is in B1
c to calculate x2, we can copute: A22*X2=B2-A23*X3
call prod(N2,N3,1,A23,B3,A23X3)
do i=1,N2
B2(i)=B2(i)-A23X3(i)
enddo
call DGETRS(TRANSN,N2,1,A22C,N2,IPIVZ2,B2,N2,INFO)
c now, the solution X2 is in B2
end
C======================
C======================
SUBROUTINE prod(M1,M2,M3,A,B,C)
double precision A,B,C
integer M1,M2,M3
dimension A(M1,M2),B(M2,M3),C(M1,M3)
c let's comput C=A*B
do i=1,M1
do j=1,M3
C(i,j)=0
do k=1,M2
C(i,j)=C(i,j)+A(i,k)*B(k,j)
enddo
enddo
enddo
end
C======================
C======================
SUBROUTINE transp(M1,M2,A,B)
double precision A,B
integer M1,M2
dimension A(M1,M2),B(M2,M1)
c let's transpose A and put the values in B
do i=1,M2
do j=1,M1
B(i,j)=A(j,i)
enddo
enddo
end
