Added code for JAD storage format

src/serial/aux/ibsrch.f

@@ -23,3 +23,4 @@
 
       return
       end
+

src/serial/dp/dcrjd.f

@@ -123,7 +123,7 @@ C
  10       CONTINUE
 
           IP2(1) = 0
-c$$$          write(0,*) 'Calling DGBLOCK first'
+
           CALL DGBLOCK(M,IA2,IP1,IAN2(PIA),IAN2(PNG), AUX, LAUX*2)
           
           PJA = PIA + 3*(IAN2(PNG)+1)
@@ -138,7 +138,6 @@ C
           ENDIF
 
           LJA = LIAN2-PJA
-c$$$          write(0,*) 'Into DGINDEX: ',lja,pja,lian2
           CALL DGINDEX(M,IAN2(PNG),AR,IA1,IA2,ARN,IAN1,IAN2(PIA), 
      +         IAN2(PJA), INFON, LARN,LIAN1,
      +         LJA,IP1, AUX, LAUX*2, SIZE_REQ,IERROR)

src/serial/jad/djadmv.f

@@ -96,8 +96,6 @@ c
         K   = IA(2,IPG)
         NPG = JA(K+1)-JA(K)
 
-c$$$         WRITE(0,*)NPG
-
         IF (NPG.EQ.4) THEN
           IPX = IA(1,IPG)                                                 
           Y0  = ZERO

src/serial/psb_dspgtrow.f90

@@ -71,10 +71,8 @@ subroutine psb_dspgtrow(irw,a,b,info,append,iren,lrw)
      call coo_dspgtrow(irw_,a,b,append_,iren_,lrw_)
      
   else if (a%fida == 'JAD') then 
-     info=135
+     call jad_dspgtrow(irw_,a,b,append_,iren_,lrw_)
-     ch_err=a%fida(1:3)
+
-     call psb_errpush(info,name,a_err=ch_err)
-     goto 9999
   else
      info=136
      ch_err=a%fida(1:3)
@@ -101,7 +99,7 @@ contains
     use psb_spmat_type
     use psb_const_mod
     implicit none
-    
+
     type(psb_dspmat_type), intent(in)     :: a
     integer                               :: irw
     type(psb_dspmat_type), intent(inout)  :: b
@@ -109,62 +107,102 @@ contains
     integer, pointer                      :: iren(:)
     integer                               :: lrw
 
-    integer :: idx,i,j ,nr,nz,nzb
+    integer :: idx,i,j ,nr,nz,nzb, row_idx
+    integer, pointer :: indices(:)
 
-    if (a%pl(1) /= 0) then
-      write(0,*) 'Fatal error in SPGTROW: do not feed a permuted mat so far!',&
-           & a%pl(1)
-      idx = -1 
-    else
-      idx = irw
-    endif
-!!$    write(0,*) 'csr_gtrow: ',irw,lrw,a%pl(1),idx
-    if (idx<0) then 
-      write(0,*) ' spgtrow Error : idx no good ',idx
-      return
-    end if
-    nr = lrw - irw + 1 
-    nz = a%ia2(idx+nr) - a%ia2(idx)
     if (append) then 
-      nzb = b%infoa(psb_nnz_)
+       nzb = b%infoa(psb_nnz_)
     else
-      nzb = 0
+       nzb = 0
     endif
