subroutine rscsf3 (lftblk, array, right, nrwblk, rgtblk,
     *               beta, phi, pivot, minblk, remblk, nparts, blaws)
c
      double precision lftblk(nrwblk,1), array(nrwblk,2*nrwblk,1),
     *                 right(nrwblk,nrwblk,1), rgtblk(nrwblk,1),
     *                 beta(1), phi(1), blaws(nrwblk,1)
      integer nrwblk, pivot(1), minblk, remblk, nparts
c    ***----------------------------------------------------------***
c    *   The following notation is used in the comments:            *
c    *                                                              *
c    *      B_a, B_b <=> lftblk, rgtblk                             *
c    *      V_k, W_k <=> array(1,1,k), array(1,nrwblk+1,k)          *
c    *           R_k <=> right(1,1,k)                               *
c    *          beta <=> beta(1) (<=> phi_0)                        *
c    *         phi_k <=> phi(k*nrwblk+1), k >= 1                    *
c    *                                                              *
c    *   In addition, the affix '/'' designates that the            *
c    *   vector was (or at least could have been) transformed       *
c    *   at the first/second level of the forward solve.            *
c    *   (In order to be consistent with other variants of          *
c    *   this algorithm, all transformations to beta occur          *
c    *   at the third level of the forward solve.)                  *
c    ***----------------------------------------------------------***
         integer base, top, info
c
         if (nparts .gt. 1) then
c    ***----------------------------------------------------------***
c    *      If there is more than one partition, the third-level    *
c    *      system is of the form                                   *
c    *                                                              *
c    *       [B_a''  B_b        [y_a  = [beta''                     *
c    *         T''   W_nbloks ]  y_b]    phi_base_nparts''],        *
c    *                                                              *
c    *      where B_a'', T'', and phi_base_nparts'' are as          *
c    *      described in rscf3.f and rscsf2.f, and                  *
c    *                                                              *
c    *         beta'' = beta' + B_a' phi_base_1''                   *
c    *                = beta' + (B_a + B_a V_1') phi_base_1''       *
c    *                = beta + B_a phi_1'                           *
c    *                       + (B_a + B_a V_1') phi_base_1''        *
c    *                = beta + B_a (phi_1' + phi_base_1''           *
c    *                                     + V_1' phi_base_1'')     *
c    *                                                              *
c    *      The right-hand side of the reduced nrwblkxnrwblk        *
c    *      system for y_a (see rscf3.f) is then                    *
c    *                                                              *
c    *      beta''' = beta'' - B_b W_nbloks^-1 phi_base_nparts''.   *
c    ***----------------------------------------------------------***
            call partx(minblk,remblk,1,base,top)
c    ***----------------------------------------------------------***
c    *      beta'' <- beta + B_a (phi_1' + phi_base_1''             *
c    *                                   + V_1' phi_base_1'')       *
c    ***----------------------------------------------------------***
            call DCOPY(nrwblk,phi(base*nrwblk+1),1,blaws,1)
            call DGEMV('N',nrwblk,nrwblk,1.d0,array,nrwblk,
     *                 phi(base*nrwblk+1),1,1.d0,blaws,1)
            call DAXPY(nrwblk,1.d0,phi(nrwblk+1),1,blaws,1)
            call DGEMV('N',nrwblk,nrwblk,1.d0,lftblk,nrwblk,
     *                 blaws,1,1.d0,beta,1)
c    ***----------------------------------------------------------***
c    *      beta''' <- beta'' - B_b W_nbloks^-1 phi_base_nparts''   *
c    *                                                              *
c    *  [ Note:  The LU factorization and pivot indices             *
c    *           for W_nbloks were stored in R_nbloks and           *
c    *           pivot(1..nrwblk), respectively, during the         *
c    *           third level of the factorization.               ]  *
c    ***----------------------------------------------------------***
            call partx(minblk,remblk,nparts,base,top)
            call DCOPY(nrwblk,phi(base*nrwblk+1),1,blaws,1)
            call DGETRS('N',nrwblk,1,right(1,1,base),nrwblk,
     *                  pivot,blaws,nrwblk,info)
            call DGEMV('N',nrwblk,nrwblk,-1.d0,rgtblk,nrwblk,
     *                 blaws,1,1.d0,beta,1)
         else
c    ***----------------------------------------------------------***
c    *      If there is only one partition, the third-level         *
c    *      system is of the form                                   *
c    *                                                              *
c    *             [B_a'  B_b        [y_a  = [beta'                 *
c    *              S_1'  W_nbloks ]  y_b]    phi_base_1'],         *
c    *                                                              *
c    *      where B_a', S_1' are as described in rscf3.f, and       *
c    *                                                              *
c    *         beta' = beta + B_a phi_1'                            *
c    *                                                              *
c    *      The right-hand side of the reduced nrwblkxnrwblk        *
c    *      system for y_a (see rscf3.f) is then                    *
c    *                                                              *
c    *         beta''' = beta' - B_b W_nbloks^-1 phi_base_1'.       *
c    ***----------------------------------------------------------***
            call partx(minblk,remblk,1,base,top)
c    ***----------------------------------------------------------***
c    *         beta' <- beta + B_a phi_1'                           *
c    ***----------------------------------------------------------***
            call DGEMV('N',nrwblk,nrwblk,1.d0,lftblk,nrwblk,
     *                 phi(nrwblk+1),1,1.d0,beta,1)
c    ***----------------------------------------------------------***
c    *         beta''' <- beta' - B_b W_nbloks^-1 phi_base_1'       *
c    ***----------------------------------------------------------***
            call DCOPY(nrwblk,phi(base*nrwblk+1),1,blaws,1)
            call DGETRS('N',nrwblk,1,right(1,1,base),nrwblk,
     *                  pivot,blaws,nrwblk,info)
            call DGEMV('N',nrwblk,nrwblk,-1.d0,rgtblk,nrwblk,
     *                 blaws,1,1.d0,beta,1)
         endif
      return
      end