* this file contains some utility routines for some sparse
* operations (extraction / insertion)
*
* I will comment a little more in the future
* (bruno december 2004)
*

      subroutine  set_perm_id(p, n)
      implicit none
      integer p(*), n
      integer i
      do i = 1, n
         p(i) = i
      enddo
      end ! subroutine set_perm_id


      subroutine insert_row(ka, A_it, A_mnel, A_icol, A_R, A_I, 
     $                      kc, C_it, C_mnel, C_icol, C_R, C_I,
     $                      j, p, nj, ib, B_m, B_it, B_R, B_I,
     $                      B_is_scalar, nelmax, ierr)

*      ka, kc : pointent sur les debuts de ligne de A et C

      implicit none
      integer ka, A_it, A_mnel, A_icol(*), kc, C_it, C_mnel, C_icol(*),
     $        j(*), p(*), nj, it, ib, nelmax, ierr, B_m, B_it
      logical B_is_scalar
      double precision A_R(*), A_i(*), C_R(*), C_I(*),
     $                 B_R(B_m,*), B_I(B_m,*)

*     local vars
      integer i, j1, j2, jp, k, kamax
      double precision Bij_R, Bij_I
*     external functions and subroutines
      external insert_j1j2

      if (B_is_scalar) then
         Bij_R = B_R(1,1)
         if (B_it .eq. 1) Bij_I = B_I(1,1)
      endif

*     indice max pour la ligne de A
      kamax = ka + A_mnel - 1
      j1  = 1

      k = 1
*     loop on k from 1 to nj but some values must be skiped
 100     jp = j(p(k))
         if ( k .lt. nj ) then
            if ( j(p(k+1)) .eq. jp ) then
               k = k+1
               goto 100
            endif
         endif

         j2 = jp - 1
         call insert_j1j2(j1, j2, A_it, A_icol, A_R, A_I, ka, 
     $                    kamax, C_it, C_mnel, C_icol, C_R, C_I, kc, 
     $                    nelmax, ierr)
         if (ierr .ne. 0) return
                  
       ! insertion de B ...
         if (kc .gt. nelmax) then ! test if there is enough memory
            ierr = -1
            return
         endif

         if ( .not. B_is_scalar ) then
            Bij_R = B_R(ib,p(k))
            if (B_it .eq. 1) Bij_I = B_I(ib,p(k))
         endif

         if (C_it .eq. 0) then        ! real case
            if ( Bij_R .ne. 0 ) then
               C_icol(kc) = jp
               C_mnel = C_mnel + 1
               C_R(kc) = Bij_R
               kc = kc + 1
             endif

         elseif (B_it .eq. 0) then        ! complex case but B is real
            if ( Bij_R .ne. 0.d0 ) then
               C_R(kc) = Bij_R
               C_I(kc) = 0.d0
               C_icol(kc) = jp
               C_mnel = C_mnel + 1
               kc = kc + 1
            endif

         else                             ! complex case , B complex
            if ( Bij_R.ne.0.d0 .or. Bij_I.ne.0.d0 ) then
               C_R(kc) =  Bij_R
               C_I(kc) =  Bij_I
               C_icol(kc) = jp
               C_mnel = C_mnel + 1
               kc = kc + 1
            endif
         endif

         j1 = jp + 1
         k = k + 1
         if (k .le. nj) goto 100
*      endloop on k
 
      j2 = A_icol(kamax)
      call insert_j1j2(j1, j2, A_it, A_icol, A_R, A_I, ka, 
     $                 kamax, C_it, C_mnel, C_icol, C_R, C_I, kc, 
     $                 nelmax, ierr)

      end ! subroutine insert_row
         

      subroutine copy_sprow(i1,i2, ka, A_it, A_mnel, A_icol, A_R, A_I,
     $                             kc, C_it, C_mnel, C_icol, C_R, C_I,
     $                             nelmax, ierr)
*
*     recopie les lignes i1  i2 de A sur C
*         
      implicit none
      integer i1, i2, ka, A_it, A_mnel(*), A_icol(*),
     $                kc, C_it, C_mnel(*), C_icol(*),
     $        it, nelmax, ierr
      double precision A_R(*), A_I(*), C_R(*), C_I(*)