-    if (min(size(b%ia1),size(b%ia2),size(b%aspk)) < nzb+nz) then 
+
-      call psb_spreall(b,nzb+nz,iret)
+    if (a%pl(1) /= 0) then
-    endif
+
-    b%fida='COO'
+       nr = lrw - irw + 1 
-!!$    write(0,*) 'csr_gtrow: ',out_,b%fida,nzb
+       allocate(indices(nr))
-    if (associated(iren)) then 
+       do i=1,nr
-      k=0
+          indices(i)=a%pl(irw+i-1)
-      do i=irw,lrw
+          nz=nz+a%ia2(indices(i)+1)-a%ia2(indices(i))
-        do j=a%ia2(i),a%ia2(i+1)-1
+       end do
-          k             = k + 1
+
-          b%aspk(nzb+k) = a%aspk(j)
+       if (min(size(b%ia1),size(b%ia2),size(b%aspk)) < nzb+nz) then 
-          b%ia1(nzb+k)  = iren(i)
+          call psb_spreall(b,nzb+nz,iret)
-          b%ia2(nzb+k)  = iren(a%ia1(j))
+       endif
-        end do
+
-      enddo
+       k=0
+       if(associated(iren)) then
+          do i=1,nr
+             row_idx=indices(i)
+             do j=a%ia2(row_idx),a%ia2(row_idx+1)-1
+                k             = k + 1
+                b%aspk(nzb+k) = a%aspk(j)
+                b%ia1(nzb+k)  = iren(row_idx)
+                b%ia2(nzb+k)  = iren(a%ia1(j))
+             end do
+          end do
+       else
+          do i=1,nr
+             row_idx=indices(i)
+             do j=a%ia2(row_idx),a%ia2(row_idx+1)-1
+                k             = k + 1
+                b%aspk(nzb+k) = a%aspk(j)
+                b%ia1(nzb+k)  = row_idx
+                b%ia2(nzb+k)  = a%ia1(j)
+             end do
+          end do
+       end if
+
+       b%infoa(psb_nnz_) = nzb+k
+       b%m = b%m+nr
+       b%k = max(b%k,a%k)
+
     else
-      k=0
+       idx = irw
-!!$      write(0,*) 'csr_gtrow: ilp',irw,lrw
+
-      do i=irw,lrw
+       if (idx<0) then 
-!!$        write(0,*) 'csr_gtrow: jlp',a%ia2(i),a%ia2(i+1)-1
+          write(0,*) ' spgtrow Error : idx no good ',idx
-        do j=a%ia2(i),a%ia2(i+1)-1
+          return
-          k             = k + 1
+       end if
-          b%aspk(nzb+k) = a%aspk(j)
+       nr = lrw - irw + 1 
-          b%ia1(nzb+k)  = i
+       nz = a%ia2(idx+nr) - a%ia2(idx)
-          b%ia2(nzb+k)  = a%ia1(j)
+
+       if (min(size(b%ia1),size(b%ia2),size(b%aspk)) < nzb+nz) then 
+          call psb_spreall(b,nzb+nz,iret)
+       endif
+       b%fida='COO'
+
+       if (associated(iren)) then 
+          k=0
+          do i=irw,lrw
+             do j=a%ia2(i),a%ia2(i+1)-1
+                k             = k + 1
+                b%aspk(nzb+k) = a%aspk(j)
+                b%ia1(nzb+k)  = iren(i)
+                b%ia2(nzb+k)  = iren(a%ia1(j))
+             end do
+          enddo
+       else
+          k=0
+
+          do i=irw,lrw
+
+             do j=a%ia2(i),a%ia2(i+1)-1
+                k             = k + 1
+                b%aspk(nzb+k) = a%aspk(j)
+                b%ia1(nzb+k)  = i
+                b%ia2(nzb+k)  = a%ia1(j)
 !!$          write(0,*) 'csr_gtrow: in:',a%aspk(j),i,a%ia1(j)
-        end do
+             end do
-      enddo
+          enddo
-    end if
+       end if
-    b%infoa(psb_nnz_) = nzb+nz
+       b%infoa(psb_nnz_) = nzb+nz
-    if (a%pr(1) /= 0) then
+       if (a%pr(1) /= 0) then
-      write(0,*) 'Feeling lazy today, Right Permutation will have to wait'
+          write(0,*) 'Feeling lazy today, Right Permutation will have to wait'
+       endif
+       b%m = b%m+lrw-irw+1
+       b%k = max(b%k,a%k)
+
     endif
-    b%m = b%m+lrw-irw+1
-    b%k = max(b%k,a%k)
 
   end subroutine csr_dspgtrow
 
@@ -309,5 +347,124 @@ contains
     b%k = max(b%k,a%k)
   end subroutine coo_dspgtrow
 
