1c----------------------------------------------------------------------c
2       subroutine amub (nrow,ncol,job,a,ja,ia,b,jb,ib,
3     *                  c,jc,ic,nzmax,iw,ierr)
4      double precision a(*), b(*), c(*)
5      integer ja(*),jb(*),jc(*),ia(nrow+1),ib(*),ic(*),iw(ncol)
6c-----------------------------------------------------------------------
7c performs the matrix by matrix product C = A B
8c-----------------------------------------------------------------------
9c on entry:
10c ---------
11c nrow  = integer. The row dimension of A = row dimension of C
12c ncol  = integer. The column dimension of B = column dimension of C
13c job   = integer. Job indicator. When job = 0, only the structure
14c                  (i.e. the arrays jc, ic) is computed and the
15c                  real values are ignored.
16c
17c a,
18c ja,
19c ia   = Matrix A in compressed sparse row format.
20c
21c b,
22c jb,
23c ib    =  Matrix B in compressed sparse row format.
24c
25c nzmax = integer. The  length of the arrays c and jc.
26c         amub will stop if the result matrix C  has a number
27c         of elements that exceeds exceeds nzmax. See ierr.
28c
29c on return:
30c----------
31c c,
32c jc,
33c ic    = resulting matrix C in compressed sparse row sparse format.
34c
35c ierr  = integer. serving as error message.
36c         ierr = 0 means normal return,
37c         ierr .gt. 0 means that amub stopped while computing the
38c         i-th row  of C with i=ierr, because the number
39c         of elements in C exceeds nzmax.
40c
41c work arrays:
42c------------
43c iw    = integer work array of length equal to the number of
44c         columns in B.
45c Note:
46c-------
47c   The row dimension of B is not needed. However there is no checking
48c   on the condition that ncol(A) = nrow(B).
49c
50c-----------------------------------------------------------------------
51      double precision scal
52      logical values
53      values = (job .ne. 0)
54      len = 0
55      ic(1) = 1
56      ierr = 0
57c     initialize array iw.
58      do 1 j=1, ncol
59         iw(j) = 0
60 1    continue
61c
62      do 500 ii=1, nrow
63c     row i
64         do 200 ka=ia(ii), ia(ii+1)-1
65            if (values) scal = a(ka)
66            jj   = ja(ka)
67            do 100 kb=ib(jj),ib(jj+1)-1
68               jcol = jb(kb)
69               jpos = iw(jcol)
70               if (jpos .eq. 0) then
71                  len = len+1
72                  if (len .gt. nzmax) then
73                     ierr = ii
74                     return
75                  endif
76                  jc(len) = jcol
77                  iw(jcol)= len
78                  if (values) c(len)  = scal*b(kb)
79               else
80                  if (values) c(jpos) = c(jpos) + scal*b(kb)
81               endif
82 100     continue
83 200     continue
84         do 201 k=ic(ii), len
85            iw(jc(k)) = 0
86 201     continue
87         ic(ii+1) = len+1
88 500  continue
89      return
90c-------------end-of-amub-----------------------------------------------
91c-----------------------------------------------------------------------
92      end
93c-----------------------------------------------------------------------
94      subroutine amudia (nrow,job, a, ja, ia, diag, b, jb, ib)
95      double precision a(*), b(*), diag(*)
96      integer ja(*),jb(*), ia(nrow+1),ib(nrow+1)
97c-----------------------------------------------------------------------
98c performs the matrix by matrix product B = A * Diag  (in place)
99c-----------------------------------------------------------------------
100c on entry:
101c ---------
102c nrow	= integer. The row dimension of A
103c
104c job   = integer. job indicator. Job=0 means get array b only
105c         job = 1 means get b, and the integer arrays ib, jb.
106c
107c a,
108c ja,
109c ia   = Matrix A in compressed sparse row format.
110c
111c diag = diagonal matrix stored as a vector dig(1:n)
112c
113c on return:
114c----------
115c
116c b,
117c jb,
118c ib	= resulting matrix B in compressed sparse row sparse format.
119c
120c Notes:
121c-------
122c 1)        The column dimension of A is not needed.
123c 2)        algorithm in place (B can take the place of A).
124c-----------------------------------------------------------------
125      do 1 ii=1,nrow
126c
127c     scale each element
128c
129         k1 = ia(ii)
130         k2 = ia(ii+1)-1
131         do 2 k=k1, k2
132            b(k) = a(k)*diag(ja(k))
133 2       continue
134 1    continue
135c
136      if (job .eq. 0) return
137c
138      do 3 ii=1, nrow+1
139         ib(ii) = ia(ii)
140 3    continue
141      do 31 k=ia(1), ia(nrow+1) -1
142         jb(k) = ja(k)
143 31   continue
144      return
145c-----------------------------------------------------------------------
146c-----------end-of-amudiag----------------------------------------------
147      end
148c----------------------------------------------------------------------c
149      subroutine amux (n, x, y, a,ja,ia)
150      double precision  x(*), y(*), a(*)
151      integer n, ja(*), ia(*)
152c-----------------------------------------------------------------------
153c         A times a vector
154c-----------------------------------------------------------------------
155c multiplies a matrix by a vector using the dot product form
156c Matrix A is stored in compressed sparse row storage.
157c
158c on entry:
159c----------
160c n     = row dimension of A
161c x     = real array of length equal to the column dimension of
162c         the A matrix.
163c a, ja,
164c    ia = input matrix in compressed sparse row format.
165c
166c on return:
167c-----------
168c y     = real array of length n, containing the product y=Ax
169c
170c-----------------------------------------------------------------------
171c local variables
172c
173      double precision t
174      integer i, k
175c-----------------------------------------------------------------------
176      do 100 i = 1,n
177c
178c     compute the inner product of row i with vector x
179c
180         t = 0.0d0
181         do 99 k=ia(i), ia(i+1)-1
182            t = t + a(k)*x(ja(k))
183 99      continue
184c
185c     store result in y(i)
186c
187         y(i) = t
188 100  continue
189c
190      return
191      end
192c-----------------------------------------------------------------------
193      subroutine aplb (nrow,ncol,job,a,ja,ia,b,jb,ib,
194     *     c,jc,ic,nzmax,iw,ierr)
195      double precision a(*), b(*), c(*)
196      integer ja(*),jb(*),jc(*),ia(nrow+1),ib(nrow+1),ic(nrow+1),
197     *     iw(ncol)
198c-----------------------------------------------------------------------
199c performs the matrix sum  C = A+B.
200c-----------------------------------------------------------------------
201c on entry:
202c ---------
203c nrow	= integer. The row dimension of A and B
204c ncol  = integer. The column dimension of A and B.
205c job   = integer. Job indicator. When job = 0, only the structure
206c                  (i.e. the arrays jc, ic) is computed and the
207c                  real values are ignored.
208c
209c a,
210c ja,
211c ia   = Matrix A in compressed sparse row format.
212c
213c b,
214c jb,
215c ib	=  Matrix B in compressed sparse row format.
216c
217c nzmax	= integer. The  length of the arrays c and jc.
218c         amub will stop if the result matrix C  has a number
219c         of elements that exceeds exceeds nzmax. See ierr.
220c
221c on return:
222c----------
223c c,
224c jc,
225c ic	= resulting matrix C in compressed sparse row sparse format.
226c
227c ierr	= integer. serving as error message.
228c         ierr = 0 means normal return,
229c         ierr .gt. 0 means that amub stopped while computing the
230c         i-th row  of C with i=ierr, because the number
231c         of elements in C exceeds nzmax.
232c
233c work arrays:
234c------------
235c iw	= integer work array of length equal to the number of
236c         columns in A.
237c
238c-----------------------------------------------------------------------
239      logical values
240      values = (job .ne. 0)
241      ierr = 0
242      len = 0
243      ic(1) = 1
244      do 1 j=1, ncol
245         iw(j) = 0
246 1    continue
247c
248      do 500 ii=1, nrow
249c     row i
250         do 200 ka=ia(ii), ia(ii+1)-1
251            len = len+1
252            jcol    = ja(ka)
253            if (len .gt. nzmax) goto 999
254            jc(len) = jcol
255            if (values) c(len)  = a(ka)
256            iw(jcol)= len
257 200     continue
258c
259         do 300 kb=ib(ii),ib(ii+1)-1
260            jcol = jb(kb)
261            jpos = iw(jcol)
262            if (jpos .eq. 0) then
263               len = len+1
264               if (len .gt. nzmax) goto 999
265               jc(len) = jcol
266               if (values) c(len)  = b(kb)
267               iw(jcol)= len
268            else
269               if (values) c(jpos) = c(jpos) + b(kb)
270            endif
271 300     continue
272         do 301 k=ic(ii), len
273            iw(jc(k)) = 0
274 301     continue
275         ic(ii+1) = len+1
276 500  continue
277      return
278 999  ierr = ii
279      return
280c------------end of aplb -----------------------------------------------
281c-----------------------------------------------------------------------
282      end
283      subroutine csrmsr (n,a,ja,ia,ao,jao,wk,iwk,nnzao,ierr)
284      double precision a(*),ao(*),wk(n)
285      integer ia(n+1),ja(*),jao(*),iwk(n+1),nnzao,ierr
286c-----------------------------------------------------------------------
287c Compressed Sparse Row   to      Modified - Sparse Row
288c                                 Sparse row with separate main diagonal
289c-----------------------------------------------------------------------
290c converts a general sparse matrix a, ja, ia into
291c a compressed matrix using a separated diagonal (referred to as
292c the bell-labs format as it is used by bell labs semi conductor
293c group. We refer to it here as the modified sparse row format.
294c Note: this has been coded in such a way that one can overwrite
295c the output matrix onto the input matrix if desired by a call of
296c the form
297c
298c     call csrmsr (n, a, ja, ia, a, ja, wk,iwk)
299c
300c In case ao, jao, are different from a, ja, then one can
301c use ao, jao as the work arrays in the calling sequence:
302c
303c     call csrmsr (n, a, ja, ia, ao, jao, ao,jao)
304c
305c-----------------------------------------------------------------------
306c
307c on entry :
308c---------
309c a, ja, ia = matrix in csr format. note that the
310c	     algorithm is in place: ao, jao can be the same
311c            as a, ja, in which case it will be overwritten on it
312c            upon return.
313c nnzao = the number of non-zero entries in ao, jao
314c
315c on return :
316c-----------
317c
318c ao, jao  = sparse matrix in modified sparse row storage format:
319c	   +  ao(1:n) contains the diagonal of the matrix.
320c	   +  ao(n+2:nnz) contains the nondiagonal elements of the
321c             matrix, stored rowwise.
322c	   +  jao(n+2:nnz) : their column indices
323c	   +  jao(1:n+1) contains the pointer array for the nondiagonal
324c             elements in ao(n+2:nnz) and jao(n+2:nnz).
325c             i.e., for i .le. n+1 jao(i) points to beginning of row i
326c	      in arrays ao, jao.
327c	       here nnz = number of nonzero elements+1
328c ierr:
329c    = -1 : length of ao, jao < itpr
330c
331c work arrays:
332c------------
333c wk	= real work array of length n
334c iwk   = integer work array of length n+1
335c
336c notes:
337c-------
338c        Algorithm is in place.  i.e. both:
339c
340c          call csrmsr (n, a, ja, ia, ao, jao, ao,jao)
341c          (in which  ao, jao, are different from a, ja)
342c           and
343c          call csrmsr (n, a, ja, ia, a, ja, wk,iwk)
344c          (in which  wk, jwk, are different from a, ja)
345c        are OK.
346c--------
347c coded by Y. Saad Sep. 1989. Rechecked Feb 27, 1990.
348c Modified by Pin Ng on June 11, 2002 to provide warning when
349c    iptr > nnzao+1
350c-----------------------------------------------------------------------
351      icount = 0
352c
353c store away diagonal elements and count nonzero diagonal elements.
354c
355      do 1 i=1,n
356         wk(i) = 0.0d0
357         iwk(i+1) = ia(i+1)-ia(i)
358         do 2 k=ia(i),ia(i+1)-1
359            if (ja(k) .eq. i) then
360               wk(i) = a(k)
361               icount = icount + 1
362               iwk(i+1) = iwk(i+1)-1
363            endif
364 2       continue
365 1    continue
366c
367c compute total length
368c
369      iptr = n + ia(n+1) - icount
370      if (iptr .gt. nnzao+1) then
371         ierr = -1
372         return
373      endif
374c
375c     copy backwards (to avoid collisions)
376c
377      do 500 ii=n,1,-1
378         do 100 k=ia(ii+1)-1,ia(ii),-1
379            j = ja(k)
380            if (j .ne. ii) then
381               ao(iptr) = a(k)
382               jao(iptr) = j
383               iptr = iptr-1
384            endif
385 100     continue
386 500  continue
387c
388c compute pointer values and copy wk(*)
389c
390      jao(1) = n+2
391      do 600 i=1,n
392         ao(i) = wk(i)
393         jao(i+1) = jao(i)+iwk(i+1)
394 600  continue
395      return
396c------------ end of subroutine csrmsr ---------------------------------
397c-----------------------------------------------------------------------
398      end
399