*     local vars
      integer nbe, i, n
*     external functions and subroutines
      external icopy, unsfdcopy

      if (i1 .gt. i2) return  ! a mettre au niveau de l'appelant ?

      nbe = 0
      do i = i1, i2
         nbe = nbe + A_mnel(i)
      enddo
      if (kc + nbe .gt. nelmax) then
         ierr = -1
         return
      endif

      n = i2 - i1 + 1
      call icopy(n,   A_mnel(i1), 1, C_mnel(i1), 1)
      call icopy(nbe, A_icol(ka), 1, C_icol(kc), 1)
      if (C_it .ge. 0) then
         call unsfdcopy(nbe, A_R(ka), 1, C_R(kc), 1)
         if (C_it .eq. 1) then     ! C is complex
            if (A_it .eq. 1) then   ! A is complex
               call unsfdcopy(nbe, A_I(ka), 1, C_I(kc), 1)
            else                   ! A is real
               call dset(nbe, 0.d0, C_I(kc), 1)
            endif
         endif
      endif

      ! mise  jour de ka et kc
      ka = ka + nbe
      kc = kc + nbe

      end ! subroutine copy_sprow


      subroutine insert_j1j2(j1, j2, A_it, A_icol, A_R, A_I, ka, kamax, 
     $                       C_it, C_mnel, C_icol, C_R, C_I, kc, 
     $                       nelmax, ierr)
*
*     insere les coefs de la ligne A dans C dont les indices de
*     colonnes sont compris entre j1 et j2
      
      implicit none
      integer j1, j2, A_it, A_icol(*), ka, kamax, C_it, C_mnel, 
     $        C_icol(*), kc, nelmax, ierr
      double precision A_R(*), A_I(*), C_R(*), C_I(*)

*     local var
      integer j

      if ( ka .gt. kamax ) return

*     goto the first element with column index >= j1
      do while ( A_icol(ka) .lt. j1 ) 
         ka = ka + 1
         if (ka .gt. kamax) return
      enddo

      j = A_icol(ka)
      do while (j .le. j2)   ! add a new element to C 
         if (kc .gt. nelmax) then  ! but test if there is enough memory before
            ierr = -1
            return
         endif
         C_icol(kc) = j
         C_R(kc) = A_R(ka)
         if ( C_it .eq. 1 ) then ! C is complex
            if ( A_it .eq. 0 ) then ! A is real
               C_I(kc) = 0.d0
            else                ! A is complex
               C_I(kc) = A_I(ka)
            endif
         endif
         kc = kc + 1          ! number of the next element of C
         C_mnel = C_mnel + 1  ! increment the number of nnz of the row for C
         ka = ka + 1
         if ( ka .gt. kamax ) return
         j = A_icol(ka)
      enddo
            
      end ! subroutine insert_j1j2


      subroutine copy_fullrow2sprow(i, kc, C_it, C_mnel, C_icol, C_R,
     $                              C_I, B_m, B_n, B_it, B_R, B_I,
     $                              B_is_scalar, nelmax, ierr)
      implicit none
      integer i, kc, C_it, C_mnel, C_icol(*), B_m, B_n, B_it, nelmax, 
     $        ierr
      logical B_is_scalar
      double precision C_R(*), C_I(*), B_R(B_m,*), B_I(B_m,*)
      integer j
      double precision Bij_R, Bij_I

      if (B_is_scalar) then
         Bij_R = B_R(1,1)
         if (B_it .eq. 1) Bij_I = B_I(1,1)
      endif

      do j = 1, B_n

         if (kc .gt. nelmax) then ! test if enough memory  
            ierr = -1
            return
         endif

         if ( .not. B_is_scalar ) then
            Bij_R = B_R(i,j)
            if (B_it .eq. 1) Bij_I = B_I(i,j)
         endif

         if (C_it .eq. 0) then      ! real case (ie A and B are real)
            if ( Bij_R .ne. 0.d0 ) then
               C_icol(kc) = j
               C_R(kc) = Bij_R
               C_mnel = C_mnel + 1
               kc = kc + 1
            endif

         else ! C_it = 1 this is the complex case (but B may be real or complex)
            if (B_it .eq. 0) then   ! B is real
               if ( Bij_R .ne. 0.d0 ) then
                  C_icol(kc) = j
                  C_R(kc) = Bij_R
                  C_I(kc) = 0.d0
                  C_mnel = C_mnel + 1
                  kc = kc + 1
               endif
            else                    ! B is complex
               if ( Bij_R.ne.0.d0 .or. Bij_I.ne.0.d0 ) then
                  C_icol(kc) = j
                  C_R(kc) = Bij_R
                  C_I(kc) = Bij_I
                  C_mnel = C_mnel + 1
                  kc = kc + 1
               endif
            endif
         endif

      enddo

      end  ! subroutine copy_fullrow2sprow

      subroutine isorti(ix, p, n)

      implicit none
      integer n, ix(n), p(n)
      
      integer i, j, kmax, temp
      call set_perm_id(p, n)
      kmax = 1
      do i = 2,n
         if (ix(i) .lt. ix(kmax)) kmax = i
      enddo
      p(1) = kmax
      p(kmax) = 1

      do i = 3, n
         j = i
         do while ( ix(p(j-1)) .gt. ix(p(j)) )
            temp = p(j)
            p(j) = p(j-1)
            p(j-1) = temp
            j = j-1
         enddo
      enddo

      end ! subroutine isorti
      
                              
      integer function dicho_search(jj, val, n)
