subroutine rscsf2 (array, right, prodx1, nrwblk, phi,
     *                   pivot, minblk, remblk, nparts, blaws)
c
      double precision array(nrwblk,2*nrwblk,1),
     *                 right(nrwblk,nrwblk,1), phi(1),
     *                 prodx1(nrwblk,nrwblk,1), blaws(nrwblk,1)
      integer nrwblk, pivot(1), minblk, remblk, 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    *         phi_k <=> phi(k*nrwblk+1)                            *
c    *                                                              *
c    *   In addition, the affix '/'' designates that the            *
c    *   vector was/is transformed at the first/second level        *
c    *   of the forward solve.                                      *
c    ***----------------------------------------------------------***
         integer kpart, base, basep, top, info
c    ***----------------------------------------------------------***
c    *      If there is only one partition, nothing needs to be     *
c    *      done at the second level of the forward solve.          *
c    ***----------------------------------------------------------***
         if (nparts .eq. 1) return
c    ***----------------------------------------------------------***
c    *      Forward elimination starts at the second-last           *
c    *      block-row of the second-level array.                    *
c    ***----------------------------------------------------------***
         call partx(minblk,remblk,nparts-1,base,top)
c    ***----------------------------------------------------------***
c    *      phi_base'' <- (W_base' - S_nparts-1')^-1 phi_base'      *
c    ***----------------------------------------------------------***
         call DGETRS('N',nrwblk,1,array(1,nrwblk+1,base),
     *               nrwblk,pivot(base*nrwblk+1),
     *               phi(base*nrwblk+1),nrwblk,info)
c    ***----------------------------------------------------------***
c    *      Forward elimination now proceeds sequentially           *
c    *      from the third-last block row to the top of the         *
c    *      second-level array.                                     *
c    ***----------------------------------------------------------***
         do 10 kpart = nparts-2, 1, -1
            basep = base
            call partx(minblk,remblk,kpart,base,top)
c    ***----------------------------------------------------------***
c    *       phi_base'' <- phi_base' + W_base' phi_basep''          *
c    *                                                              *
c    *  [ Note:  A copy of W_base' was stored in R_base-1 during    *
c    *           the first-level factorization.                  ]  *
c    ***----------------------------------------------------------***
            call DGEMV('N',nrwblk,nrwblk,
     *               1.d0,right(1,1,base-1),nrwblk,
     *               phi(basep*nrwblk+1),1,1.d0,phi(base*nrwblk+1),1)
c    ***----------------------------------------------------------***
c    *      phi_base'' <- (W_base'(I + S_kpart+1'')                 *
c    *                               - S_kpart')^-1 phi_base''      *
c    ***----------------------------------------------------------***
            call DGETRS('N',nrwblk,1,array(1,nrwblk+1,base),
     *                  nrwblk,pivot(base*nrwblk+1),
     *                  phi(base*nrwblk+1),nrwblk,info)
  10     continue
c    ***----------------------------------------------------------***
c    *   phi_base_nparts'' is now computed.  First, phi_base_1''    *
c    *   is stored in a temporary vector (temp <- phi_base_1'').    *
c    ***----------------------------------------------------------***
         call partx(minblk,remblk,1,base,top)
         call DCOPY(nrwblk,phi(base*nrwblk+1),1,blaws,1)
c    ***----------------------------------------------------------***
c    *      The temporary vector is then processed sequentially     *
c    *      from the second block-row to the bottom of the          *
c    *      second-level array.                                     *
c    ***----------------------------------------------------------***
         do 20 kpart = 2, nparts-1
            call partx(minblk,remblk,kpart,base,top)
c    ***----------------------------------------------------------***
c    *         temp <- phi_base'' - S_kpart'' temp                  *
c    ***----------------------------------------------------------***
            call DCOPY(nrwblk,phi(base*nrwblk+1),1,blaws(1,2),1)
            call DGEMV('N',nrwblk,nrwblk,
     *                 -1.d0,prodx1(1,1,kpart),nrwblk,
     *                 blaws,1,1.d0,blaws(1,2),1)
            call DCOPY(nrwblk,blaws(1,2),1,blaws,1)
  20     continue
c    ***----------------------------------------------------------***
c    *   phi_base_nparts'' <- phi_base_nparts' - S_nparts' temp     *
c    ***----------------------------------------------------------***
         call partx(minblk,remblk,nparts,base,top)
         call DGEMV('N',nrwblk,nrwblk,
     *              -1.d0,prodx1(1,1,nparts),nrwblk,
     *              blaws,1,1.d0,phi(base*nrwblk+1),1)
      return
      end