I'm trying to solve the following linear (6x6) system that contains a row
of zeros in it using LU decomposition and LU back substitution.
coefficient matrix =
0.07192 -0.05664 -0.22443 -159.36051 0.00000 -0.01644
0.00000 -0.00000 -0.00000 -0.00000 0.00000 -0.00000
-0.14818 0.11671 0.46242 328.35612 0.00000 0.03387
0.02470 -0.01945 -0.07708 -54.73203 0.00000 -0.00565
-0.00002 0.00002 0.00007 0.04627 0.00000 0.00000
-0.00012 0.00009 0.00037 0.26397 1.00000 0.00003
right hand side vector =
-0.04441696
0.00000000
0.08007247
-0.00272736
-0.00002423
-0.00008325
My idea is to "eliminate" the zero-row together with the corresponding
column and then solve the "reduced" system. My problem is that I do manage to pick out the non-zero row(s), but I cannot pick the corresponding column(s). This my piece of code.
I don't know if there is an algorithm or implementable method that can handle this case of a linear system. I would appreciate any suggestions.
Thank you.
of zeros in it using LU decomposition and LU back substitution.
coefficient matrix =
0.07192 -0.05664 -0.22443 -159.36051 0.00000 -0.01644
0.00000 -0.00000 -0.00000 -0.00000 0.00000 -0.00000
-0.14818 0.11671 0.46242 328.35612 0.00000 0.03387
0.02470 -0.01945 -0.07708 -54.73203 0.00000 -0.00565
-0.00002 0.00002 0.00007 0.04627 0.00000 0.00000
-0.00012 0.00009 0.00037 0.26397 1.00000 0.00003
right hand side vector =
-0.04441696
0.00000000
0.08007247
-0.00272736
-0.00002423
-0.00008325
My idea is to "eliminate" the zero-row together with the corresponding
column and then solve the "reduced" system. My problem is that I do manage to pick out the non-zero row(s), but I cannot pick the corresponding column(s). This my piece of code.
Code:
program driver_lusolver
use nrtype
use nrutil, only : ludcmp, lubksb
implicit none
integer, parameter :: mrows = 6, ncols = 6 ! rows & columns
integer :: mlog, nlog ! log dimemsions of matrix to be 4md from a
real(sp) :: a(mrows, ncols), b(mrows)
integer, dimension(mrows) :: zero_row_indx
integer :: i, j, kk, zeros_in_row, jj, ll
logical :: any_nonzero
real(sp) :: d
integer(i4b) :: indx(mrows)
real(sp), allocatable :: a_prime(:, :) , b_prime(:)
! matrix A
open(unit=7, file='matrix.dat', status = 'old', action='read')
read(7,*) ((a(i,j), j=1,ncols), i=1,mrows)
close(7)
! matrix b
open(unit=11, file='bb.dat', status = 'old', action='read')
read(11,*) (b(i), i=1,mrows)
close(11)
zero_row_indx = 0 ! array for storing zero-row indices
do i = 1 , mrows
if ( all(abs(a(i,:)) < 1.0e-12) ) then
zero_row_indx(i) = i
end if
end do
kk = count( mask= zero_row_indx > 0 )
! write(*,*) "there are", kk, "zero rows"
! write(*,"(i2)") zero_row_indx
mlog = mrows - kk ! compute logical dimensions of new arrays
nlog = ncols - kk
allocate(a_prime(mlog, nlog), b_prime(mlog)) !create storage for new arrays
a_prime = 0.000 ; b_prime = 0.00 !initialize new matrix & vector
jj = 1 ! initialize row counter of new matrix & vecor
do i = 1, mrows
if ( any(abs(a(i,:)) > 1.0e-12) ) then
a_prime(jj,:) = a(i,:) ! picking row
b_prime(jj) = b(i)
jj = jj + 1
end if
end do
! Initialise 'indx' ( a vector of length mrows)
! containing row permutations
indx(1:mrows) = 0 ; d = -1.0_sp
! Solve by lu decomposition
call ludcmp(a_prime, indx, d)
call lubksb(a_prime, indx, b)
! print solutions x(n)
! write (*,203)
! write (*,201) (b(i),i=1,mlog)
! expected soln: b = [3.000000 1.000000 -2.000000 1.000000]^T
200 format (' LUdecomp and LUbksb: Matrix A and vector')
201 format (6f13.5)
202 format (/,' Matrix A after decomposition')
203 format (/,' Solutions x(n)')
end program driver_lusolver
I don't know if there is an algorithm or implementable method that can handle this case of a linear system. I would appreciate any suggestions.
Thank you.