*
*     PURPOSE
*        val(1..n) being a increasing array, this
*        routines by the mean of a dichotomic search computes
*        the smallest indice of jj in val, ie the smallest 
*        integer k such that val(k) = jj or return 0 if jj 
*        is not in val.
*
*     EXAMPLE
*        if val = [1 4 5 5 8 8 12] we must
*        get  k = 3 if jj = 5
*             k = 5 if jj = 8
*             k = 0 if jj = 3
*             etc... 
*
*     code brings by Bruno to speed up sparse extraction

      implicit none
      integer jj, n, val(n)

      integer k, k1, k2

      if ( n .lt. 1 ) then ! void search array 
         dicho_search = 0
         return
      endif

      if ( val(1) .le. jj  .and.  jj .le. val(n) ) then
*        find k such that jj = val(k) by a dicho search
         k1 = 1
         k2 = n
         do while ( k2 - k1 .gt. 1 )
            k = (k1 + k2)/2
            if ( jj .le. val(k) ) then
               k2 = k
            else
               k1 = k
            endif
         enddo
*        here we know that val(k1) <= jj <= val(k2)  with k2 = k1 + 1 
         if (jj .eq. val(k1)) then
            dicho_search = k1
         else if (jj .eq. val(k2)) then
            dicho_search = k2
         else
            dicho_search = 0
         endif
      else
         dicho_search = 0
      endif
      
      end 

      integer function dicho_search_bis(jj, val, p, n)
*
*     PURPOSE
*        same operation than the previus function but here
*        the array val is in increasing order through the
*        the permutation p, ie we have :
*            val(p(1)) <= val(p(2)) <= ... <= val(p(n))
*
*     code brings by Bruno to speed up sparse extraction

      implicit none
      integer jj, n, val(n), p(n)

      integer k, k1, k2

      if ( n .lt. 1 ) then ! void search array 
         dicho_search_bis = 0
         return
      endif

      if ( val(p(1)) .le. jj  .and.  jj .le. val(p(n)) ) then
*        find the smallest k such that jj = val(p(k)) by a dicho search
         k1 = 1
         k2 = n
         do while ( k2 - k1 .gt. 1 )
            k = (k1 + k2)/2
            if ( jj .le. val(p(k)) ) then
               k2 = k
            else
               k1 = k
            endif
         enddo
*        here we know that val(p(k1)) <= jj <= val(p(k2))  with k2 = k1 + 1 
         if (jj .eq. val(p(k1))) then
            dicho_search_bis = k1
         else if (jj .eq. val(p(k2))) then
            dicho_search_bis = k2
         else
            dicho_search_bis = 0
         endif
      else
         dicho_search_bis = 0
      endif
      
      end 

      logical function is_in_order(val, n)
