Problem with fortran and lapack

[Telemachus]

Hi there. I'm trying to solve a linear system which I have constructed with the aim of learning how to do this in fortran. The idea is to solve:

$A\vec{x}=\vec{b}$

The thing is that when I call the subroutine degsv from lapack to solve the linear system, and then write the solution, it only prints the vector b, and I can't find what I am doing wrong. It should solve the system by calling dgesv, and then print x in b.

This is the code:

      program matrix

!... This program constructs a matrix A, a vector b, and then solves
!... the system Ax=b using the lapack package.

    
      implicit none      integer mxpts
      parameter (mxpts=5000)      
      integer i,j,rows,cols,k,aux,pivot(3),ok

    
      real*8 A(mxpts,mxpts),b(mxpts),c,x(mxpts)

      open (unit=10,file='matrix.dat',status='unknown')
      open (unit=20,file='rhs.dat',status='unknown')

!... Inputs
    
      rows=3
      cols=rows !square matrix

!... Coefficients
!... Initialization      do i=1,rows
        b=0.0d0
        do j=1,cols
          A(i,j)=0.0d0
       enddo
      enddo

!... Matrix preparation

      do i=1,rows
        do j=1,cols
          k=i+j
          aux=mod(k,3)
          c=1.0d0
          if(aux.eq.1) c=-1.0d0
!          print*,k,aux,c
          A(i,j)=(i+j)*j*c
        print 200,i,j,A(i,j)
       enddo
      enddo

!... Vector preparation for the right hand side of A*x=b.

      do i=1,rows
        b(i)=i
      print*,i,b(i)
      enddo
    200   format(4x,"A(",i1,",",i1,")=",1pg15.8)  !,f10.8)

!... 4x: 4 espacios en blanco, i1 imprime entero y le asigna 1 espacio
!... 1pg10.5 imprime real, le asigna 10 espacios donde 5 son para
!... decimales.  

!... Print Matrix

        do i=1,rows
                do j=1,cols
                    write (10,*)A(i,j)
                enddo
        enddo

!      do i=1,rows
!        print*,( A(i,j), j=1,cols)
!      enddo

!... Find the solution x=A-1*b

        call dgesv (rows, 1, A, 3, pivot, b, 3, ok)
c
c parameters in the order as they appear in the function call
c    order of matrix A, number of right hand sides (b), matrix A,
c    leading dimension of A, array that records pivoting,
c    result vector b on entry, x on exit, leading dimension of b
c    return value   
c  
c print the vector x

   do i=1,rows
      print*,b(i)
           write (20,*) b(i)
   end do

       stop    
       end

<<Moderator's note: please use code tags around your code>>

 

[Mark44]

Are all of your parameters to dgesv() properly formed? Here's some documentation for this function: http://www.math.utah.edu/software/lapack/lapack-d/dgesv.html

What value do you get in the last parameter, which you're calling ok? A successful call to dgesv() should set this to 0.

 

[Telemachus]

Hello Mark. Thank you very much for your feedback. I'm not really sure about that, because its the first time I'm attempting to use lapack. I have set those parameters while following this example: http://physics.oregonstate.edu/~landaur/nacphy/lapack/codes/linear-f.html , which works well, but uses sgesv (I have assumed that for dgesv should be the same, and I had took a look at the lapack documentation, but as I said, I'm not an expert).

I have tried as you said to put that parameter equal to zero, but I get this error when I run the program:

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0 0x7F9354A59E08
#1 0x7F9354A58F90
#2 0x7F93546AA4AF
#3 0x7F9354F11F18
#4 0x400F2D in MAIN__ at matrix.f:?
Segmentation fault (core dumped)

That parameter is associated with an output from lapack.

 

[eys_physics]

Here are my sugggestions:

1) The Lapack routine dgesv solves the linear system AX=B, where A and B are matrices. So, if you want to solve for the case where the right-hand side is a vector you should declare B as a two-dimensional array of size (n,1) where n is the number of rows.
2) Before you call dgesv you need to compute the LU factorization of A by using the routine DGETRF. This routine produces an array which contains the pivot indices. In your case it is the variable called pivot. But, I suppose that for the moment this array only contains zeros because you don't call DGETRF.
3) As Mark44 said please print out and post the value of the variable "ok".

 

[Mark44]

I have tried as you said to put that parameter equal to zero

No, the 6th parameter is an output variable. You don't set it -- the function sets it, with 0 indicating a successful completion. If the value is negative or positive, the value gives useful information, as explained in the documentation page I linked to.

3) As Mark44 said please print out and post the value of the variable "ok".

From the page in the link:
SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
Of these six parameters, A and B are input/output, which means they must be set when the subroutine is called, but can be changed by the subroutine.
IPIV and INFO are output parameters, which means they can be anything when the subroutine is called, but get set by the subroutine.

 

[DrClaude]

1) The Lapack routine dgesv solves the linear system AX=B, where A and B are matrices. So, if you want to solve for the case where the right-hand side is a vector you should declare B as a two-dimensional array of size (n,1) where n is the number of rows.

This is not necessary. Since Fortran is "pass by reference," the address of the first element is passed to the routine, and the shape of b doesn't matter.

Here are my sugggestions:
2) Before you call dgesv you need to compute the LU factorization of A by using the routine DGETRF. This routine produces an array which contains the pivot indices. In your case it is the variable called pivot. But, I suppose that for the moment this array only contains zeros because you don't call DGETRF.

That's incorrect. dgsev calls dgetrf itself on line 167.

 

[DrClaude]

integer mxpts
parameter (mxpts=5000)
integer i,j,rows,cols,k,aux,pivot(3),ok

There is a first problem here: pivot has to be the same size as the actual matrix A, so you should declare it as size mxpts.

call dgesv (rows, 1, A, 3, pivot, b, 3, ok)

The dimensions here are all wrong.
SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )

The first 3 parameters are fine, but LDA = 3 is incorrect. LDA is the distance in memory between the the first element in the column and the first element in the column next to it. In other words, it is the first dimension of A, so mxpts.

As I said above, IPIV is of size N, which is not how you declared pivot.

As for LDA, LDB should be the first dimension of B, so also mxpts.

This is all explained in the comments at the top of the LAPACK routine. You can also check this information on netlib.

 

[Telemachus]

I've found the problem, the thing seems to be that I have specified the matrix to be from a size which doesn't correspond to the size I was trying to solve. I have changed the value of the parameter mxpts, and then used mxpts as the input to call dgesv.

Thank you all.

 

[Telemachus]

There is a first problem here: pivot has to be the same size as the actual matrix A, so you should declare it as size mxpts.
The dimensions here are all wrong.
SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )

The first 3 parameters are fine, but LDA = 3 is incorrect. LDA is the distance in memory between the the first element in the column and the first element in the column next to it. In other words, it is the first dimension of A, so mxpts.

As I said above, IPIV is of size N, which is not how you declared pivot.

As for LDA, LDB should be the first dimension of B, so also mxpts.

This is all explained in the comments at the top of the LAPACK routine. You can also check this information on netlib.

You are totally right. Thank you very much.

 

[Telemachus]

I don't know if I should start another topic for this. The thing is that fortran after running dgesv to solve $A\vec{x}=b$ gives me this message:user@user-desktop:~/Desktop/radtransf$ ./matrix.x
rmax= 10.000 npoints= 5000 Dr= 2.00000E-03
Rows,hrows= 10000 5000
Info= 0
Note: The following floating-point exceptions are signalling: IEEE_DENORMAL

What could be the problem?

 

[DrClaude]

That seems to be some kind of underflow. What compiler are you using?

Since you get back info = 0, I wouldn't worry too much. You can always check the result: Does Ax give you back b?

 

[Telemachus]

