*** empty log message ***

src/modules/psb_methd_mod.f90

@@ -63,7 +63,6 @@ Module psb_methd_mod
       use psb_serial_mod
       use psb_descriptor_type
       Use psb_prec_type
-!!$  parameters 
       Type(psb_dspmat_type), Intent(in)  :: a
       Type(psb_desc_type), Intent(in)    :: desc_a
       type(psb_dprec_type), intent(in)   :: prec

src/modules/psb_psblas_mod.f90

@@ -230,4 +230,21 @@ module psb_psblas_mod
   end interface
 
 
+  interface psb_gelp
+     subroutine psb_dgelp(trans,iperm,x,desc_a,info)
+       use psb_descriptor_type
+       type(psb_desc_type), intent(in)      ::  desc_a
+       real(kind(1.d0)), intent(inout)      ::  x(:,:)
+       integer, intent(inout)               ::  iperm(:),info
+       character, intent(in)                ::  trans
+     end subroutine psb_dgelp
+     subroutine psb_dgelpv(trans,iperm,x,desc_a,info)
+       use psb_descriptor_type
+       type(psb_desc_type), intent(in)      ::  desc_a
+       real(kind(1.d0)), intent(inout)      ::  x(:)
+       integer, intent(inout)               ::  iperm(:),info
+       character, intent(in)                ::  trans
+     end subroutine psb_dgelpv
+  end interface
+     
 end module psb_psblas_mod

src/prec/psb_dbldaggrmat.f90

@@ -185,7 +185,6 @@ contains
       enddo
     end if
 
-
     call psb_fixcoo(b,info)
     if(info /= 0) then
       call psb_errpush(4010,name,a_err='fixcoo')

src/prec/psb_dprecbld.f90

@@ -11,7 +11,7 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd)
   Implicit None
 
   integer, intent(out)                       :: info
-  type(psb_dspmat_type), intent(in), target  :: a
+  type(psb_dspmat_type), target              :: a
   type(psb_dprec_type),intent(inout)         :: p
   type(psb_desc_type), intent(in)            :: desc_a
   character, intent(in), optional            :: upd
@@ -125,7 +125,7 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd)
            call psb_errpush(info,name)
            goto 9999
         end if
-        call  psb_dgelp('n',a%Pl,p%baseprecv(1)%d,desc_a,info)
+        call  psb_gelp('n',a%pl,p%baseprecv(1)%d,desc_a,info)
         if(info /= 0) then
            info=4010
            ch_err='psb_dgelp'
@@ -239,6 +239,9 @@ subroutine psb_dprecbld(a,p,desc_a,info,upd)
     call psb_check_def(p%baseprecv(2)%iprcparm(jac_sweeps_),'Jacobi sweeps',&
          & 1,is_legal_jac_sweeps)
 
+    call blacs_barrier(icontxt,'All') ! to be removed
+    write(0,'(i2," Calling mlprecbld")')me
+    call blacs_barrier(icontxt,'All') ! to be removed
     call psb_mlprec_bld(a,desc_a,p%baseprecv(2),info)
     if(info /= 0) then
        info=4010
@@ -497,6 +500,8 @@ subroutine psb_mlprec_bld(a,desc_a,p,info)
      end subroutine psb_dbldaggrmat
   end interface
 
+  integer :: icontxt, nprow, npcol, me, mycol
+  
   name='psb_mlprec_bld'
   info=0
   call psb_erractionsave(err_act)

src/serial/psb_dneigh.f90

@@ -16,23 +16,9 @@ subroutine psb_dneigh(a,idx,neigh,n,info,lev)
   integer, pointer                  :: neigh(:) ! the neighbours
   integer, optional                 :: lev ! level of neighbours to find
 
-  interface psb_spgtrow
-     subroutine psb_dspgtrow(irw,a,b,info,append,iren,lrw)
-       use psb_spmat_type
-       type(psb_dspmat_type), intent(in) :: a
-       integer, intent(in)       :: irw
-       type(psb_dspmat_type), intent(inout)    :: b
-       logical, intent(in), optional :: append
-       integer, intent(in), target, optional :: iren(:)
-       integer, intent(in), optional :: lrw
-       integer, intent(out)  :: info
-     end subroutine psb_dspgtrow
-  end interface
-
-  type(psb_dspmat_type) :: atmp
   integer, pointer :: tmpn(:)
   integer :: level, dim, i, j, k, r, c, brow,&
