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 lspisp(ma,na,nela,inda,i,ni,j,nj,
     $     mb,nb,nelb,indb,mr,nr,nelr,indr,ptrb,ierr)
c     extract a submatrix from a sparse boolean matrix
c!
      integer inda(*),indr(*),i(*),j(*),ptrb(*)
      integer indb(*),mb,nelb
      integer ma,na,ni,nj,mr,nr,nela,nelr,ierr
      logical allrow,allcol
      integer findl
      external findl
c
      mr=ni
      nr=nj
      ierr=0
      nelmx=nelr
      allrow=ni.lt.0
      allcol=nj.lt.0
      if(allrow) then 
         mr=ma
         ni=mr
      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
         nj=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
         call icopy(mb+nelb,indb,1,indr,1)
         nelr=nelb
         return
      elseif(allcol) then
c     a(i,:)=b
         jr=1
         jb=1
         ja=1
         call sz2ptr(indb,mb,ptrb)
         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)
                  if(jr+indr(l).gt.nelmx) then
                     ierr=1
                     return
                  endif
                  jr=jr+indr(l)
                  ja=ja+indr(l)
               else
                  indr(l)=0
               endif
            else
c     all this line is replaced by corresponding b line
               jb=ptrb(ii)
               indr(l)=indb(ii)
               if(jr+indr(l).gt.nelmx) then
                  ierr=1
                  return
               endif
               call icopy(indr(l),indb(mb+jb),1,indr(mr+jr),1)
               jr=jr+indr(l)
               if(l.le.ma) then
                  ja=ja+inda(l)
               endif
            endif
 20      continue
         nelr=jr-1
      elseif(allrow) then
c     a(:,j)=b
         jr=1
         ja=0
         jb=0
         do 35 l=1,ma
            ja1=1
            nbl=indb(l)
            nal=inda(l)
            indr(l)=0
            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
                     if(jr+1.gt.nelmx) then
                        ierr=1
                        return
                     endif
                     indr(l)=indr(l)+1
                     indr(mr+jr)=k
                     jr=jr+1
                     ja1=ja1+1
                  endif
               else
c     the  a(l,k) element is replaced by b(l,jj) element if non zero
                  jb1=findl(jj,indb(mb+jb+1),nbl)
                  if(jb1.ne.0) then
c     *           b(l,jj) element if non zero
                     if(jr+1.gt.nelmx) then
                        ierr=1
                        return
                     endif
                     indr(l)=indr(l)+1
                     indr(mr+jr)=k
                     jr=jr+1
                  endif
                  if(ja1.le.nal.and.inda(ma+ja+ja1).eq.k) ja1=ja1+1
               endif
 31         continue
            ja=ja+nal
            jb=jb+nbl
 35     continue
        nelr=jr-1
        return
      else
c     a(i,j)=b
         jr=1
         ja=0
         call sz2ptr(indb,mb,ptrb)
         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
               jb=ptrb(ii)-1
               jb1=1
               nbl=indb(ii)
               do 42 k=1,nr
                  jj=findl(k,j,nj)
                  if(jj.eq.0) then
c     *           insert a(l,k) element insert  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
                        if(jr+1.gt.nelmx) then
                           ierr=1
                           return
                        endif
                        indr(l)=indr(l)+1
                        indr(mr+jr)=k
                        jr=jr+1
                        ja1=ja1+1
                     endif
                  else
c     *           replace a(l,k) element  by b(ii,jj) element if non zero
                     jb1=findl(jj,indb(mb+jb+1),nbl)
                     if(jb1.ne.0) then
c     *               b(l,jj) element if non zero
                        if(jr+1.gt.nelmx) then
                           ierr=1
                           return
                        endif
                        indr(l)=indr(l)+1
                        indr(mr+jr)=k
                        jr=jr+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