I'm using gfortran. I'll try to do Ax to get b. Now one of the $x_i$ is wrong, and I know that because it is set to zero when constructing the matrix, its a boundary condition, but when lapack does its work it seems to change its value to some small but non zero number.

 

[DrClaude]

I'm using fortran.

From what I could find on the web, there was a problem with gfortran throwing floating point exceptions. I really wouldn't worry too much.

I'll try to do Ax to get b. Now one of the $x_i$ is wrong, and I know that because it is set to zero when constructing the matrix, its a boundary condition, but when lapack does its work it seems to change its value to some small but non zero number.

You have to have some tolerance for numerical errors. What kind of numbers are you getting instead of 0?

 

[Telemachus]

4.6338717745555646E-015

With typical values beeing like: E-04 to E-06 the algorithm doesn't seems to be working. I'm trying to solve an integrodifferential equation, and I get nusty things.

BTW: Have you ever worked with the $S_N$ approximation to the Boltzmann equation?

After checking, its clearly not working.

This is the code, i'll attach it with the input if you want to check.

        program radtransf                            !Nov 29 2016
 
 
!... This program calculates the solution for the radiative transfer
!... equation in the Sn approximation, with n=2. The program starts
!... by building a Matrix for the discretized system of coupled
!... differential equations, obtained by the Diamond difference scheme.!... Version 0.0.1

        implicit none
 
        real*8 Dr,rmax
        integer npoints
        common/bkinput/ Dr,rmax,npoints!... Input file
        open (unit=10,file='radtrans.inp',status='old')!... call input
        call readinput

!... Matrix generation and resolution
        call matrix!... Close files

        close (unit=10)
        close (unit=17)
        close (unit=18)
        close (unit=19)
        close (unit=20)
 
        stop
        end
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        subroutine readinput
        implicit none
 
!... Input parameters

 
        real*8 Dr,rmax
        integer npoints
        common/bkinput/ Dr,rmax,npoints
 
        real*8 rzero, one, two, five
        data rzero,one,two,five/0.0d0,1.0d0,2.0d0,5.0d0/
 
 
        namelist /griddata/Dr,npoints,rmax
     
        Dr=rzero
        npoints=0
        rmax=rzero
 
!... read griddata

        read(10,griddata)
        if (rmax.eq.rzero) rmax=Dr*npoints
        if (Dr.eq.0) Dr=rmax/npoints
        write(*,111) rmax,npoints,Dr
111     format(/5x,'rmax=',f9.3,5x,'npoints=',i8,5x,'Dr=',1pg12.5)        return
        end

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      subroutine matrix

!... This program constructs a matrix A, a vector b, and then solves
!... the system Ax=b using the lapack package.

 
      implicit none      real*8 Dr,rmax
      integer npoints
      common/bkinput/ Dr,rmax,npoints
 
      integer mxpts
      parameter (mxpts=10000)
      integer i,j,rows,hrows,cols,k,aux,pivot(mxpts),info

 
      real*8 A(mxpts,mxpts),b(mxpts),v(mxpts),c,x,mu1,mu2      real*8 one, two, four, hundred
      data one, two, four, hundred/ 1.0d0, 2.0d0, 4.0d0, 100.0d0/!... Output files

!      open (unit=16,file='rhs.dat',status='unknown')
      open (unit=17,file='matrix.dat',status='unknown')
      open (unit=18,file='Psi1.dat',status='unknown')
      open (unit=19,file='Psi2.dat',status='unknown')
      open (unit=20,file='check.dat',status='unknown')

!... Inputs
 
      rows=2*npoints
      hrows=npoints
      cols=rows !square matrix

      print*,'Rows,hrows=',rows,hrows

!... Coefficients
!... Initialization

        mu1=dsqrt(1.0d0/3.0d0)
        mu2=-1.0d0*dsqrt(1.0d0/3.0d0)

      do i=1,rows
        b=0.0d0
        do j=1,cols
          A(i,j)=0.0d0
       enddo
      enddo