-       & elem_pt, ii, n1, col_idx, ne, err_act
+       & elem_pt, ii, n1, col_idx, ne, err_act, nn, nidx
   character(len=20)                 :: name, ch_err
 
   name='psb_dneigh'
@@ -52,105 +38,245 @@ subroutine psb_dneigh(a,idx,neigh,n,info,lev)
     level=1
   end if
 
-  if ((a%fida /= 'CSR')) then 
+  call psb_dneigh1l(a,idx,neigh,n)
-
-     call psb_spall(atmp,1,info)
-     if(info.ne.0) then
-        call psb_errpush(4010,name)
-        goto 9999
-     end if
-     call psb_spgtrow(idx,a,atmp,info)
-     if(info.ne.0) then
-        call psb_errpush(4010,name)
-        goto 9999
-     end if
-     dim=atmp%infoa(psb_nnz_)
-     allocate(tmpn(dim))
-     tmpn(1:dim)=atmp%ia2(1:dim)
   
   if(level.eq.2) then
-        do i=1,dim
+     allocate(tmpn(1))
-           if ((1<=tmpn(i)).and.(tmpn(i)<=a%m).and.(tmpn(i).ne.idx)) then
+     n1=n
-              call psb_spgtrow(tmpn(i),a,atmp,info,append=.true.)
+     do i=1,n1
-              if(info.ne.0) then
+        nidx=neigh(i)
-                 call psb_errpush(4010,name)
+        if((nidx.ne.idx).and.(nidx.gt.0).and.(nidx.le.a%m)) then
-                 goto 9999
+           call psb_dneigh1l(a,nidx,tmpn,nn)
-              end if
+           if((n+nn).gt.size(neigh)) call psb_realloc(n+nn,neigh,info)
+           neigh(n+1:n+nn)=tmpn(1:nn)
+           n=n+nn
         end if
      end do
   end if
   
+  call psb_erractionrestore(err_act)
+  return
   
-     dim=atmp%infoa(psb_nnz_)
+9999 continue
-     if(dim.gt.size(neigh)) &
+  call psb_erractionrestore(err_act)
-          & call psb_realloc(dim,neigh,info)
+  if (err_act.eq.act_abort) then
-     if(info.ne.0) then
+     call psb_error()
-        call psb_errpush(4010,name)
+     return
-        goto 9999
-     end if
-     call psb_spfree(atmp,info)
-     if(info.ne.0) then
-        call psb_errpush(4010,name)
-        goto 9999
   end if
-     call psb_nullify_sp(atmp)
+  return
-     deallocate(tmpn)
 
-  else if(a%fida.eq.'CSR') then
      
-    dim=0
+contains
-    if(level.eq.1) dim=(a%ia2(idx+1)-a%ia2(idx))
+
-    if(dim >size(neigh)) call psb_realloc(dim,neigh,info)
+  subroutine psb_dneigh1l(a,idx,neigh,n)
-    if(info.ne.izero) then
+    
-       info=4010
+    use psb_realloc_mod
-       ch_err='psrealloc'
+    use psb_const_mod
-       call psb_errpush(info,name,a_err=ch_err)
+    use psb_spmat_type
-       goto 9999
+    implicit none
+    
+    
+    type(psb_dspmat_type), intent(in) :: a   ! the sparse matrix
+    integer, intent(in)               :: idx ! the index whose neighbours we want to find
+    integer, intent(out)              :: n   ! the number of neighbours and the info
+    integer, pointer                  :: neigh(:) ! the neighbours
+
+
+    select case(a%fida(1:3))
+    case('CSR')
+       call csr_dneigh1l(a,idx,neigh,n)
+    case('COO')
+       call coo_dneigh1l(a,idx,neigh,n)
+    case('JAD')
+       call jad_dneigh1l(a,idx,neigh,n)
+    end select
+
+  end subroutine psb_dneigh1l
+
+  subroutine csr_dneigh1l(a,idx,neigh,n)
+    
+    use psb_realloc_mod
+    use psb_const_mod
+    use psb_spmat_type
+    implicit none
+    
+    
+    type(psb_dspmat_type), intent(in) :: a   ! the sparse matrix
+    integer, intent(in)               :: idx ! the index whose neighbours we want to find
+    integer, intent(out)              :: n   ! the number of neighbours and the info
+    integer, pointer                  :: neigh(:) ! the neighbours
+    
+    integer :: dim, i, iidx
+
+    if(a%pl(1).ne.0) then
+       iidx=a%pl(idx)
+    else
+       iidx=idx
     end if
 