+
+
+
+  subroutine jad_dspgtrow(irw,a,b,append,iren,lrw)
+
+    type(psb_dspmat_type), intent(in)     :: a
+    integer                               :: irw
+    type(psb_dspmat_type), intent(inout)  :: b
+    logical, intent(in)                   :: append
+    integer, pointer                      :: iren(:)
+    integer                               :: lrw
+
+    integer, pointer                      :: ia1(:), ia2(:), ia3(:),&
+         & ja(:), ka(:), indices(:), blks(:)
+    integer  :: png, pia, pja, ipx, blk, rb, row, k_pt, npg, col, ng
+
+
+    png = a%ia2(1) ! points to the number of blocks
+    pia = a%ia2(2) ! points to the beginning of ia(3,png)
+    pja = a%ia2(3) ! points to the beginning of ja(:)
+    
+    ng  =  a%ia2(png)              ! the number of blocks
+    ja  => a%ia2(pja:)             ! the array containing the pointers to ka and aspk
+    ka  => a%ia1(:)                ! the array containing the column indices
+    ia1 => a%ia2(pia:pja-1:3)      ! the array containing the first row index of each block
+    ia2 => a%ia2(pia+1:pja-1:3)    ! the array containing a pointer to the pos. in ja to the first jad column
+    ia3 => a%ia2(pia+2:pja-1:3)    ! the array containing a pointer to the pos. in ja to the first csr column
+
+
+    if (append) then 
+       nzb = b%infoa(psb_nnz_)
+    else
+       nzb = 0
+    endif
+
+    if (a%pl(1) /= 0) then
+
+       nr = lrw - irw + 1 
+       allocate(indices(nr),blks(nr))
+       nz = 0
+
+       do i=1,nr
+          indices(i)=a%pl(irw+i-1)
+          j=0
+          blkfnd: do
+             j=j+1
+             if(ia1(j).eq.indices(i)) then
+                blks(i)=j
+                nz=nz+ia3(j)-ia2(j)
+                ipx = ia1(j)         ! the first row index of the block
+                rb  = indices(i)-ipx   ! the row offset within the block
+                row = ia3(j)+rb
+                nz  = nz+ja(row+1)-ja(row)
+                exit blkfnd
+             else if(ia1(j).gt.indices(i)) then
+                blks(i)=j-1
+                nz=nz+ia3(j-1)-ia2(j-1)
+                ipx = ia1(j-1)         ! the first row index of the block
+                rb  = indices(i)-ipx   ! the row offset within the block
+                row = ia3(j-1)+rb
+                nz  = nz+ja(row+1)-ja(row)
+                exit blkfnd
+             end if
+          end do blkfnd
+       end do
+       
+       if (size(b%ia1) < nzb+nz) then 
+          call psb_spreall(b,nzb+nz,iret)
+       endif
+
+       k=0
+       ! cycle over rows
+       do i=1,nr
+
+          ! find which block the row belongs to
+          blk = blks(i)
+
+          ! extract first part of the row from the jad block
+          ipx = ia1(blk)             ! the first row index of the block
+          k_pt= ia2(blk)             ! the pointer to the beginning of a column in ja
+          rb  = indices(i)-ipx       ! the row offset within the block
+          npg = ja(k_pt+1)-ja(k_pt)  ! the number of rows in the block
+
+          do  col = ia2(blk), ia3(blk)-1 
+             k=k+1
+             b%aspk(nzb+k) = a%aspk(ja(col)+rb)
+             b%ia1(nzb+k)  = irw+i-1
+             b%ia2(nzb+k)  = ka(ja(col)+rb)
+          end do
+
+          ! extract second part of the row from the csr tail
+          row=ia3(blk)+rb
+          do j=ja(row), ja(row+1)-1
+             k=k+1
+             b%aspk(nzb+k) = a%aspk(j)
+             b%ia1(nzb+k)  = irw+i-1
+             b%ia2(nzb+k)  = ka(j)
+          end do
+       end do
+
+       b%infoa(psb_nnz_) = nzb+k
+       b%m = b%m+lrw-irw+1
+       b%k = max(b%k,a%k)
+       b%fida='COO'
+       
+    else
+       ! There might be some problems
+       info=134
+       ch_err=a%fida(1:3)
+       call psb_errpush(info,name,a_err=ch_err)
+    end if
+
+    
+
+  end subroutine jad_dspgtrow
+
+
+
+
 end subroutine psb_dspgtrow
 

