subroutine rscsb3 (array, right, prodx1, prodx2, nrwblk,
     *                   nbloks, beta, phi, pivot, nparts, blaws)
c
      double precision array(nrwblk,2*nrwblk,1),
     *             right(nrwblk,nrwblk,1), prodx1(nrwblk,nrwblk,1),
     *             prodx2(nrwblk,1), beta(1), phi(1), blaws(nrwblk,1)
      integer nrwblk, nbloks, pivot(1), nparts
c    ***----------------------------------------------------------***
c    *   The following notation is used in the comments:            *
c    *                                                              *
c    *      V_k, W_k <=> array(1,1,k), array(1,nrwblk+1,k)          *
c    *           R_k <=> right(1,1,k)                               *
c    *           S_k <=> prodx1(1,1,k)                              *
c    *           T   <=> prodx2(1,1)                                *
c    *          beta <=> beta(1) (<=> phi_0)                        *
c    *         phi_k <=> phi(k*nrwblk+1), k >= 1                    *
c    *                                                              *
c    *   Since [beta; phi] is overwritten with the solution,        *
c    *                                                              *
c    *      y_a <=> y_0      <=> beta                               *
c    *              y_k      <=> phi_k,  k = 1, nbloks-1            *
c    *      y_b <=> y_nbloks <=> phi(nbloks*nrwblk+1)               *
c    *                                                              *
c    *   In addition, the affix ''' designates that the solution    *
c    *   vector is obtained at the third level of the back-solve.   *
c    ***----------------------------------------------------------***
         integer info
c    ***----------------------------------------------------------***
c    *   First y_a''' is obtained by solving the reduced system     *
c    *                                                              *
c    *      [Ba'' - Bb W_nbloks^-1 T''] y_a''' = beta'''  or        *
c    *      [Ba' - Bb W_nbloks^-1 S_1'] y_a''' = beta'''            *
c    *                                                              *
c    *   if nparts > 1, or nparts = 1, respectively.                *
c    *                                                              *
c    *   In either case, the factorization and pivot indices        *
c    *   for the reduced system are stored in W_nbloks and          *
c    *   pivot_nbloks, respectively.                                *
c    *                                                              *
c    *   Note that y_a''' _is not_ y_a.  Rescaling results in       *
c    *   a change of variable which is undone at levels 2 and 1     *
c    *   of the back-solve.                                         *
c    ***----------------------------------------------------------***
         call DGETRS('N',nrwblk,1,array(1,nrwblk+1,nbloks),nrwblk,
     *               pivot(nbloks*nrwblk+1),beta(1),nrwblk,info)
c    ***----------------------------------------------------------***
c    *   Next, y_b''' is obtained from                              *
c    *                                                              *
c    *       W_nbloks y_b''' = phi_nbloks'' - T'' y_a'''  or        *
c    *       W_nbloks y_b''' = phi_nbloks' - S_1' y_a'''            *
c    *                                                              *
c    *   if nparts > 1, or nparts = 1, respectively.                *
c    *                                                              *
c    *   Note that y_b''' _is_ y_b.  Rescaling at levels 1 and 2    *
c    *   of the factorization did not change y_b.                   *
c    ***----------------------------------------------------------***
         if (nparts .gt. 1) then
            call DGEMV('N',nrwblk,nrwblk,-1.d0,prodx2,nrwblk,
     *                 beta,1,1.d0,phi(nbloks*nrwblk+1),1)
         else
            call DGEMV('N',nrwblk,nrwblk,-1.d0,prodx1,nrwblk,
     *                 beta,1,1.d0,phi(nbloks*nrwblk+1),1)
         endif
         call DGETRS('N',nrwblk,1,right(1,1,nbloks),nrwblk,
     *               pivot,phi(nbloks*nrwblk+1),nrwblk,info)
      return
      end