Python3, 475 bytes

def C(p,L):
 q,s=[([],[],p)],[]
 for a,b,p in q:
  if len(a)==L and(a,b+p)not in s:yield a,b+p;s+=[(a,b+p)]
  if len(a)<L and p:
   q+=[(a,b+[p[0]],p[1:])]
   if len(V:=a+p[0])<=L:q+=[(V,b,p[1:])]
def f(n):
 W=[[i]*i for i in range(1,n+1)for _ in range(i)]
 n-=1
 while(n:=n+1):
  if sum(map(len,W))%n:continue
  q=[([a],b)for a,b in C(W,n)]
  while q:
   (a,p),*q=sorted(q,key=lambda x:sum(map(len,x[1]))) 
   if[]==p and len(a)>1:return a
   q+=[(a+[A],P)for A,P in C(p,n)]

Try it online!

Wolfram Language (Mathematica), 148 bytes

finds random matrix with only horizontal tiles.
It is very slow for n>5 but this is code-golf...

(While@!MatrixQ[Join@@@(w=TakeList[d=RandomSample@Flatten[Table[#~Table~#,#]&/@Range@#,1],r=RandomChoice@IntegerPartitions[#(#+1)/2][[2;;-2]]])];w)&

Try it online!

here is also a very quick random generator for test cases up to 6

Wolfram Language (Mathematica), 97 bytes

(While@!MatrixQ[Join@@@(w=Partition[RandomSample@Flatten[Table[#~Table~#,#]&/@Range@#,1],3])];w)&

Try it online!

Python 2, 238 ... 186 bytes

-2 bytes thanks to @Arnauld!

-5 bytes thanks to @ovs!!

-6 bytes thanks to @Jonathan Allan!!

n=input()
p=w=~n*n*(~n-n)/6
e=range(n+1)
def f(a,b,c,d=[]):1>a>exit(c);1>b>f(a-1,w,c+[d]);g=n;exec"e[g]-=1;g<=b>-1<e[g]>f(a,b-g,c,d+[g]*g);e[g]+=1;g-=1;"*n
while 1:1>p%w>f(p/w,w,[]);w-=1

Try it online!

Notice that the above code only used horizontal tiles. To prove it's correct...

C++ (gcc), 4355 bytes

This code is an explicit construction using only horizontal tiles for \$n>6\$. It can produce correct solutions (at least) up to \$n=99\$.

#include <bits/stdc++.h>
using namespace std;
void dim(int x,int&p,int&q)
{
	int u[]={x,x+1,x+x+1},o=0;
	while(u[o]%2) ++o; u[o]/=2;
	o=0; while(u[o]%3) ++o; u[o]/=3;
	sort(u,u+3);if(x%3==1)swap(u[0],u[1]);
	p=u[0];q=u[1]*u[2];
}
int o[100][99999],c[99999],g[99999];
void brute(int n,int& v,int& u,bool t)
{
	for(int i=1;i<=n;++i) c[i]=i;
	if(t)
	{
		u=n*(n+1)*(2*n+1)/6,v=2;
		while(u%v) ++v;
		if(n>=5) v=u/v;
		u/=v; c[2]-=n==6;
	}
	for(int i=1;i<=v;++i)
	{
		int s=u,cnt=0;
		for(int j=n;j>=1;--j) while(s>=j&&c[j])
		{
			s-=j, --c[j]; for(int k=j;k--;) o[i][++cnt]=j;
		}
	}
	c[2]+=t&&n==6;
	for(int i=u;i>=u-1&&i;--i)
		for(int j=1;j<=v;++j) if(!o[j][i])
		{
			int x=0; while(!c[x]) ++x;
			--c[x]; int l=x; --j;
			while(l--) o[++j][i]=x;
		}
}
int s0[99999],s1[99999];
//find a subset of a with sum b
//guaranteed b is half of a's sum
//guaranteed a is consecutive
pair<vector<int>,vector<int>> solve2(vector<int> a,int b)
{
	if(!a.size()) return make_pair(a,a);
	int w=0,s=0,as=a.size();
	for(int i=0;i<as;++i)
		s0[i+1]=s0[i]+a[i],
		s1[i+1]=s1[i]+a[as-1-i];
	while(w<a.size()&&s+a[w]<=b)
		s+=a[w++];
	assert(w!=a.size()&&w>0);
	for(int l=0;l<w;++l)
	{
		int r=w-1-l;
		int p=s0[l]+s1[r];
		if(p>b) continue;
		int q=b-p;
		if(!(a[l]<=q&&q<=a[as-r-1]))
			continue;
		//first l, last r, q.
		vector<int> A,B;
		for(int j=0;j<l;++j) A.push_back(a[j]);
		A.push_back(q);
		for(int j=as-r;j<as;++j) A.push_back(a[j]);
		for(int j=l;j<as-r;++j)
		{
			if(a[j]==q) q=-1;
			else B.push_back(a[j]);
		}
		return make_pair(A,B);
	}
	assert(0);
}
int main()
{
	int n,p,q;
	cin>>n;
	assert(n>=3);
	if(n<=6) brute(n,p,q,1);
	else
	{
		if(n%3==1)
		{
			int x=n%3,pp,qq;
			dim(x,p,q);
			brute(x,p,q,0);
			while(n!=x)
			{
				//x+1...x+3
				int xx=x+3; dim(xx,pp,qq);
				assert(pp-p==2&&qq-q==x*3+6);
				for(int i=x+1;i<=x+3;++i) c[i]=i;
				int mr=x/6;
				for(int i=1;i<=p;++i)
				{
					int cnt=q;
					vector<int> rv;
					if(i<=mr) rv=vector<int>{x+2,x+2,x+2};
					else rv=vector<int>{x+1,x+2,x+3};
					for(int j:rv)
					{
						assert(c[j]); --c[j];
						for(int k=j;k--;) o[i][++cnt]=j;
					}
				}
				vector<int> rv;
				for(int i=x+1;i<=x+3;++i)
					for(int j=c[i];j;--j) rv.push_back(i);
				pair<vector<int>,vector<int>> s=solve2(rv,qq);
				for(int i=p+1;i<=pp;++i)
				{
					int cnt=0;
					for(auto j:(i==pp)?s.first:s.second)
					{
						assert(c[j]); --c[j];
						for(int k=j;k--;) o[i][++cnt]=j;
					}
				}
				p=pp; q=qq; x=xx;
			}
		}
		else
		{
			int x=n%6,pp,qq;
			dim(x,p,q);
			brute(x,p,q,0);
			while(n!=x)
			{
				int xx=x+6; dim(xx,pp,qq);
				if(n%6==0)
					assert(pp-p==1&&qq-q==24*x+90); //[1 2 3 4 5 6]*3+[2 3 5 5 6 6]
				else if(n%6==2)
					assert(pp-p==2&&qq-q==12*x+39); //[1 2 3 4 5 6]+[1 1 2 3 5 6]
				else if(n%6==3)
					assert(pp-p==2&&qq-q==12*x+45); //[1 2 3 4 5 6]+[1 2 4 5 6 6]
				else if(n%6==5)
					assert(pp-p==1&&qq-q==24*x+78); //[1 2 3 4 5 6]*3+[1 1 2 2 4 5]
				else assert(0);
				vector<int> uv;
				if(n%6==0) uv=vector<int>{2,3,5,5,6,6};
				else if(n%6==2) uv=vector<int>{1,1,2,3,5,6};
				else if(n%6==3) uv=vector<int>{1,2,4,5,6,6};
				else uv=vector<int>{1,1,2,2,4,5};
				for(int j=1;j<4/(pp-p);++j)
					for(int k=1;k<=6;++k) uv.push_back(k);
				for(int i=x+1;i<=x+6;++i) c[i]=i;
				for(int i=1;i<=p;++i)
				{
					int cnt=q;
					for(auto j_:uv)
					{
						int j=j_+x;
						assert(c[j]); --c[j];
						for(int k=j;k--;) o[i][++cnt]=j;
					}
				}
				vector<int> rv;
				for(int i=x+1;i<=x+6;++i)
					for(int j=c[i];j;--j) rv.push_back(i);
				if(pp-p==2)
				{
					pair<vector<int>,vector<int>> s=solve2(rv,qq);
					for(int i=p+1;i<=pp;++i)
					{
						int cnt=0;
						for(auto j:(i==pp)?s.first:s.second)
						{
							assert(c[j]); --c[j];
							for(int k=j;k--;) o[i][++cnt]=j;
						}
					}
				}
				else
				{
					int i=pp,cnt=0;
					for(auto j:rv)
					{
						assert(c[j]); --c[j];
						for(int k=j;k--;) o[i][++cnt]=j;
					}
				}
				p=pp; q=qq; x=xx;
			}
		}
	}
	cerr<<p<<","<<q<<"\n";
	for(int i=1;i<=p;++i,cout<<"\n")
		for(int j=1;j<=q;++j)
			cout<<setw(2)<<o[i][j]<<" ";
	if(n>6)
	{
		for(int i=1;i<=n;++i) c[i]=i;
		for(int i=1;i<=p;++i)
			for(int j=1;j<=q;++j)
			{
				int u=o[i][j];
				assert(u>=1&&u<=n);
				for(int k=2;k<=u;++k)
					assert(o[i][++j]==u);
				--c[u];
			}
		for(int i=1;i<=n;++i) assert(!c[i]);
	}
}

Try it online!

How does the construction work?

It's an incremental construction. Consider \$n \bmod 6\$, we can have these values for height and width of rectangles:

• \$n/6\times (n+1)(2n+1)~(n\bmod 6=0)\$

• \$(2n+1)/3\times n(n+1)/2~(n\bmod 6=1)\$

• \$(n+1)/3\times n(2n+1)/2~(n\bmod 6=2)\$

• \$n/3\times (n+1)(2n+1)/2~(n\bmod 6=3)\$

• \$(2n+1)/3\times (n+1)n/2~(n\bmod 6=4)\$

• \$(n+1)/6\times n(2n+1)~(n\bmod 6=5)\$

(dimensions might be \$1\$ for \$n\leq 6\$ so these small cases are handled manually)

So the main idea of my construction is:

• We construct rectangles with the heights and widths as in above list.

• If \$n\bmod 3 \neq 1\$, construct the solution for \$n-6\$ recursively, add in \$n-5,n-4\cdots n\$. The height of the rectangle will only increase 1 or 2.

• If \$n\bmod 3=1\$, construct the solution for \$n-3\$ recursively and add in \$n-2,n-1,n\$. The height of the rectangle will only increase 2.

• We first carefully assign the new numbers to the added columns, and then put rest of the numbers into one or two added rows.

The rest of the job is some careful casework to pick the numbers. These details are left as an exercise to the readers (Be ready for some long and tedious casework!). If you complete all the details, this should become a formal proof for the existence of the solutions (and using only horizontal tiles!).

JavaScript (ES6), 188 bytes

A shorter version inspired by @newbie's answer.

f=(n,w=2)=>(g=(h,a,m,r=[])=>h%1||r[w]?1:r[w-1]?--h*g(h,a,M=[...m,r]):a.every((_,j,[...a])=>a[j]++>j++||g(h,a,m,[...r,...Array(j).fill(j)])))(n*(~n-n)*~n/6/w,Array(n).fill(0),[])?f(n,w+1):M

Try it online!

Commented

f = (                           // f is a recursive function taking:
  n, w = 2                      //   n = input, w = width of matrix
) => (                          //
g = (                           // g is a recursive function taking:
  h, a,                         //   h = height of matrix, a[] = array of counters
  m, r = []                     //   m[] = matrix, r[] = current row
) =>                            // (g returns 0 for success or 1 for failure)
  h % 1 ||                      // if h is not an integer or
  r[w] ?                        // the length of r[] is w + 1 (i.e. r[] is too long):
    1                           //   abort
  :                             // else:
    r[w - 1] ?                  //   if the length of r[] is w:
      --h *                     //     decrement h and force success if h = 0
      g(                        //     do a recursive call with:
        h, a,                   //       h and a[] unchanged
        M = [...m, r]           //       a new matrix M[] with r[] appended
      )                         //     end of recursive call
    :                           //   else:
      a.every((_, j, [...a]) => //     for each entry at position j in a[]:
        a[j]++ > j++ ||         //       unless a[j] is greater than j,
        g(                      //       do a recursive call with:
          h, a, m,              //         h, a[] and m[] unchanged
          [ ...r,               //         j added j times to the current row
            ...Array(j).fill(j) //         NB1: both j and a[j] where incremented above
          ]                     //         NB2: a[] is a local copy defined in this loop
        )                       //       end of recursive call
      )                         //     end of every()
)(                              // initial call to g with:
  n * (~n - n) * ~n / 6 / w,    //   h = n(n+1)(2n+1) / 6 / w
  Array(n).fill(0),             //   a[] initialized to n 0's
  []                            //   an empty matrix
) ? f(n, w + 1) : M             // return M[] on success, or try again with w + 1

JavaScript (ES6),  297  287 bytes

A brute-force search, which always try to put the biggest available rectangles first.

f=(n,i=2,k=n*(~n-n)*~n/6,A=n=>n?[0,...A(n-1)]:[])=>k%i||!(g=(m,a,x,y=m.findIndex(r=>r.some(v=>!v*~++x,x=-1)))=>~y?a.some((v,j)=>[0,1].some(r=>v<(o=n-j)&o<=(r?i-y:k/i-x)&&g(M=m.map(r=>[...r]),b=[...a],b[(h=p=>p--?h(p,M[y+r*p][x+!r*p]=o):j)(o)]++))):1)(A(i).map(_=>A(k/i)),A(n))?f(n,i+1):M

Try it online!