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

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)
integer jp(neqs)
data au/-2.,-2.,-3.,-1.,0.,0.,4./
data b/0.,1.,0.,0.,0./
data jp/0,1,2,3,7/


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 b/1.,1.,1.,1.,1.,1./
data jp/0,0,2,5,9,14/


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 b/1.,1.,1.,1.,1.,1./
data jp/0,0,0,0,0,1/
Print*,au


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

• Where did these matrices come from? How did you copy the entries to create the au and ad vectors? What are you comparing with to know that your solution is wrong? – Bill Greene Jan 19 '16 at 13:31
• I working on a discrete element method with spring network. The first matrix was from a randomly generated sparse matrix on Matlab and the second one is an actual stiffness matrix obtained between two node elements. Both of them have been verified on Matlab. Could you please check and tell me the results obtained in your code using 'datri' and ' dasol' – Chaitanya Krishna Jan 19 '16 at 17:05

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.

• Sorry, I couldn't understand exactly. I tried to check everything in matlab as double precision but I do not get the same answers. Could you explain me what you you get with feap? – Chaitanya Krishna Jan 19 '16 at 17:58
• In adddition to this, I could check my results with hand calculations and doesn't really match with the answers. I feel the problem is due to some errors as you have mentioned in your answer. – Chaitanya Krishna Jan 20 '16 at 4:01

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.