numerical accuracy problem

[angel]

I have a numerical accuracy problem in the general band and tridiagonal routines for a linear system, in particular lower accuracy than the sequential code.

I try to solve a linear system of equations Ax = b. The matrix A is a tridiagonal matrix with variable coefficients. I have tried to use some routines, like as:
- PDDBTRF/S: LU, no pivoting, general band.
- PDDTTRF/S: LU, no pivoting, tridiagonal.
- PDGBTRF/S: LU, partial pivoting, general band.

In all cases, I have obtained a solution with an absolute error of order of 10^-3.
I have prepared next program to obtain a vector solution with all terms setted to 1.

Here I attach the code for the PDDBTRF/S case, to run with only one processor:

PROGRAM PDDBTRF_EXAMPLE

*
* This is an example of using PDGBTRF and PDGBTRS.
* A matrix of size 9x9 is distributed on a 1x1 process
* grid, factored, and solved in parallel.
*
*
INTEGER ICTXT, INFO, MYCOL, MYROW, NPCOL, NPROW
INTEGER BWL, BWU, N, NB, LAF, LWORK
DOUBLE PRECISION deltax, deltat, c, xmin, xmax, tmp

PARAMETER (BWL = 1, BWU = 1, NB = 1599, N = 1599)
PARAMETER (LAF=NB*(bwl+bwu)+6*max(bwl,bwu)*max(bwl,bwu))
! LAF=(NB+BWU)*(BWL+BWU)+6*(BWL+BWU)*(BWL+2*BWU))
PARAMETER (LWORK = LAF)

INTEGER DESCA(7), DESCB(7), IPIV(NB), FILL_IN
DOUBLE PRECISION A(1+BWL+BWU,NB),B(NB),AF(LAF),WORK(LWORK)
DOUBLE PRECISION xi(0:NB+1), ux(0:NB+1)
*
* INITIALIZE THE PROCESS GRID
*
NPROW = 1
NPCOL = 1
CALL SL_INIT( ICTXT, NPROW, NPCOL )
CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )

*
* DISTRIBUTE THE MATRIX ON THE PROCESS GRID
* Initialize the array descriptors for the matrices A and B
*
DESCA( 1 ) = 501 ! descriptor type
DESCA( 2 ) = ICTXT ! BLACS process grid handle
DESCA( 3 ) = N ! number of rows in A
DESCA( 4 ) = NB ! Blocking factor of the distribution
DESCA( 5 ) = 0 ! size of block rows
DESCA( 6 ) = 1+BWL+BWU ! leading dimension of A
DESCA( 7 ) = 0 ! process row for 1st row of A

DESCB( 1 ) = 502 ! descriptor type
DESCB( 2 ) = ICTXT ! BLACS process grid handle
DESCB( 3 ) = N ! number of rows in B
DESCB( 4 ) = NB ! Blocking factor of the distribution
DESCB( 5 ) = 0 ! size of block rows
DESCB( 6 ) = NB ! leading dimension of B
DESCB( 7 ) = 0 ! process row for 1st row of B

