      program DiagUpdt
c
c     Example code. Given two CSC matrices A and B in canonical form, a scalar
c     w, compute
c                 C = diag(A) + w B
c     in the CSC representation.
c
c     It is assumed that A, B and C have no missing diagonal entries, i.e.,
c     missing entries are represented with zero values.
c
      implicit none
      integer, parameter :: NA=8,NNZA=18, NB=NA, NNZB=20,
     &                      ND=NA, NC=NA, NNZC=NNZB
      integer ia(NA+1),ja(NNZA), ib(NB+1),jb(NNZB),
     &        ic(NA+1),jc(NNZB)
      REAL :: w, va(NNZA), vb(NNZB), vd(NA), vc(NNZB)
      integer i,j
      DATA w/2.0/
      DATA ia /1,5,8,10,12,15,17,18,19/
      DATA ja
     1 /1,  3,    6,7,
     2    2,3,  5,
     3      3,        8,
     4        4,    7,
     5          5,6,7,
     6            6,  8,
     7              7,
     8                8/
      DATA va
     1 /7.0,      1.0,           2.0, 7.0,
     2       -4.0,8.0,      2.0,
     3            1.0,                     5.0,
     4                 7.0,      9.0,
     5                      5.0, 1.0, 5.0,
     6                          -1.0,      5.0,
     7                               11.0,
     8                                     5.0/
      DATA ib /1,5,8,10,12,15,17,19,21/
      DATA jb
     1 /1,  3,    6,7,
     2    2,3,  5,
     3      3,        8,
     4        4,    7,
     5          5,6,7,
     6            6,  8,
     7        4,    7,
     8      3,        8/
      DATA vb
     1 /7.0,      1.0,           2.0, 7.0,
     2       -4.0,8.0,      2.0,
     3            1.0,                     5.0,
     4                 7.0,      9.0,
     5                      5.0, 1.0, 5.0,
     6                          -1.0,      5.0,
     7                 2.0,          11.0,
     8            5.0,                     5.0/
c
c     collect diagonal(A)
c
      do i=1,NA
         do j=ia(i), ia(i+1)-1
            if(ja(j).eq.i)then   !diagonal entry located
               vd(i) = va(j)
            end if
         end do
      end do
c
c     calculate C = diag(A) + w B
c
      do i=1,NB
         ic(i) = ib(i)
         do j=ib(i), ib(i+1)-1
            jc(j) = jb(j)
            if(jb(j).eq.i)then
               vc(j) = vd(i) + w*vb(j)
            else
               vc(j) = w*vb(j)
            endif
         end do
      end do
      ic(NB+1) = NNZC + 1
c
c     print result matrix C
c
      do i=1,NC
         do j = ic(i), ic(i+1)-1
            write(*,'(1x,2i3,2x,F8.2)')i,jc(j),vc(j)
         end do
      end do
      end


