solving linear system with zero rows(picking rows and columns) 3

[dibas]

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.

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.

 

[ArkM]

Are you sure that you have "a row of zeroes"? If so, no solution of such a system. Look at minus sign in some -0.00000s: it's a very small negative value!..

 

[dibas]

Those numbers are practically zeros, even when printed to 16 decimal places. They are just floating point artifacts.

 

[gummibaer]

Think about your system as linear equations. Then you have one trivial equation that says 0 = 0 for every row of zeroes. If you have kk of these trivial equations, then you have (6 - kk) equations to define the 6 unknown entities --> no unambiguous solution. If you had the same number of zero columns, then you once again have a solvable system.

But why bother to code a new solver of linear equations ? As this is a standard problem there are a lot of libraries around. Even if I did not use it as yet, I heard a lot about LAPACK. See this link:

http://www.netlib.org/lapack/

Norbert

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[ArkM]

That's a problem: no zero columns in the matrix - no solution...

 

[dibas]

Yes gummibaer, I have tried a countless number of times to get LAPACK onto my system but failed. My PC is running on Linux Ubuntu 11.10 & uses the Gnome desktop (&terminal). My compiler is "gfortran".
Even though my professors don not fully approve of using "black-box" solvers, I always wish to compare my results with what LAPACK would give. Starting from having the "lapack-3.4.1.tgz" file, is it possible for you to provide me with practical steps of installing the package? In one of my attempts, I did some steps to "build" and "test" but the package was not
properly installed such that I could use it. So I then deleted everything.

 

[mikrom]

is it possible for you to provide me with practical steps of installing the package?

Didn't this help you a little bit?

http://www.tek-tips.com/viewthread.cfm?qid=1678628

 

[salgerman]

It's good to want to compare the results from your own fortran programs with algorithms implemented by you to the results from third party libraries or programs; but LAPACK is proving difficult to install, you can always compare to results from other programs like octave, freemat, python/numpy/scipy, etc...all which, I believe are available for installation to your Ubuntu with a single click.

 

[dibas]

mikrom,I'm sorry the example program you gave me (

http://www.tek-tips.com/viewthread.cfm?qid=1678628

) couldn't work because I had failed to properly
install LAPACK on my system. These were the error messages from the two options you gave me:
(1) /tmp/cczJi5YX.o: In function `MAIN__':
example.f95.text+0xe1): undefined reference to `dgesv_'
collect2: ld returned 1 exit status

(2) /tmp/ccKYV7Q0.o: In function `MAIN__':
example.f95.text+0xe1): undefined reference to `dgesv_'
collect2: ld returned 1 exit status

But I think the liblapack.dll is a Windows-based library.


Anyway, I tried to define a lower-order square matrix with from the main
one such that the zero row and the corresponding columns were left out.
The decision to do that does not conflict with the theoretical background of my problem. I'm modifying the solving steps to see if the solving part
works fine.

 

[mikrom]

But I think the liblapack.dll is a Windows-based library.

Yes, I tried it only at windows, because I'm not experienced Ubuntu-Linux user and I don't have it installed now.

But it seems, that there is Lapack package for Ubuntu available:

http://packages.ubuntu.com/precise/liblapack-dev

The Ubuntu's APT packages could be installed simply using the command apt-get.
Because it seems that liblapack-dev is dependent on other packages libblas-dev and liblapack3gf I would try something like this to install them all:

sudo apt-get install liblapack-dev libblas-dev liblapack3gf

 

[gummibaer]

one such that the zero row and the corresponding columns were left out.

What would be the corresponding column ? I do not know of any relationship of rows to columns. Imagine your problem as a system of linear equations again. Then the rows correspond to the sequence how you noted the equations one after the other. So it is clear you could exchange any to rows - together with the vector of course - without any significant change to the equations and their solutions. Similar with the columns. They correspond to the sequence in which the parts of the equations are written down. Without any significant change to the system you can exchange any two columns. So it should be clear that you cannot lay your hand on a column and say this corresponds in any conceivable way to any row.

Hope I could clarify this. You can do whatever you want with your system of equations, with five valid equations your system of six unknown variables does not have a solution. Howgh.

The thing you can do is to ignore that one of the variables is unknown, say you assume x(1) as known. Solving the remaining 5 by 5 system gives you a set of five equations of the form

x = a * x(1) + b

to define x(2) to x(6) when you know a value of x(1).

Only this and nothing more.

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[gummibaer]

Somehow your equations intrigue me.

For certain, having five equations for six unknown entities is not sufficient for an unambiguous solution and that is a fact.

Buuuut....

Looking at your matrix one can see that all the values in the fifth column have a completely different magnitude then the rest. Roughly there is a factor of about 1000 or more between the coefficients in the fith column and the others and the vector elements of your right hand side(well except the last one, but still...). So I would assume that x(5) is smaller than the other x by the factor of 1000 or more.

So maybe you could get an estimate of x(1) to x(4) and x(6) if you approximate x(5) to be zero and thus ignore the fifth column. But, mind you, this is not because this column corresponds to the zero-row in any way but because of the numeric properties of these values. And this solution can only be an approximation to the full solution.