!... Build matrix A

      do i=1,hrows
        do j=1,cols
        if(i.eq.1 .and. i.ne.j) then
          A(i,j)=0.0d0
        else if(i.eq.1 .and. i.eq.j) then
          A(i,j)=1.0d0
        else if(i.ne.1 .and. i.eq.j) then
          A(i,j)=two*mu1
        else if(i.ne. 1 .and. i-1 .eq. j) then
          A(i,j)=-two*mu1
        else if(i.ne.1 .and. i+j.eq.rows+1) then
          A(i,j)=-hundred*Dr
        else if(i.ne.1) then
          A(i,j)=0.0d0
        endif
       enddo
      enddo

      do i=hrows+1,rows
        do j=1,cols
        if(i.eq.hrows+1 .and. i.ne.j) then
          A(i,j)=0.0d0
        else if(i.eq.hrows+1 .and. i.eq.j) then
          A(i,j)=1.0d0
        else if(i.ne.hrows+1 .and. i.eq.j) then
          A(i,j)=-two*mu2
        else if(i.ne.hrows+1 .and. i-1.eq.j) then
          A(i,j)=two*mu2
        else if(i.ne.hrows+1 .and. i+j.eq.rows+1) then
          A(i,j)=-hundred*Dr
        else if(i.ne.hrows+1) then
          A(i,j)=0.0d0
        endif
       enddo
      enddo

!        do i=1,rows
!           do j=1,cols
!                print 200,i,j,A(i,j)
!           enddo
!       enddo

!... Vector preparation for the right hand side of Ax=b.

        do i=1,rows
                 b(i)=0.01d0*Dr
        enddo

       b(1)=0.0d0
       b(hrows+1)=0.0d0!    do i=1,rows
!      print*,b(i)
!    end do!200   format(4x,"A(",i5,",",i5,")=",1pg15.8)  !,f10.8)!... Print Matrix

        do i=1,rows
                do j=1,cols
                    write (17,200)i,j,A(i,j)
                enddo
        enddo

200   format(4x,"A(",i5,",",i5,")=",1pg15.8)  !,f10.8)

!      do i=1,rows
!        print*,( A(i,j), j=1,cols)
!      enddo

!... Find the solution x=A-1*b

        call dgesv (rows, 1, A, mxpts, pivot, b, mxpts, info)

        print*,"Info=",info    do j=1,hrows+1
           x=(j-1)*Dr
           write(18,*)x,b(j)
    end do

 
    do i=hrows,rows
           x=(i-hrows)*Dr
           write(19,*)x,b(hrows+1+rows-i)
    end do

!... Verification for Ax=b

      do i=1,rows
        do j=1,cols
          v(i)=A(i,j)*b(j)
        write(20,*)v(i)
       enddo
      enddo

!      do i=1,rows
!        b(i)=0.01d0*Dr
!        write(20,*)v(i),b(i)
!      enddo

       return
       end

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ok, now I see that I was thinking too much symbolically.

I think now the matrix product is ok:

!... Verification for Ax=b
!... Initialization

      do i=1,rows
        v(i)=0.0d0
      enddo

      rsum=0.0d0
      sumaux=0.0d0

!... Verification for Ax=b

      do i=1,rows
        do j=1,cols
          sumaux=A(i,j)*b(j)
          rsum=rsum+sumaux
          v(i)=rsum
       enddo
       write(20,*)v(i)
      enddo

Anyway, only the first element of the array gives the correct value.

 

[eys_physics]

Please, notice that the matrix A is changed by dgesv. So, to do the test you need to make a copy of A before calling that subroutine.

 

[Mark44]

Please, notice that the matrix A is changed by dgesv. So, to do the test you need to make a copy of A before calling that subroutine.

@Telemachus, here's what I said in post #5:

Of these six parameters, A and B are input/output, which means they must be set when the subroutine is called, but can be changed by the subroutine.