src/serial/psb_dspinfo.f90

@@ -20,7 +20,8 @@ subroutine psb_dspinfo(ireq,a,ires,info,iaux)
   integer, intent(out)              :: ires, info
   integer, intent(in), optional     :: iaux
 
-  integer :: i,j,k,ip,jp,nr,irw,nz, err_act
+  integer :: i,j,k,ip,jp,nr,irw,nz, err_act, row, ipx, pia, pja, rb,idx
+  integer, pointer :: ia1(:), ia2(:), ia3(:), ja(:)
   character(len=20)                 :: name, ch_err
 
   name='psb_dspinfo'
@@ -29,17 +30,14 @@ subroutine psb_dspinfo(ireq,a,ires,info,iaux)
 
 
   if (ireq == psb_nztotreq_) then 
+     ! The number of nonzeroes
      if (a%fida == 'CSR') then 
         nr   = a%m
         ires = a%ia2(nr+1)-1
      else if ((a%fida == 'COO').or.(a%fida == 'COI')) then 
         ires = a%infoa(psb_nnz_)
      else if (a%fida == 'JAD') then 
-        ires=-1
+        ires = a%infoa(psb_nnz_)
-        info=135
-        ch_err=a%fida(1:3)
-        call psb_errpush(info,name,a_err=ch_err)
-        goto 9999
      else
         ires=-1
         info=136
@@ -49,6 +47,7 @@ subroutine psb_dspinfo(ireq,a,ires,info,iaux)
      end if
 
   else if (ireq == psb_nzrowreq_) then 
+     ! The number of nonzeroes in row iaux
      if (.not.present(iaux)) then 
         write(0,*) 'Need IAUX when ireq=nzrowreq'
         ires=-1
@@ -92,11 +91,35 @@ subroutine psb_dspinfo(ireq,a,ires,info,iaux)
 !!$        if (a%ia1(i) == irw) ires = ires + 1
 !!$      enddo
      else if (a%fida == 'JAD') then 
-        ires=-1
+        pia = a%ia2(2) ! points to the beginning of ia(3,png)
-        info=135
+        pja = a%ia2(3) ! points to the beginning of ja(:)
-        ch_err=a%fida(1:3)
+        ja  => a%ia2(pja:)             ! the array containing the pointers to ka and aspk
-        call psb_errpush(info,name,a_err=ch_err)
+        ia1 => a%ia2(pia:pja-1:3)      ! the array containing the first row index of each block
-        goto 9999
+        ia2 => a%ia2(pia+1:pja-1:3)    ! the array containing a pointer to the pos. in ja to the first jad column
+        ia3 => a%ia2(pia+2:pja-1:3)    ! the array containing a pointer to the pos. in ja to the first csr column
+
+        idx=a%pl(irw)
+        j=0
+        nz=0
+        blkfnd: do
+           j=j+1
+           if(ia1(j).eq.idx) then
+              nz=nz+ia3(j)-ia2(j)
+              ipx = ia1(j)         ! the first row index of the block
+              rb  = idx-ipx        ! the row offset within the block
+              row = ia3(j)+rb
+              nz  = nz+ja(row+1)-ja(row)
+              exit blkfnd
+           else if(ia1(j).gt.idx) then
+              nz=nz+ia3(j-1)-ia2(j-1)
+              ipx = ia1(j-1)         ! the first row index of the block
+              rb  = idx-ipx          ! the row offset within the block
+              row = ia3(j-1)+rb
+              nz  = nz+ja(row+1)-ja(row)
+              exit blkfnd
+           end if
+        end do blkfnd
+        ires=nz
      else
         ires=-1
         info=136