ilu0.f 4.73 KB
Newer Older
Kirill Terekhov's avatar
Kirill Terekhov committed
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
c-----------------------------------------------------------------------
c This routine solves the system (LU) y = x
c-----------------------------------------------------------------------
        SUBROUTINE prevec0( iprevec, dummy, x, y, iwork, dwork)
        INTEGER   iprevec(*), dummy, iwork(*), n
        REAL*8    x(*),y(*),dwork(*)
        n = iprevec(1)
        call lusol(n, x, y, dwork(1), iwork(n+2), iwork(1))
        return
        end

	subroutine ilu0(n, a, ja, ia, alu, jlu, ju, iw, ierr)
	implicit real*8 (a-h,o-z)
	real*8 a(*), alu(*)
        integer ja(*), ia(*), ju(*), jlu(*), iw(*)
c input:
c---------
c n       = dimension of matrix
c a, ja,
c ia      = original matrix in compressed sparse row storage.
c
c output:   
c alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
c           the L and U factors together. The diagonal (stored in
c           alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix
c           contains the i-th row of L (excluding the diagonal entry=1)
c           followed by the i-th row of U.
c
c ju	  = pointer to the diagonal elements in alu, jlu.
c
c ierr	  = integer indicating error code on return
c	     ierr = 0 --> normal return
c	     ierr = k --> code encountered a zero pivot at step k.
c work arrays:
c-------------
c iw	    = integer work array of length n.
c-----------------------------------------------------------------------
      integer maxbuf, nbuf
      parameter (maxbuf=1000) ! max number of non-zeros in a row
      integer jbuf(maxbuf)
      real*8  abuf(maxbuf)
c-----------------------------------------------------------------------
        ju0 = n+2
        jlu(1) = ju0
c
c initialize work vector to zero's
c
	do 31 i=1, n
           iw(i) = 0
 31     continue
c
c main loop
c
	do 500 ii = 1, n
           js = ju0
c
c sort column indexes
c
           nbuf=0
           do j=ia(ii),ia(ii+1)-1
            nbuf=nbuf+1
            abuf(nbuf)=a(j)
            jbuf(nbuf)=ja(j)
           end do
           if (nbuf.gt.maxbuf) stop 'ilu0: maxbuf too small'
           call dsortilu (jbuf, abuf, nbuf, 2, ierr)

c generating row number ii of L and U.
c
           do 100 j=1,nbuf  ! j=ia(ii),ia(ii+1)-1
c
c     copy row ii of a, ja, ia into row ii of alu, jlu (L/U) matrix.
c
              jcol = jbuf(j)  ! ja(j)
              if (jcol .eq. ii) then
                 alu(ii) = abuf(j)  ! a(j)
                 iw(jcol) = ii
                 ju(ii)  = ju0
              else
                 alu(ju0) = abuf(j)  ! a(j)
                 jlu(ju0) = jbuf(j)  ! ja(j)
                 iw(jcol) = ju0
                 ju0 = ju0+1
              endif
 100       continue
           jlu(ii+1) = ju0
           jf = ju0-1
           jm = ju(ii)-1
c
c     exit if diagonal element is reached.
c
           do 150 j=js, jm
              jrow = jlu(j)
              tl = alu(j)*alu(jrow)
              alu(j) = tl
c
c     perform  linear combination
c
              do 140 jj = ju(jrow), jlu(jrow+1)-1
                 jw = iw(jlu(jj))
                 if (jw .ne. 0) alu(jw) = alu(jw) - tl*alu(jj)
 140          continue
 150       continue
c
c     invert  and store diagonal element.
c
           if (alu(ii) .eq. 0.0d0) goto 600
           alu(ii) = 1.0d0/alu(ii)
c
c     reset pointer iw to zero
c
           iw(ii) = 0
           do 201 i = js, jf
 201          iw(jlu(i)) = 0
 500       continue
           ierr = 0
           return
c
c     zero pivot :
c
 600       ierr = ii
c
           return
           end

        subroutine lusol(n, y, x, alu, jlu, ju)
        real*8 x(n), y(n), alu(*)
        integer n, jlu(*), ju(*)
c-----------------------------------------------------------------------
c
c This routine solves the system (LU) x = y,
c given an LU decomposition of a matrix stored in (alu, jlu, ju)
c modified sparse row format
c
c-----------------------------------------------------------------------
c on entry:
c n   = dimension of system
c y   = the right-hand-side vector
c alu, jlu, ju
c     = the LU matrix as provided from the ILU routines.
c
c on return
c x   = solution of LU x = y.
c-----------------------------------------------------------------------
c
c Note: routine is in place: call lusol (n, x, x, alu, jlu, ju)
c       will solve the system with rhs x and overwrite the result on x .
c
c-----------------------------------------------------------------------
c local variables
c
        integer i,k
c
c forward solve
c
        do 40 i = 1, n
           x(i) = y(i)
           do 41 k=jlu(i),ju(i)-1
              x(i) = x(i) - alu(k)* x(jlu(k))
 41        continue
 40     continue
c
c     backward solve.
c
        do 90 i = n, 1, -1
           do 91 k=ju(i),jlu(i+1)-1
              x(i) = x(i) - alu(k)*x(jlu(k))
 91        continue
           x(i) = alu(i)*x(i)
 90     continue
c
        return
        end