*
*     PURPOSE
*        test if the array val(1..n) is in increasing order
*
*     code brings by Bruno to speed up sparse extraction
*
      implicit none
      integer n, val(n), i
      
      is_in_order = .true.

      do i = 2, n
         if ( val(i) .lt. val(i-1) ) then
            is_in_order = .false.
            return
         endif
      enddo

      end


      subroutine insert_in_order(icol, kleft, kright, ind, it, B_R, B_I, 
     $                           val_re, val_im)
*
*    expliquer ....
*
      implicit none
      integer icol(*), kleft, kright, ind, it
      double precision B_R(*), B_I(*), val_re, val_im

      integer k

      k = kright
      if (k .gt. kleft) then  ! a caution due to the fortran evaluation of .and.
         do while (k .gt. kleft  .and.  icol(k-1) .gt. ind)
            icol(k) = icol(k-1)
            if (it .ge. 0) B_R(k) = B_R(k-1)
            if (it .eq. 1) B_I(k) = B_I(k-1)
            k = k-1
         enddo
      endif
      icol(k) = ind
      if (it .ge. 0) B_R(k) = val_re
      if (it .eq. 1) B_I(k) = val_im

      end


      SUBROUTINE QSORTI (X, IND, N)
      INTEGER N, X(N), IND(N)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C
C   THIS SUBROUTINE USES AN ORDER N*LOG(N) QUICK SORT TO
C SORT AN INTEGER ARRAY X INTO INCREASING ORDER.  THE ALGOR-
C ITHM IS AS FOLLOWS.  IND IS INITIALIZED TO THE ORDERED
C SEQUENCE OF INDICES 1,...,N, AND ALL INTERCHANGES ARE
C APPLIED TO IND.  X IS DIVIDED INTO TWO PORTIONS BY PICKING
C A CENTRAL ELEMENT T.  THE FIRST AND LAST ELEMENTS ARE COM-
C PARED WITH T, AND INTERCHANGES ARE APPLIED AS NECESSARY SO
C THAT THE THREE VALUES ARE IN ASCENDING ORDER.  INTER-
C CHANGES ARE THEN APPLIED SO THAT ALL ELEMENTS GREATER THAN
C T ARE IN THE UPPER PORTION OF THE ARRAY AND ALL ELEMENTS
C LESS THAN T ARE IN THE LOWER PORTION.  THE UPPER AND LOWER
C INDICES OF ONE OF THE PORTIONS ARE SAVED IN LOCAL ARRAYS,
C AND THE PROCESS IS REPEATED ITERATIVELY ON THE OTHER
C PORTION.  WHEN A PORTION IS COMPLETELY SORTED, THE PROCESS
C BEGINS AGAIN BY RETRIEVING THE INDICES BOUNDING ANOTHER
C UNSORTED PORTION.
C
C INPUT PARAMETERS -   N - LENGTH OF THE ARRAY X.
C
C                      X - VECTOR OF LENGTH N TO BE SORTED.
C
C                    IND - VECTOR OF LENGTH .GE. N.
C
C N AND X ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETER - IND - SEQUENCE OF INDICES 1,...,N
C                          PERMUTED IN THE SAME FASHION AS X
C                          WOULD BE.  THUS, THE ORDERING ON
C                          X IS DEFINED BY Y(I) = X(IND(I)).
C
C INTRINSIC FUNCTIONS CALLED BY QSORTI - IFIX, FLOAT
C
C***********************************************************
C
C NOTE -- IU AND IL MUST BE DIMENSIONED .GE. LOG(N) WHERE
C         LOG HAS BASE 2.
C
C***********************************************************
C
      INTEGER IU(21), IL(21)
      INTEGER M, I, J, K, L, IJ, IT, ITT, INDX, T
      REAL    R