When you're using the lapack routines (or, in fact, the routines in any library), it's very important to read and understand the documenation for the routines.

 

[Telemachus]

Thanks. And sorry for not giving the proper attention to those important recommendations, I am too anxious trying to make this work, so sorry for that, and thank you all for your patience, I'll try to not commit those mistakes again, and to be more lucid in the future.

I have rebuilt the matrix after calling dgesv. Still not giving the correct right hand side after matrix multiplication. Could this be a consequence of the matrix A being ill conditioned? I haven't proven it is actually ill conditioned, but I think it could be the case.

 

[Telemachus]

Hi again. I'm having additional problems with this code. When I try to build a bigger matrix, so I can get a finer grid, I have an upper bound. When I try to set the parameter mxpts>20000, I get this error when I try to compile it::~/Desktop/radtransf$ gfortran radtrans.f -llapack -o radtrans.x
/tmp/ccNJ2Cfs.o: In function `matrix_':
radtrans.f:(.text+0x198): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
radtrans.f:(.text+0x1a3): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
radtrans.f:(.text+0x494): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
radtrans.f:(.text+0x6d0): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
radtrans.f:(.text+0x77f): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
radtrans.f:(.text+0x922): relocation truncated to fit: R_X86_64_32 against `.bss'
radtrans.f:(.text+0x9eb): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
radtrans.f:(.text+0xacc): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
radtrans.f:(.text+0xda8): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
radtrans.f:(.text+0xfe4): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
radtrans.f:(.text+0x1098): additional relocation overflows omitted from the output

If the parameter (which is used to store the arrays for matrix and vectors) is smaller than 20000, then the program compiles without any problem.

 

[eys_physics]

Maybe, you don't have enough RAM memory. How much memory do you have on your computer?
With the dimension 20 000, you will have 4*10^8 matrix elements, where each of them is a double-precision number. And with several arrays/matrices it gets even worse.

 

[Telemachus]

The compilller knows that? I'm short of ram, 3.8gb.

 

[Telemachus]

Now I'm going to explain in more detail what the program does, so perhaps if I'm committing some mistake in the derivation of the difference scheme, you can advise me.

I am trying to solve the stationary Boltzmann transport equation in slab geometry with isotropic scattering.

It reads:

