Solve $AX=B$ where $A$ is a skyline matrix

Question

Solve a matrix equation of the type $AX=B$, where $A$ is an $n \times n$ symmetric matrix stored in the form of symmetric skyline matrix.

With the solution given by Bill and some more research on factorization, I used the prototype code provided by S.H Lee along with the same test code:

    parameter (neqs = 5)
    real*8 au(7), ad(neqs), b(neqs), energy
    integer jp(neqs)
    data au/-2.,-2.,-3.,-1.,0.,0.,4./
    data ad/2.,3.,5.,10.,10./
    data b/0.,1.,0.,0.,0./
    data jp/0,1,2,3,7/
    CALL REDUCl(neqs,jp,ad,au)
    CALL SOLVE1(neqs,jp,ad,au,b)

The results are satisfactory:

solution 636. 619. 292. 74. 34.

The subroutines REDUC1 and SOLVE1 are as follows:

SUBROUTINE  REDUCl  (N,IP,D,A)
IMPLICIT  REAL(8) (A-H,O-Z)
DIMENSION  IP(*),D(*),A(*)
!FUNCTION  :  CROUT  DECOMPOSITION  A  =  LDU
!REDUCTION  COLUMN  BY  COLUMN
DO  11  K=2,N
K1=K-1
LK=IP(K)-K1
KH=IP(K1)-LK+1
S=D(K)-A(LK+KH)*A(LK+KH)/D(KH)
DO  22  J=KH+1,K1
J1=J-1
LJ=IP(J)-J1
JH=IP(J1)-LJ+1
IF  (KH.GT.JH)  JH=KH
T=A(LK+J)
DO  33  M=JH,J1
33  T=T-A(LJ+M)*A(LK+M)
A(LK+J)=T
22  S=S-T*T/D(J)
DO  44 J=KH,K1
L=LK+J
44  A(L)=A(L)/D(J)
11  D(K)=S
RETURN
END

SUBROUTINE  SOLVE1  (N,IP,D,A,B)
IMPLICIT  REAL(8)  (A-H,O-Z)
DIMENSION  IP(*),D(*),A(*),B(*)
!FUNCTION  :  SOLVE  FOR  X,  LDUx  =  b  where  U  =  transpose  of  L
!1.  FORWARD  SUBSTITUTION  :  Lz  =  b
DO  11  J=2,N
J1=J-1
LJ=IP(J)-J1
JH=IP(J1)-LJ+1
T=B(J)
DO  22  M=JH,J1
22  T=T-A(LJ+M)*B(M)
11  B(J)=T
!2.  DIVIDING  BY  DIAGONAL  ELEMENTS  :  Dy  =  2
DO  33  K=1,N
33  B(K)=B(K)/D(K)
!BACKWARD  SUBSTITUTION  :  UX  =  y
DO  44  K=N,2,-1
K1=K-1
LK=IP(K)-K1
KH=IP(K1)-LK+1
T=B(K)
DO  44  J=KH,K1
44  B(J)=B(J)-T*A(LK+J)
RETURN
END

However, now I am stuck in a new problem! this works fine when my matrix is like above but fails (I get wrong answers) for the below two cases:
Case 1:

    data au/0.6801,0.3846,0.4820,1.4074,1.2731,0.7047,0.1078,0.6171,0.9063,0.5895,0.,0.5260,0.3439,0.8178/
    data ad/0.3818,0.4524,0.,1.4605,1.7593,1.1887/
    data b/1.,1.,1.,1.,1.,1./
    data jp/0,0,2,5,9,14/
    CALL REDUCl(neqs,jp,ad,au)
    CALL SOLVE1(neqs,jp,ad,au,b)

Case 2:

!!C  Solve A*U = b
!!C        [1    0    0     0     0            0         ]
!!C        [0    1    0     0     0            0         ]
!!C  A =   [0    0    1     0     0            0         ]
!!         [0    0    0     1     0            0         ]
!!C        [0    0    0           1.6667E10   -0.4167E10 ]
!!C        [0    0    0    0     -0.4167E10    0.1041e10 ]
data au/-0.4166667E10/
data ad/1.,1.,1.,1.,1.6667E10,0.1041667E10/
data b/1.,1.,1.,1.,1.,1./
data jp/0,0,0,0,0,1/
Print*,au
Print*,ad
CALL REDUCl(neqs,jp,ad,au)
CALL SOLVE1(neqs,jp,ad,au,b)

Where am I going wrong? Is it due to some allocation or use of precision?

Answer 1

The matrix in your case 2 example has a much higher condition number than the one in the original 5-equation example. Accordingly, errors in the definition of the matrix terms result in larger errors in the solution.

I believe that if you output the matrices from MATLAB with 16 digits of precision, the fortran codes will yield the same solution as MATLAB. That is, the entries in the au and ad vectors in your fortran code need to be defined with 16 digits of precision instead of the (roughly) 6 digits of precision in the posted code.

In addition, you should be careful to define floating point constants correctly, e.g. 1.d0 instead of 1.e0. This was an oversight in the original 5-equation example that caused no great problem in that case because of the low condition number.

Answer 2

To expand on Bill Greene's answer, Fortran uses different syntax for single and double precision floating point numbers. 1.e0 is a single precision representation of the number '1.0', whereas 1.d0 is double precision representation of the number '1.0'. If you have a rather straight forward problem, and an associated low condition number, you will not notice a difference because the problem can be solved in single precision.

As problems become more difficult, they require higher levels of precision - that's why it's good practice to set your variables to double precision values; however, there is a caveat if you mix single and double precision numbers. Fortran will reduce to the lowest common denominator, which means that multiplying a double precision number with a single precision number will result in a single precision number. This won't be an issue for simple problems, but as problems get more complex and require consistent usage of double precision, this will represent itself as a bug in your code.

---

Copying from this answer:

4.5.1 Double Precision Exponent.

The form of a double precision exponent is the letter D followed by an optionally signed integer constant. A double precision exponent denotes a power of ten. Note that the form and interpretation of a double precision exponent are identical to those of a real exponent, except that the letter D is used instead of the letter E.