C
C LOCAL PARAMETERS -
C
C IU,IL =  TEMPORARY STORAGE FOR THE UPPER AND LOWER
C            INDICES OF PORTIONS OF THE ARRAY X
C M =      INDEX FOR IU AND IL
C I,J =    LOWER AND UPPER INDICES OF A PORTION OF X
C K,L =    INDICES IN THE RANGE I,...,J
C IJ =     RANDOMLY CHOSEN INDEX BETWEEN I AND J
C IT,ITT = TEMPORARY STORAGE FOR INTERCHANGES IN IND
C INDX =   TEMPORARY INDEX FOR X
C R =      PSEUDO RANDOM NUMBER FOR GENERATING IJ
C T =      CENTRAL ELEMENT OF X
C
      IF (N .LE. 0) RETURN
C
C INITIALIZE IND, M, I, J, AND R
C
      DO 1 I = 1,N
    1   IND(I) = I
      M = 1
      I = 1
      J = N
      R = .375
C
C TOP OF LOOP
C
    2 IF (I .GE. J) GO TO 10
      IF (R .GT. .5898437) GO TO 3
      R = R + .0390625
      GO TO 4
    3 R = R - .21875
C
C INITIALIZE K
C
    4 K = I
C
C SELECT A CENTRAL ELEMENT OF X AND SAVE IT IN T
C
      IJ = I + IFIX(R*FLOAT(J-I))
      IT = IND(IJ)
      T = X(IT)
C
C IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T,
C   INTERCHANGE IT WITH T
C
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 5
      IND(IJ) = INDX
      IND(I) = IT
      IT = INDX
      T = X(IT)
C
C INITIALIZE L
C
    5 L = J
C
C IF THE LAST ELEMENT OF THE ARRAY IS LESS THAN T,
C   INTERCHANGE IT WITH T
C
      INDX = IND(J)
      IF (X(INDX) .GE. T) GO TO 7
      IND(IJ) = INDX
      IND(J) = IT
      IT = INDX
      T = X(IT)
C
C IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T,
C   INTERCHANGE IT WITH T
C
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 7
      IND(IJ) = INDX
      IND(I) = IT
      IT = INDX
      T = X(IT)
      GO TO 7
C
C INTERCHANGE ELEMENTS K AND L
C
    6 ITT = IND(L)
      IND(L) = IND(K)
      IND(K) = ITT
C
C FIND AN ELEMENT IN THE UPPER PART OF THE ARRAY WHICH IS
C   NOT LARGER THAN T
C
    7 L = L - 1
      INDX = IND(L)
      IF (X(INDX) .GT. T) GO TO 7
C
C FIND AN ELEMENT IN THE LOWER PART OF THE ARRAY WHCIH IS
C   NOT SMALLER THAN T
C
    8 K = K + 1
      INDX = IND(K)
      IF (X(INDX) .LT. T) GO TO 8
C
C IF K .LE. L, INTERCHANGE ELEMENTS K AND L
C
      IF (K .LE. L) GO TO 6
C
C SAVE THE UPPER AND LOWER SUBSCRIPTS OF THE PORTION OF THE
C   ARRAY YET TO BE SORTED
C
      IF (L-I .LE. J-K) GO TO 9
      IL(M) = I
      IU(M) = L
      I = K
      M = M + 1
      GO TO 11
C
    9 IL(M) = K
      IU(M) = J
      J = L
      M = M + 1
      GO TO 11
C
C BEGIN AGAIN ON ANOTHER UNSORTED PORTION OF THE ARRAY
C
   10 M = M - 1
      IF (M .EQ. 0) RETURN
      I = IL(M)
      J = IU(M)
C
   11 IF (J-I .GE. 11) GO TO 4
      IF (I .EQ. 1) GO TO 2
      I = I - 1
C
C SORT ELEMENTS I+1,...,J.  NOTE THAT 1 .LE. I .LT. J AND
C   J-I .LT. 11.
C
   12 I = I + 1
      IF (I .EQ. J) GO TO 10
      INDX = IND(I+1)
      T = X(INDX)
      IT = INDX
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 12
      K = I
C
   13 IND(K+1) = IND(K)
      K = K - 1
      INDX = IND(K)
      IF (T .LT. X(INDX)) GO TO 13
      IND(K+1) = IT
      GO TO 12
      END