$\displaystyle \mu \frac{d \psi(x,\mu)}{dx}+\sigma_T \psi(x,\mu)= \sigma_s \int_{-1}^{1}\psi(x,\mu')d\mu'+Q(x,\mu)$

I'm being quiet general, and then I'll go through the steps into the particular case I am trying to solve. So, this equation model the stationary angular flow of particles $\psi(x,\mu)$ at point $x$ in the direction of $\mu$. $\sigma_T$ is the total cross section: $\sigma_T=\sigma_s+\sigma_a$, $\sigma_s$ is the scattering cross section, and $\sigma_a$ is the absorption cross section. This is a very general equation that models many types of transport phenomena. I'm interested in particular in radiative transfer problems, but its used in plasma physics, hydrodynamics, nuclear reactor design, etc, etc. As you can see its an integrodifferential equation, and is not that trivial to solve. I am trying to reproduce a case in this paper (the first in the section of numerical results): http://www.sciencedirect.com/science/article/pii/0021999187901707

This equation in this regime asymptotically behaves as a diffusion equation, which I have solved analytically (its very simple), and the numerical result is show in the paper, with which it agrees, so I in advance know what the solution should be.

Now, this equation has to be discretized in both variables it depends on (even in one dimension one has scattering from all directions). To solve the integral on the right hand side it is used a Gauss-Legendre quadrature (known as the $S_N$ approximation, I have started with N=2, the paper does it for N=8).

In general one has:

$\displaystyle \mu_m \frac{d \psi_m(x)}{dx}+\sigma_T \psi_m(x)=\sigma_s \sum_{m'} w_{m'} \psi_m' +Q(x)$ (1)

This is a set of coupled differential equations (the scattering term in the rhs couples all directions in the flux).

To do the spatial discretization, one defines the cell average flux:

$\displaystyle \psi_{m,j}=\frac{1}{\Delta x_j} \int_{x_{j-1/2}}^{x_{j+1/2}} \psi_m(x)dx$

Then, integrating the eq. (1) in the j cell and dividing by $\Delta x_j=x_{j+1/2}-x_{j-1/2}$

One arrives to:

$\mu_m \left ( \psi_{m,j+1/2}-\psi_{m,j-1/2} \right )+\Delta x_j \sigma_{T,j} \psi_{m,j}=\sigma_{s,j} \Delta x_j \sum_{m'} \psi_{m',j}w_{m'}+Q_j \Delta x_j$

At this point one defines the discretization scheme by connecting the cell edges (all half integers in the subindexes of x's), with the cell centers (all integer subindexes in x's).

I am using the Diamond difference scheme:

$\psi_{m,j}=\frac{1}{2} \left ( \psi_{m,j+1/2}+\psi_{m,j-1/2} \right )$

From this equation I eliminate $\psi_{m,j+1/2}$ for the equations in which $\mu>0$

And $\psi_{m,j-1/2}$ for the equations in which $\mu<0$

The case I want to solve is with $\sigma_T=\sigma_s=100$ independent of x, and $Q=0.01$, with boundary conditions:

$\psi_m(0)=0$ for $\mu>0$
$\psi_m(10)=0$ for $\mu<0$

When I construct the matrix I drop the half integer indexes, that could be a mistake, but I'm not sure, I just thought that I could redefine the half integers, and everything would be the same, but now I'm not that sure.

For $S_2$, with only two spacing I get this system of equations:

$2\mu_1 \psi_{1,i}-2\mu_1\psi_{1,i-1}-100\Delta x_i \psi_{2,i}=0.01\Delta x_i$

$-2\mu_2 \psi_{2,i}+2\mu_2\psi_{2,i+1}-100\Delta x_i \psi_{1,i}=0.01\Delta x_i$

For this case I build this matrix, the first and third row indicate the boundary condition:

$\begin{bmatrix} 1 & 0 & 0 & 0 \\ -2\mu_1 & 2\mu_1 & -100\Delta x_1 & 0 \\ 0 & 0 & 1 & 0 \\ -100 \Delta x_1 & 0 & 2mu_2 & -2\mu_2 \end{bmatrix} \left[ \begin{array}{c} \psi_{1,1} \\ \psi_{1,2} \\ \psi_{2,2} \\ \psi_{2,1} \end{array} \right] = 0.01\times \left[ \begin{array}{c} 0 \\ \Delta x_1 \\ 0 \\ \Delta x_1 \end{array} \right]$

Then I did the same for 4 spatial intervals, and then I've generalized for the routine I wrote in the program, which works well for this two test cases.

 

[DrClaude]

The memory problem seems to be related to a common block. These have to be allocated as a single chunck of memory, so there might be a size limit. In any case, I don't think that solving for a 20000x20000 matrix can be done in any reasonable amount of time.

I'm not sure I understand the structure of the matrix you obtain. Is it tridiagonal?

 

[Telemachus]

No, it isn't, because I have terms that couples all equations. For the 4 spacing case this is the matrix I have obtained and from which I did the general case for any spacing:

$\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -2\mu_1 & 2\mu_1 & 0 & 0 & 0 & 0 & -100\Delta x_2 & 0 \\ 0 & -2\mu_1 & 2\mu_1 & 0 & 0 & -100\Delta x_3 & 0 & 0 \\ 0 & 0 & -2\mu_1 & 2\mu_1 & -100\Delta x_4 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & -100 \Delta x_3 & 0 & 2\mu_2 & -2\mu_2 & 0 & 0 \\ 0 & -100 \Delta x_2 & 0 & 0 & 0 & 2\mu_2 & -2 \mu_2 & 0 \\ -100 \Delta x_1 & 0 & 0 & 0 & 0 & 0 & 2\mu_2 & -2\mu_2 \end{bmatrix} \left[ \begin{array}{c} \psi_{1,1} \\ \psi_{1,2} \\ \psi_{1,3} \\ \psi_{1,4} \\ \psi_{2,4} \\ \psi_{2,3} \\ \psi_{2,2} \\ \psi_{2,1} \end{array} \right] = 0.01\times \left[ \begin{array}{c} 0 \\ \Delta x_2 \\ \Delta x_3 \\ \Delta x_4 \\ 0 \\ \Delta x_3 \\ \Delta x_2 \\ \Delta x_1 \end{array} \right]$

In the case I'm working with all $\Delta x_i$ are equal, I'm using a uniform grid.

I have seen methods in the bibliography that obtain a Lower triangular matrix by using an iteration scheme in a way that the upper lower terms goes to the right hand side of the equation (those terms comes from the integral, they solve the ingegral recursively by updating the $\psi$). Perhaps I'll try with that, but I have also read that the convergence for the high scattering case (the one I'm dealing with belongs to that category) is slow (this can be explained on physical grounds), so I would have complications in that sense too. The paper anyways gets the answer, and is from 1987, so I must be doing something wrong.

The problem could be at this point:

This is the equation I have to solve:
(1) $\mu_m \left ( \psi_{m,j+1/2}-\psi_{m,j-1/2} \right )+\Delta x_j 100 \psi_{m,j}=100 \Delta x_j \sum_{m'} \psi_{m',j}w_{m'}+ 0.01 \Delta x_j$I am using the Diamond difference scheme:

$\psi_{m,j}=\frac{1}{2} \left ( \psi_{m,j+1/2}+\psi_{m,j-1/2} \right )$

So what I do is to eliminate $\psi_{m,j+1/2}$ when $\mu_m>0$ (in this case the boundary condition is $\psi(0)=0$, so I need to go in the forward direction, and I need the previous step to update):

$\psi_{m,j+1/2}=2\psi_{m,j}-\psi_{m,j-1/2}$

And I eliminate $\psi_{m,j-1/2}$ when $\mu_m<0$ (in this case I go in the backward direction, and the boundary condition is $\psi(10)=0$):

$\psi_{m,j-1/2}=2\psi_{m,j}-\psi_{m,j+1/2}$

So, depending on the sign of $\mu$ I use this conditions in equation (1), so, for example, for $S_2$ with $\mu>0$ I get:

$2\mu_1 \psi_{1,j}-2\mu_1\psi_{1,j-1/2}-100\Delta x_i\psi_{2,i}=0.01 \Delta x_i$

And for $\mu<0$ (this always comes in pairs with opposite sign) I get:

$-2\mu_2 \psi_{2,j}+2\mu_2\psi_{1,j+1/2}-100\Delta x_i \psi_{1,i}=0.01 \Delta x_i$

You must have in mind that in this case $w_{m'}=1$, all weights are equal, so I have canceled the terms on the right hand side with those on the lhs, and $\mu_1=-\mu_2=\sqrt{\frac{1}{3}}$

Now, when I build the matrix, I dropped the subindexes $\psi_{m,j-1/2}$ and $\psi_{m,j+1/2}$ and just called them $\psi_{m,i-1}$ and $\psi_{m,i+1}$. Now I think that could be a mistake, but I'm not sure on how to fix it. I think that the matrix should take the information from the half integer x's which corresponds to the boundaries of the cell's, but should only operate on the center of the cells somehow (but I'm not totally sure of this).

This is a cell, the center and the boundaries:
__________________________
|[COLOR=#black]......[/COLOR]|[COLOR=#black]......[/COLOR]|
$x_{1/2}$[COLOR=#black]...[/COLOR]$x_{1}$[COLOR=#black].....[/COLOR] $x_{3/2}$