+    dim=a%ia2(iidx+1)-a%ia2(iidx)
+    if(dim.gt.size(neigh)) call psb_realloc(dim,neigh,info)
+
     n=0
-    if(level.eq.1) then
+    do i=a%ia2(iidx), a%ia2(iidx+1)-1
-      do i=a%ia2(idx), a%ia2(idx+1)-1
        n=n+1
        neigh(n)=a%ia1(i)
     end do
     
-    else
+  end subroutine csr_dneigh1l
 
-      do i=a%ia2(idx), a%ia2(idx+1)-1
 
-        j=a%ia1(i)
+  subroutine coo_dneigh1l(a,idx,neigh,n)
-        if ((1<=j).and.(j<=a%m).and.(j.ne.idx)) then
     
-           dim=dim+ a%ia2(j+1)-a%ia2(j)
+    use psb_realloc_mod
-           if(dim.gt.size(neigh)) call psb_realloc(dim,neigh,info)
+    use psb_const_mod
-           if(info.ne.izero) then
+    use psb_spmat_type
-              info=4010
+    implicit none
-              ch_err='psrealloc'
+    
-              call psb_errpush(info,name,a_err=ch_err)
+    
-              goto 9999
+    type(psb_dspmat_type), intent(in) :: a   ! the sparse matrix
+    integer, intent(in)               :: idx ! the index whose neighbours we want to find
+    integer, intent(out)              :: n   ! the number of neighbours and the info
+    integer, pointer                  :: neigh(:) ! the neighbours
+    
+    integer :: dim, i, iidx, ip, nza
+
+    if(a%pl(1).ne.0) then
+       iidx=a%pl(idx)
+    else
+       iidx=idx
     end if
 
-          do k=a%ia2(j), a%ia2(j+1)-1
+    nza=a%infoa(psb_nnz_)
-            n=n+1
+
-            neigh(n)=a%ia1(k)
+    if (a%infoa(psb_srtd_) == psb_isrtdcoo_) then 
+       call ibsrch(ip,iidx,nza,a%ia1)
+       if (ip /= -1) then 
+          ! bring ip backward to the beginning of the row
+          do 
+             if (ip < 2) exit
+             if (a%ia1(ip-1) == iidx) then  
+                ip = ip -1 
+             else 
+                exit
+             end if
           end do
        end if
+       
+       dim=0
+       do
+          if(a%ia1(ip).eq.iidx) then
+             dim=dim+1
+             if(dim.gt.size(neigh)) call psb_realloc(dim*3/2,neigh,info)
+             neigh(dim)=a%ia2(ip)
+             ip=ip+1
+          else
+             exit
+          end if
        end do
 
+    else
+
+       dim=0
+       do i=1,nza
+          if(a%ia1(i).eq.iidx) then
+             dim=dim+1
+             if(dim.gt.size(neigh)) call psb_realloc(dim*3/2,neigh,info)
+             neigh(dim)=a%ia2(ip)
           end if
+       end do
+    end if
+
+    n=dim
+
+  end subroutine coo_dneigh1l
+    
+
+
+  subroutine jad_dneigh1l(a,idx,neigh,n)
     
