1c $Id: blassm.f,v 1.1 2008-04-11 06:01:06 geuzaine Exp $
2c----------------------------------------------------------------------c
3c                          S P A R S K I T                             c
4c----------------------------------------------------------------------c
5c        BASIC LINEAR ALGEBRA FOR SPARSE MATRICES. BLASSM MODULE       c
6c----------------------------------------------------------------------c
7c amub   :   computes     C = A*B                                      c
8c aplb   :   computes     C = A+B                                      c
9c aplb1  :   computes     C = A+B  [Sorted version: A, B, C sorted]    c
10c aplsb  :   computes     C = A + s B                                  c
11c aplsb1 :   computes     C = A+sB  [Sorted version: A, B, C sorted]   c
12c apmbt  :   Computes     C = A +/- transp(B)                          c
13c aplsbt :   Computes     C = A + s * transp(B)                        c
14c diamua :   Computes     C = Diag * A                                 c
15c amudia :   Computes     C = A* Diag                                  c
16c aplsca :   Computes     A:= A + s I    (s = scalar)                  c
17c apldia :   Computes     C = A + Diag.                                c
18c----------------------------------------------------------------------c
19c Note: this module still incomplete.                                  c
20c----------------------------------------------------------------------c
21       subroutine amub (nrow,ncol,job,a,ja,ia,b,jb,ib,
22     *                  c,jc,ic,nzmax,iw,ierr)
23      real*8 a(*), b(*), c(*)
24      integer ja(*),jb(*),jc(*),ia(nrow+1),ib(*),ic(*),iw(ncol)
25c-----------------------------------------------------------------------
26c performs the matrix by matrix product C = A B
27c-----------------------------------------------------------------------
28c on entry:
29c ---------
30c nrow  = integer. The row dimension of A = row dimension of C
31c ncol  = integer. The column dimension of B = column dimension of C
32c job   = integer. Job indicator. When job = 0, only the structure
33c                  (i.e. the arrays jc, ic) is computed and the
34c                  real values are ignored.
35c
36c a,
37c ja,
38c ia   = Matrix A in compressed sparse row format.
39c
40c b,
41c jb,
42c ib    =  Matrix B in compressed sparse row format.
43c
44c nzmax = integer. The  length of the arrays c and jc.
45c         amub will stop if the result matrix C  has a number
46c         of elements that exceeds exceeds nzmax. See ierr.
47c
48c on return:
49c----------
50c c,
51c jc,
52c ic    = resulting matrix C in compressed sparse row sparse format.
53c
54c ierr  = integer. serving as error message.
55c         ierr = 0 means normal return,
56c         ierr .gt. 0 means that amub stopped while computing the
57c         i-th row  of C with i=ierr, because the number
58c         of elements in C exceeds nzmax.
59c
60c work arrays:
61c------------
62c iw    = integer work array of length equal to the number of
63c         columns in A.
64c Note:
65c-------
66c   The row dimension of B is not needed. However there is no checking
67c   on the condition that ncol(A) = nrow(B).
68c
69c-----------------------------------------------------------------------
70      real*8 scal
71      logical values
72      values = (job .ne. 0)
73      len = 0
74      ic(1) = 1
75      ierr = 0
76c     initialize array iw.
77      do 1 j=1, ncol
78         iw(j) = 0
79 1    continue
80c
81      do 500 ii=1, nrow
82c     row i
83         do 200 ka=ia(ii), ia(ii+1)-1
84	    if (values) scal = a(ka)
85	    jj   = ja(ka)
86	    do 100 kb=ib(jj),ib(jj+1)-1
87               jcol = jb(kb)
88               jpos = iw(jcol)
89               if (jpos .eq. 0) then
90                  len = len+1
91                  if (len .gt. nzmax) then
92                     ierr = ii
93                     return
94                  endif
95                  jc(len) = jcol
96                  iw(jcol)= len
97                  if (values) c(len)  = scal*b(kb)
98               else
99                  if (values) c(jpos) = c(jpos) + scal*b(kb)
100               endif
101 100	    continue
102 200     continue
103         do 201 k=ic(ii), len
104	    iw(jc(k)) = 0
105 201     continue
106         ic(ii+1) = len+1
107 500  continue
108      return
109c-------------end-of-amub-----------------------------------------------
110c-----------------------------------------------------------------------
111      end
112c-----------------------------------------------------------------------
113      subroutine aplb (nrow,ncol,job,a,ja,ia,b,jb,ib,
114     *     c,jc,ic,nzmax,iw,ierr)
115      real*8 a(*), b(*), c(*)
116      integer ja(*),jb(*),jc(*),ia(nrow+1),ib(nrow+1),ic(nrow+1),
117     *     iw(ncol)
118c-----------------------------------------------------------------------
119c performs the matrix sum  C = A+B.
120c-----------------------------------------------------------------------
121c on entry:
122c ---------
123c nrow	= integer. The row dimension of A and B
124c ncol  = integer. The column dimension of A and B.
125c job   = integer. Job indicator. When job = 0, only the structure
126c                  (i.e. the arrays jc, ic) is computed and the
127c                  real values are ignored.
128c
129c a,
130c ja,
131c ia   = Matrix A in compressed sparse row format.
132c
133c b,
134c jb,
135c ib	=  Matrix B in compressed sparse row format.
136c
137c nzmax	= integer. The  length of the arrays c and jc.
138c         amub will stop if the result matrix C  has a number
139c         of elements that exceeds exceeds nzmax. See ierr.
140c
141c on return:
142c----------
143c c,
144c jc,
145c ic	= resulting matrix C in compressed sparse row sparse format.
146c
147c ierr	= integer. serving as error message.
148c         ierr = 0 means normal return,
149c         ierr .gt. 0 means that amub stopped while computing the
150c         i-th row  of C with i=ierr, because the number
151c         of elements in C exceeds nzmax.
152c
153c work arrays:
154c------------
155c iw	= integer work array of length equal to the number of
156c         columns in A.
157c
158c-----------------------------------------------------------------------
159      logical values
160      values = (job .ne. 0)
161      ierr = 0
162      len = 0
163      ic(1) = 1
164      do 1 j=1, ncol
165         iw(j) = 0
166 1    continue
167c
168      do 500 ii=1, nrow
169c     row i
170         do 200 ka=ia(ii), ia(ii+1)-1
171            len = len+1
172            jcol    = ja(ka)
173            if (len .gt. nzmax) goto 999
174            jc(len) = jcol
175            if (values) c(len)  = a(ka)
176            iw(jcol)= len
177 200     continue
178c
179         do 300 kb=ib(ii),ib(ii+1)-1
180            jcol = jb(kb)
181            jpos = iw(jcol)
182            if (jpos .eq. 0) then
183               len = len+1
184               if (len .gt. nzmax) goto 999
185               jc(len) = jcol
186               if (values) c(len)  = b(kb)
187               iw(jcol)= len
188            else
189               if (values) c(jpos) = c(jpos) + b(kb)
190            endif
191 300     continue
192         do 301 k=ic(ii), len
193	    iw(jc(k)) = 0
194 301     continue
195         ic(ii+1) = len+1
196 500  continue
197      return
198 999  ierr = ii
199      return
200c------------end of aplb -----------------------------------------------
201c-----------------------------------------------------------------------
202      end
203c-----------------------------------------------------------------------
204      subroutine aplb1(nrow,ncol,job,a,ja,ia,b,jb,ib,c,jc,ic,nzmax,ierr)
205      real*8 a(*), b(*), c(*)
206      integer ja(*),jb(*),jc(*),ia(nrow+1),ib(nrow+1),ic(nrow+1)
207c-----------------------------------------------------------------------
208c performs the matrix sum  C = A+B for matrices in sorted CSR format.
209c the difference with aplb  is that the resulting matrix is such that
210c the elements of each row are sorted with increasing column indices in
211c each row, provided the original matrices are sorted in the same way.
212c-----------------------------------------------------------------------
213c on entry:
214c ---------
215c nrow	= integer. The row dimension of A and B
216c ncol  = integer. The column dimension of A and B.
217c job   = integer. Job indicator. When job = 0, only the structure
218c                  (i.e. the arrays jc, ic) is computed and the
219c                  real values are ignored.
220c
221c a,
222c ja,
223c ia   = Matrix A in compressed sparse row format with entries sorted
224c
225c b,
226c jb,
227c ib	=  Matrix B in compressed sparse row format with entries sorted
228c        ascendly in each row
229c
230c nzmax	= integer. The  length of the arrays c and jc.
231c         amub will stop if the result matrix C  has a number
232c         of elements that exceeds exceeds nzmax. See ierr.
233c
234c on return:
235c----------
236c c,
237c jc,
238c ic	= resulting matrix C in compressed sparse row sparse format
239c         with entries sorted ascendly in each row.
240c
241c ierr	= integer. serving as error message.
242c         ierr = 0 means normal return,
243c         ierr .gt. 0 means that amub stopped while computing the
244c         i-th row  of C with i=ierr, because the number
245c         of elements in C exceeds nzmax.
246c
247c Notes:
248c-------
249c     this will not work if any of the two input matrices is not sorted
250c-----------------------------------------------------------------------
251      logical values
252      values = (job .ne. 0)
253      ierr = 0
254      kc = 1
255      ic(1) = kc
256c
257      do 6 i=1, nrow
258         ka = ia(i)
259         kb = ib(i)
260         kamax = ia(i+1)-1
261         kbmax = ib(i+1)-1
262 5       continue
263         if (ka .le. kamax) then
264            j1 = ja(ka)
265         else
266            j1 = ncol+1
267         endif
268         if (kb .le. kbmax) then
269            j2 = jb(kb)
270         else
271            j2 = ncol+1
272         endif
273c
274c     three cases
275c
276         if (j1 .eq. j2) then
277            if (values) c(kc) = a(ka)+b(kb)
278            jc(kc) = j1
279            ka = ka+1
280            kb = kb+1
281            kc = kc+1
282         else if (j1 .lt. j2) then
283            jc(kc) = j1
284            if (values) c(kc) = a(ka)
285            ka = ka+1
286            kc = kc+1
287         else if (j1 .gt. j2) then
288            jc(kc) = j2
289            if (values) c(kc) = b(kb)
290            kb = kb+1
291            kc = kc+1
292         endif
293         if (kc .gt. nzmax) goto 999
294         if (ka .le. kamax .or. kb .le. kbmax) goto 5
295         ic(i+1) = kc
296 6    continue
297      return
298 999  ierr = i
299      return
300c------------end-of-aplb1-----------------------------------------------
301c-----------------------------------------------------------------------
302      end
303c-----------------------------------------------------------------------
304      subroutine aplsb (nrow,ncol,a,ja,ia,s,b,jb,ib,c,jc,ic,
305     *     nzmax,ierr)
306      real*8 a(*), b(*), c(*), s
307      integer ja(*),jb(*),jc(*),ia(nrow+1),ib(nrow+1),ic(nrow+1)
308c-----------------------------------------------------------------------
309c performs the operation C = A+s B for matrices in sorted CSR format.
310c the difference with aplsb is that the resulting matrix is such that
311c the elements of each row are sorted with increasing column indices in
312c each row, provided the original matrices are sorted in the same way.
313c-----------------------------------------------------------------------
314c on entry:
315c ---------
316c nrow	= integer. The row dimension of A and B
317c ncol  = integer. The column dimension of A and B.
318c
319c a,
320c ja,
321c ia   = Matrix A in compressed sparse row format with entries sorted
322c
323c s	= real. scalar factor for B.
324c
325c b,
326c jb,
327c ib	=  Matrix B in compressed sparse row format with entries sorted
328c        ascendly in each row
329c
330c nzmax	= integer. The  length of the arrays c and jc.
331c         amub will stop if the result matrix C  has a number
332c         of elements that exceeds exceeds nzmax. See ierr.
333c
334c on return:
335c----------
336c c,
337c jc,
338c ic	= resulting matrix C in compressed sparse row sparse format
339c         with entries sorted ascendly in each row.
340c
341c ierr	= integer. serving as error message.
342c         ierr = 0 means normal return,
343c         ierr .gt. 0 means that amub stopped while computing the
344c         i-th row  of C with i=ierr, because the number
345c         of elements in C exceeds nzmax.
346c
347c Notes:
348c-------
349c     this will not work if any of the two input matrices is not sorted
350c-----------------------------------------------------------------------
351      ierr = 0
352      kc = 1
353      ic(1) = kc
354c
355c     the following loop does a merge of two sparse rows + adds  them.
356c
357      do 6 i=1, nrow
358         ka = ia(i)
359         kb = ib(i)
360         kamax = ia(i+1)-1
361         kbmax = ib(i+1)-1
362 5       continue
363c
364c     this is a while  -- do loop --
365c
366         if (ka .le. kamax .or. kb .le. kbmax) then
367c
368            if (ka .le. kamax) then
369               j1 = ja(ka)
370            else
371c     take j1 large enough  that always j2 .lt. j1
372               j1 = ncol+1
373            endif
374            if (kb .le. kbmax) then
375               j2 = jb(kb)
376            else
377c     similarly take j2 large enough  that always j1 .lt. j2
378               j2 = ncol+1
379            endif
380c
381c     three cases
382c
383            if (j1 .eq. j2) then
384               c(kc) = a(ka)+s*b(kb)
385               jc(kc) = j1
386               ka = ka+1
387               kb = kb+1
388               kc = kc+1
389            else if (j1 .lt. j2) then
390               jc(kc) = j1
391               c(kc) = a(ka)
392               ka = ka+1
393               kc = kc+1
394            else if (j1 .gt. j2) then
395               jc(kc) = j2
396               c(kc) = s*b(kb)
397               kb = kb+1
398               kc = kc+1
399            endif
400            if (kc .gt. nzmax) goto 999
401            goto 5
402c
403c     end while loop
404c
405         endif
406         ic(i+1) = kc
407 6    continue
408      return
409 999  ierr = i
410      return
411c------------end-of-aplsb ---------------------------------------------
412c-----------------------------------------------------------------------
413      end
414c-----------------------------------------------------------------------
415      subroutine aplsb1 (nrow,ncol,a,ja,ia,s,b,jb,ib,c,jc,ic,
416     *     nzmax,ierr)
417      real*8 a(*), b(*), c(*), s
418      integer ja(*),jb(*),jc(*),ia(nrow+1),ib(nrow+1),ic(nrow+1)
419c-----------------------------------------------------------------------
420c performs the operation C = A+s B for matrices in sorted CSR format.
421c the difference with aplsb is that the resulting matrix is such that
422c the elements of each row are sorted with increasing column indices in
423c each row, provided the original matrices are sorted in the same way.
424c-----------------------------------------------------------------------
425c on entry:
426c ---------
427c nrow	= integer. The row dimension of A and B
428c ncol  = integer. The column dimension of A and B.
429c
430c a,
431c ja,
432c ia   = Matrix A in compressed sparse row format with entries sorted
433c
434c s	= real. scalar factor for B.
435c
436c b,
437c jb,
438c ib	=  Matrix B in compressed sparse row format with entries sorted
439c        ascendly in each row
440c
441c nzmax	= integer. The  length of the arrays c and jc.
442c         amub will stop if the result matrix C  has a number
443c         of elements that exceeds exceeds nzmax. See ierr.
444c
445c on return:
446c----------
447c c,
448c jc,
449c ic	= resulting matrix C in compressed sparse row sparse format
450c         with entries sorted ascendly in each row.
451c
452c ierr	= integer. serving as error message.
453c         ierr = 0 means normal return,
454c         ierr .gt. 0 means that amub stopped while computing the
455c         i-th row  of C with i=ierr, because the number
456c         of elements in C exceeds nzmax.
457c
458c Notes:
459c-------
460c     this will not work if any of the two input matrices is not sorted
461c-----------------------------------------------------------------------
462      ierr = 0
463      kc = 1
464      ic(1) = kc
465c
466c     the following loop does a merge of two sparse rows + adds  them.
467c
468      do 6 i=1, nrow
469         ka = ia(i)
470         kb = ib(i)
471         kamax = ia(i+1)-1
472         kbmax = ib(i+1)-1
473 5       continue
474c
475c     this is a while  -- do loop --
476c
477         if (ka .le. kamax .or. kb .le. kbmax) then
478c
479            if (ka .le. kamax) then
480               j1 = ja(ka)
481            else
482c     take j1 large enough  that always j2 .lt. j1
483               j1 = ncol+1
484            endif
485            if (kb .le. kbmax) then
486               j2 = jb(kb)
487            else
488c     similarly take j2 large enough  that always j1 .lt. j2
489               j2 = ncol+1
490            endif
491c
492c     three cases
493c
494            if (j1 .eq. j2) then
495               c(kc) = a(ka)+s*b(kb)
496               jc(kc) = j1
497               ka = ka+1
498               kb = kb+1
499               kc = kc+1
500            else if (j1 .lt. j2) then
501               jc(kc) = j1
502               c(kc) = a(ka)
503               ka = ka+1
504               kc = kc+1
505            else if (j1 .gt. j2) then
506               jc(kc) = j2
507               c(kc) = s*b(kb)
508               kb = kb+1
509               kc = kc+1
510            endif
511            if (kc .gt. nzmax) goto 999
512            goto 5
513c
514c     end while loop
515c
516         endif
517         ic(i+1) = kc
518 6    continue
519      return
520 999  ierr = i
521      return
522c------------end-of-aplsb1 ---------------------------------------------
523c-----------------------------------------------------------------------
524      end
525c-----------------------------------------------------------------------
526      subroutine apmbt (nrow,ncol,job,a,ja,ia,b,jb,ib,
527     *     c,jc,ic,nzmax,iw,ierr)
528      real*8 a(*), b(*), c(*)
529      integer ja(*),jb(*),jc(*),ia(nrow+1),ib(ncol+1),ic(*),iw(*)
530c-----------------------------------------------------------------------
531c performs the matrix sum  C = A + transp(B) or C = A - transp(B)
532c-----------------------------------------------------------------------
533c on entry:
534c ---------
535c nrow	= integer. The row dimension of A and transp(B)
536c ncol  = integer. The column dimension of A. Also the row
537c                  dimension of B.
538c
539c job	= integer. if job = -1, apmbt will compute C= A - transp(B)
540c         (structure + values)
541c         if (job .eq. 1)  it will compute C=A+transp(A)
542c         (structure+ values)
543c         if (job .eq. 0) it will compute the structure of
544c         C= A+/-transp(B) only (ignoring all real values).
545c         any other value of job will be treated as  job=1
546c a,
547c ja,
548c ia    = Matrix A in compressed sparse row format.
549c
550c b,
551c jb,
552c ib	=  Matrix B in compressed sparse row format.
553c
554c nzmax	= integer. The  length of the arrays c, jc, and ic.
555c         amub will stop if the result matrix C  has a number
556c         of elements that exceeds exceeds nzmax. See ierr.
557c
558c on return:
559c----------
560c c,
561c jc,
562c ic	= resulting matrix C in compressed sparse row format.
563c
564c ierr	= integer. serving as error message.
565c         ierr = 0 means normal return.
566c         ierr = -1 means that nzmax was .lt. either the number of
567c         nonzero elements of A or the number of nonzero elements in B.
568c         ierr .gt. 0 means that amub stopped while computing the
569c         i-th row  of C with i=ierr, because the number
570c         of elements in C exceeds nzmax.
571c
572c work arrays:
573c------------
574c iw	= integer work array of length at least max(ncol,nrow)
575c
576c Notes:
577c------- It is important to note that here all of three arrays c, ic,
578c        and jc are assumed to be of length nnz(c). This is because
579c        the matrix is internally converted in coordinate format.
580c
581c-----------------------------------------------------------------------
582      logical values
583      values = (job .ne. 0)
584c
585      ierr = 0
586      do 1 j=1, ncol
587         iw(j) = 0
588 1    continue
589c
590      nnza = ia(nrow+1)-1
591      nnzb = ib(ncol+1)-1
592      len = nnzb
593      if (nzmax .lt. nnzb .or. nzmax .lt. nnza) then
594         ierr = -1
595         return
596      endif
597c
598c trasnpose matrix b into c
599c
600      ljob = 0
601      if (values) ljob = 1
602      ipos = 1
603      call csrcsc (ncol,ljob,ipos,b,jb,ib,c,jc,ic)
604c-----------------------------------------------------------------------
605      if (job .eq. -1) then
606         do 2 k=1,len
607	    c(k) = -c(k)
608 2       continue
609      endif
610c
611c--------------- main loop --------------------------------------------
612c
613      do 500 ii=1, nrow
614         do 200 k = ic(ii),ic(ii+1)-1
615            iw(jc(k)) = k
616 200     continue
617c-----------------------------------------------------------------------
618         do 300 ka = ia(ii), ia(ii+1)-1
619            jcol = ja(ka)
620            jpos = iw(jcol)
621            if (jpos .eq. 0) then
622c
623c     if fill-in append in coordinate format to matrix.
624c
625               len = len+1
626               if (len .gt. nzmax) goto 999
627               jc(len) = jcol
628
629               ic(len) = ii
630               if (values) c(len)  = a(ka)
631            else
632c     else do addition.
633               if (values) c(jpos) = c(jpos) + a(ka)
634            endif
635 300     continue
636         do 301 k=ic(ii), ic(ii+1)-1
637	    iw(jc(k)) = 0
638 301     continue
639 500  continue
640c
641c     convert first part of matrix (without fill-ins) into coo format
642c
643      ljob = 2
644      if (values) ljob = 3
645      do 501 i=1, nrow+1
646         iw(i) = ic(i)
647 501  continue
648      call csrcoo (nrow,ljob,nnzb,c,jc,iw,nnzb,c,ic,jc,ierr)
649c
650c     convert the whole thing back to csr format.
651c
652      ljob = 0
653      if (values) ljob = 1
654      call coicsr (nrow,len,ljob,c,jc,ic,iw)
655      return
656 999  ierr = ii
657      return
658c--------end-of-apmbt---------------------------------------------------
659c-----------------------------------------------------------------------
660      end
661c-----------------------------------------------------------------------
662      subroutine aplsbt(nrow,ncol,a,ja,ia,s,b,jb,ib,
663     *     c,jc,ic,nzmax,iw,ierr)
664      real*8 a(*), b(*), c(*), s
665      integer ja(*),jb(*),jc(*),ia(nrow+1),ib(ncol+1),ic(*),iw(*)
666c-----------------------------------------------------------------------
667c performs the matrix sum  C = A + transp(B).
668c-----------------------------------------------------------------------
669c on entry:
670c ---------
671c nrow	= integer. The row dimension of A and transp(B)
672c ncol  = integer. The column dimension of A. Also the row
673c                  dimension of B.
674c
675c a,
676c ja,
677c ia    = Matrix A in compressed sparse row format.
678c
679c s	= real. scalar factor for B.
680c
681c
682c b,
683c jb,
684c ib	=  Matrix B in compressed sparse row format.
685c
686c nzmax	= integer. The  length of the arrays c, jc, and ic.
687c         amub will stop if the result matrix C  has a number
688c         of elements that exceeds exceeds nzmax. See ierr.
689c
690c on return:
691c----------
692c c,
693c jc,
694c ic	= resulting matrix C in compressed sparse row format.
695c
696c ierr	= integer. serving as error message.
697c         ierr = 0 means normal return.
698c         ierr = -1 means that nzmax was .lt. either the number of
699c         nonzero elements of A or the number of nonzero elements in B.
700c         ierr .gt. 0 means that amub stopped while computing the
701c         i-th row  of C with i=ierr, because the number
702c         of elements in C exceeds nzmax.
703c
704c work arrays:
705c------------
706c iw	= integer work array of length at least max(nrow,ncol)
707c
708c Notes:
709c------- It is important to note that here all of three arrays c, ic,
710c        and jc are assumed to be of length nnz(c). This is because
711c        the matrix is internally converted in coordinate format.
712c
713c-----------------------------------------------------------------------
714      ierr = 0
715      do 1 j=1, ncol
716         iw(j) = 0
717 1    continue
718c
719      nnza = ia(nrow+1)-1
720      nnzb = ib(ncol+1)-1
721      len = nnzb
722      if (nzmax .lt. nnzb .or. nzmax .lt. nnza) then
723         ierr = -1
724         return
725      endif
726c
727c     transpose matrix b into c
728c
729      ljob = 1
730      ipos = 1
731      call csrcsc (ncol,ljob,ipos,b,jb,ib,c,jc,ic)
732      do 2 k=1,len
733 2       c(k) = c(k)*s
734c
735c     main loop. add rows from ii = 1 to nrow.
736c
737         do 500 ii=1, nrow
738c     iw is used as a system to recognize whether there
739c     was a nonzero element in c.
740            do 200 k = ic(ii),ic(ii+1)-1
741               iw(jc(k)) = k
742 200        continue
743c
744            do 300 ka = ia(ii), ia(ii+1)-1
745               jcol = ja(ka)
746               jpos = iw(jcol)
747           if (jpos .eq. 0) then
748c
749c     if fill-in append in coordinate format to matrix.
750c
751              len = len+1
752              if (len .gt. nzmax) goto 999
753              jc(len) = jcol
754              ic(len) = ii
755              c(len)  = a(ka)
756           else
757c     else do addition.
758              c(jpos) = c(jpos) + a(ka)
759           endif
760 300    continue
761        do 301 k=ic(ii), ic(ii+1)-1
762           iw(jc(k)) = 0
763 301    continue
764 500  continue
765c
766c     convert first part of matrix (without fill-ins) into coo format
767c
768      ljob = 3
769      do 501 i=1, nrow+1
770         iw(i) = ic(i)
771 501  continue
772      call csrcoo (nrow,ljob,nnzb,c,jc,iw,nnzb,c,ic,jc,ierr)
773c
774c     convert the whole thing back to csr format.
775c
776      ljob = 1
777      call coicsr (nrow,len,ljob,c,jc,ic,iw)
778      return
779 999  ierr = ii
780      return
781c--------end-of-aplsbt--------------------------------------------------
782c-----------------------------------------------------------------------
783      end
784c-----------------------------------------------------------------------
785      subroutine diamua (nrow,job, a, ja, ia, diag, b, jb, ib)
786      real*8 a(*), b(*), diag(nrow), scal
787      integer ja(*),jb(*), ia(nrow+1),ib(nrow+1)
788c-----------------------------------------------------------------------
789c performs the matrix by matrix product B = Diag * A  (in place)
790c-----------------------------------------------------------------------
791c on entry:
792c ---------
793c nrow	= integer. The row dimension of A
794c
795c job   = integer. job indicator. Job=0 means get array b only
796c         job = 1 means get b, and the integer arrays ib, jb.
797c
798c a,
799c ja,
800c ia   = Matrix A in compressed sparse row format.
801c
802c diag = diagonal matrix stored as a vector dig(1:n)
803c
804c on return:
805c----------
806c
807c b,
808c jb,
809c ib	= resulting matrix B in compressed sparse row sparse format.
810c
811c Notes:
812c-------
813c 1)        The column dimension of A is not needed.
814c 2)        algorithm in place (B can take the place of A).
815c           in this case use job=0.
816c-----------------------------------------------------------------
817      do 1 ii=1,nrow
818c
819c     normalize each row
820c
821         k1 = ia(ii)
822         k2 = ia(ii+1)-1
823         scal = diag(ii)
824         do 2 k=k1, k2
825            b(k) = a(k)*scal
826 2       continue
827 1    continue
828c
829      if (job .eq. 0) return
830c
831      do 3 ii=1, nrow+1
832         ib(ii) = ia(ii)
833 3    continue
834      do 31 k=ia(1), ia(nrow+1) -1
835         jb(k) = ja(k)
836 31   continue
837      return
838c----------end-of-diamua------------------------------------------------
839c-----------------------------------------------------------------------
840      end
841c-----------------------------------------------------------------------
842      subroutine amudia (nrow,job, a, ja, ia, diag, b, jb, ib)
843      real*8 a(*), b(*), diag(nrow)
844      integer ja(*),jb(*), ia(nrow+1),ib(nrow+1)
845c-----------------------------------------------------------------------
846c performs the matrix by matrix product B = A * Diag  (in place)
847c-----------------------------------------------------------------------
848c on entry:
849c ---------
850c nrow	= integer. The row dimension of A
851c
852c job   = integer. job indicator. Job=0 means get array b only
853c         job = 1 means get b, and the integer arrays ib, jb.
854c
855c a,
856c ja,
857c ia   = Matrix A in compressed sparse row format.
858c
859c diag = diagonal matrix stored as a vector dig(1:n)
860c
861c on return:
862c----------
863c
864c b,
865c jb,
866c ib	= resulting matrix B in compressed sparse row sparse format.
867c
868c Notes:
869c-------
870c 1)        The column dimension of A is not needed.
871c 2)        algorithm in place (B can take the place of A).
872c-----------------------------------------------------------------
873      do 1 ii=1,nrow
874c
875c     scale each element
876c
877         k1 = ia(ii)
878         k2 = ia(ii+1)-1
879         do 2 k=k1, k2
880            b(k) = a(k)*diag(ja(k))
881 2       continue
882 1    continue
883c
884      if (job .eq. 0) return
885c
886      do 3 ii=1, nrow+1
887         ib(ii) = ia(ii)
888 3    continue
889      do 31 k=ia(1), ia(nrow+1) -1
890         jb(k) = ja(k)
891 31   continue
892      return
893c-----------------------------------------------------------------------
894c-----------end-of-amudiag----------------------------------------------
895      end
896c-----------------------------------------------------------------------
897      subroutine aplsca (nrow, a, ja, ia, scal,iw)
898      real*8 a(*), scal
899      integer ja(*), ia(nrow+1),iw(*)
900c-----------------------------------------------------------------------
901c Adds a scalar to the diagonal entries of a sparse matrix A :=A + s I
902c-----------------------------------------------------------------------
903c on entry:
904c ---------
905c nrow	= integer. The row dimension of A
906c
907c a,
908c ja,
909c ia    = Matrix A in compressed sparse row format.
910c
911c scal  = real. scalar to add to the diagonal entries.
912c
913c on return:
914c----------
915c
916c a,
917c ja,
918c ia	= matrix A with diagonal elements shifted (or created).
919c
920c iw    = integer work array of length n. On return iw will
921c         contain  the positions of the diagonal entries in the
922c         output matrix. (i.e., a(iw(k)), ja(iw(k)), k=1,...n,
923c         are the values/column indices of the diagonal elements
924c         of the output matrix. ).
925c
926c Notes:
927c-------
928c     The column dimension of A is not needed.
929c     important: the matrix a may be expanded slightly to allow for
930c     additions of nonzero elements to previously nonexisting diagonals.
931c     The is no checking as to whether there is enough space appended
932c     to the arrays a and ja. if not sure allow for n additional
933c     elemnts.
934c     coded by Y. Saad. Latest version July, 19, 1990
935c-----------------------------------------------------------------------
936      logical test
937c
938      call diapos (nrow,ja,ia,iw)
939      icount = 0
940      do 1 j=1, nrow
941         if (iw(j) .eq. 0) then
942            icount = icount+1
943         else
944            a(iw(j)) = a(iw(j)) + scal
945         endif
946 1    continue
947c
948c     if no diagonal elements to insert in data structure return.
949c
950      if (icount .eq. 0) return
951c
952c shift the nonzero elements if needed, to allow for created
953c diagonal elements.
954c
955      ko = ia(nrow+1)+icount
956c
957c     copy rows backward
958c
959      do 5 ii=nrow, 1, -1
960c
961c     go through  row ii
962c
963         k1 = ia(ii)
964         k2 = ia(ii+1)-1
965         ia(ii+1) = ko
966         test = (iw(ii) .eq. 0)
967         do 4 k = k2,k1,-1
968            j = ja(k)
969            if (test .and. (j .lt. ii)) then
970               test = .false.
971               ko = ko - 1
972               a(ko) = scal
973               ja(ko) = ii
974               iw(ii) = ko
975            endif
976            ko = ko-1
977            a(ko) = a(k)
978            ja(ko) = j
979 4       continue
980c     diagonal element has not been added yet.
981         if (test) then
982            ko = ko-1
983            a(ko) = scal
984            ja(ko) = ii
985            iw(ii) = ko
986         endif
987 5    continue
988      ia(1) = ko
989      return
990c-----------------------------------------------------------------------
991c----------end-of-aplsca------------------------------------------------
992      end
993c-----------------------------------------------------------------------
994      subroutine apldia (nrow, job, a, ja, ia, diag, b, jb, ib, iw)
995      real*8 a(*), b(*), diag(nrow)
996      integer ja(*),jb(*), ia(nrow+1),ib(nrow+1), iw(*)
997c-----------------------------------------------------------------------
998c Adds a diagonal matrix to a general sparse matrix:  B = A + Diag
999c-----------------------------------------------------------------------
1000c on entry:
1001c ---------
1002c nrow	= integer. The row dimension of A
1003c
1004c job   = integer. job indicator. Job=0 means get array b only
1005c         (i.e. assume that a has already been copied into array b,
1006c         or that algorithm is used in place. ) For all practical
1007c         purposes enter job=0 for an in-place call and job=1 otherwise
1008c
1009c         Note: in case there are missing diagonal elements in A,
1010c         then the option job =0 will be ignored, since the algorithm
1011c         must modify the data structure (i.e. jb, ib) in this
1012c         situation.
1013c
1014c a,
1015c ja,
1016c ia   = Matrix A in compressed sparse row format.
1017c
1018c diag = diagonal matrix stored as a vector dig(1:n)
1019c
1020c on return:
1021c----------
1022c
1023c b,
1024c jb,
1025c ib	= resulting matrix B in compressed sparse row sparse format.
1026c
1027c
1028c iw    = integer work array of length n. On return iw will
1029c         contain  the positions of the diagonal entries in the
1030c         output matrix. (i.e., a(iw(k)), ja(iw(k)), k=1,...n,
1031c         are the values/column indices of the diagonal elements
1032c         of the output matrix. ).
1033c
1034c Notes:
1035c-------
1036c 1)        The column dimension of A is not needed.
1037c 2)        algorithm in place (b, jb, ib, can be the same as
1038c           a, ja, ia, on entry). See comments for parameter job.
1039c
1040c coded by Y. Saad. Latest version July, 19, 1990
1041c-----------------------------------------------------------------
1042      logical test
1043c
1044c     copy integer arrays into b's data structure if required
1045c
1046      if (job .ne. 0) then
1047         nnz = ia(nrow+1)-1
1048         do 2  k=1, nnz
1049            jb(k) = ja(k)
1050            b(k)  = a(k)
1051 2       continue
1052         do 3 k=1, nrow+1
1053            ib(k) = ia(k)
1054 3       continue
1055      endif
1056c
1057c     get positions of diagonal elements in data structure.
1058c
1059      call diapos (nrow,ja,ia,iw)
1060c
1061c     count number of holes in diagonal and add diag(*) elements to
1062c     valid diagonal entries.
1063c
1064      icount = 0
1065      do 1 j=1, nrow
1066         if (iw(j) .eq. 0) then
1067            icount = icount+1
1068         else
1069            b(iw(j)) = a(iw(j)) + diag(j)
1070         endif
1071 1    continue
1072c
1073c     if no diagonal elements to insert return
1074c
1075      if (icount .eq. 0) return
1076c
1077c     shift the nonzero elements if needed, to allow for created
1078c     diagonal elements.
1079c
1080      ko = ib(nrow+1)+icount
1081c
1082c     copy rows backward
1083c
1084      do 5 ii=nrow, 1, -1
1085c
1086c     go through  row ii
1087c
1088         k1 = ib(ii)
1089         k2 = ib(ii+1)-1
1090         ib(ii+1) = ko
1091         test = (iw(ii) .eq. 0)
1092         do 4 k = k2,k1,-1
1093            j = jb(k)
1094            if (test .and. (j .lt. ii)) then
1095               test = .false.
1096               ko = ko - 1
1097               b(ko) = diag(ii)
1098               jb(ko) = ii
1099               iw(ii) = ko
1100            endif
1101            ko = ko-1
1102            b(ko) = a(k)
1103            jb(ko) = j
1104 4       continue
1105c     diagonal element has not been added yet.
1106         if (test) then
1107            ko = ko-1
1108            b(ko) =  diag(ii)
1109            jb(ko) = ii
1110            iw(ii) = ko
1111         endif
1112 5    continue
1113      ib(1) = ko
1114      return
1115c-----------------------------------------------------------------------
1116c------------end-of-apldiag---------------------------------------------
1117      end
1118