*
* Generate matrices A and B and distribute them to the process grid
*
* P1 P2 P3 X B
*
* / 1 1 . . \ / 1 \ / 2 \
* | 1 1 1 . . | | 1 | | 3 |
* | 1 1 . 2 . | | 1 | | 4 |
* | 1 . 2 2 . | | 1 | | 5 |
* | . 2 2 2 . | | 1 | = | 6 |
* | . 2 2 . 3 | | 1 | | 7 |
* | . 2 . 3 3 | | 1 | | 8 |
* | . . 3 3 3 | | 1 | | 9 |
* \ . . 3 3 / \ 1 / \ 6 /
*
FILL_IN = 0
deltat = 0.001
deltax = 0.01
c = 1.0
xmin = -8.0
xmax = 8.0
tmp = xmin
DO I=0,NB+1
xi(I) = tmp
tmp = tmp + deltax
ux(I) = cos((xi(I)*3.141592654*0.5)/xmax)
$ *cos((xi(I)*3.141592654*0.5)/xmax)
END DO
! Only one task
IF (MYCOL == 0) THEN
! Matrix A
A(2+FILL_IN,1) = 1.d0 + (g(xi(1))*(deltat/(deltax*deltax)))
. + c*0.5*deltat - 0.5*deltat*ux(1) - deltat*0.5*ux(1)*ux(1)
A(3+FILL_IN,1) = -(g(xi(1))*0.5*deltat)/(deltax*deltax)
write(*,*) ' ', A(2+FILL_IN,1),' ',A(3+FILL_IN,1)
DO I = 2, NB-1
A(1+FILL_IN,I) = -(g(xi(I))*0.5*deltat)/(deltax*deltax)
A(2+FILL_IN,I) = 1.d0 + (g(xi(I))*(deltat/(deltax*deltax)))
. + c*0.5*deltat - 0.5*deltat*ux(I) - deltat*0.5*ux(I)*ux(I)
A(3+FILL_IN,I) = -(g(xi(I))*0.5*deltat)/(deltax*deltax)
write(*,*) A(1+FILL_IN,I),' ',A(2+FILL_IN,I),' ',
. A(3+FILL_IN,I)
IF (ABS(A(2+FILL_IN,I)) < ABS(A(1+FILL_IN,I))+
. ABS(A(3+FILL_IN,I))) THEN
write (*,*) 'No diagonal dominante!!'
END IF
END DO
A(1+FILL_IN,NB) = -(g(xi(NB))*0.5*deltat)/(deltax*deltax)
A(2+FILL_IN,NB) = 1.d0 + (g(xi(NB))*(deltat/(deltax*deltax)))
. + c*0.5*deltat - 0.5*deltat*ux(NB) - deltat*0.5*ux(NB)*ux(NB)
write(*,*) A(1+FILL_IN,NB),' ',A(2+FILL_IN,NB)
! RHS
B(1) = A(2+FILL_IN,1) + A(3+FILL_IN,1)
DO I = 2, NB-1
B(I) = A(1+FILL_IN,I) + A(2+FILL_IN,I) + A(3+FILL_IN,I)
write(*,300) 'A: ', A(1+FILL_IN,I),' ',A(2+FILL_IN,I),' ',
. A(3+FILL_IN,I),' B = ', B(I)
300 format((a2,1(f13.8,2x),a1,1(f13.8,2x),a1,1(f13.8,2x),
. a5,1(f12.8,2x)))
END DO
B(NB) = A(1+FILL_IN,NB) + A(2+FILL_IN,NB)
END IF

*
* Perform LU factorization
*
! CALL PDGBTRF( N, BWL, BWU, A(1,1), 1, DESCA, IPIV,
! $ AF, LAF, WORK, LWORK, INFO )
CALL PDDBTRF( N, BWL, BWU, A(1,1), 1, DESCA, AF, LAF, WORK,
$ LWORK, INFO )


IF (INFO/=0) THEN
write(*,*) 'Info flag from PDDBTRF = ',INFO, ', Col = ',MYCOL
GOTO 100
END IF

*
* Solve using the LU factorization from PDGBTRF
*
! CALL PDGBTRS('N', N, BWL, BWU, 1, A(1,1), 1, DESCA, IPIV,
! $ B(1), 1, DESCB, AF, LAF, WORK, LWORK, INFO)
CALL PDDBTRS( 'N', N, BWL, BWU, 1, A(1,1), 1, DESCA, B, 1,
$ DESCB, AF, LAF, WORK, LWORK, INFO )


IF (INFO/=0) THEN
write(*,*) 'Info flag from PDDBTRS = ',INFO, ', Col = ',MYCOL
GOTO 100
END IF

DO I=1,NB
write(*,200) 'X(',I,') = ',B(I)
END DO
200 format((a2,I4,a4,1(f10.8,2x)))
*
* Release the process grid and exit BLACS
*
100 CALL BLACS_GRIDEXIT( ICTXT )
CALL BLACS_EXIT( 0 )

STOP
END

FUNCTION g(x)
DOUBLE PRECISION x
g = x*x + 1.d0;
END FUNCTION

FUNCTION f(x)
DOUBLE PRECISION x
f = 6.d0;
END FUNCTION


Thank you in advance.

Angel Rodriguez