+    use psb_realloc_mod
+    use psb_const_mod
+    use psb_spmat_type
+    implicit none
+    
+    
+    type(psb_dspmat_type), intent(in) :: a   ! the sparse matrix
+    integer, intent(in)               :: idx ! the index whose neighbours we want to find
+    integer, intent(out)              :: n   ! the number of neighbours and the info
+    integer, pointer                  :: neigh(:) ! the neighbours
+    
+    integer :: dim, i, iidx, ip, nza
+    integer, pointer                      :: ia1(:), ia2(:), ia3(:),&
+         & ja(:), ka(:)
+    integer  :: png, pia, pja, ipx, blk, rb, row, k_pt, npg, col, ng
+
+    if(a%pl(1).ne.0) then
+       iidx=a%pl(idx)
+    else
+       iidx=idx
     end if
 
-  call psb_erractionrestore(err_act)
+    nza=a%infoa(psb_nnz_)
-  return
 
-9999 continue
+    png = a%ia2(1) ! points to the number of blocks
-  call psb_erractionrestore(err_act)
+    pia = a%ia2(2) ! points to the beginning of ia(3,png)
-  if (err_act.eq.act_abort) then
+    pja = a%ia2(3) ! points to the beginning of ja(:)
-     call psb_error()
+    
-     return
+    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
+
+
+    i=0
+    dim=0
+    blkfnd: do
+       i=i+1
+       if(ia1(i).eq.iidx) then
+          blk=i
+          dim=dim+ia3(i)-ia2(i)
+          ipx = ia1(i)         ! the first row index of the block
+          rb  = iidx-ipx   ! the row offset within the block
+          row = ia3(i)+rb
+          dim  = dim+ja(row+1)-ja(row)
+          exit blkfnd
+       else if(ia1(i).gt.iidx) then
+          blk=i-1
+          dim=dim+ia3(i-1)-ia2(i-1)
+          ipx = ia1(i-1)         ! the first row index of the block
+          rb  = iidx-ipx   ! the row offset within the block
+          row = ia3(i-1)+rb
+          dim  = dim+ja(row+1)-ja(row)
+          exit blkfnd
        end if
