How to form groups before randomizing the treatment assignment? Consider the following example. You want to assess the effect of a new class program on educational outcomes relative to the old program. You have recruited $N$ subjects for a randomized control trial and have obtained a measure of their prior skills (e.g., through a test conducted at the registration) for each of them. 
Your task is to provide an experimental design that consists of a protocol for 1. sorting the subjects into $C$ classrooms and 2. assign the resulting classrooms to either the new or the old program.
What is the optimal experimental design?
My first solution was to randomize subjects into the $C$ classes (forming classes of size $s=N/C$) and then randomly assign half of the classes to the new program and the remainder to the old program. 
Can I do better? How to use the baseline skill measure to improve my solution?
 A: A reasonable approach here is to use block randomisation, where you create non-randomised blocks of subjects, grouping together like subjects (e.g., by prior skill variables), and then you create random groups by randomly allocating people from the blocks into the different treatment groups.  This gives you randomised treatment groups, but it also reduces colinearity between the treatment group and the prior skill variable, which has later advantages when you come to analyse your model via regression methods.
For simplicity, suppose you have $N = sC$ subjects, where $s \in \mathbb{N}$ (i.e., suppose your number of subjects is an exact multiple of the desired group size), and suppose that each subject has a covariate $z_i$ representing their prior skill.  In this case, you would form $s$ blocks of $C$ subjects by ordering the subjects by their score and then taking consecutive blocks over those scores.  So you would have one group of $C$ subjects with the lowest scores, another group of $C$ subjects with the next lowest scores, etc.  You would then allocated each block of $C$ subjects to the $C$ random treatment groups (via a random permutation).  For the simple case where you have an exact multiple of the desired group size, you can implement this in R as follows:
#Generate mock data containing ID and Score for each subject
N    <- 400;
DATA <- data.frame(ID = 1:N, Score = ceiling(runif(N)*100));

#Allocated blocks to subjects
S    <- 40;
C    <- N/S;
BBB  <- rep(1:S, each = C);
RRR  <- rank(DATA$Score, ties = 'first');
DATA$Block <- BBB[RRR];

#Randomise into treatments
set.seed(12345);
TTT  <- rep(0, N);
for (s in 1:S) { TTT[((s-1)*C+1):(s*C)] <- order(runif(C)); }
DATA$Treat <- TTT[RRR];

A: You want to balance prior skills across classes, because they may influence the outcome and confound the results. This is analogous to clinical trials where covariates that are known to influence the prognosis need to be controlled for. This is done at the level of assignment through stratified randomization (https://www.statisticshowto.datasciencecentral.com/stratified-randomization/).
In your case, that would consist in dividing the sample of baseline measures in $k$ quantiles. Then, each of the $k$ subgroups of students is randomly but evenly assigned to the $C$ classes. A good value of $k$ will depend on $N$ and $C$.
A: What is nice is that regardless of how you allocate the students to the classes, as long as you randomly assign the classes to treatment, your effect estimate will be free of confounding. Ideally, though, you want the class composition to be the same on average between the treatment groups (i.e., for the distribution of baseline skill to be the same across treated classrooms as across control classrooms). Here's one way you could do this:
Perform a match on students with respect to baseline skill. That means for each student, find a student with similar baseline skill and make them a pair. You could try to minimize some global measure of pairwise imbalance to ensure no one is paired with someone too different from them. For example, if you were to use a greedy algorithm, the first students to be paired would be paired with students close to them, while the last students to be paired will have to be paired with whoever is left, who might not be so close. An optimal matching algorithm might instead yield a set of pairs for which the distance between each member is small on average, and large distance are penalized. One member of each pair will eventually go to treatment and the other to control. This type of matching (i.e., before treatment has been assigned) is called non-bipartite matching. There are a few software packages that can do this. One I found is the nbpMatching R package. I believe the designmatch R package can do this as well.
Once you have your matched pairs, we want to assign the pairs into classroom blocks. It doesn't really matter how this is done because the treatment effect will be unbiased regardless, but to reduce the within-classroom treatment effect variability you could create strata of the pairs based on average pairwise baseline skill and then make each stratum (or randomly slices block of a stratum) into a classroom-block. This way, each classroom-block will contain pairs roughly homogeneous in baseline skill level. Within each classroom block, randomly assign one member of each pair to treatment and the other to control. The treated students become their own classroom and the control students become their own classroom (i.e., each classroom-block splits into two classrooms, one treated and the other control).
The value of this design is that it makes extensive use of pairing and stratification. Students are cross-classified both into student-pairs and classrooms. Because the students within pairs are similar to each other, conditioning on pair membership in the analysis will greatly improve the precision of your effect estimate. Likewise, because classrooms within classroom-blocks are similar to each other, conditioning on classroom-block membership will improve the precision of the estimate. Allowing the treatment effect to vary across classroom-blocks also allows you to assess the extent to which the average classroom-level baseline skill might change the treatment effect.
Your final analysis could be a multilevel model with student performance as the outcome, a random effect for student-pairs, fixed or random effects for classrooms, fixed effects for classroom-blocks, and a treatment by classroom-block interaction. The effect of treatment overall could be computed as the average of the classroom-block-specific treatment effects, and you can use the classroom-block-specific treatment effect to assess the degree to which the treatment effect varies across classrooms with different levels of average baseline skill. If that seems like too much, you can run a GEE or regression model of the outcome on treatment alone (or with some fixed effects) that accounts for the pairing and nesting within classrooms and classroom-blocks using cluster-robust standard errors. No matter what analysis you perform, as long as you are comparing the treatment group means, you effect estimate will be unbiased. All this design stuff is just a matter of improving precision and adding interpretational benefits.
A: You can do better by randomizing treatment at the student level and then forming the classes deterministically. Consider the following procedure. 
1) Match the N students based on baseline skill and any other important covariates you want to include, e.g. gender, race, etc. This will form N/2 student pairs. References on how to do this are below.
2) Randomly assign one student in each pair to treatment A and the other to treatment B. 
The randomization is now done. Your next steps are to form the classrooms. The following forms classrooms balanced on baseline skill.
3) Create C/2 class pairs labeled (1-A, 1-B), (2-A, 2B), ..., (C/2-A, C/2-B).
4) Sort the pairs of students by the mean baseline skill score within each pair. Don't sort by individual skill scores or any individual variables. Only sort by the pairs' mean values within each pair.
5) Place the highest ranked pair of students in class 1 with the student randomized to A going into 1-A and the student randomized to B going into 1-B. Place the next highest rank pair into class 2. Repeat until all C/2 class pairs have 1 student in them. Then start again at class 1 for the next highest ranked pair of students.
You now have C/2 classes that are fairly well balanced between each other on baseline skill. However, that isn't essential to your inference. What is important is that you have treatment A and B very well balanced on the student level and your randomization has occured at the student level. 
If you randomize the classes to A and B, you have a cluster randomized trial and should analyze it as such. Here you have a student randomized trial where the intervention is delivered in clusters. This makes it easier to justify a wider range of analysis methods. 
I really like Noah's answer in this thread. Our ideas are very similar. I recommend reading his post too.
References:
1) Optimal multivariate matching before randomization: 
https://www.ncbi.nlm.nih.gov/pubmed/15054030
2) Optimal Nonbipartite Matching and Its Statistical Applications: 
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3501247/
3) nbpMatching: 
https://cran.r-project.org/web/packages/nbpMatching/nbpMatching.pdf
4) Tutorials: 
http://biostat.mc.vanderbilt.edu/MatchedRandomization
