subroutine spextr(A_m, A_n, A_nel, A_mnel, A_icol, A_R, A_I, $ B_m, B_n, B_nel, B_mnel, B_icol, B_R, B_I, $ it, i, ni, j, nj, nel_max, ptr, p, ierr) * * PURPOSE * perform a sparse matrix extraction, ie, computes the * sparse matrix B such that : B = A(i,j) where i and * j are vectors of indices of length ni and nj respectively. * * depending upon the argument it, A may be (and so B) : * a/ a boolean sparse (it = -1). In this case the arrays * A_R, A_I, B_R, B_I are not used. * b/ a real sparse (it = 0). The arrays A_I and B_I are not used. * * c/ a complex sparse (it = 1) * * CAUTION * nel_max is the maximum non zeros elements authorized * for B => if this limit is not enough the error * indicator ierr is set to -1 (else 0) * * AUTHOR * Bruno Pincon * implicit none integer A_m, A_n, A_nel, A_mnel(*), A_icol(*) integer B_m, B_n, B_nel, B_mnel(*), B_icol(*) integer it, i(*), ni, j(*), nj, nel_max, ptr(*), p(*), ierr double precision A_R(*), A_I(*), B_R(*), B_I(*) * internal vars integer k, kk, ka, kb, ib, ia, ja, ptrb logical allrow, allcol, j_in_order * external functions and subroutines called integer dicho_search, dicho_search_bis logical is_in_order external dicho_search, dicho_search_bis, is_in_order, qsorti, $ insert_in_order, icopy, unsfdcopy, sz2ptr ierr = 0 allrow = ni.lt.0 ! ni < 0 for a A(:,j) extraction allcol = nj.lt.0 ! nj < 0 for a A(i,:) extraction if (allrow) then ni = A_m endif if (allcol) then nj = A_n else j_in_order = is_in_order(j, nj) if ( .not. j_in_order ) call qsorti(j, p, nj) endif B_m = ni B_n = nj call sz2ptr(A_mnel, A_m, ptr) kb = 1 do ib = 1, B_m B_mnel(ib) = 0 if (allrow) then ia = ib else ia = i(ib) endif if ( A_mnel(ia) .gt. 0 ) then if (allcol) then if ( kb+A_mnel(ia) .ge. nel_max ) then ! not enough memory ierr = -1 return endif B_mnel(ib) = A_mnel(ia) call icopy(A_mnel(ia), A_icol(ptr(ia)), 1, B_icol(kb), 1) if ( it .ge. 0 ) $ call unsfdcopy(A_mnel(ia), A_R(ptr(ia)), 1, B_R(kb), 1) if ( it .eq. 1 ) $ call unsfdcopy(A_mnel(ia), A_I(ptr(ia)), 1, B_I(kb), 1) kb = kb + A_mnel(ia) elseif ( nj .gt. A_mnel(ia) .and. j_in_order ) then ! algorithm with loop on the row ia of A do ka = ptr(ia), ptr(ia+1)-1 ja = A_icol(ka) k = dicho_search(ja, j, nj) if ( k .ne. 0 ) then ! we have found (the smallest) k such that ja = j(k) 100 if (kb .gt. nel_max) then ! not enough memory ierr = -1 return endif B_mnel(ib) = B_mnel(ib) + 1 B_icol(kb) = k if ( it .ge. 0 ) B_R(kb) = A_R(ka) if ( it .eq. 1 ) B_I(kb) = A_I(ka) kb = kb + 1 if ( k .lt. nj) then ! verify if j(k+1) = ja k = k + 1 if ( j(k) .eq. ja ) goto 100 endif endif enddo elseif ( nj .gt. 2*A_mnel(ia) .and. .not.j_in_order ) then ! algorithm with loop on the row ia of A but the array ! j is not in order : nevertheless we may work with j(p(.)) where ! p is the permutation which reorders the array j ptrb = kb do ka = ptr(ia), ptr(ia+1)-1 ja = A_icol(ka) k = dicho_search_bis(ja, j, p, nj) if ( k .ne. 0 ) then ! we have found (the smallest) k such that ja = j(p(k)) 200 if (kb .gt. nel_max) then ! not enough memory ierr = -1 return endif B_mnel(ib) = B_mnel(ib) + 1 call insert_in_order( B_icol, ptrb, kb, p(k), it, $ B_R, B_I, A_R(ka), A_I(ka) ) kb = kb + 1 if ( k .lt. nj) then ! verify if j(p(k+1)) = ja k = k + 1 if ( j(p(k)) .eq. ja ) goto 200 endif endif enddo else ! algorithm with loop on j(1),...,j(nj) do k = 1, nj kk = dicho_search(j(k), A_icol(ptr(ia)), A_mnel(ia)) if ( kk .ne. 0 ) then ! we have found kk such that j(k) = A_icol(ptr(ia)-1+kk) ka = kk + ptr(ia) - 1 if (kb .gt. nel_max) then ! not enough memory ierr = -1 return endif B_mnel(ib) = B_mnel(ib) + 1 B_icol(kb) = k if ( it .ge. 0 ) B_R(kb) = A_R(ka) if ( it .eq. 1 ) B_I(kb) = A_I(ka) kb = kb + 1 endif enddo endif endif enddo B_nel = kb - 1 end