c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
c Copyright (C) INRIA
c 
c This file must be used under the terms of the CeCILL.
c This source file is licensed as described in the file COPYING, which
c you should have received as part of this distribution.  The terms
c are also available at    
c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
      subroutine lspis(ma,na,nela,inda,i,ni,j,nj,
     $     mb,nb,b,mr,nr,nelr,indr,ierr)
c     extract a submatrix from a sparse matrix
c!
      integer inda(*),indr(*),i(*),j(*)
      integer mb,nb
      integer ma,na,ni,nj,mr,nr,nela,nelr,ierr
      integer b(mb,*)
      logical allrow,allcol
      integer findl
      external findl
c
      nelmx=nelr
      ierr=0
      mr=ni
      nr=nj
      allrow=ni.lt.0
      allcol=nj.lt.0
      if(allrow) then 
         mr=ma
      else
         mi=0
         do 01 kk=1,ni
            mi=max(mi,i(kk))
 01      continue
         mr=max(ma,mi)
      endif
      if(allcol) then 
         nr=na
      else
         mj=0
         do 02 kk=1,nj
            mj=max(mj,j(kk))
 02      continue
         nr=max(na,mj)
      endif
      if (allrow.and.allcol) then
c     a(:,:)=b
         if(mb*nb.gt.nelmx) then
            ierr=1
            return
         endif
         call iset(mb,0,indr,1)
         jr=1
         if(mb*nb.ne.0) then
            do 04 l=1,mb
               do 03 k=1,nb
                  if(b(l,k).ne.0) then
                     indr(mb+jr)=k
                     indr(l)=indr(l)+1
                     jr=jr+1
                  endif
 03            continue
 04         continue
         endif
         nelr=jr-1
         return
      elseif(allcol) then
c     a(i,:)=b
         jr=1
         ja=1
         do 20 l=1,mr
            indr(l)=0
            ii=findl(l,i,ni)
            if(ii.eq.0) then
c     this line is not modified
               if(l.le.ma) then
                  indr(l)=inda(l)
                  call icopy(indr(l),inda(ma+ja),1,indr(mr+jr),1)
                  jr=jr+indr(l)
                  ja=ja+indr(l)
               else
                  indr(l)=0
               endif
            else
c     all this line is replaced by corresponding b line
               if(nb.gt.0) then
                  if(jr+nb.gt.nelmx) then
                     ierr=1
                     return
                  endif
                  do 19 k1=1,nb
                     if(b(ii,k1).ne.0) then
                        indr(mr+jr)=k1
                        indr(l)=indr(l)+1
                        jr=jr+1
                     endif
 19               continue
               endif
               if(l.le.ma) then
                  ja=ja+inda(l)
               endif
            endif
 20      continue
         nelr=jr-1
         if(nb.eq.0) mr=mr-ni
      elseif(allrow) then
c     a(:,j)=b
         jr=1
         ja=0
         do 35 l=1,ma
            ja1=1
            nal=inda(l)
            indr(l)=0
            kr=1
            do 31 k=1,nr
               jj=findl(k,j,nj)
               if(jj.eq.0) then
c     the  a(l,k) element is not modified insert it in r if non zero
                  if(ja1.le.nal.and.inda(ma+ja+ja1).eq.k) then
c     *              a(l,k) element is non zero
                     indr(l)=indr(l)+1
                     indr(mr+jr)=kr
                     jr=jr+1
                     ja1=ja1+1
                  endif
                  kr=kr+1
               else
c     the  a(l,k) element is replaced by b(l,jj) element 
                  if(nb.gt.0) then
                     if(b(l,jj).ne.0) then
                        if(jr+1.gt.nelmx) then
                           ierr=1
                           return
                        endif
                        indr(l)=indr(l)+1
                        indr(mr+jr)=kr
                        jr=jr+1
                     endif
                     kr=kr+1
                  endif
               endif
               if(ja1.le.nal.and.inda(ma+ja+ja1).eq.k) ja1=ja1+1
 31         continue
            ja=ja+nal
 35     continue
        nelr=jr-1
        if(nb.eq.0) nr=nr-nj
        return
      else
c     a(i,j)=b
         jr=1
         ja=0
         do 45 l=1,mr
            ja1=1
            if(l.le.ma) then 
               nal=inda(l)
            else
               nal=0
            endif
            indr(l)=0
            ii=findl(l,i,ni)
            if(ii.eq.0) then
c     *     the a(l,:) is not modified
               if(l.le.ma) then
                  indr(l)=inda(l)
                  call icopy(indr(l),inda(ma+ja+ja1),1,indr(mr+jr),1)
                  jr=jr+indr(l)
                  ja1=ja1+indr(l)
               else
                  indr(l)=0
               endif
            else
               kr=1
               do 42 k=1,nr
                  jj=findl(k,j,nj)
                  if(jj.eq.0) then
c     *           insert a(l,k) element  in r if non zero
                     if(ja1.le.nal.and.inda(ma+ja+ja1).eq.k) then
c     *              a(l,k) element is non zero
                        indr(l)=indr(l)+1
                        indr(mr+jr)=kr
                        jr=jr+1
                        ja1=ja1+1
                     endif
                     kr=kr+1
                  else
c     *           replace a(l,k) element  by b(ii,jj) element if non zero
                     if(nb.gt.0) then
                        if(b(ii,jj).ne.0) then
                           if(jr+1.gt.nelmx) then
                              ierr=1
                              return
                           endif
                           indr(l)=indr(l)+1
                           indr(mr+jr)=kr
                           jr=jr+1
                        endif
                        kr=kr+1
                     endif
                     if(ja1.le.nal.and.inda(ma+ja+ja1).eq.k) ja1=ja1+1
                  endif
 42            continue
            endif
            ja=ja+nal
 45      continue
         nelr=jr-1
      endif
      end