-  return
+    end do blkfnd
+    
+    if(dim.gt.size(neigh)) call psb_realloc(dim,neigh,info)
+    
+    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  = iidx-ipx       ! the row offset within the block
+    npg = ja(k_pt+1)-ja(k_pt)  ! the number of rows in the block
+    
+    k=0
+    do  col = ia2(blk), ia3(blk)-1 
+       k=k+1
+       neigh(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
+       neigh(k)  = ka(j)
+    end do
+    
+    n=k
+
+  end subroutine jad_dneigh1l
 
 end subroutine psb_dneigh
+  

src/serial/psb_dspgtdiag.f90

@@ -68,6 +68,8 @@ subroutine psb_dspgtdiag(a,d,info)
  else if (a%fida == 'JAD') then 
 
     rng=min(a%m,a%k)
+    nrb=16
+    write(0,*)'in spgtdiag'
     do i=1, rng, nrb
        irb=min(i+nrb-1,rng)
        call psb_spgtrow(i,a,tmpa,info,lrw=irb)
@@ -86,6 +88,7 @@ subroutine psb_dspgtdiag(a,d,info)
        enddo
        
     end do
+    write(0,*)'leaving spgtdiag'
 
  end if
 

src/serial/psb_dspgtrow.f90

@@ -429,21 +429,38 @@ contains
           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
 
+          if(associated(iren))then
+             do  col = ia2(blk), ia3(blk)-1 
+                k=k+1
+                b%aspk(nzb+k) = a%aspk(ja(col)+rb)
+                b%ia1(nzb+k)  = iren(irw+i-1)
+                b%ia2(nzb+k)  = iren(ka(ja(col)+rb))
+             end do
+          else
              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
-
+          end if
           ! extract second part of the row from the csr tail
           row=ia3(blk)+rb
+          if(associated(iren))then
+             do j=ja(row), ja(row+1)-1
+                k=k+1
+                b%aspk(nzb+k) = a%aspk(j)
+                b%ia1(nzb+k)  = iren(irw+i-1)
+                b%ia2(nzb+k)  = iren(ka(j))
+             end do
+          else
              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 if
        end do
 
        b%infoa(psb_nnz_) = nzb+k

test/Fileread/RUNS/rtst.inp

@@ -1,13 +1,13 @@
 11                  Number of inputs
-kivap001.mtx      This (and others) from: http://math.nist.gov/MatrixMarket/ 
+bcsstk35.mtx      This (and others) from: http://math.nist.gov/MatrixMarket/ 
 NONE
-BICGSTAB 
+CGS        
 ILU                 !!!! Actually, it's IPREC below. Should take this line out.
 CSR
 2                   IPART: Partition method
 1                   ISTOPC
 00800               ITMAX
--1                  ITRACE
+6                   ITRACE
 7                   IPREC       0:NONE  1:DIAGSC  2:ILU  3: AS 4: RAS 5,6: variants
 1                   ML
 1.d-6               EPS

test/Fileread/df_sample.f90

@@ -64,8 +64,8 @@ program df_sample
   integer            :: internal, m,ii,nnzero
   real(kind(1.d0)) :: mpi_wtime, t1, t2, tprec, r_amax, b_amax,&
        &scale,resmx,resmxp
-  integer :: nrhs, nrow, n_row, dim, nv
+  integer :: nrhs, nrow, n_row, dim, nv, ne
-  integer, pointer :: ivg(:), ipv(:)
+  integer, pointer :: ivg(:), ipv(:), neigh(:)
 
   external mpi_wtime
   
@@ -260,6 +260,9 @@ program df_sample
   else if (cmethd.eq.'CGS') then
     call  psb_cgs(a,pre,b_col,x_col,eps,desc_a,info,& 
        & itmax,iter,err,itrace)     
+  else if (cmethd.eq.'CG') then
+    call  psb_cg(a,pre,b_col,x_col,eps,desc_a,info,& 
+       & itmax,iter,err,itrace)     
   else if (cmethd.eq.'BICGSTABL') then
     call  psb_bicgstabl(a,pre,b_col,x_col,eps,desc_a,info,& 
        & itmax,iter,err,ierr,itrace,ml)     
@@ -297,12 +300,12 @@ program df_sample
     if (amroot) then
       write(0,'(" ")')
       write(0,'("Saving x on file")')
-      write(20,*) 'matrix: ',mtrx_file
+!!$      write(20,*) 'matrix: ',mtrx_file
-      write(20,*) 'computed solution on ',nprow,' processors.'
+!!$      write(20,*) 'computed solution on ',nprow,' processors.'
-      write(20,*) 'iterations to convergence: ',iter
+!!$      write(20,*) 'iterations to convergence: ',iter
-      write(20,*) 'error indicator (infinity norm) on exit:', &
+!!$      write(20,*) 'error indicator (infinity norm) on exit:', &
-           & ' ||r||/(||a||||x||+||b||) = ',err
+!!$           & ' ||r||/(||a||||x||+||b||) = ',err
-      write(20,*) 'max residual = ',resmx, resmxp
+!!$      write(20,*) 'max residual = ',resmx, resmxp
       do i=1,m_problem
         write(20,998) i,x_col_glob(i),r_col_glob(i),b_col_glob(i)
       enddo

test/Fileread/getp.f90

@@ -104,6 +104,7 @@ CONTAINS
           write(*,'("Preconditioner       : ",i)')iprec
           if(iprec.gt.2) write(*,'("Overlapping levels   : ",i)')novr
           write(*,'("Iterative method     : ",a)')cmethd
+          write(*,'("Storage format       : ",a3)')afmt(1:3)
           write(*,'(" ")')
       else
         CALL PR_USAGE(0)

test/Fileread/mmio.f90

@@ -149,7 +149,6 @@ contains
         a%ia1(1:nzr)  = ia2_loc(1:nzr)
         tmp(1:nzr)    = ia1_loc(1:nzr)
       else
-        write(0,*) 'After DESYM: ',nzr,ia2_loc(1:10)
         do i=1,nzr
           a%aspk(i) = as_loc(i)
           a%ia1(i)  = ia2_loc(i)