subroutine couple (lftblk, nrwtop, array, nrwblk, nbloks,
     *                   rgtblk, nrwbot, b, nparts, work)
c
      double precision lftblk(1), array(1), rgtblk(1), b(1), work(1)
      integer nrwtop, nrwblk, nbloks, nrwbot, nparts
c    ***----------------------------------------------------------***
c    *   'couple' inputs an ABD system in uncoupled-boundary        *
c    *   right-block-general format:                                *
c    *                                                              *
c    *          lftblk                                b_0t          *
c    *          array(,,1)                            b_1           *
c    *               array(,,2)                       b_2           *
c    *                    .                            .            *
c    *                      .                          .            *
c    *                        .                        .            *
c    *                           array(,,nbloks)      b_nbloks      *
c    *                                    rgtblk      b_0b          *
c    *                                                              *
c    *   lftblk and rgtblk are nrwtopxnrwblk and nrwbotxnrwblk,     *
c    *   respectively, array(,,k) is nrwblkx2*nrwblk, array(,,k)    *
c    *   and array(,,k+1) overlap by nrwblk columns, b_k refers     *
c    *   to the sub-vector of b ranging from k*nrwblk+1-nrwbot to   *
c    *   (k+1)*nrwblk-nrwbot, and b_0t/b_0b hold the first/last     *
c    *   nrwtop/nrwbot elements of b.                               *
c    *                                                              *
c    *   A sequence of row permutations are used to transform       *
c    *   this system so that on output it is in coupled-boundary    *
c    *   right-block-general format:                                *
c    *                                                              *
c    *          lftblk                    rgtblk      b_0           *
c    *          array(,,1)                            b_1           *
c    *               array(,,2)                       b_2           *
c    *                    .                            .            *
c    *                      .                          .            *
c    *                        .                        .            *
c    *                           array(,,nbloks)      b_nbloks      *
c    *                                                              *
c    *   Here lftblk is nrwblkxnrwblk (nrwbot rows of zero have     *
c    *   been appended to the bottom), rgtblk is nrwblkxnrwblk      *
c    *   (nrwtop rows of zero have been inserted at the top),       *
c    *   array is unchanged, b_k, k = 1, nbloks is unchanged,       *
c    *   and b_0 is the catenation of b_0t and b_0b.                *
c    *                                                              *
c    *   After transforming the initial system, subsequent right-   *
c    *   hand sides may be transformed by calling cplrhs directly.  *
c    *                                                              *
c    *   PARAMETERS:                                                *
c    *                                                              *
c    *    on entry                                                  *
c    *                                                              *
c    *      lftblk, nrwtop, array, nrwblk, nbloks,                  *
c    *      rgtblk, nrwbot, b                                       *
c    *             The ABD system in uncoupled-boundary right-      *
c    *             block general format as described above.         *
c    *                                                              *
c    *      nparts [integer]                                        *
c    *             The requested number of partitions.  If nparts   *
c    *             processors are available, each partition may     *
c    *             be processed concurrently with the others when   *
c    *             transforming the right-hand-side b.              *
c    *                                                              *
c    *        work [double precision(nparts*nrwblk**2)]             *
c    *             Work space to hold a temporary copy of the       *
c    *             bottom nrwbot slice of b in each partition.      *
c    *             This workspace is required to be larger than     *
c    *             its optimal size of nparts*nrwbot in order to    *
c    *             avoid processor contention for cached memory     *
c    *             in the parallel implementation.                  *
c    *                                                              *
c    *   on return                                                  *
c    *                                                              *
c    *      lftblk, nrwtop, array, nrwblk, nbloks,                  *
c    *      rgtblk, nrwbot, b                                       *
c    *             The ABD system in coupled-boundary right-block   *
c    *             general rotated format as described above.       *
c    *                                                              *
c    *      nparts [integer]                                        *
c    *             Normally unchanged.  If, however, the            *
c    *             requested number of partitions would             *
c    *             result in fewer than 2 nrwblk slices of b        *
c    *             per partition (i.e. if nbloks < 2*nparts),       *
c    *             the subroutine automatically resets nparts       *
c    *             to 1 and uses non-partitioned coupling.          *
c    ***----------------------------------------------------------***
         call cpllhs(lftblk,nrwtop,nrwblk,rgtblk,nrwbot)
         call cplrhs(nrwtop,nrwblk,nbloks,nrwbot,b,nparts,work)
      return
      end