Norbert.

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[salgerman]

I wonder where these variables come from.

Allow me to play Devil's Advocate here, for a moment, and challenge the validity of the choice of variables....is it 100% certain that these variables are true state variables? Or could it be that one of them was mistaken by one and included in the equations when if should have been left out?

Other than that, there is nothing wrong in multiplying one of the rows by 100, or 1000 or whatever in order for it to have similar order of magnitude as the rest in order to make it easy on the solver...otherwise, it may declare it singular or something.

Yet, I don't know, of top of my head, what to do about the very small column.

 

[dibas]

Yes, looking at the coefficient matrix, there equations are intriguing for sure. Just some brief background of the problem. This is part of an evolution problem which uses implicit time-stepping (integrator). So this
involves a Newton-Raphson step at some point. It is in this step that we
then encounter this "linear" system - actually the system is solved at
each iteration for "Newton corrections(dx)" to the state variables,
not the state variables themselves.

According to the theory of the problem, the state variable x2 should not
evolve, which means the Jacobian matrix will have zeros in row2 (hence we
can also set dx ~ 0). It is only included for generalization purposes -
should it be allowed in some instance to evolve?


Also, the last two state variables are special in the sense that x5 is zero
on initially (at t=0) but will then build up in future times. So to avoid
some of the Jacobian elements being undefined because of this, I set them
to zero and gave the value 1 to the diagonal elements so that the Jacobian
matrix (in J.dx = -F) could be 'invertible'. So this brings in the
different orders of magnitude. I'm yet to experiment with this and see if
the one or some other number will be an appropriate choice.

 

[salgerman]

What? ... Hold on a second!

After reading what you have said above, I am getting impression that you are confusing variable values with variable coefficients. I mean, a system of equations should look something like:

a11x1 + a12x2 + a13x3 + a14x4 + a15x5 + a16x6 = b1
.
.
.
a61x1 + a62x2 + a63x3 + a64x4 + a65x5 + a66x6 = b6

So, just because x5 is initially zero, that does not meant that you need to set a15, a25, a35, a45, a55 and a56 to zeroes in your coefficient matrix...the product of ai5x5 will zero out naturally if x5 is zero, YOU do not need to worry about that...a45 or a55 or whatever, could take on any values...you know what I am saying?

So, where do the coefficients in your matrix (that you show in the original posting) come from? Do they have a genuine origin, or did you make some of them up like the zeroes and the 1 that you talked about?

Then again, if you know that x2 is not going to evolve and if your matrix is actually solving for dx's and not x's and, hence, you know that dx2 is zero, THEN, you can actually take the second column out of the picture.

I don't remember why the second row is full of zeroes, but if it is like for some reason and taking second column out brings you back a square system...go for it

That is, of course, if those zeroes truly represent the matrix coefficient and NOT the variable values!!!

Hope it is me who is confused about this and not you

 

[gummibaer]

... so line 2 and column 5 are added to the set artificially so to speak ? Then, if you leave them off, it would be a 5 by 5 system ?

Norbert

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[dibas]

...I am getting impression that you are confusing variable values with variable coefficients." Oh, I agree salgerman. My language use was very
sloppy (or even inaccurate) and hence the confusion between matrix
coefficients and variable. Also, the coefficient matrix I printed out in
my original post should have been:

0.07192 -0.05664 -0.22443 -159.36051 -0.01644 0.00000
0.00000 -0.00000 -0.00000 -0.00000 -0.00000 0.00000
-0.14818 0.11671 0.46242 328.35612 0.03387 0.00000
0.02470 -0.01945 -0.07708 -54.73203 -0.00565 0.00000
-0.00002 0.00002 0.00007 0.04627 0.00000 0.00000
-0.00012 0.00009 0.00037 0.26397 0.00003 1.00000

I'm sorry for that. So when referring to x5 (or ai5), it should
have been x6 (and ai6). I picked this up because I know it is the
sixth variable in my problem that behaves the way I was describing. I
think I mixed up my printing statements from the code.
Yes gummibaer, I want to then reduce the system to a 5x5 because of the
zero row (row2).

Will follow the rest of your details salgerman. It looks I'll get something from there.

 

[mikrom]

Maybe I missed something, but IMHO:
If you have given 5 equations of 6 variables, how want you to reduce it to 5 equations of 5 variables?
It's simply undedetermined system and you cannot solve it using Gauss Elimination Method or derived methods like LU decomposition, because of singularity of the matrix.

But you can use an approximation method, for example least-norm method.
Because you suggested me to install Lapack I tried to use the subroutine DGELS to solve your system and got this result:

$ dibas
 lwork =          198
 info =            0
 Solution was succesfully computed:
x[1]= -33964.12600132240913808300
x[2]=     -0.00000000000000000000
x[3]= -12768.14057842434522172000
x[4]=      4.57661538381135368780
x[5]= -18639.50908847808022983400
x[6]=     -0.00041115380137691132

Here is the source code:

[COLOR=#a020f0]Program[/color] LinearEquations
  [COLOR=#0000ff]! solving the matrix equation A*x=b using LAPACK[/color]
  [COLOR=#2e8b57][b]Implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]

  [COLOR=#0000ff]! declarations[/color]
  [COLOR=#2e8b57][b]integer[/b][/color] :: i, m, n, lda, ldb, lwork, lwmax, info
  [COLOR=#2e8b57][b]parameter[/b][/color](m [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]6[/color], n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]6[/color], lwmax [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]300[/color])
  [COLOR=#0000ff]! lda = max(1, m)[/color]
  [COLOR=#0000ff]! ldb = max(1, m, n)[/color]
  [COLOR=#2e8b57][b]parameter[/b][/color](lda [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]6[/color], ldb [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]6[/color])
  [COLOR=#2e8b57][b]double precision[/b][/color] :: A(lda,n), b(ldb), work(lwmax)
  
  [COLOR=#0000ff]! matrix A[/color]
  A([COLOR=#ff00ff]1[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]0.07192[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.05664[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.22443[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]159.36051[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.01644[/color],[COLOR=#ff00ff]0.00000[/color][COLOR=#804040][b]/[/b][/color])
  A([COLOR=#ff00ff]2[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]0.00000[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.00000[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.00000[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.00000[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.00000[/color],[COLOR=#ff00ff]0.00000[/color][COLOR=#804040][b]/[/b][/color])  
  A([COLOR=#ff00ff]2[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/-[/b][/color][COLOR=#ff00ff]0.14818[/color],[COLOR=#ff00ff]0.11671[/color],[COLOR=#ff00ff]0.46242[/color],[COLOR=#ff00ff]328.35612[/color],[COLOR=#ff00ff]0.03387[/color],[COLOR=#ff00ff]0.00000[/color][COLOR=#804040][b]/[/b][/color])   
  A([COLOR=#ff00ff]3[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]0.02470[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.01945[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.07708[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]54.73203[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.00565[/color],[COLOR=#ff00ff]0.00000[/color][COLOR=#804040][b]/[/b][/color])    
  A([COLOR=#ff00ff]4[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/-[/b][/color][COLOR=#ff00ff]0.00002[/color],[COLOR=#ff00ff]0.00002[/color],[COLOR=#ff00ff]0.00007[/color],[COLOR=#ff00ff]0.04627[/color],[COLOR=#ff00ff]0.00000[/color],[COLOR=#ff00ff]0.00000[/color][COLOR=#804040][b]/[/b][/color])
  A([COLOR=#ff00ff]5[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/-[/b][/color][COLOR=#ff00ff]0.00012[/color],[COLOR=#ff00ff]0.00009[/color],[COLOR=#ff00ff]0.00037[/color],[COLOR=#ff00ff]0.26397[/color],[COLOR=#ff00ff]0.00003[/color],[COLOR=#ff00ff]1.00000[/color][COLOR=#804040][b]/[/b][/color])

  [COLOR=#0000ff]! vector b[/color]

  b(:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/-[/b][/color][COLOR=#ff00ff]0.04441696[/color],[COLOR=#ff00ff]0.00000000[/color],[COLOR=#ff00ff]0.08007247[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.00272736[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.00002423[/color],[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]0.00008325[/color][COLOR=#804040][b]/[/b][/color])

  [COLOR=#0000ff]! workspace query: calculates the optimal size of the WORK array[/color]
  lwork [COLOR=#804040][b]=[/b][/color] [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]
  [COLOR=#008080]call[/color] DGELS( [COLOR=#ff00ff]'No transpose'[/color], m, n, [COLOR=#ff00ff]1[/color], A, lda, b, ldb, work, lwork, info)
  [COLOR=#0000ff]! compute workspace size[/color]
  lwork [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]min[/color](lwmax, [COLOR=#008080]int[/color](work([COLOR=#ff00ff]1[/color])))
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'lwork = '[/color], lwork

  [COLOR=#0000ff]! solve[/color]
  [COLOR=#008080]call[/color] DGELS([COLOR=#ff00ff]'No transpose'[/color], m, n, [COLOR=#ff00ff]1[/color], A, lda, b, ldb, work, lwork, info)

  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'info = '[/color], info
  [COLOR=#804040][b]if[/b][/color] (info [COLOR=#804040][b].eq.[/b][/color] [COLOR=#ff00ff]0[/color]) [COLOR=#804040][b]then[/b][/color]
    [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Solution was succesfully computed:'[/color]
    [COLOR=#0000ff]! print the solution x[/color]
    [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], [COLOR=#ff00ff]6[/color]
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]9[/color]) i, b(i)
    [COLOR=#804040][b]end do[/b][/color]  
  [COLOR=#804040][b]else[/b][/color]
    [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'* Error computing solution!'[/color]
  [COLOR=#804040][b]end if[/b][/color]
  
[COLOR=#6a5acd]9[/color] [COLOR=#804040][b]format[/b][/color]([COLOR=#ff00ff]'x['[/color], i1, [COLOR=#ff00ff]']= '[/color], [COLOR=#008080]f27.20[/color])  
[COLOR=#a020f0]end program[/color] LinearEquations