subroutine rscsf1 (array, right, nrwblk, phi,
     *                   pivot, minblk, remblk, nparts, blaws)
c
      double precision array(nrwblk,2*nrwblk,1),
     *         right(nrwblk,nrwblk,1), phi(1), blaws(nrwblk,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    *         phi_k <=> phi(k*nrwblk+1)                            *
c    *                                                              *
c    *   In addition, the affix ' designates that the vector is     *
c    *   transformed at the first level of the forward solve.       *
c    ***----------------------------------------------------------***
         integer kpart, kblok, base, base1, top, info
c    ***----------------------------------------------------------***
c    *      Each loop 30 iteration is independent and could         *
c    *      execute concurrently with the others.                   *
c    ***----------------------------------------------------------***
C$DOACROSS SHARE (array, right, nrwblk, phi,
C$&               pivot, minblk, remblk, nparts, blaws),
C$&        LOCAL (kpart, kblok, base, base1, top, info)
         do 30 kpart = 1, nparts
c    ***----------------------------------------------------------***
c    *         Forward elimination starts at the second-last        *
c    *         block-row of each partition.                         *
c    ***----------------------------------------------------------***
            call partx(minblk,remblk,kpart,base,top)
            base1 = base - 1
c    ***----------------------------------------------------------***
c    *         phi_base1' <- (W_base1 - V_base1)^-1 phi_base1       *
c    ***----------------------------------------------------------***
            call DGETRS('N',nrwblk,1,array(1,nrwblk+1,base1),
     *                  nrwblk,pivot(base1*nrwblk+1),
     *                  phi(base1*nrwblk+1),nrwblk,info)
c    ***----------------------------------------------------------***
c    *         Each partition is now processed sequentially         *
c    *         from the third-last block row to the top.            *
c    ***----------------------------------------------------------***
            do 10 kblok = base-2, top, -1
c    ***----------------------------------------------------------***
c    *        phi_kblok' <- phi_kblok + R_kblok phi_kblok+1'        *
c    ***----------------------------------------------------------***
               call DGEMV('N',nrwblk,nrwblk,
     *                    1.d0,right(1,1,kblok),nrwblk,
     *                    phi((kblok+1)*nrwblk+1),1,
     *                    1.d0,phi(kblok*nrwblk+1),1)
c    ***----------------------------------------------------------***
c    *        phi_kblok' <- (W_kblok(I + V_kblok+1')                *
c    *                               - V_kblok)^-1 phi_kblok'       *
c    ***----------------------------------------------------------***
               call DGETRS('N',nrwblk,1,array(1,nrwblk+1,kblok),
     *                     nrwblk,pivot(kblok*nrwblk+1),
     *                     phi(kblok*nrwblk+1),nrwblk,info)
  10        continue
c    ***----------------------------------------------------------***
c    *      phi_base is now computed.  First, phi_top is stored     *
c    *      in a temporary vector (temp <- phi_top).                *
c    ***----------------------------------------------------------***
            call DCOPY(nrwblk,phi(top*nrwblk+1),1,blaws(1,1,kpart),1)
c    ***----------------------------------------------------------***
c    *      The temporary vector is then processed sequentially     *
c    *      from the second block-row to the bottom.                *
c    ***----------------------------------------------------------***
            do 20 kblok = top+1, base-1
c    ***----------------------------------------------------------***
c    *            temp <- phi_kblok' - V_kblok' temp                *
c    ***----------------------------------------------------------***
               call DCOPY(nrwblk,phi(kblok*nrwblk+1),1,
     *                           blaws(1,2,kpart),1)
               call DGEMV('N',nrwblk,nrwblk,
     *                    -1.d0,array(1,1,kblok),nrwblk,
     *                    blaws(1,1,kpart),1,1.d0,blaws(1,2,kpart),1)
               call DCOPY(nrwblk,blaws(1,2,kpart),1,
     *                           blaws(1,1,kpart),1)
  20        continue
c    ***----------------------------------------------------------***
c    *         phi_base' <- phi_base - V_base temp                  *
c    ***----------------------------------------------------------***
            call DGEMV('N',nrwblk,nrwblk,
     *                 -1.d0,array(1,1,base),nrwblk,
     *                 blaws(1,1,kpart),1,1.d0,phi(base*nrwblk+1),1)
  30     continue
c    ***----------------------------------------------------------***
c    *      Finally, phi_base is updated in each partition by       *
c    *      computing a vector sum across partition boundaries.     *
c    *      (This could be done concurrently.)                      *
c    ***----------------------------------------------------------***
         do 40 kpart = 1, nparts - 1
c    ***----------------------------------------------------------***
c    *         phi_base' <- phi_base' + R_base phi_base+1'          *
c    ***----------------------------------------------------------***
            call partx(minblk,remblk,kpart,base,top)
            call DGEMV('N',nrwblk,nrwblk,
     *                 1.d0,right(1,1,base),nrwblk,
     *                 phi((base+1)*nrwblk+1),1,
     *                 1.d0,phi(base*nrwblk+1),1)
  40     continue
      return
      end