1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
*DECK MSRLUS
SUBROUTINE MSRLUS(LFORW,N,LC,IM,MCU,JU,ILUDF,ILUCF,XIN,XOUT)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Solve LU y = x where LU is stored in the "L\U-I" form in MSR format.
* Can be use "in-place" i.e. XOUT=XIN.
*
*Copyright:
* Copyright (C) 2002 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): R. Le Tellier
*
*Parameters: input
* LFORW flag set to .false. to transpose the coefficient matrix.
* N order of the system.
* LC dimension of IM.
* IM elements of row i in MCU(IM(i):IM(i+1)-1)
* MCU
* JU MCU(JU(i):IM(i+1)) corresponds to U.
* MCU(IM(i)+1:JU(i)-1) correspond to L.
* ILUDF diagonal elements of U (inversed diagonal).
* ILUCF non-diagonal elements of L and U.
* XIN input vector x.
*
*Parameters: output
* XOUT output vector y.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*---
* SUBROUTINE ARGUMENTS
*---
INTEGER N,LC,IM(N+1),MCU(LC),JU(N)
REAL ILUDF(N),ILUCF(LC)
DOUBLE PRECISION XIN(N),XOUT(N)
LOGICAL LFORW
*---
* LOCAL VARIABLES
*---
INTEGER I,J,IJ
IF(LFORW) THEN
*---
* FORWARD SOLVE
*---
DO I=1,N
XOUT(I)=XIN(I)
DO IJ=IM(I)+1,JU(I)-1
J=MCU(IJ)
IF ((J.GT.0).AND.(J.LE.N)) XOUT(I)=XOUT(I)
1 -ILUCF(IJ)*XOUT(J)
ENDDO
ENDDO
*---
* BACKWARD SOLVE
*---
DO I=N,1,-1
DO IJ=JU(I),IM(I+1)
J=MCU(IJ)
IF ((J.GT.0).AND.(J.LE.N)) XOUT(I)=XOUT(I)
1 -ILUCF(IJ)*XOUT(J)
ENDDO
XOUT(I)=ILUDF(I)*XOUT(I)
ENDDO
ELSE
*---
* FORWARD SOLVE (TRANSPOSED LINEAR SYSTEM)
*---
DO I=1,N
XOUT(I)=XIN(I)
ENDDO
DO I=1,N
DO IJ=JU(I),IM(I+1)
J=MCU(IJ)
IF ((J.GT.0).AND.(J.LE.N)) XOUT(J)=XOUT(J)
1 -ILUDF(I)*ILUCF(IJ)*XOUT(I)
ENDDO
ENDDO
*---
* BACKWARD SOLVE (TRANSPOSED LINEAR SYSTEM)
*---
DO I=N,1,-1
DO IJ=IM(I)+1,JU(I)-1
J=MCU(IJ)
IF ((J.GT.0).AND.(J.LE.N)) XOUT(J)=XOUT(J)
1 -ILUCF(IJ)*XOUT(I)
ENDDO
ENDDO
ENDIF
*
RETURN
END
|
